BICPL: smooth_polygons.c Source File

00001 /* ---------------------------------------------------------------------------- 00002 @COPYRIGHT : 00003 Copyright 1993,1994,1995 David MacDonald, 00004 McConnell Brain Imaging Centre, 00005 Montreal Neurological Institute, McGill University. 00006 Permission to use, copy, modify, and distribute this 00007 software and its documentation for any purpose and without 00008 fee is hereby granted, provided that the above copyright 00009 notice appear in all copies. The author and McGill University 00010 make no representations about the suitability of this 00011 software for any purpose. It is provided "as is" without 00012 express or implied warranty. 00013 ---------------------------------------------------------------------------- */ 00014 00015 #include <volume_io/internal_volume_io.h> 00016 #include <bicpl/geom.h> 00017 00018 #ifndef lint 00019 static char rcsid[] = "$Header: /software/source//libraries/bicpl/Geometry/smooth_polygons.c,v 1.12 2000/02/06 15:30:18 stever Exp $"; 00020 #endif 00021 00022 #define CHECK_INTERVAL 1.0 00023 00024 private void smooth_points( 00025 polygons_struct *polygons, 00026 Point current_points[], 00027 Point new_points[], 00028 Smallest_int point_done[], 00029 Real max_dist_from_original, 00030 Real fraction_to_move, 00031 Real normal_ratio, 00032 BOOLEAN range_flag, 00033 volume_struct *volume, 00034 int min_value, 00035 int max_value, 00036 Real *avg_moved, 00037 Real *max_moved ); 00038 private Real update_point_position( 00039 polygons_struct *polygons, 00040 int poly, 00041 int vertex_index, 00042 int point_index, 00043 Point current_points[], 00044 Point *new_point, 00045 Real max_dist_from_original, 00046 Real fraction_to_move, 00047 Real normal_ratio, 00048 BOOLEAN range_flag, 00049 volume_struct *volume, 00050 int min_value, 00051 int max_value ); 00052 private BOOLEAN point_inside_range( 00053 volume_struct *volume, 00054 Point *point, 00055 int min_value, 00056 int max_value ); 00057 00058 /* ----------------------------- MNI Header ----------------------------------- 00059 @NAME : smooth_polygon 00060 @INPUT : polygons 00061 max_dist_from_original 00062 fraction_to_move 00063 stop_threshold 00064 normal_ratio 00065 range_flag 00066 volume 00067 min_value 00068 max_value 00069 @OUTPUT : 00070 @RETURNS : 00071 @DESCRIPTION: Smooths the polygons by moving vertices towards the centroid 00072 of their neighbours. 00073 @METHOD : 00074 @GLOBALS : 00075 @CALLS : 00076 @CREATED : 1993 David MacDonald 00077 @MODIFIED : 00078 ---------------------------------------------------------------------------- */ 00079 00080 public void smooth_polygon( 00081 polygons_struct *polygons, 00082 Real max_dist_from_original, 00083 Real fraction_to_move, 00084 Real stop_threshold, 00085 Real normal_ratio, 00086 BOOLEAN range_flag, 00087 volume_struct *volume, 00088 int min_value, 00089 int max_value ) 00090 { 00091 Real avg_moved, max_moved; 00092 int i, iteration; 00093 Point *new_points, *tmp, *current_points; 00094 Smallest_int *point_done; 00095 Real next_check_time; 00096 00097 if( polygons->n_points <= 0 ) 00098 return; 00099 00100 ALLOC( new_points, polygons->n_points ); 00101 ALLOC( current_points, polygons->n_points ); 00102 ALLOC( point_done, polygons->n_points ); 00103 00104 check_polygons_neighbours_computed( polygons ); 00105 00106 for_less( i, 0, polygons->n_points ) 00107 { 00108 current_points[i] = polygons->points[i]; 00109 } 00110 00111 iteration = 0; 00112 next_check_time = current_realtime_seconds() + CHECK_INTERVAL; 00113 do 00114 { 00115 smooth_points( polygons, current_points, new_points, point_done, 00116 max_dist_from_original, fraction_to_move, 00117 normal_ratio, range_flag, 00118 volume, min_value, max_value, 00119 &avg_moved, &max_moved ); 00120 00121 tmp = current_points; 00122 current_points = new_points; 00123 new_points = tmp; 00124 00125 ++iteration; 00126 00127 print( "Iteration %d -- avg distance %g max distance %g\n", 00128 iteration, avg_moved, max_moved ); 00129 if( current_realtime_seconds() > next_check_time ) 00130 { 00131 next_check_time = current_realtime_seconds() + CHECK_INTERVAL; 00132 00133 if( file_exists("interrupt") ) 00134 { 00135 print( "Interrupting as requested\n" ); 00136 remove_file( "interrupt" ); 00137 break; 00138 } 00139 } 00140 } 00141 while( max_moved > stop_threshold ); 00142 00143 for_less( i, 0, polygons->n_points ) 00144 polygons->points[i] = current_points[i]; 00145 00146 FREE( new_points ); 00147 FREE( current_points ); 00148 FREE( point_done ); 00149 } 00150 00151 /* ----------------------------- MNI Header ----------------------------------- 00152 @NAME : smooth_points 00153 @INPUT : polygons 00154 current_points[] 00155 new_points 00156 point_done 00157 max_dist_from_original 00158 fraction_to_move 00159 normal_ratio 00160 range_flag 00161 volume 00162 min_value 00163 max_value 00164 @OUTPUT : avg_moved 00165 max_moved 00166 @RETURNS : 00167 @DESCRIPTION: 00168 @METHOD : 00169 @GLOBALS : 00170 @CALLS : 00171 @CREATED : 1993 David MacDonald 00172 @MODIFIED : 00173 ---------------------------------------------------------------------------- */ 00174 00175 private void smooth_points( 00176 polygons_struct *polygons, 00177 Point current_points[], 00178 Point new_points[], 00179 Smallest_int point_done[], 00180 Real max_dist_from_original, 00181 Real fraction_to_move, 00182 Real normal_ratio, 00183 BOOLEAN range_flag, 00184 volume_struct *volume, 00185 int min_value, 00186 int max_value, 00187 Real *avg_moved, 00188 Real *max_moved ) 00189 { 00190 int p, vertex_index, point_index, poly, size; 00191 progress_struct progress; 00192 Real moved; 00193 00194 for_less( p, 0, polygons->n_points ) 00195 { 00196 new_points[p] = current_points[p]; 00197 point_done[p] = FALSE; 00198 } 00199 00200 *avg_moved = 0.0; 00201 *max_moved = 0.0; 00202 00203 initialize_progress_report( &progress, TRUE, polygons->n_items, 00204 "Averaging Points" ); 00205 00206 for_less( poly, 0, polygons->n_items ) 00207 { 00208 size = GET_OBJECT_SIZE( *polygons, poly ); 00209 00210 for_less( vertex_index, 0, size ) 00211 { 00212 point_index = polygons->indices[ 00213 POINT_INDEX(polygons->end_indices,poly,vertex_index)]; 00214 00215 if( !point_done[point_index] ) 00216 { 00217 point_done[point_index] = TRUE; 00218 00219 moved = update_point_position( polygons, poly, vertex_index, 00220 point_index, current_points, &new_points[point_index], 00221 max_dist_from_original, fraction_to_move, normal_ratio, 00222 range_flag, volume, min_value, max_value ); 00223 00224 *avg_moved += moved; 00225 if( moved > *max_moved ) 00226 *max_moved = moved; 00227 } 00228 } 00229 00230 update_progress_report( &progress, poly+1 ); 00231 } 00232 00233 terminate_progress_report( &progress ); 00234 00235 *avg_moved /= (Real) polygons->n_points; 00236 } 00237 00238 #define MAX_NEIGHBOURS 100 00239 00240 /* ----------------------------- MNI Header ----------------------------------- 00241 @NAME : update_point_position 00242 @INPUT : 00243 @OUTPUT : 00244 @RETURNS : 00245 @DESCRIPTION: 00246 @METHOD : 00247 @GLOBALS : 00248 @CALLS : 00249 @CREATED : 1993 David MacDonald 00250 @MODIFIED : 00251 ---------------------------------------------------------------------------- */ 00252 00253 private Real update_point_position( 00254 polygons_struct *polygons, 00255 int poly, 00256 int vertex_index, 00257 int point_index, 00258 Point current_points[], 00259 Point *new_point, 00260 Real max_dist_from_original, 00261 Real fraction_to_move, 00262 Real normal_ratio, 00263 BOOLEAN range_flag, 00264 volume_struct *volume, 00265 int min_value, 00266 int max_value ) 00267 { 00268 int i, halves, n_neighbours; 00269 int neighbours[MAX_NEIGHBOURS]; 00270 Point neigh_points[MAX_NEIGHBOURS]; 00271 Point new_pos, current_scaled, centroid, destination; 00272 Real len, ratio, movement, sin_of_angle; 00273 Vector diff, delta, normal, unit_delta, unit_normal, cross; 00274 BOOLEAN interior_point; 00275 00276 n_neighbours = get_neighbours_of_point( polygons, poly, vertex_index, 00277 neighbours, MAX_NEIGHBOURS, 00278 &interior_point ); 00279 00280 if( interior_point && n_neighbours > 2 ) 00281 { 00282 for_less( i, 0, n_neighbours ) 00283 neigh_points[i] = current_points[neighbours[i]]; 00284 00285 get_points_centroid( n_neighbours, neigh_points, 00286 &centroid ); 00287 00288 destination = centroid; 00289 00290 if( normal_ratio > 0.0 ) 00291 { 00292 find_polygon_normal( n_neighbours, neigh_points, 00293 &normal ); 00294 00295 if( !null_Vector(&normal) ) 00296 { 00297 NORMALIZE_VECTOR( unit_normal, normal ); 00298 SUB_POINTS( delta, centroid, current_points[point_index] ); 00299 NORMALIZE_VECTOR( unit_delta, delta ); 00300 CROSS_VECTORS( cross, unit_delta, unit_normal ); 00301 sin_of_angle = MAGNITUDE( cross ); 00302 ratio = DOT_VECTORS( delta, normal ) / 00303 DOT_VECTORS( normal, normal ); 00304 SCALE_VECTOR( normal, normal, -ratio * normal_ratio * 00305 sin_of_angle); 00306 ADD_POINT_VECTOR( destination, centroid, normal ); 00307 } 00308 } 00309 00310 SCALE_POINT( new_pos, destination, fraction_to_move ); 00311 SCALE_POINT( current_scaled, current_points[point_index], 00312 1.0-fraction_to_move ); 00313 ADD_POINTS( new_pos, new_pos, current_scaled ); 00314 00315 SUB_POINTS( diff, new_pos, polygons->points[point_index] ); 00316 len = MAGNITUDE( diff ); 00317 00318 if( len > max_dist_from_original ) 00319 { 00320 ratio = max_dist_from_original / len; 00321 SCALE_VECTOR( diff, diff, ratio ); 00322 ADD_POINT_VECTOR( new_pos, polygons->points[point_index], diff ); 00323 } 00324 00325 if( range_flag ) 00326 { 00327 #define MAX_CUTS_IN_HALF 4 00328 for_less( halves, 0, MAX_CUTS_IN_HALF ) 00329 { 00330 if( point_inside_range( volume, &new_pos, 00331 min_value, max_value ) ) 00332 { 00333 break; 00334 } 00335 00336 ADD_POINTS( new_pos, new_pos, current_points[point_index] ); 00337 SCALE_POINT( new_pos, new_pos, 0.5 ); 00338 } 00339 00340 if( halves == MAX_CUTS_IN_HALF ) 00341 new_pos = current_points[point_index]; 00342 } 00343 00344 *new_point = new_pos; 00345 00346 SUB_POINTS( diff, new_pos, current_points[point_index] ); 00347 00348 movement = MAGNITUDE( diff ); 00349 } 00350 else 00351 movement = 0.0; 00352 00353 return( movement ); 00354 } 00355 00356 /* ----------------------------- MNI Header ----------------------------------- 00357 @NAME : point_inside_range 00358 @INPUT : volume 00359 point 00360 min_value 00361 max_value 00362 @OUTPUT : 00363 @RETURNS : TRUE or FALSE 00364 @DESCRIPTION: Tests if the volume value at the point is within the specified 00365 value range. 00366 @METHOD : 00367 @GLOBALS : 00368 @CALLS : 00369 @CREATED : 1993 David MacDonald 00370 @MODIFIED : 00371 ---------------------------------------------------------------------------- */ 00372 00373 private BOOLEAN point_inside_range( 00374 volume_struct *volume, 00375 Point *point, 00376 int min_value, 00377 int max_value ) 00378 { 00379 int vx, vy, vz, val; 00380 Real x, y, z; 00381 00382 convert_3D_world_to_voxel( volume, 00383 (Real) Point_x(*point), (Real) Point_y(*point), 00384 (Real) Point_z(*point), 00385 &x, &y, &z ); 00386 00387 vx = ROUND( x ); 00388 vy = ROUND( y ); 00389 vz = ROUND( z ); 00390 00391 val = (int) get_volume_voxel_value( volume, vx, vy, vz, 0, 0 ); 00392 00393 return( min_value <= val && val <= max_value ); 00394 }