BICPL: scan_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.h> 00017 00018 #ifndef lint 00019 static char rcsid[] = "$Header: /software/source//libraries/bicpl/Volumes/scan_polygons.c,v 1.14 2001/08/11 20:09:08 stever Exp $"; 00020 #endif 00021 00022 #define MAX_TEMP_STORAGE 1000 00023 00024 00040 private void recursive_scan_polygon_to_voxels( 00041 int size, 00042 Point points[], 00043 Volume label_volume, 00044 int label, 00045 int min_voxel[], 00046 int max_voxel[] ) 00047 { 00048 int n_left, n_right, dim, max_dim, save, pos; 00049 Real slice_pos; 00050 Point right_vertices[MAX_TEMP_STORAGE]; 00051 Point left_vertices[MAX_TEMP_STORAGE]; 00052 Vector normal; 00053 00054 max_dim = 0; 00055 for_less( dim, 1, N_DIMENSIONS ) 00056 { 00057 if( max_voxel[dim] - min_voxel[dim] > 00058 max_voxel[max_dim] - min_voxel[max_dim] ) 00059 max_dim = dim; 00060 } 00061 00062 if( min_voxel[max_dim] == max_voxel[max_dim] ) 00063 { 00064 set_volume_label_data( label_volume, min_voxel, label ); 00065 return; 00066 } 00067 00068 pos = (min_voxel[max_dim] + max_voxel[max_dim]) / 2; 00069 slice_pos = (Real) pos + 0.5; 00070 00071 fill_Vector( normal, 0.0, 0.0, 0.0 ); 00072 Vector_coord( normal, max_dim ) = -1.0f; 00073 00074 split_polygon_with_plane( size, points, slice_pos, 00075 &normal, &n_left, left_vertices, 00076 &n_right, right_vertices ); 00077 00078 if( n_left > 0 ) 00079 { 00080 save = max_voxel[max_dim]; 00081 max_voxel[max_dim] = pos; 00082 recursive_scan_polygon_to_voxels( n_left, left_vertices, 00083 label_volume, label, 00084 min_voxel, max_voxel ); 00085 max_voxel[max_dim] = save; 00086 } 00087 00088 if( n_right > 0 ) 00089 { 00090 save = min_voxel[max_dim]; 00091 min_voxel[max_dim] = pos+1; 00092 recursive_scan_polygon_to_voxels( n_right, right_vertices, 00093 label_volume, label, 00094 min_voxel, max_voxel ); 00095 min_voxel[max_dim] = save; 00096 } 00097 } 00098 00099 /* ----------------------------- MNI Header ----------------------------------- 00100 @NAME : scan_a_polygon 00101 @INPUT : size 00102 vertices 00103 voxels 00104 output_vertices 00105 polygons 00106 volume 00107 label_volume 00108 sizes 00109 label 00110 @OUTPUT : 00111 @RETURNS : 00112 @DESCRIPTION: Scans a polygon into the volume. 00113 @METHOD : 00114 @GLOBALS : 00115 @CALLS : 00116 @CREATED : Jul. 23, 1997 David MacDonald 00117 @MODIFIED : 00118 ---------------------------------------------------------------------------- */ 00119 00134 private void scan_a_polygon( 00135 int size, 00136 Point vertices[], 00137 Point voxels[], 00138 int n_output_vertices, 00139 Point output_vertices[], 00140 Volume volume, 00141 Volume label_volume, 00142 int sizes[], 00143 int label ) 00144 { 00145 int vertex, dim, n_clip; 00146 Real voxel[N_DIMENSIONS]; 00147 Real min_voxel[N_DIMENSIONS], max_voxel[N_DIMENSIONS]; 00148 int min_iv[N_DIMENSIONS], max_iv[N_DIMENSIONS]; 00149 00150 min_voxel[X] = 0.0; /*--- to avoid warnings */ 00151 max_voxel[X] = 0.0; 00152 00153 for_less( vertex, 0, size ) 00154 { 00155 convert_world_to_voxel( volume, 00156 RPoint_x(vertices[vertex]), 00157 RPoint_y(vertices[vertex]), 00158 RPoint_z(vertices[vertex]), 00159 voxel ); 00160 00161 fill_Point( voxels[vertex], voxel[X], voxel[Y], voxel[Z] ); 00162 00163 if( vertex == 0 ) 00164 { 00165 for_less( dim, 0, N_DIMENSIONS ) 00166 { 00167 min_voxel[dim] = voxel[dim]; 00168 max_voxel[dim] = voxel[dim]; 00169 } 00170 } 00171 else 00172 { 00173 for_less( dim, 0, N_DIMENSIONS ) 00174 { 00175 if( voxel[dim] < min_voxel[dim] ) 00176 min_voxel[dim] = voxel[dim]; 00177 if( voxel[dim] > max_voxel[dim] ) 00178 max_voxel[dim] = voxel[dim]; 00179 } 00180 } 00181 } 00182 00183 for_less( dim, 0, N_DIMENSIONS ) 00184 { 00185 min_iv[dim] = ROUND( min_voxel[dim] ); 00186 if( min_iv[dim] < 0 ) 00187 min_iv[dim] = 0; 00188 00189 max_iv[dim] = ROUND( max_voxel[dim] ); 00190 if( max_iv[dim] >= sizes[dim] ) 00191 max_iv[dim] = sizes[dim]-1; 00192 } 00193 00194 n_clip = clip_polygon_against_box( size, voxels, 00195 (Real) min_iv[X] - 0.5, 00196 (Real) max_iv[X] + 0.5, 00197 (Real) min_iv[Y] - 0.5, 00198 (Real) max_iv[Y] + 0.5, 00199 (Real) min_iv[Z] - 0.5, 00200 (Real) max_iv[Z] + 0.5, 00201 n_output_vertices, 00202 output_vertices ); 00203 00204 if( n_clip > 0 ) 00205 { 00206 recursive_scan_polygon_to_voxels( n_clip, output_vertices, 00207 label_volume, label, 00208 min_iv, max_iv ); 00209 } 00210 } 00211 00212 /* ----------------------------- MNI Header ----------------------------------- 00213 @NAME : scan_polygons_to_voxels 00214 @INPUT : polygons 00215 volume 00216 label_volume 00217 label 00218 max_size 00219 @OUTPUT : 00220 @RETURNS : 00221 @DESCRIPTION: Subdivides polygons until they are smaller than max_size 00222 size, then labels the voxel corresponding to the centre of 00223 each fragment. Only polygons with 3 or 4 vertices will be 00224 scanned, should be rewritten to handle more. 00225 @METHOD : 00226 @GLOBALS : 00227 @CALLS : 00228 @CREATED : 1993 David MacDonald 00229 @MODIFIED : Jul. 27, 1997 D. MacDonald - broke into 2 parts 00230 ---------------------------------------------------------------------------- */ 00231 00240 public void scan_polygons_to_voxels( 00241 polygons_struct *polygons, 00242 Volume volume, 00243 Volume label_volume, 00244 int label, 00245 Real max_distance ) 00246 { 00247 int vertex, poly, size, point_index, max_size; 00248 Point *vertices, *voxels, *output_vertices; 00249 int sizes[N_DIMENSIONS]; 00250 00251 get_volume_sizes( label_volume, sizes ); 00252 00253 max_size = 0; 00254 for_less( poly, 0, polygons->n_items ) 00255 { 00256 size = GET_OBJECT_SIZE( *polygons, poly ); 00257 max_size = MAX( max_size, size ); 00258 } 00259 00260 ALLOC( vertices, max_size ); 00261 ALLOC( voxels, max_size ); 00262 /* Possible bug here: a triangle can be cut into a 7-gon after 00263 * clipping against a box. What is true upper limit on output_vertices? 00264 */ 00265 ALLOC( output_vertices, 2 * max_size ); 00266 00267 for_less( poly, 0, polygons->n_items ) 00268 { 00269 size = GET_OBJECT_SIZE( *polygons, poly ); 00270 00271 for_less( vertex, 0, size ) 00272 { 00273 point_index = polygons->indices[POINT_INDEX( 00274 polygons->end_indices,poly,vertex)]; 00275 vertices[vertex] = polygons->points[point_index]; 00276 } 00277 00278 scan_a_polygon( size, vertices, voxels, 00279 2*max_size, output_vertices, 00280 volume, label_volume, sizes, label ); 00281 00282 } 00283 00284 FREE( vertices ); 00285 FREE( voxels ); 00286 FREE( output_vertices ); 00287 } 00288 00289 /* ----------------------------- MNI Header ----------------------------------- 00290 @NAME : scan_quadmesh_to_voxels 00291 @INPUT : quadmesh 00292 volume 00293 label_volume 00294 label 00295 max_distance - ignored 00296 @OUTPUT : 00297 @RETURNS : 00298 @DESCRIPTION: Scans a quadmesh into the volume. 00299 @METHOD : 00300 @GLOBALS : 00301 @CALLS : 00302 @CREATED : Jul. 23, 1997 David MacDonald 00303 @MODIFIED : 00304 ---------------------------------------------------------------------------- */ 00305 00314 public void scan_quadmesh_to_voxels( 00315 quadmesh_struct *quadmesh, 00316 Volume volume, 00317 Volume label_volume, 00318 int label, 00319 Real max_distance ) 00320 { 00321 int i, j, m, n; 00322 Point vertices[4], voxels[4], output_vertices[4]; 00323 int sizes[N_DIMENSIONS]; 00324 00325 get_volume_sizes( label_volume, sizes ); 00326 00327 get_quadmesh_n_objects( quadmesh, &m, &n ); 00328 00329 for_less( i, 0, m ) 00330 { 00331 for_less( j, 0, n ) 00332 { 00333 get_quadmesh_patch( quadmesh, i, j, vertices ); 00334 00335 scan_a_polygon( 4, vertices, voxels, 00336 4, output_vertices, 00337 volume, label_volume, sizes, label ); 00338 00339 } 00340 } 00341 }