BICPL: fill_volume.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/vols.h> 00017 00018 #ifndef lint 00019 static char rcsid[] = "$Header: /software/source//libraries/bicpl/Volumes/fill_volume.c,v 1.10 2000/02/06 15:30:54 stever Exp $"; 00020 #endif 00021 00022 typedef struct 00023 { 00024 int x, y, z; 00025 } xyz_struct; 00026 00027 /* ----------------------------- MNI Header ----------------------------------- 00028 @NAME : fill_connected_voxels 00029 @INPUT : volume 00030 label_volume 00031 connectivity 00032 voxel 00033 min_label_threshold 00034 max_label_threshold 00035 desired_label 00036 min_threshold 00037 max_threshold 00038 @OUTPUT : 00039 @RETURNS : TRUE if changed 00040 @DESCRIPTION: Performs a 3D fill from the specified voxel. 00041 @METHOD : 00042 @GLOBALS : 00043 @CALLS : 00044 @CREATED : 1993 David MacDonald 00045 @MODIFIED : 00046 ---------------------------------------------------------------------------- */ 00047 00048 public BOOLEAN fill_connected_voxels( 00049 Volume volume, 00050 Volume label_volume, 00051 Neighbour_types connectivity, 00052 int voxel[], 00053 int min_label_threshold, 00054 int max_label_threshold, 00055 int desired_label, 00056 Real min_threshold, 00057 Real max_threshold, 00058 int range_changed[2][N_DIMENSIONS] ) 00059 { 00060 int dir, n_dirs, *dx, *dy, *dz, dim; 00061 int x, y, z, tx, ty, tz; 00062 int sizes[N_DIMENSIONS]; 00063 int voxel_index[MAX_DIMENSIONS]; 00064 xyz_struct entry; 00065 QUEUE_STRUCT( xyz_struct ) queue; 00066 bitlist_3d_struct checked_flags, change_flags; 00067 BOOLEAN first; 00068 00069 if( !should_change_this_one( volume, label_volume, voxel, 00070 min_threshold, max_threshold, 00071 min_label_threshold, max_label_threshold, 00072 desired_label ) ) 00073 return( FALSE ); 00074 00075 n_dirs = get_3D_neighbour_directions( connectivity, &dx, &dy, &dz ); 00076 00077 get_volume_sizes( volume, sizes ); 00078 00079 create_bitlist_3d( sizes[X], sizes[Y], sizes[Z], &checked_flags ); 00080 create_bitlist_3d( sizes[X], sizes[Y], sizes[Z], &change_flags ); 00081 00082 INITIALIZE_QUEUE( queue ); 00083 00084 set_bitlist_bit_3d( &checked_flags, voxel[X], voxel[Y], voxel[Z], TRUE ); 00085 set_bitlist_bit_3d( &change_flags, voxel[X], voxel[Y], voxel[Z], TRUE ); 00086 00087 entry.x = voxel[X]; 00088 entry.y = voxel[Y]; 00089 entry.z = voxel[Z]; 00090 INSERT_IN_QUEUE( queue, entry ); 00091 00092 while( !IS_QUEUE_EMPTY( queue ) ) 00093 { 00094 REMOVE_FROM_QUEUE( queue, entry ); 00095 00096 x = entry.x; 00097 y = entry.y; 00098 z = entry.z; 00099 00100 for_less( dir, 0, n_dirs ) 00101 { 00102 tx = x + dx[dir]; 00103 ty = y + dy[dir]; 00104 tz = z + dz[dir]; 00105 00106 if( tx >= 0 && tx < sizes[X] && 00107 ty >= 0 && ty < sizes[Y] && 00108 tz >= 0 && tz < sizes[Z] && 00109 !get_bitlist_bit_3d( &checked_flags, tx, ty, tz ) ) 00110 { 00111 set_bitlist_bit_3d( &checked_flags, tx, ty, tz, TRUE ); 00112 00113 voxel_index[X] = tx; 00114 voxel_index[Y] = ty; 00115 voxel_index[Z] = tz; 00116 if( should_change_this_one( volume, label_volume, voxel_index, 00117 min_threshold, max_threshold, 00118 min_label_threshold, max_label_threshold, 00119 desired_label ) ) 00120 { 00121 set_bitlist_bit_3d( &change_flags, tx, ty, tz, TRUE ); 00122 entry.x = tx; 00123 entry.y = ty; 00124 entry.z = tz; 00125 INSERT_IN_QUEUE( queue, entry ); 00126 } 00127 } 00128 } 00129 } 00130 00131 first = TRUE; 00132 00133 for_less( voxel_index[X], 0, sizes[X] ) 00134 for_less( voxel_index[Y], 0, sizes[Y] ) 00135 for_less( voxel_index[Z], 0, sizes[Z] ) 00136 { 00137 if( get_bitlist_bit_3d( &change_flags, voxel_index[X], voxel_index[Y], 00138 voxel_index[Z] ) ) 00139 { 00140 set_volume_label_data( label_volume, voxel_index, desired_label ); 00141 00142 for_less( dim, 0, N_DIMENSIONS ) 00143 { 00144 if( first || voxel_index[dim] < range_changed[0][dim] ) 00145 range_changed[0][dim] = voxel_index[dim]; 00146 if( first || voxel_index[dim] > range_changed[1][dim] ) 00147 range_changed[1][dim] = voxel_index[dim]; 00148 } 00149 00150 first = FALSE; 00151 } 00152 } 00153 00154 delete_bitlist_3d( &checked_flags ); 00155 delete_bitlist_3d( &change_flags ); 00156 00157 DELETE_QUEUE( queue ); 00158 00159 return( TRUE ); 00160 } 00161 00162 private int Dx4[4] = { 1, 0, -1, 0 }; 00163 private int Dy4[4] = { 0, 1, 0, -1 }; 00164 00165 private int Dx8[8] = { 1, 1, 0, -1, -1, -1, 0, 1 }; 00166 private int Dy8[8] = { 0, 1, 1, 1, 0, -1, -1, -1 }; 00167 00168 /* ----------------------------- MNI Header ----------------------------------- 00169 @NAME : get_neighbour_directions 00170 @INPUT : connectivity 00171 @OUTPUT : dx 00172 dy 00173 @RETURNS : number of directions 00174 @DESCRIPTION: Passes back the neighbour directions of a 2D grid. 00175 @METHOD : 00176 @GLOBALS : 00177 @CALLS : 00178 @CREATED : 1993 David MacDonald 00179 @MODIFIED : 00180 ---------------------------------------------------------------------------- */ 00181 00182 public int get_neighbour_directions( 00183 Neighbour_types connectivity, 00184 int *dx[], 00185 int *dy[] ) 00186 { 00187 int n_dirs; 00188 00189 switch( connectivity ) 00190 { 00191 case FOUR_NEIGHBOURS: 00192 *dx = Dx4; 00193 *dy = Dy4; 00194 n_dirs = 4; 00195 break; 00196 00197 case EIGHT_NEIGHBOURS: 00198 *dx = Dx8; 00199 *dy = Dy8; 00200 n_dirs = 8; 00201 break; 00202 } 00203 00204 return( n_dirs ); 00205 } 00206 00207 private int Dx6[6] = { 1, 0, 0, -1, 0, 0 }; 00208 private int Dy6[6] = { 0, 1, 0, 0, -1, 0 }; 00209 private int Dz6[6] = { 0, 0, 1, 0, 0, -1 }; 00210 00211 private int Dx26[26]; 00212 private int Dy26[26]; 00213 private int Dz26[26]; 00214 00215 /* ----------------------------- MNI Header ----------------------------------- 00216 @NAME : create_3D_neighbours 00217 @INPUT : 00218 @OUTPUT : 00219 @RETURNS : 00220 @DESCRIPTION: Initializes the 3D neighbours. 00221 @METHOD : 00222 @GLOBALS : 00223 @CALLS : 00224 @CREATED : 1993 David MacDonald 00225 @MODIFIED : 00226 ---------------------------------------------------------------------------- */ 00227 00228 private void create_3D_neighbours( void ) 00229 { 00230 int x, y, z, n; 00231 00232 n = 0; 00233 for_inclusive( x, -1, 1 ) 00234 for_inclusive( y, -1, 1 ) 00235 for_inclusive( z, -1, 1 ) 00236 { 00237 if( x != 0 || y != 0 || z != 0 ) 00238 { 00239 Dx26[n] = x; 00240 Dy26[n] = y; 00241 Dz26[n] = z; 00242 ++n; 00243 } 00244 } 00245 00246 if( n != 26 ) 00247 { 00248 HANDLE_INTERNAL_ERROR( "create_3D_neighbours" ); 00249 } 00250 } 00251 00252 /* ----------------------------- MNI Header ----------------------------------- 00253 @NAME : get_3D_neighbour_directions 00254 @INPUT : connectivity 00255 @OUTPUT : dx 00256 dy 00257 dz 00258 @RETURNS : number of directions 00259 @DESCRIPTION: Passes back the neighbour directions of a 3D grid. 00260 @METHOD : 00261 @GLOBALS : 00262 @CALLS : 00263 @CREATED : 1993 David MacDonald 00264 @MODIFIED : 00265 ---------------------------------------------------------------------------- */ 00266 00267 public int get_3D_neighbour_directions( 00268 Neighbour_types connectivity, 00269 int *dx[], 00270 int *dy[], 00271 int *dz[] ) 00272 { 00273 static BOOLEAN first = TRUE; 00274 int n_dirs; 00275 00276 if( first ) 00277 { 00278 first = FALSE; 00279 create_3D_neighbours(); 00280 } 00281 00282 switch( connectivity ) 00283 { 00284 case FOUR_NEIGHBOURS: 00285 *dx = Dx6; 00286 *dy = Dy6; 00287 *dz = Dz6; 00288 n_dirs = 6; 00289 break; 00290 00291 case EIGHT_NEIGHBOURS: 00292 *dx = Dx26; 00293 *dy = Dy26; 00294 *dz = Dz26; 00295 n_dirs = 26; 00296 break; 00297 } 00298 00299 return( n_dirs ); 00300 } 00301 00302 /* ----------------------------- MNI Header ----------------------------------- 00303 @NAME : should_change_this_one 00304 @INPUT : volume 00305 label_volume 00306 voxel 00307 min_threshold 00308 max_threshold 00309 label_min_threshold 00310 label_max_threshold 00311 desired_label 00312 @OUTPUT : 00313 @RETURNS : TRUE if should change this 00314 @DESCRIPTION: 00315 @METHOD : 00316 @GLOBALS : 00317 @CALLS : 00318 @CREATED : 1993 David MacDonald 00319 @MODIFIED : 00320 ---------------------------------------------------------------------------- */ 00321 00322 public BOOLEAN should_change_this_one( 00323 Volume volume, 00324 Volume label_volume, 00325 int voxel[], 00326 Real min_threshold, 00327 Real max_threshold, 00328 int label_min_threshold, 00329 int label_max_threshold, 00330 int desired_label ) 00331 { 00332 int label; 00333 Real value; 00334 00335 label = get_volume_label_data( label_volume, voxel ); 00336 00337 if( label == desired_label ) 00338 return( FALSE ); 00339 00340 if( label_min_threshold <= label_max_threshold && 00341 (label < label_min_threshold || label > label_max_threshold) ) 00342 return( FALSE ); 00343 00344 if( min_threshold > max_threshold ) 00345 return( TRUE ); 00346 00347 value = get_volume_real_value( volume, voxel[X], voxel[Y], voxel[Z], 0, 0 ); 00348 00349 return( min_threshold <= value && value <= max_threshold ); 00350 }