BICPL: box_filter.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 #include <bicpl/numerical.h> 00018 00019 #ifndef lint 00020 static char rcsid[] = "$Header: /software/source//libraries/bicpl/Volumes/box_filter.c,v 1.19 2000/02/06 15:30:52 stever Exp $"; 00021 #endif 00022 00023 #define DEBUG 00024 #undef DEBUG 00025 00026 #ifdef DEBUG 00027 private Real get_correct_amount( 00028 Volume volume, 00029 int x, 00030 int y, 00031 int z, 00032 Real x_width, 00033 Real y_width, 00034 Real z_width ); 00035 #endif 00036 00037 #define INITIALIZE_SAMPLE( half_width, receding, advancing, left_weight, right_weight ) \ 00038 { \ 00039 Real min_pos, max_pos; \ 00040 \ 00041 min_pos = -(half_width); \ 00042 max_pos = (half_width); \ 00043 \ 00044 (receding) = ROUND( min_pos ); \ 00045 (advancing) = ROUND( max_pos ); \ 00046 \ 00047 (left_weight) = (Real) (receding) + 0.5 - min_pos; \ 00048 if( (left_weight) >= 1.0 ) \ 00049 --(advancing); \ 00050 (right_weight) = 1.0 - (left_weight); \ 00051 } 00052 00053 #define GET_FIRST_SAMPLE( advancing, left_weight, size, data, sample ) \ 00054 { \ 00055 int _I; \ 00056 \ 00057 (sample) = 0.0; \ 00058 for_inclusive( _I, 0, MIN( (size)-1, (advancing)-1) ) \ 00059 (sample) += (data); \ 00060 \ 00061 if( (advancing) < (size) ) \ 00062 (sample) += (data) * (left_weight); \ 00063 } 00064 00065 #define GET_NEXT_SAMPLE( receding, advancing, left_weight, right_weight, size, data, sample ) \ 00066 { \ 00067 int _I; \ 00068 \ 00069 if( (advancing) < (size) ) \ 00070 { \ 00071 _I = (advancing); \ 00072 (sample) += (data) * (right_weight); \ 00073 } \ 00074 if( (advancing) < (size)-1 ) \ 00075 { \ 00076 _I = (advancing) + 1; \ 00077 (sample) += (data) * (left_weight); \ 00078 } \ 00079 \ 00080 if( (receding) >= 0 ) \ 00081 { \ 00082 _I = (receding); \ 00083 (sample) -= (data) * (left_weight); \ 00084 } \ 00085 if( (receding) >= -1 ) \ 00086 { \ 00087 _I = (receding) + 1; \ 00088 (sample) -= (data) * (right_weight); \ 00089 } \ 00090 } 00091 00092 /* ----------------------------- MNI Header ----------------------------------- 00093 @NAME : create_box_filtered_volume 00094 @INPUT : volume - volume to filter 00095 nc_data_type - if NC_UNSPECIFIED, this and next 3 args ignored 00096 sign_flag 00097 real_min_value 00098 real_max_value 00099 x_width - full width of box filter 00100 y_width 00101 z_width 00102 @OUTPUT : 00103 @RETURNS : filtered volume 00104 @DESCRIPTION: Filters a volume with a box filter, creating a new volume with 00105 the same number of samples per dimension. 00106 @METHOD : 00107 @GLOBALS : 00108 @CALLS : 00109 @CREATED : 1993 David MacDonald 00110 @MODIFIED : 00111 ---------------------------------------------------------------------------- */ 00112 00113 public Volume create_box_filtered_volume( 00114 Volume volume, 00115 nc_type nc_data_type, 00116 BOOLEAN sign_flag, 00117 Real real_min_value, 00118 Real real_max_value, 00119 Real x_width, 00120 Real y_width, 00121 Real z_width ) 00122 { 00123 Volume resampled_volume; 00124 int x, y, z; 00125 int sizes[MAX_DIMENSIONS], nx_required; 00126 #ifdef DEBUG 00127 Real correct_voxel; 00128 #endif 00129 Real total_volume, value, sum, sample; 00130 Real *volume_row, ***volume_cache; 00131 int n_alloced, max_alloced; 00132 int x_init_recede_index, y_init_recede_index; 00133 int z_init_recede_index; 00134 int x_init_advance_index, y_init_advance_index; 00135 int z_init_advance_index; 00136 int x_recede_index, y_recede_index, z_recede_index; 00137 int x_advance_index, y_advance_index, z_advance_index; 00138 Real x_left_weight, x_right_weight; 00139 Real y_left_weight, y_right_weight; 00140 Real z_left_weight, z_right_weight; 00141 float **slice, *row; 00142 progress_struct progress; 00143 00144 if( get_volume_n_dimensions(volume) != 3 ) 00145 { 00146 handle_internal_error( 00147 "create_box_filtered_volume: volume must be 3D.\n" ); 00148 } 00149 00150 x_width = FABS( x_width ); 00151 y_width = FABS( y_width ); 00152 z_width = FABS( z_width ); 00153 00154 if( x_width < 1.0 ) 00155 x_width = 1.0; 00156 if( y_width < 1.0 ) 00157 y_width = 1.0; 00158 if( z_width < 1.0 ) 00159 z_width = 1.0; 00160 00161 get_volume_sizes( volume, sizes ); 00162 00163 resampled_volume = copy_volume_definition( volume, nc_data_type, sign_flag, 00164 real_min_value, real_max_value ); 00165 00166 total_volume = x_width * y_width * z_width; 00167 00168 x_width /= 2.0; 00169 y_width /= 2.0; 00170 z_width /= 2.0; 00171 00172 INITIALIZE_SAMPLE( x_width, x_init_recede_index, x_init_advance_index, 00173 x_left_weight, x_right_weight ) 00174 INITIALIZE_SAMPLE( y_width, y_init_recede_index, y_init_advance_index, 00175 y_left_weight, y_right_weight ) 00176 INITIALIZE_SAMPLE( z_width, z_init_recede_index, z_init_advance_index, 00177 z_left_weight, z_right_weight ) 00178 00179 initialize_progress_report( &progress, FALSE, sizes[X] * sizes[Y], 00180 "Box Filtering" ); 00181 00182 ALLOC2D( slice, sizes[Y], sizes[Z] ); 00183 ALLOC( row, sizes[Z] ); 00184 00185 nx_required = MAX( sizes[X], x_init_advance_index + 1 ); 00186 ALLOC( volume_row, nx_required ); 00187 00188 for_less( y, 0, sizes[Y] ) 00189 { 00190 for_less( z, 0, sizes[Z] ) 00191 { 00192 get_volume_value_hyperslab_3d( volume, 0, y, z, nx_required, 1, 1, 00193 volume_row ); 00194 GET_FIRST_SAMPLE( x_init_advance_index, 00195 x_left_weight, sizes[X], 00196 volume_row[_I], sample ); 00197 slice[y][z] = (float) sample; 00198 } 00199 } 00200 00201 FREE( volume_row ); 00202 00203 ALLOC( volume_cache, sizes[X] ); 00204 for_less( x, 0, sizes[X] ) 00205 volume_cache[x] = NULL; 00206 n_alloced = 0; 00207 ALLOC2D( volume_cache[0], sizes[Y], sizes[Z] ); 00208 get_volume_value_hyperslab_3d( volume, 0, 0, 0, 1, sizes[Y], sizes[Z], 00209 &volume_cache[0][0][0] ); 00210 ++n_alloced; 00211 if( x_init_recede_index >= 0 ) 00212 { 00213 ALLOC2D( volume_cache[1], sizes[Y], sizes[Z] ); 00214 get_volume_value_hyperslab_3d( volume, 1, 0, 0, 1, sizes[Y], sizes[Z], 00215 &volume_cache[1][0][0] ); 00216 ++n_alloced; 00217 } 00218 00219 if( x_init_advance_index < sizes[X] ) 00220 { 00221 ALLOC2D( volume_cache[x_init_advance_index], sizes[Y], sizes[Z] ); 00222 get_volume_value_hyperslab_3d( volume, x_init_advance_index, 0, 0, 00223 1, sizes[Y], sizes[Z], 00224 &volume_cache[x_init_advance_index][0][0] ); 00225 ++n_alloced; 00226 } 00227 00228 if( x_init_advance_index+1 < sizes[X] ) 00229 { 00230 ALLOC2D( volume_cache[x_init_advance_index+1], sizes[Y], sizes[Z] ); 00231 get_volume_value_hyperslab_3d( volume, x_init_advance_index+1, 0, 0, 00232 1, sizes[Y], sizes[Z], 00233 &volume_cache[x_init_advance_index+1][0][0] ); 00234 ++n_alloced; 00235 } 00236 00237 if( x_init_recede_index == x_init_advance_index ) 00238 max_alloced = 2; 00239 else if( x_init_recede_index+1 == x_init_advance_index ) 00240 max_alloced = 3; 00241 else 00242 max_alloced = 4; 00243 00244 x_recede_index = x_init_recede_index; 00245 x_advance_index = x_init_advance_index; 00246 00247 for_less( x, 0, sizes[X] ) 00248 { 00249 for_less( z, 0, sizes[Z] ) 00250 { 00251 GET_FIRST_SAMPLE( y_init_advance_index, y_left_weight, sizes[Y], 00252 (Real) slice[_I][z], sample ) 00253 row[z] = (float) sample; 00254 } 00255 00256 y_recede_index = y_init_recede_index; 00257 y_advance_index = y_init_advance_index; 00258 00259 for_less( y, 0, sizes[Y] ) 00260 { 00261 GET_FIRST_SAMPLE( z_init_advance_index, z_left_weight, 00262 sizes[Z], (Real) row[_I], sum ) 00263 z_recede_index = z_init_recede_index; 00264 z_advance_index = z_init_advance_index; 00265 00266 for_less( z, 0, sizes[Z] ) 00267 { 00268 #ifdef DEBUG 00269 correct_voxel = get_correct_amount( volume, x, y, z, 00270 x_width, y_width, z_width ); 00271 00272 if( !numerically_close( sum, correct_voxel, 0.001 ) ) 00273 handle_internal_error( "Dang" ); 00274 #endif 00275 00276 value = sum / total_volume; 00277 set_volume_real_value( resampled_volume, x, y, z, 0, 0, value ); 00278 00279 if( z == sizes[Z]-1 ) 00280 continue; 00281 00282 GET_NEXT_SAMPLE( z_recede_index, z_advance_index, 00283 z_left_weight, z_right_weight, 00284 sizes[Z], (Real) row[_I], sum ) 00285 00286 ++z_recede_index; 00287 ++z_advance_index; 00288 } 00289 00290 if( y == sizes[Y]-1 ) 00291 continue; 00292 00293 for_less( z, 0, sizes[Z] ) 00294 { 00295 sample = (Real) row[z]; 00296 GET_NEXT_SAMPLE( y_recede_index, y_advance_index, 00297 y_left_weight, y_right_weight, 00298 sizes[Y], (Real) slice[_I][z], sample ) 00299 row[z] = (float) sample; 00300 } 00301 00302 ++y_recede_index; 00303 ++y_advance_index; 00304 00305 update_progress_report( &progress, x * sizes[Y] + y + 1 ); 00306 } 00307 00308 if( x == sizes[X] - 1 ) 00309 continue; 00310 00311 for_less( y, 0, sizes[Y] ) 00312 { 00313 for_less( z, 0, sizes[Z] ) 00314 { 00315 sample = (Real) slice[y][z]; 00316 GET_NEXT_SAMPLE( x_recede_index, x_advance_index, 00317 x_left_weight, x_right_weight, 00318 sizes[X], volume_cache[_I][y][z], 00319 sample ) 00320 slice[y][z] = (float) sample; 00321 } 00322 } 00323 00324 ++x_recede_index; 00325 ++x_advance_index; 00326 00327 if( x_recede_index+1 >= 0 && x_recede_index+1 < sizes[X] && 00328 volume_cache[x_recede_index+1] == NULL ) 00329 { 00330 if( n_alloced < max_alloced ) 00331 { 00332 ALLOC2D( volume_cache[x_recede_index+1], sizes[Y], sizes[Z] ); 00333 ++n_alloced; 00334 } 00335 else 00336 { 00337 volume_cache[x_recede_index+1] = volume_cache[x_recede_index-1]; 00338 volume_cache[x_recede_index-1] = NULL; 00339 } 00340 00341 get_volume_value_hyperslab_3d( volume, x_recede_index+1, 0, 0, 00342 1, sizes[Y], sizes[Z], 00343 &volume_cache[x_recede_index+1][0][0] ); 00344 } 00345 00346 if( x_advance_index+1 >= 0 && x_advance_index+1 < sizes[X] && 00347 volume_cache[x_advance_index+1] == NULL ) 00348 { 00349 if( x_advance_index - 1 == x_recede_index || 00350 x_advance_index - 1 == x_recede_index+1 ) 00351 { 00352 ALLOC2D( volume_cache[x_advance_index+1], sizes[Y], sizes[Z] ); 00353 ++n_alloced; 00354 } 00355 else 00356 { 00357 volume_cache[x_advance_index+1] = 00358 volume_cache[x_advance_index-1]; 00359 volume_cache[x_advance_index-1] = NULL; 00360 } 00361 00362 get_volume_value_hyperslab_3d( volume, x_advance_index+1, 0, 0, 00363 1, sizes[Y], sizes[Z], 00364 &volume_cache[x_advance_index+1][0][0] ); 00365 } 00366 } 00367 00368 terminate_progress_report( &progress ); 00369 00370 for_less( x, 0, sizes[X] ) 00371 { 00372 if( volume_cache[x] != NULL ) 00373 { 00374 FREE2D( volume_cache[x] ); 00375 --n_alloced; 00376 } 00377 } 00378 00379 FREE( volume_cache ); 00380 00381 FREE2D( slice ); 00382 FREE( row ); 00383 00384 return( resampled_volume ); 00385 } 00386 00387 #ifdef DEBUG 00388 private void get_voxel_range( 00389 int size, 00390 Real min_pos, 00391 Real max_pos, 00392 int *min_voxel, 00393 int *max_voxel, 00394 Real *start_weight, 00395 Real *end_weight ) 00396 { 00397 if( min_pos < -0.5 ) 00398 min_pos = -0.5; 00399 if( max_pos > (Real) size - 0.5 ) 00400 max_pos = (Real) size - 0.5; 00401 00402 *min_voxel = ROUND( min_pos ); 00403 *max_voxel = ROUND( max_pos ); 00404 00405 if( (Real) (*max_voxel) - 0.5 == max_pos ) 00406 --(*max_voxel); 00407 00408 if( *min_voxel == *max_voxel ) 00409 *start_weight = max_pos - min_pos; 00410 else 00411 { 00412 *start_weight = ((Real) (*min_voxel) + 0.5) - min_pos; 00413 *end_weight = max_pos - ((Real) (*max_voxel) - 0.5); 00414 } 00415 } 00416 00417 private Real get_amount_in_box( 00418 Volume volume, 00419 int x_min_voxel, 00420 int x_max_voxel, 00421 int y_min_voxel, 00422 int y_max_voxel, 00423 int z_min_voxel, 00424 int z_max_voxel, 00425 Real x_weight_start, 00426 Real x_weight_end, 00427 Real y_weight_start, 00428 Real y_weight_end, 00429 Real z_weight_start, 00430 Real z_weight_end ) 00431 { 00432 int x, y, z; 00433 Real sum, z_sum, x_sum, value; 00434 00435 sum = 0.0; 00436 00437 for_inclusive( z, z_min_voxel, z_max_voxel ) 00438 { 00439 z_sum = 0.0; 00440 00441 for_inclusive( x, x_min_voxel, x_max_voxel ) 00442 { 00443 x_sum = 0.0; 00444 00445 for_inclusive( y, y_min_voxel, y_max_voxel ) 00446 { 00447 value = get_volume_real_value( volume, x, y, z, 0, 0 ); 00448 00449 if( y == y_min_voxel ) 00450 value *= y_weight_start; 00451 else if( y == y_max_voxel ) 00452 value *= y_weight_end; 00453 00454 x_sum += value; 00455 } 00456 00457 if( x == x_min_voxel ) 00458 x_sum *= x_weight_start; 00459 else if( x == x_max_voxel ) 00460 x_sum *= x_weight_end; 00461 00462 z_sum += x_sum; 00463 } 00464 00465 if( z == z_min_voxel ) 00466 z_sum *= z_weight_start; 00467 else if( z == z_max_voxel ) 00468 z_sum *= z_weight_end; 00469 00470 sum += z_sum; 00471 } 00472 00473 return( sum ); 00474 } 00475 00476 private Real get_correct_amount( 00477 Volume volume, 00478 int x, 00479 int y, 00480 int z, 00481 Real x_width, 00482 Real y_width, 00483 Real z_width ) 00484 { 00485 int sizes[MAX_DIMENSIONS]; 00486 Real x_start_weight, x_end_weight; 00487 Real y_start_weight, y_end_weight; 00488 Real z_start_weight, z_end_weight; 00489 Real sum; 00490 int x_min_voxel, x_max_voxel; 00491 int y_min_voxel, y_max_voxel; 00492 int z_min_voxel, z_max_voxel; 00493 00494 get_volume_sizes( volume, sizes ); 00495 00496 get_voxel_range( sizes[X], (Real) x - x_width, (Real) x + x_width, 00497 &x_min_voxel, &x_max_voxel, &x_start_weight, &x_end_weight); 00498 get_voxel_range( sizes[Y], (Real) y - y_width, (Real) y + y_width, 00499 &y_min_voxel, &y_max_voxel, &y_start_weight, &y_end_weight); 00500 get_voxel_range( sizes[Z], (Real) z - z_width, (Real) z + z_width, 00501 &z_min_voxel, &z_max_voxel, &z_start_weight, &z_end_weight); 00502 00503 sum = get_amount_in_box( volume, 00504 x_min_voxel, x_max_voxel, 00505 y_min_voxel, y_max_voxel, 00506 z_min_voxel, z_max_voxel, 00507 x_start_weight, x_end_weight, 00508 y_start_weight, y_end_weight, 00509 z_start_weight, z_end_weight ); 00510 00511 return( sum ); 00512 } 00513 #endif