BICPL: dilate.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/dilate.c,v 1.13 2000/02/06 15:30:54 stever Exp $"; 00020 #endif 00021 00022 typedef enum { NOT_INVOLVED, INSIDE_REGION, CANDIDATE } 00023 Voxel_classes; 00024 00025 /* ----------------------------- MNI Header ----------------------------------- 00026 @NAME : dilate_voxels_3d 00027 @INPUT : volume 00028 label_volume 00029 min_inside_label 00030 max_inside_label 00031 min_inside_value 00032 max_inside_value 00033 min_outside_label 00034 max_outside_label 00035 min_outside_value 00036 max_outside_value 00037 new_label 00038 connectivity 00039 @OUTPUT : 00040 @RETURNS : 00041 @DESCRIPTION: Dilates the label volume. 00042 @METHOD : 00043 @GLOBALS : 00044 @CALLS : 00045 @CREATED : 1993 David MacDonald 00046 @MODIFIED : 00047 ---------------------------------------------------------------------------- */ 00048 00049 public int dilate_voxels_3d( 00050 Volume volume, 00051 Volume label_volume, 00052 Real min_inside_label, 00053 Real max_inside_label, 00054 Real min_inside_value, 00055 Real max_inside_value, 00056 Real min_outside_label, 00057 Real max_outside_label, 00058 Real min_outside_value, 00059 Real max_outside_value, 00060 Real new_label, 00061 Neighbour_types connectivity, 00062 int range_changed[2][N_DIMENSIONS] ) 00063 { 00064 int n_changed; 00065 int x, y, z, delta_x, tx, ty, tz; 00066 int sizes[N_DIMENSIONS]; 00067 int dir, n_dirs, *dx, *dy, *dz; 00068 Real value, label, *value_row, *label_row; 00069 Smallest_int **voxel_classes[3], **swap; 00070 progress_struct progress; 00071 Voxel_classes voxel_class; 00072 BOOLEAN use_label_volume, use_volume, at_end, at_edge_y; 00073 BOOLEAN inside_specified, outside_specified; 00074 BOOLEAN inside, outside; 00075 00076 use_label_volume = (min_inside_label <= max_inside_label || 00077 min_outside_label <= max_outside_label); 00078 use_volume = (min_inside_value <= max_inside_value || 00079 min_outside_value <= max_outside_value); 00080 00081 inside_specified = (min_inside_label <= max_inside_label || 00082 min_inside_value <= max_inside_value); 00083 outside_specified = (min_outside_label <= max_outside_label || 00084 min_outside_value <= max_outside_value); 00085 00086 if( !inside_specified && !outside_specified ) 00087 { 00088 min_inside_label = new_label; 00089 max_inside_label = new_label; 00090 inside_specified = TRUE; 00091 use_label_volume = TRUE; 00092 } 00093 00094 n_dirs = get_3D_neighbour_directions( connectivity, &dx, &dy, &dz ); 00095 00096 get_volume_sizes( label_volume, sizes ); 00097 00098 for_less( x, 0, 3 ) 00099 ALLOC2D( voxel_classes[x], sizes[Y]+2, sizes[Z]+2 ); 00100 00101 for_less( x, 0, 1 ) 00102 { 00103 for_less( y, 0, sizes[Y] + 2 ) 00104 { 00105 for_less( z, 0, sizes[Z] + 2 ) 00106 voxel_classes[x][y][z] = (Smallest_int) NOT_INVOLVED; 00107 } 00108 } 00109 00110 ALLOC( value_row, sizes[Z] ); 00111 ALLOC( label_row, sizes[Z] ); 00112 00113 initialize_progress_report( &progress, FALSE, sizes[X], 00114 "Expanding labeled voxels" ); 00115 00116 n_changed = 0; 00117 00118 for_less( x, 0, sizes[X] ) 00119 { 00120 for_less( delta_x, (x == 0) ? 0 : 1, 2 ) 00121 { 00122 at_end = (x + delta_x == sizes[X]); 00123 00124 for_less( y, -1, sizes[Y] + 1 ) 00125 { 00126 at_edge_y = (y == -1 || y == sizes[Y]); 00127 voxel_classes[delta_x+1][y+1][0] = (Smallest_int) NOT_INVOLVED; 00128 voxel_classes[delta_x+1][y+1][sizes[Z]+1] = NOT_INVOLVED; 00129 00130 if( !at_edge_y && !at_end ) 00131 { 00132 if( use_label_volume ) 00133 { 00134 get_volume_value_hyperslab_3d( label_volume, 00135 x + delta_x, y, 0, 00136 1, 1, sizes[Z], 00137 label_row ); 00138 } 00139 00140 if( use_volume ) 00141 { 00142 get_volume_value_hyperslab_3d( volume, 00143 x + delta_x, y, 0, 00144 1, 1, sizes[Z], 00145 value_row ); 00146 } 00147 } 00148 00149 for_less( z, 0, sizes[Z] ) 00150 { 00151 if( at_edge_y || at_end ) 00152 { 00153 voxel_class = NOT_INVOLVED; 00154 } 00155 else 00156 { 00157 if( use_label_volume ) 00158 label = label_row[z]; 00159 00160 if( use_volume ) 00161 value = value_row[z]; 00162 00163 inside = (min_inside_label > max_inside_label || 00164 (min_inside_label <= label && 00165 label <= max_inside_label)) && 00166 (min_inside_value > max_inside_value || 00167 (min_inside_value <= value && 00168 value <= max_inside_value)); 00169 00170 outside = (min_outside_label > max_outside_label || 00171 (min_outside_label <= label && 00172 label <= max_outside_label)) && 00173 (min_outside_value > max_outside_value || 00174 (min_outside_value <= value && 00175 value <= max_outside_value)); 00176 00177 if( inside_specified ) 00178 { 00179 if( inside ) 00180 voxel_class = INSIDE_REGION; 00181 else if( outside ) 00182 voxel_class = CANDIDATE; 00183 else 00184 voxel_class = NOT_INVOLVED; 00185 } 00186 else 00187 { 00188 if( outside ) 00189 voxel_class = CANDIDATE; 00190 else 00191 voxel_class = INSIDE_REGION; 00192 } 00193 } 00194 00195 voxel_classes[delta_x+1][y+1][z+1] = 00196 (Smallest_int) voxel_class; 00197 } 00198 } 00199 } 00200 00201 for_less( y, 0, sizes[Y] ) 00202 { 00203 for_less( z, 0, sizes[Z] ) 00204 { 00205 if( voxel_classes[1][y+1][z+1] == CANDIDATE ) 00206 { 00207 for_less( dir, 0, n_dirs ) 00208 { 00209 tx = 0 + dx[dir] + 1; 00210 ty = y + dy[dir] + 1; 00211 tz = z + dz[dir] + 1; 00212 00213 if( voxel_classes[tx][ty][tz] == INSIDE_REGION ) 00214 { 00215 set_volume_real_value( label_volume, x, y, z, 0, 0, 00216 new_label ); 00217 00218 if( n_changed == 0 || x < range_changed[0][X] ) 00219 range_changed[0][X] = x; 00220 if( n_changed == 0 || x > range_changed[1][X] ) 00221 range_changed[1][X] = x; 00222 00223 if( n_changed == 0 || y < range_changed[0][Y] ) 00224 range_changed[0][Y] = y; 00225 if( n_changed == 0 || y > range_changed[1][Y] ) 00226 range_changed[1][Y] = y; 00227 00228 if( n_changed == 0 || z < range_changed[0][Z] ) 00229 range_changed[0][Z] = z; 00230 if( n_changed == 0 || z > range_changed[1][Z] ) 00231 range_changed[1][Z] = z; 00232 00233 ++n_changed; 00234 break; 00235 } 00236 } 00237 } 00238 } 00239 } 00240 00241 swap = voxel_classes[0]; 00242 voxel_classes[0] = voxel_classes[1]; 00243 voxel_classes[1] = voxel_classes[2]; 00244 voxel_classes[2] = swap; 00245 00246 update_progress_report( &progress, x + 1 ); 00247 } 00248 00249 terminate_progress_report( &progress ); 00250 00251 for_less( x, 0, 3 ) 00252 FREE2D( voxel_classes[x] ); 00253 00254 FREE( value_row ); 00255 FREE( label_row ); 00256 00257 return( n_changed ); 00258 }