BICPL: find_in_direction.c Source File

00001 #include <volume_io/internal_volume_io.h> 00002 #include <bicpl/deform.h> 00003 00004 private BOOLEAN voxel_might_contain_boundary( 00005 Volume volume, 00006 bitlist_3d_struct *done_bits, 00007 bitlist_3d_struct *surface_bits, 00008 int degrees_continuity, 00009 int voxel_indices[], 00010 boundary_definition_struct *boundary_def ); 00011 00012 private BOOLEAN check_voxel_for_isovalue( 00013 voxel_coef_struct *lookup, 00014 Volume volume, 00015 Volume label_volume, 00016 int degrees_continuity, 00017 int voxel_indices[], 00018 Real line_origin[], 00019 Real line_direction[], 00020 BOOLEAN searching_inwards, 00021 Real min_line_t, 00022 Real max_line_t, 00023 Real isovalue, 00024 Normal_directions normal_direction, 00025 BOOLEAN already_found, 00026 Real *dist_found ); 00027 00028 private BOOLEAN check_voxel_for_boundary( 00029 voxel_coef_struct *lookup, 00030 Volume volume, 00031 Volume label_volume, 00032 int degrees_continuity, 00033 int voxel_indices[], 00034 Real line_origin[], 00035 Real line_direction[], 00036 BOOLEAN *inside_surface, 00037 BOOLEAN searching_inwards, 00038 Real min_line_t, 00039 Real max_line_t, 00040 boundary_definition_struct *boundary_def, 00041 BOOLEAN already_found, 00042 Real *dist_found ); 00043 00044 #define GET_RAY_POINT( result, origin, direction, dist ) \ 00045 { \ 00046 (result)[X] = (origin)[X] + (dist) * (direction)[X]; \ 00047 (result)[Y] = (origin)[Y] + (dist) * (direction)[Y]; \ 00048 (result)[Z] = (origin)[Z] + (dist) * (direction)[Z]; \ 00049 } 00050 00051 private void clip_line_to_volume( 00052 Volume volume, 00053 int degrees_continuity, 00054 Real line_origin[], 00055 Real line_direction[], 00056 Real *t_min, 00057 Real *t_max ) 00058 { 00059 int sizes[N_DIMENSIONS]; 00060 Point origin; 00061 Vector direction; 00062 00063 get_volume_sizes( volume, sizes ); 00064 00065 fill_Point( origin, line_origin[X], line_origin[Y], line_origin[Z] ); 00066 fill_Point( direction, 00067 line_direction[X], line_direction[Y], line_direction[Z] ); 00068 00069 (void) clip_line_to_box( &origin, &direction, 00070 0.0 + (Real) degrees_continuity * 0.5, 00071 (Real) sizes[X]-1.0 - (Real)degrees_continuity*0.5, 00072 0.0 + (Real) degrees_continuity * 0.5, 00073 (Real) sizes[Y]-1.0 - (Real)degrees_continuity*0.5, 00074 0.0 + (Real) degrees_continuity * 0.5, 00075 (Real) sizes[Z]-1.0 - (Real)degrees_continuity*0.5, 00076 t_min, t_max ); 00077 } 00078 00079 private void set_up_to_search_ray( 00080 Volume volume, 00081 int degrees_continuity, 00082 Real ray_origin[], 00083 Real unit_pos_dir[], 00084 Real unit_neg_dir[], 00085 Real model_dist, 00086 Real start_dist, 00087 Real end_dist, 00088 Real origin[], 00089 Real direction[], 00090 Real *current_distance, 00091 Real *stop_distance, 00092 Real next_distance[N_DIMENSIONS], 00093 Real delta_distance[N_DIMENSIONS], 00094 int voxel_index[N_DIMENSIONS], 00095 int delta_voxel[N_DIMENSIONS], 00096 Real *next_closest_distance ) 00097 { 00098 BOOLEAN found_exit; 00099 Real voxel_exit, t_min, t_max, voxel_offset; 00100 Real direct; 00101 Real start_voxel[N_DIMENSIONS], model_point[N_DIMENSIONS]; 00102 Real *dir; 00103 00104 if( model_dist + start_dist < 0.0 || model_dist + end_dist < 0.0 ) 00105 dir = &unit_neg_dir[0]; 00106 else 00107 dir = &unit_pos_dir[0]; 00108 00109 GET_RAY_POINT( model_point, ray_origin, dir, model_dist ); 00110 00111 convert_world_to_voxel( volume, 00112 model_point[X], model_point[Y], model_point[Z], 00113 origin ); 00114 00115 if( end_dist >= start_dist ) 00116 { 00117 GET_RAY_POINT( direction, model_point, dir, 1.0 ); 00118 } 00119 else 00120 { 00121 GET_RAY_POINT( direction, model_point, dir, -1.0 ); 00122 } 00123 00124 convert_world_to_voxel( volume, 00125 direction[X], 00126 direction[Y], 00127 direction[Z], direction ); 00128 direction[X] -= origin[X]; 00129 direction[Y] -= origin[Y]; 00130 direction[Z] -= origin[Z]; 00131 00132 clip_line_to_volume( volume, degrees_continuity, 00133 origin, direction, &t_min, &t_max ); 00134 00135 *current_distance = FABS(start_dist); 00136 if( t_min > *current_distance ) 00137 *current_distance = t_min; 00138 *stop_distance = FABS(end_dist); 00139 if( t_max < *stop_distance ) 00140 *stop_distance = t_max; 00141 00142 GET_RAY_POINT( start_voxel, origin, direction, *current_distance ); 00143 00144 found_exit = FALSE; 00145 00146 if( degrees_continuity % 2 == 1 ) 00147 voxel_offset = 0.5; 00148 else 00149 voxel_offset = 0.0; 00150 00151 /*--- do x loop */ 00152 00153 direct = direction[X]; 00154 00155 voxel_index[X] = (int) (start_voxel[X]+voxel_offset); 00156 voxel_exit = (Real) voxel_index[X] - voxel_offset; 00157 00158 if( direct == 0.0 ) 00159 { 00160 delta_voxel[X] = 1; 00161 delta_distance[X] = 0.0; 00162 next_distance[X] = *stop_distance; 00163 } 00164 else 00165 { 00166 direct = 1.0 / direct; 00167 if( direct < 0.0 ) 00168 { 00169 delta_voxel[X] = -1; 00170 delta_distance[X] = -direct; 00171 } 00172 else 00173 { 00174 voxel_exit += 1.0; 00175 delta_voxel[X] = 1; 00176 delta_distance[X] = direct; 00177 } 00178 00179 next_distance[X] = (voxel_exit - origin[X]) * direct; 00180 00181 *next_closest_distance = next_distance[X]; 00182 found_exit = TRUE; 00183 } 00184 /*--- do y loop */ 00185 00186 direct = direction[Y]; 00187 00188 voxel_index[Y] = (int) (start_voxel[Y]+voxel_offset); 00189 voxel_exit = (Real) voxel_index[Y] - voxel_offset; 00190 00191 if( direct == 0.0 ) 00192 { 00193 delta_voxel[Y] = 1; 00194 delta_distance[Y] = 0.0; 00195 next_distance[Y] = *stop_distance; 00196 } 00197 else 00198 { 00199 direct = 1.0 / direct; 00200 if( direct < 0.0 ) 00201 { 00202 delta_voxel[Y] = -1; 00203 delta_distance[Y] = -direct; 00204 } 00205 else 00206 { 00207 voxel_exit += 1.0; 00208 delta_voxel[Y] = 1; 00209 delta_distance[Y] = direct; 00210 } 00211 00212 next_distance[Y] = (voxel_exit - origin[Y]) * direct; 00213 00214 if( !found_exit || next_distance[Y] < *next_closest_distance ) 00215 { 00216 *next_closest_distance = next_distance[Y]; 00217 found_exit = TRUE; 00218 } 00219 } 00220 /*--- do z loop */ 00221 00222 direct = direction[Z]; 00223 00224 voxel_index[Z] = (int) (start_voxel[Z]+voxel_offset); 00225 voxel_exit = (Real) voxel_index[Z] - voxel_offset; 00226 00227 if( direct == 0.0 ) 00228 { 00229 delta_voxel[Z] = 1; 00230 delta_distance[Z] = 0.0; 00231 next_distance[Z] = *stop_distance; 00232 } 00233 else 00234 { 00235 direct = 1.0 / direct; 00236 if( direct < 0.0 ) 00237 { 00238 delta_voxel[Z] = -1; 00239 delta_distance[Z] = -direct; 00240 } 00241 else 00242 { 00243 voxel_exit += 1.0; 00244 delta_voxel[Z] = 1; 00245 delta_distance[Z] = direct; 00246 } 00247 00248 next_distance[Z] = (voxel_exit - origin[Z]) * direct; 00249 00250 if( !found_exit || next_distance[Z] < *next_closest_distance ) 00251 { 00252 *next_closest_distance = next_distance[Z]; 00253 found_exit = TRUE; 00254 } 00255 } 00256 } 00257 00258 #ifdef DEBUGGING 00259 static int count = 0; 00260 #endif 00261 00262 public BOOLEAN find_boundary_in_direction( 00263 Volume volume, 00264 Volume label_volume, 00265 voxel_coef_struct *lookup, 00266 bitlist_3d_struct *done_bits, 00267 bitlist_3d_struct *surface_bits, 00268 Real model_dist, 00269 Point *ray_origin, 00270 Vector *unit_pos_dir, 00271 Vector *unit_neg_dir, 00272 Real max_outwards_search_distance, 00273 Real max_inwards_search_distance, 00274 int degrees_continuity, 00275 boundary_definition_struct *boundary_def, 00276 Real *boundary_distance ) 00277 { 00278 private const Real TOLERANCE = 0.0001; 00279 BOOLEAN found, done0, done1, second_leg0, second_leg1, isovalue; 00280 Real next_closest0, next_closest1; 00281 int voxel_index0[N_DIMENSIONS], voxel_index1[N_DIMENSIONS]; 00282 int delta_voxel0[N_DIMENSIONS], delta_voxel1[N_DIMENSIONS]; 00283 Real current_distance0, current_distance1; 00284 Real next_distance0[N_DIMENSIONS], next_distance1[N_DIMENSIONS]; 00285 Real stop_distance0, stop_distance1; 00286 Real delta_distance0[N_DIMENSIONS], delta_distance1[N_DIMENSIONS]; 00287 Real max_dist, found_dist; 00288 Real ray_origin_real[N_DIMENSIONS]; 00289 Real unit_pos[N_DIMENSIONS], unit_neg[N_DIMENSIONS]; 00290 BOOLEAN inside_surface0, inside_surface1, do_first; 00291 Real origin0[N_DIMENSIONS], origin1[N_DIMENSIONS]; 00292 Real direction0[N_DIMENSIONS], direction1[N_DIMENSIONS]; 00293 00294 #ifdef DEBUGGING 00295 { 00296 static int stop_at = -2; 00297 00298 if( stop_at == -2 ) 00299 { 00300 if( getenv("STOP") == NULL || 00301 sscanf( getenv("STOP"), "%d", &stop_at ) != 1) 00302 stop_at = -1; 00303 } 00304 00305 ++count; 00306 00307 if( count == stop_at ) 00308 { 00309 (void) printf( "Stop: %d\n", stop_at ); 00310 (void) getchar(); 00311 } 00312 } 00313 #endif 00314 00315 ray_origin_real[X] = (Real) Point_x(*ray_origin); 00316 unit_pos[X] = (Real) Vector_x(*unit_pos_dir); 00317 unit_neg[X] = (Real) Vector_x(*unit_neg_dir); 00318 00319 ray_origin_real[Y] = (Real) Point_y(*ray_origin); 00320 unit_pos[Y] = (Real) Vector_y(*unit_pos_dir); 00321 unit_neg[Y] = (Real) Vector_y(*unit_neg_dir); 00322 00323 ray_origin_real[Z] = (Real) Point_z(*ray_origin); 00324 unit_pos[Z] = (Real) Vector_z(*unit_pos_dir); 00325 unit_neg[Z] = (Real) Vector_z(*unit_neg_dir); 00326 00327 second_leg0 = FALSE; 00328 second_leg1 = FALSE; 00329 00330 if( model_dist >= 0.0 ) 00331 { 00332 set_up_to_search_ray( volume, degrees_continuity, 00333 ray_origin_real, unit_pos, unit_neg, 00334 model_dist, 0.0, max_outwards_search_distance, 00335 origin0, direction0, &current_distance0, 00336 &stop_distance0, next_distance0, 00337 delta_distance0, voxel_index0, delta_voxel0, 00338 &next_closest0 ); 00339 } 00340 else 00341 { 00342 set_up_to_search_ray( volume, degrees_continuity, 00343 ray_origin_real, unit_pos, unit_neg, 00344 model_dist, 00345 0.0, 00346 MIN( max_outwards_search_distance, -model_dist), 00347 origin0, direction0, &current_distance0, 00348 &stop_distance0, next_distance0, 00349 delta_distance0, voxel_index0, delta_voxel0, 00350 &next_closest0 ); 00351 if( -model_dist < max_outwards_search_distance ) 00352 second_leg0 = TRUE; 00353 } 00354 00355 if( model_dist <= 0.0 ) 00356 { 00357 set_up_to_search_ray( volume, degrees_continuity, 00358 ray_origin_real, unit_pos, unit_neg, 00359 model_dist, 0.0, -max_inwards_search_distance, 00360 origin1, direction1, &current_distance1, 00361 &stop_distance1, next_distance1, 00362 delta_distance1, voxel_index1, delta_voxel1, 00363 &next_closest1 ); 00364 } 00365 else 00366 { 00367 set_up_to_search_ray( volume, degrees_continuity, 00368 ray_origin_real, unit_pos, unit_neg, 00369 model_dist, 00370 0.0, 00371 - MIN( max_inwards_search_distance, model_dist ), 00372 origin1, direction1, &current_distance1, 00373 &stop_distance1, next_distance1, 00374 delta_distance1, voxel_index1, delta_voxel1, 00375 &next_closest1 ); 00376 if( model_dist < max_inwards_search_distance ) 00377 second_leg1 = TRUE; 00378 } 00379 00380 if( boundary_def->min_isovalue != boundary_def->max_isovalue ) 00381 { 00382 Vector world_direction; 00383 Real point[N_DIMENSIONS]; 00384 Real vx, vy, vz; 00385 00386 isovalue = FALSE; 00387 00388 convert_voxel_normal_vector_to_world( volume, direction0, 00389 &vx, &vy, &vz ); 00390 fill_Vector( world_direction, vx, vy, vz ); 00391 NORMALIZE_VECTOR( world_direction, world_direction ); 00392 GET_RAY_POINT( point, origin0, direction0, current_distance0 ); 00393 inside_surface0 = is_point_inside_surface( volume, label_volume, 00394 degrees_continuity, 00395 point, &world_direction, 00396 boundary_def ); 00397 00398 convert_voxel_normal_vector_to_world( volume, direction1, 00399 &vx, &vy, &vz ); 00400 fill_Vector( world_direction, vx, vy, vz ); 00401 NORMALIZE_VECTOR( world_direction, world_direction ); 00402 SCALE_VECTOR( world_direction, world_direction, -1.0 ); 00403 00404 GET_RAY_POINT( point, origin1, direction1, current_distance1 ); 00405 00406 inside_surface1 = is_point_inside_surface( volume, label_volume, 00407 degrees_continuity, 00408 point, &world_direction, 00409 boundary_def ); 00410 } 00411 else 00412 isovalue = TRUE; 00413 00414 found = FALSE; 00415 done0 = current_distance0 >= stop_distance0 - TOLERANCE; 00416 done1 = current_distance1 >= stop_distance1 - TOLERANCE; 00417 00418 do_first = ( done1 || (!done0 && current_distance0 <= current_distance1)); 00419 00420 while( !done0 || !done1 ) 00421 { 00422 if( do_first ) 00423 { 00424 max_dist = next_closest0; 00425 if( stop_distance0 < max_dist ) 00426 max_dist = stop_distance0; 00427 00428 if( voxel_might_contain_boundary( volume, done_bits, surface_bits, 00429 degrees_continuity, voxel_index0, 00430 boundary_def ) && 00431 (isovalue && 00432 check_voxel_for_isovalue( lookup, volume, label_volume, 00433 degrees_continuity, 00434 voxel_index0, 00435 origin0, direction0, 00436 FALSE, 00437 current_distance0, max_dist, 00438 boundary_def->min_isovalue, 00439 boundary_def->normal_direction, 00440 found, &found_dist ) || 00441 !isovalue && 00442 check_voxel_for_boundary( lookup, volume, label_volume, 00443 degrees_continuity, 00444 voxel_index0, 00445 origin0, direction0, 00446 &inside_surface0, 00447 FALSE, 00448 current_distance0, max_dist, 00449 boundary_def, found, &found_dist ) ) ) 00450 { 00451 found = TRUE; 00452 *boundary_distance = found_dist; 00453 } 00454 00455 if( next_closest0 >= stop_distance0 - TOLERANCE ) 00456 { 00457 if( second_leg0 ) 00458 { 00459 second_leg0 = FALSE; 00460 set_up_to_search_ray( 00461 volume, degrees_continuity, ray_origin_real, 00462 unit_pos, unit_neg, model_dist, 00463 -model_dist, max_outwards_search_distance, 00464 origin0, direction0, &current_distance0, 00465 &stop_distance0, next_distance0, 00466 delta_distance0, voxel_index0, delta_voxel0, 00467 &next_closest0 ); 00468 done0 = current_distance0 >= stop_distance0 - TOLERANCE; 00469 } 00470 else 00471 done0 = TRUE; 00472 00473 do_first = (done1 || 00474 (!done0 && current_distance0 <= current_distance1)); 00475 } 00476 else 00477 { 00478 current_distance0 = next_closest0; 00479 do_first = (done1 || current_distance0 <= current_distance1); 00480 00481 if( next_distance0[X] <= current_distance0 ) 00482 { 00483 voxel_index0[X] += delta_voxel0[X]; 00484 next_distance0[X] += delta_distance0[X]; 00485 } 00486 if( next_distance0[Y] <= current_distance0 ) 00487 { 00488 voxel_index0[Y] += delta_voxel0[Y]; 00489 next_distance0[Y] += delta_distance0[Y]; 00490 } 00491 if( next_distance0[Z] <= current_distance0 ) 00492 { 00493 voxel_index0[Z] += delta_voxel0[Z]; 00494 next_distance0[Z] += delta_distance0[Z]; 00495 } 00496 00497 next_closest0 = next_distance0[X]; 00498 if( next_distance0[Y] < next_closest0 ) 00499 next_closest0 = next_distance0[Y]; 00500 if( next_distance0[Z] < next_closest0 ) 00501 next_closest0 = next_distance0[Z]; 00502 } 00503 } 00504 else 00505 { 00506 max_dist = next_closest1; 00507 if( stop_distance1 < max_dist ) 00508 max_dist = stop_distance1; 00509 00510 if( voxel_might_contain_boundary( volume, done_bits, surface_bits, 00511 degrees_continuity, voxel_index1, 00512 boundary_def ) && 00513 (isovalue && 00514 check_voxel_for_isovalue( lookup, volume, label_volume, 00515 degrees_continuity, 00516 voxel_index1, 00517 origin1, direction1, 00518 TRUE, 00519 current_distance1, max_dist, 00520 boundary_def->min_isovalue, 00521 boundary_def->normal_direction, 00522 found, &found_dist ) || 00523 !isovalue && 00524 check_voxel_for_boundary( lookup, volume, label_volume, 00525 degrees_continuity, 00526 voxel_index1, 00527 origin1, direction1, 00528 &inside_surface1, 00529 TRUE, 00530 current_distance1, max_dist, 00531 boundary_def, found, &found_dist ) ) ) 00532 { 00533 found = TRUE; 00534 *boundary_distance = found_dist; 00535 } 00536 00537 if( next_closest1 >= stop_distance1 - TOLERANCE ) 00538 { 00539 if( second_leg1 ) 00540 { 00541 second_leg1 = FALSE; 00542 set_up_to_search_ray( 00543 volume, degrees_continuity, ray_origin_real, 00544 unit_pos, unit_neg, model_dist, 00545 -model_dist, -max_inwards_search_distance, 00546 origin1, direction1, &current_distance1, 00547 &stop_distance1, next_distance1, 00548 delta_distance1, voxel_index1, delta_voxel1, 00549 &next_closest1 ); 00550 done1 = current_distance1 >= stop_distance1 - TOLERANCE; 00551 } 00552 else 00553 done1 = TRUE; 00554 00555 do_first = !(done0 || 00556 (!done1 && current_distance1 <= current_distance0)); 00557 } 00558 else 00559 { 00560 current_distance1 = next_closest1; 00561 do_first = !(done0 || current_distance1 <= current_distance0); 00562 00563 if( next_distance1[X] <= current_distance1 ) 00564 { 00565 voxel_index1[X] += delta_voxel1[X]; 00566 next_distance1[X] += delta_distance1[X]; 00567 } 00568 if( next_distance1[Y] <= current_distance1 ) 00569 { 00570 voxel_index1[Y] += delta_voxel1[Y]; 00571 next_distance1[Y] += delta_distance1[Y]; 00572 } 00573 if( next_distance1[Z] <= current_distance1 ) 00574 { 00575 voxel_index1[Z] += delta_voxel1[Z]; 00576 next_distance1[Z] += delta_distance1[Z]; 00577 } 00578 00579 next_closest1 = next_distance1[X]; 00580 if( next_distance1[Y] < next_closest1 ) 00581 next_closest1 = next_distance1[Y]; 00582 if( next_distance1[Z] < next_closest1 ) 00583 next_closest1 = next_distance1[Z]; 00584 } 00585 } 00586 00587 if( found && 00588 (done0 || FABS(*boundary_distance) <= current_distance0) && 00589 (done1 || FABS(*boundary_distance) <= current_distance1) ) 00590 { 00591 break; 00592 } 00593 } 00594 00595 return( found ); 00596 } 00597 00598 private void get_trilinear_gradient( 00599 Real coefs[], 00600 Real u, 00601 Real v, 00602 Real w, 00603 Real derivs[] ) 00604 { 00605 Real du00, du01, du10, du11, c00, c01, c10, c11; 00606 Real dv0, dv1, c0, c1, dw, du0, du1; 00607 00608 /*--- get the 4 differences in the u direction */ 00609 00610 du00 = coefs[4] - coefs[0]; 00611 du01 = coefs[5] - coefs[1]; 00612 du10 = coefs[6] - coefs[2]; 00613 du11 = coefs[7] - coefs[3]; 00614 00615 /*--- reduce to a 2D problem, by interpolating in the u direction */ 00616 00617 c00 = coefs[0] + u * du00; 00618 c01 = coefs[1] + u * du01; 00619 c10 = coefs[2] + u * du10; 00620 c11 = coefs[3] + u * du11; 00621 00622 /*--- get the 2 differences in the v direction for the 2D problem */ 00623 00624 dv0 = c10 - c00; 00625 dv1 = c11 - c01; 00626 /*--- reduce 2D to a 1D problem, by interpolating in the v direction */ 00627 00628 c0 = c00 + v * dv0; 00629 c1 = c01 + v * dv1; 00630 00631 /*--- get the 1 difference in the w direction for the 1D problem */ 00632 00633 dw = c1 - c0; 00634 00635 /*--- reduce the 2D u derivs to 1D */ 00636 00637 du0 = INTERPOLATE( v, du00, du10 ); 00638 du1 = INTERPOLATE( v, du01, du11 ); 00639 00640 /*--- interpolate the 1D problems in w, or for Z deriv, just use dw */ 00641 00642 derivs[X] = INTERPOLATE( w, du0, du1 ); 00643 derivs[Y] = INTERPOLATE( w, dv0, dv1 ); 00644 derivs[Z] = dw; 00645 } 00646 00647 #define MAX_DERIVS 2 00648 00649 private BOOLEAN check_voxel_for_isovalue( 00650 voxel_coef_struct *lookup, 00651 Volume volume, 00652 Volume label_volume, 00653 int degrees_continuity, 00654 int voxel_indices[], 00655 Real line_origin[], 00656 Real line_direction[], 00657 BOOLEAN searching_inwards, 00658 Real min_line_t, 00659 Real max_line_t, 00660 Real isovalue, 00661 Normal_directions normal_direction, 00662 BOOLEAN already_found, 00663 Real *dist_found ) 00664 { 00665 BOOLEAN found; 00666 int i; 00667 Real value, dot_prod, voxel[N_DIMENSIONS]; 00668 Real first_deriv[N_DIMENSIONS], *first_deriv_ptr[1]; 00669 BOOLEAN active, deriv_dir_correct; 00670 int n_boundaries; 00671 Real boundary_positions[3]; 00672 Real coefs[(2+MAX_DERIVS)*(2+MAX_DERIVS)*(2+MAX_DERIVS)]; 00673 00674 found = FALSE; 00675 00676 lookup_volume_coeficients( lookup, volume, degrees_continuity, 00677 voxel_indices[X], 00678 voxel_indices[Y], 00679 voxel_indices[Z], coefs ); 00680 00681 n_boundaries = find_voxel_line_value_intersection( coefs, 00682 degrees_continuity, 00683 voxel_indices[X], 00684 voxel_indices[Y], 00685 voxel_indices[Z], 00686 line_origin, 00687 line_direction, 00688 min_line_t, 00689 max_line_t, 00690 isovalue, 00691 boundary_positions ); 00692 00693 for_less( i, 0, n_boundaries ) 00694 { 00695 if( (!already_found || boundary_positions[i]< FABS(*dist_found)) ) 00696 { 00697 GET_RAY_POINT( voxel, line_origin, line_direction, 00698 boundary_positions[i] ); 00699 active = get_volume_voxel_activity( label_volume, voxel, FALSE ); 00700 if( active ) 00701 { 00702 if( degrees_continuity == 0 ) 00703 { 00704 get_trilinear_gradient( coefs, 00705 voxel[X] - (Real) voxel_indices[X], 00706 voxel[Y] - (Real) voxel_indices[Y], 00707 voxel[Z] - (Real) voxel_indices[Z], 00708 first_deriv ); 00709 } 00710 else 00711 { 00712 first_deriv_ptr[0] = first_deriv; 00713 (void) evaluate_volume( volume, voxel, NULL, 00714 degrees_continuity, FALSE, 00715 get_volume_real_min(volume), &value, 00716 first_deriv_ptr, NULL ); 00717 } 00718 00719 deriv_dir_correct = TRUE; 00720 if( normal_direction != ANY_DIRECTION ) 00721 { 00722 dot_prod = first_deriv[X] * line_direction[X] + 00723 first_deriv[Y] * line_direction[Y] + 00724 first_deriv[Z] * line_direction[Z]; 00725 if( searching_inwards ) 00726 dot_prod = -dot_prod; 00727 if( normal_direction == TOWARDS_LOWER && dot_prod > 0.0 || 00728 normal_direction == TOWARDS_HIGHER && dot_prod < 0.0 ) 00729 { 00730 deriv_dir_correct = FALSE; 00731 } 00732 } 00733 00734 if( deriv_dir_correct ) 00735 { 00736 if( searching_inwards ) 00737 *dist_found = -boundary_positions[i]; 00738 else 00739 *dist_found = boundary_positions[i]; 00740 found = TRUE; 00741 } 00742 } 00743 } 00744 } 00745 00746 return( found ); 00747 } 00748 00749 private BOOLEAN does_voxel_contain_value_range( 00750 Volume volume, 00751 int degrees_continuity, 00752 int voxel[], 00753 Real min_value, 00754 Real max_value ) 00755 { 00756 int dim, i, n_values, start, end, sizes[MAX_DIMENSIONS]; 00757 BOOLEAN greater, less; 00758 Real values[4*4*4]; 00759 00760 get_volume_sizes( volume, sizes ); 00761 00762 /*--- special case, to be done fast */ 00763 00764 if( degrees_continuity == 0 ) 00765 { 00766 if( voxel[0] < 0 || voxel[0] >= sizes[0]-1 || 00767 voxel[1] < 0 || voxel[1] >= sizes[1]-1 || 00768 voxel[2] < 0 || voxel[2] >= sizes[2]-1 ) 00769 return( FALSE ); 00770 00771 get_volume_value_hyperslab_3d( volume, voxel[X], voxel[Y], voxel[Z], 00772 2, 2, 2, values ); 00773 00774 if( min_value <= values[0] && values[0] <= max_value ) 00775 return( TRUE ); 00776 00777 if( values[0] < min_value ) 00778 { 00779 return( values[1] >= min_value || 00780 values[2] >= min_value || 00781 values[3] >= min_value || 00782 values[4] >= min_value || 00783 values[5] >= min_value || 00784 values[6] >= min_value || 00785 values[7] >= min_value ); 00786 } 00787 else 00788 { 00789 return( values[1] <= max_value || 00790 values[2] <= max_value || 00791 values[3] <= max_value || 00792 values[4] <= max_value || 00793 values[5] <= max_value || 00794 values[6] <= max_value || 00795 values[7] <= max_value ); 00796 } 00797 } 00798 00799 start = -(degrees_continuity + 1) / 2; 00800 end = start + degrees_continuity + 2; 00801 00802 for_less( dim, 0, N_DIMENSIONS ) 00803 { 00804 if( voxel[dim] + start < 0 || voxel[dim] + end > sizes[dim] ) 00805 return( FALSE ); 00806 } 00807 00808 get_volume_value_hyperslab_3d( volume, voxel[X], voxel[Y], voxel[Z], 00809 degrees_continuity + 2, 00810 degrees_continuity + 2, 00811 degrees_continuity + 2, values ); 00812 00813 n_values = (degrees_continuity + 2) * 00814 (degrees_continuity + 2) * 00815 (degrees_continuity + 2); 00816 00817 less = FALSE; 00818 greater = FALSE; 00819 00820 for_less( i, 0, n_values ) 00821 { 00822 if( values[i] < min_value ) 00823 { 00824 if( greater ) 00825 return( TRUE ); 00826 less = TRUE; 00827 } 00828 else if( values[i] > max_value ) 00829 { 00830 if( less ) 00831 return( TRUE ); 00832 greater = TRUE; 00833 } 00834 else 00835 return( TRUE ); 00836 } 00837 00838 return( FALSE ); 00839 } 00840 00841 private BOOLEAN voxel_might_contain_boundary( 00842 Volume volume, 00843 bitlist_3d_struct *done_bits, 00844 bitlist_3d_struct *surface_bits, 00845 int degrees_continuity, 00846 int voxel_indices[], 00847 boundary_definition_struct *boundary_def ) 00848 { 00849 BOOLEAN contains; 00850 00851 if( done_bits != NULL ) 00852 { 00853 if( !get_bitlist_bit_3d( done_bits, 00854 voxel_indices[X], 00855 voxel_indices[Y], 00856 voxel_indices[Z] ) ) 00857 { 00858 contains = does_voxel_contain_value_range( 00859 volume, degrees_continuity, 00860 voxel_indices, 00861 boundary_def->min_isovalue, 00862 boundary_def->max_isovalue ); 00863 00864 set_bitlist_bit_3d( surface_bits, 00865 voxel_indices[X], 00866 voxel_indices[Y], 00867 voxel_indices[Z], contains ); 00868 set_bitlist_bit_3d( done_bits, 00869 voxel_indices[X], 00870 voxel_indices[Y], 00871 voxel_indices[Z], TRUE ); 00872 } 00873 else 00874 { 00875 contains = get_bitlist_bit_3d( surface_bits, 00876 voxel_indices[X], 00877 voxel_indices[Y], 00878 voxel_indices[Z] ); 00879 } 00880 00881 if( !contains ) 00882 return( FALSE ); 00883 } 00884 00885 return( TRUE ); 00886 } 00887 00888 private BOOLEAN check_voxel_for_boundary( 00889 voxel_coef_struct *lookup, 00890 Volume volume, 00891 Volume label_volume, 00892 int degrees_continuity, 00893 int voxel_indices[], 00894 Real line_origin[], 00895 Real line_direction[], 00896 BOOLEAN *inside_surface, 00897 BOOLEAN searching_inwards, 00898 Real min_line_t, 00899 Real max_line_t, 00900 boundary_definition_struct *boundary_def, 00901 BOOLEAN already_found, 00902 Real *dist_found ) 00903 { 00904 BOOLEAN found, inside, prev_inside; 00905 int degrees; 00906 Real line_poly[10], value, vx, vy, vz; 00907 Real t, increment, t_min, t_max; 00908 Vector direction; 00909 Real surface[N_DIMENSIONS]; 00910 Real coefs[(2+MAX_DERIVS)*(2+MAX_DERIVS)*(2+MAX_DERIVS)]; 00911 00912 found = FALSE; 00913 00914 lookup_volume_coeficients( lookup, volume, degrees_continuity, 00915 voxel_indices[X], 00916 voxel_indices[Y], 00917 voxel_indices[Z], coefs ); 00918 00919 degrees = find_voxel_line_polynomial( coefs, degrees_continuity, 00920 voxel_indices[X], 00921 voxel_indices[Y], 00922 voxel_indices[Z], 00923 line_origin, 00924 line_direction, 00925 line_poly ); 00926 00927 if( degrees == 0 ) 00928 { 00929 if( *inside_surface ) 00930 { 00931 if( !already_found || min_line_t < FABS(*dist_found) ) 00932 { 00933 if( searching_inwards ) 00934 *dist_found = -min_line_t; 00935 else 00936 *dist_found = min_line_t; 00937 found = TRUE; 00938 } 00939 } 00940 *inside_surface = FALSE; 00941 return( found ); 00942 } 00943 00944 convert_voxel_normal_vector_to_world( volume, line_direction, 00945 &vx, &vy, &vz ); 00946 fill_Vector( direction, vx, vy, vz ); 00947 NORMALIZE_VECTOR( direction, direction ); 00948 if( searching_inwards ) 00949 SCALE_VECTOR( direction, direction, -1.0 ); 00950 00951 if( !get_range_of_polynomial( degrees, line_poly, min_line_t, max_line_t, 00952 boundary_def->min_isovalue, 00953 boundary_def->max_isovalue, 00954 boundary_def->tolerance, &t_min, &t_max ) ) 00955 { 00956 00957 GET_RAY_POINT( surface, line_origin, line_direction, 00958 (min_line_t + max_line_t) / 2.0 ); 00959 00960 inside = is_point_inside_surface( volume, label_volume, 00961 degrees_continuity, 00962 surface, &direction, 00963 boundary_def ); 00964 00965 if( inside != *inside_surface ) 00966 { 00967 if( !already_found || min_line_t < FABS(*dist_found) ) 00968 { 00969 if( searching_inwards ) 00970 *dist_found = -min_line_t; 00971 else 00972 *dist_found = min_line_t; 00973 found = TRUE; 00974 } 00975 } 00976 00977 *inside_surface = inside; 00978 00979 return( found ); 00980 } 00981 00982 found = FALSE; 00983 00984 increment = boundary_def->tolerance; 00985 00986 if( t_min > min_line_t ) 00987 inside = FALSE; 00988 else 00989 inside = *inside_surface; 00990 00991 for( t = t_min; t <= t_max; t += increment ) 00992 { 00993 prev_inside = inside; 00994 00995 GET_RAY_POINT( surface, line_origin, line_direction, t ); 00996 00997 value = evaluate_polynomial( degrees, line_poly, t ); 00998 00999 if( value < boundary_def->min_isovalue ) 01000 inside = FALSE; 01001 else if( value > boundary_def->max_isovalue ) 01002 inside = TRUE; 01003 else 01004 inside = is_point_inside_surface( volume, label_volume, 01005 degrees_continuity, 01006 surface, &direction, 01007 boundary_def ); 01008 01009 if( inside != prev_inside ) 01010 { 01011 if( !already_found || t < FABS(*dist_found) ) 01012 { 01013 if( searching_inwards ) 01014 *dist_found = -t; 01015 else 01016 *dist_found = t; 01017 found = TRUE; 01018 break; 01019 } 01020 } 01021 } 01022 01023 if( t_max != max_line_t ) 01024 *inside_surface = FALSE; 01025 else 01026 *inside_surface = inside; 01027 01028 return( found ); 01029 }