BICPL: mapping.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/mapping.c,v 1.33 2000/02/06 15:30:55 stever Exp $"; 00021 #endif 00022 00023 #define DISTANCE_THRESHOLD 1.0e-10 00024 00025 /* ----------------------------- MNI Header ----------------------------------- 00026 @NAME : clip_points 00027 @INPUT : n_dims 00028 sizes 00029 origin 00030 x_axis 00031 y_axis 00032 n_points 00033 points 00034 @OUTPUT : clipped_points 00035 @RETURNS : number of points output 00036 @DESCRIPTION: Clips the points against the 2 * n_dims edges of the volume. 00037 @METHOD : 00038 @GLOBALS : 00039 @CALLS : 00040 @CREATED : 1993 David MacDonald 00041 @MODIFIED : Mar. 6, 1995 D. MacDonald -- removed complicated recursive 00042 edge clipping, replaced with iterative edge 00043 clipping. 00044 ---------------------------------------------------------------------------- */ 00045 00046 private int clip_points( 00047 int n_dims, 00048 int sizes[], 00049 Real origin[], 00050 Real x_axis[], 00051 Real y_axis[], 00052 int n_points, 00053 Real points[][2], 00054 Real clipped_points[][2] ) 00055 { 00056 int i, c, n_input_points, n_output_points, this, lim; 00057 Real prev_dist, dist, ratio, first_dist; 00058 Real input_points[2*MAX_DIMENSIONS][2]; 00059 Real output_points[2*MAX_DIMENSIONS][2]; 00060 00061 for_less( i, 0, n_points ) 00062 { 00063 input_points[i][0] = points[i][0]; 00064 input_points[i][1] = points[i][1]; 00065 } 00066 00067 n_input_points = n_points; 00068 00069 dist = 0.0; 00070 first_dist = 0.0; 00071 00072 for_less( c, 0, n_dims ) 00073 { 00074 for_less( lim, 0, 2 ) 00075 { 00076 n_output_points = 0; 00077 if( n_input_points > 0 ) 00078 { 00079 for_less( i, 0, n_input_points+1 ) 00080 { 00081 prev_dist = dist; 00082 this = i % n_input_points; 00083 00084 if( i < n_input_points ) 00085 { 00086 if( lim == 0 ) 00087 { 00088 dist = origin[c] + input_points[i][0] * x_axis[c] + 00089 input_points[i][1] * y_axis[c] + 0.5; 00090 } 00091 else 00092 { 00093 dist = (Real) sizes[c] - 0.5 - 00094 (origin[c] + input_points[i][0] * x_axis[c] + 00095 input_points[i][1] * y_axis[c]); 00096 } 00097 } 00098 else 00099 dist = first_dist; 00100 00101 if( i > 0 ) 00102 { 00103 if( (dist < -DISTANCE_THRESHOLD && 00104 prev_dist > DISTANCE_THRESHOLD) || 00105 (dist > DISTANCE_THRESHOLD && 00106 prev_dist < -DISTANCE_THRESHOLD) ) 00107 { 00108 ratio = prev_dist / (prev_dist - dist); 00109 if( ratio > DISTANCE_THRESHOLD && 00110 ratio < 1.0 - DISTANCE_THRESHOLD ) 00111 { 00112 output_points[n_output_points][0] = 00113 input_points[i-1][0] + ratio * 00114 (input_points[this][0] - input_points[i-1][0]); 00115 output_points[n_output_points][1] = 00116 input_points[i-1][1] + ratio * 00117 (input_points[this][1] - input_points[i-1][1]); 00118 ++n_output_points; 00119 } 00120 } 00121 } 00122 else 00123 first_dist = dist; 00124 00125 if( dist >= 0.0 && i != n_input_points ) 00126 { 00127 output_points[n_output_points][0] = input_points[i][0]; 00128 output_points[n_output_points][1] = input_points[i][1]; 00129 ++n_output_points; 00130 } 00131 } 00132 } 00133 00134 for_less( i, 0, n_output_points ) 00135 { 00136 input_points[i][0] = output_points[i][0]; 00137 input_points[i][1] = output_points[i][1]; 00138 } 00139 00140 n_input_points = n_output_points; 00141 } 00142 } 00143 00144 for_less( i, 0, n_output_points ) 00145 { 00146 clipped_points[i][0] = output_points[i][0]; 00147 clipped_points[i][1] = output_points[i][1]; 00148 } 00149 00150 return( n_output_points ); 00151 } 00152 00153 /* ----------------------------- MNI Header ----------------------------------- 00154 @NAME : get_cross_section 00155 @INPUT : volume 00156 origin 00157 x_axis 00158 y_axis 00159 @OUTPUT : clipped_points 00160 @RETURNS : number of clipped points 00161 @DESCRIPTION: Gets the points defining the outline of the cross section of the 00162 plane with the volume. 00163 @METHOD : 00164 @GLOBALS : 00165 @CALLS : 00166 @CREATED : 1993 David MacDonald 00167 @MODIFIED : 00168 ---------------------------------------------------------------------------- */ 00169 00170 private int get_cross_section( 00171 Volume volume, 00172 Real origin[], 00173 Real x_axis[], 00174 Real y_axis[], 00175 Real clipped_points[][2] ) 00176 { 00177 Real points[4][2], voxel[MAX_DIMENSIONS]; 00178 int d, sizes[MAX_DIMENSIONS], n_dims, n_points; 00179 int n_limits[MAX_DIMENSIONS], lim[MAX_DIMENSIONS]; 00180 Real x_pixel, y_pixel, x_min, x_max, y_min, y_max, dx, dy; 00181 BOOLEAN first; 00182 00183 get_volume_sizes( volume, sizes ); 00184 n_dims = get_volume_n_dimensions( volume ); 00185 00186 for_less( d, 0, MAX_DIMENSIONS ) 00187 { 00188 if( d < n_dims ) 00189 n_limits[d] = 2; 00190 else 00191 { 00192 n_limits[d] = 1; 00193 sizes[d] = 1; 00194 } 00195 } 00196 00197 first = TRUE; 00198 00199 for_less( lim[0], 0, n_limits[0] ) 00200 for_less( lim[1], 0, n_limits[1] ) 00201 for_less( lim[2], 0, n_limits[2] ) 00202 for_less( lim[3], 0, n_limits[3] ) 00203 for_less( lim[4], 0, n_limits[4] ) 00204 { 00205 for_less( d, 0, n_dims ) 00206 voxel[d] = -0.5 + (Real) (sizes[d] * lim[d]); 00207 00208 map_voxel_to_pixel( get_volume_n_dimensions(volume), 00209 voxel, origin, x_axis, y_axis, &x_pixel, &y_pixel ); 00210 00211 if( first ) 00212 { 00213 x_min = x_pixel; 00214 x_max = x_pixel; 00215 y_min = y_pixel; 00216 y_max = y_pixel; 00217 first = FALSE; 00218 } 00219 else 00220 { 00221 if( x_pixel < x_min ) 00222 x_min = x_pixel; 00223 else if( x_pixel > x_max ) 00224 x_max = x_pixel; 00225 if( y_pixel < y_min ) 00226 y_min = y_pixel; 00227 else if( y_pixel > y_max ) 00228 y_max = y_pixel; 00229 } 00230 } 00231 00232 dx = x_max - x_min; 00233 dy = y_max - y_min; 00234 00235 x_min -= dx; 00236 x_max += dx; 00237 y_min -= dy; 00238 y_max += dy; 00239 00240 points[0][0] = x_min; 00241 points[0][1] = y_min; 00242 points[1][0] = x_max; 00243 points[1][1] = y_min; 00244 points[2][0] = x_max; 00245 points[2][1] = y_max; 00246 points[3][0] = x_min; 00247 points[3][1] = y_max; 00248 00249 n_points = clip_points( n_dims, sizes, origin, x_axis, y_axis, 00250 4, points, clipped_points ); 00251 00252 if( n_points == 1 || n_points == 2 ) 00253 { 00254 /* 00255 print_error( "N points = %d\n", n_points ); 00256 handle_internal_error( "clipping" ); 00257 */ 00258 n_points = 0; 00259 } 00260 00261 return( n_points ); 00262 } 00263 00264 /* ----------------------------- MNI Header ----------------------------------- 00265 @NAME : get_volume_cross_section 00266 @INPUT : volume 00267 origin 00268 x_axis 00269 y_axis 00270 @OUTPUT : clipped_voxels 00271 @RETURNS : number of voxels 00272 @DESCRIPTION: Gets the cross section of the plane with the volume, in terms 00273 of a list of voxels defining the vertices of the polygon. 00274 @METHOD : 00275 @GLOBALS : 00276 @CALLS : 00277 @CREATED : 1993 David MacDonald 00278 @MODIFIED : 00279 ---------------------------------------------------------------------------- */ 00280 00281 public int get_volume_cross_section( 00282 Volume volume, 00283 Real origin[], 00284 Real x_axis[], 00285 Real y_axis[], 00286 Real clipped_voxels[][MAX_DIMENSIONS] ) 00287 { 00288 Real clipped_points[2*MAX_DIMENSIONS][2]; 00289 int i, c, n_dims, n_points; 00290 Real real_origin[MAX_DIMENSIONS]; 00291 Real real_x_axis[MAX_DIMENSIONS], real_y_axis[MAX_DIMENSIONS]; 00292 00293 get_mapping( volume, origin, x_axis, y_axis, 00294 0.0, 0.0, 1.0, 1.0, 00295 real_origin, real_x_axis, real_y_axis ); 00296 00297 n_points = get_cross_section( volume, real_origin, real_x_axis, real_y_axis, 00298 clipped_points ); 00299 00300 n_dims = get_volume_n_dimensions( volume ); 00301 00302 for_less( i, 0, n_points ) 00303 { 00304 for_less( c, 0, n_dims ) 00305 { 00306 clipped_voxels[i][c] = real_origin[c] + 00307 clipped_points[i][0] * real_x_axis[c] + 00308 clipped_points[i][1] * real_y_axis[c]; 00309 } 00310 } 00311 00312 return( n_points ); 00313 } 00314 00315 /* ----------------------------- MNI Header ----------------------------------- 00316 @NAME : get_volume_slice_range 00317 @INPUT : volume 00318 origin 00319 x_axis 00320 y_axis 00321 @OUTPUT : x_pixel_start 00322 x_pixel_end 00323 y_pixel_start 00324 y_pixel_end 00325 @RETURNS : 00326 @DESCRIPTION: Gets the range of the cross section of the volume with the 00327 plane. 00328 @METHOD : 00329 @GLOBALS : 00330 @CALLS : 00331 @CREATED : 1993 David MacDonald 00332 @MODIFIED : 00333 ---------------------------------------------------------------------------- */ 00334 00335 private void get_volume_slice_range( 00336 Volume volume, 00337 Real origin[], 00338 Real x_axis[], 00339 Real y_axis[], 00340 Real *x_pixel_start, 00341 Real *x_pixel_end, 00342 Real *y_pixel_start, 00343 Real *y_pixel_end ) 00344 { 00345 Real clipped_points[2*MAX_DIMENSIONS][2]; 00346 int i, n_points; 00347 00348 n_points = get_cross_section( volume, origin, x_axis, y_axis, 00349 clipped_points ); 00350 00351 if( n_points == 0 ) 00352 { 00353 *x_pixel_start = 1.0; 00354 *x_pixel_end = 0.0; 00355 *y_pixel_start = 1.0; 00356 *y_pixel_end = 0.0; 00357 } 00358 else 00359 { 00360 *x_pixel_start = clipped_points[0][0]; 00361 *x_pixel_end = clipped_points[0][0]; 00362 *y_pixel_start = clipped_points[0][1]; 00363 *y_pixel_end = clipped_points[0][1]; 00364 for_less( i, 0, n_points ) 00365 { 00366 if( clipped_points[i][0] < *x_pixel_start ) 00367 *x_pixel_start = clipped_points[i][0]; 00368 else if( clipped_points[i][0] > *x_pixel_end ) 00369 *x_pixel_end = clipped_points[i][0]; 00370 if( clipped_points[i][1] < *y_pixel_start ) 00371 *y_pixel_start = clipped_points[i][1]; 00372 else if( clipped_points[i][1] > *y_pixel_end ) 00373 *y_pixel_end = clipped_points[i][1]; 00374 } 00375 } 00376 } 00377 00378 /* ----------------------------- MNI Header ----------------------------------- 00379 @NAME : get_volume_mapping_range 00380 @INPUT : volume 00381 origin 00382 x_axis 00383 y_axis 00384 x_trans 00385 y_trans 00386 x_scale 00387 y_scale 00388 @OUTPUT : x_pixel_start 00389 x_pixel_end 00390 y_pixel_start 00391 y_pixel_end 00392 @RETURNS : 00393 @DESCRIPTION: Gets the pixel limits of the cross section, based on the scale 00394 and translation parameters. 00395 @METHOD : 00396 @GLOBALS : 00397 @CALLS : 00398 @CREATED : 1993 David MacDonald 00399 @MODIFIED : 00400 ---------------------------------------------------------------------------- */ 00401 00402 public void get_volume_mapping_range( 00403 Volume volume, 00404 Real origin[], 00405 Real x_axis[], 00406 Real y_axis[], 00407 Real x_trans, 00408 Real y_trans, 00409 Real x_scale, 00410 Real y_scale, 00411 Real *x_pixel_start, 00412 Real *x_pixel_end, 00413 Real *y_pixel_start, 00414 Real *y_pixel_end ) 00415 { 00416 Real real_origin[MAX_DIMENSIONS]; 00417 Real real_x_axis[MAX_DIMENSIONS], real_y_axis[MAX_DIMENSIONS]; 00418 00419 get_mapping( volume, origin, x_axis, y_axis, 00420 x_trans, y_trans, x_scale, y_scale, 00421 real_origin, real_x_axis, real_y_axis ); 00422 00423 get_volume_slice_range( volume, real_origin, real_x_axis, real_y_axis, 00424 x_pixel_start, x_pixel_end, 00425 y_pixel_start, y_pixel_end ); 00426 } 00427 00428 /* ----------------------------- MNI Header ----------------------------------- 00429 @NAME : clip_viewport_to_volume 00430 @INPUT : volume 00431 origin 00432 x_axis 00433 y_axis 00434 @OUTPUT : x_pixel_start 00435 x_pixel_end 00436 y_pixel_start 00437 y_pixel_end 00438 @RETURNS : 00439 @DESCRIPTION: Finds the pixel range of the volume cross section. 00440 @METHOD : 00441 @GLOBALS : 00442 @CALLS : 00443 @CREATED : 1993 David MacDonald 00444 @MODIFIED : 00445 ---------------------------------------------------------------------------- */ 00446 00447 public void clip_viewport_to_volume( 00448 Volume volume, 00449 Real origin[], 00450 Real x_axis[], 00451 Real y_axis[], 00452 int *x_pixel_start, 00453 int *x_pixel_end, 00454 int *y_pixel_start, 00455 int *y_pixel_end ) 00456 { 00457 int int_x_min, int_x_max, int_y_min, int_y_max; 00458 Real x_min, x_max, y_min, y_max; 00459 00460 get_volume_slice_range( volume, origin, x_axis, y_axis, 00461 &x_min, &x_max, &y_min, &y_max ); 00462 00463 int_x_min = CEILING( x_min ); 00464 int_x_max = FLOOR( x_max ); 00465 int_y_min = CEILING( y_min ); 00466 int_y_max = FLOOR( y_max ); 00467 00468 if( int_x_min > *x_pixel_start ) 00469 *x_pixel_start = int_x_min; 00470 if( int_x_max < *x_pixel_end ) 00471 *x_pixel_end = int_x_max; 00472 00473 if( int_y_min > *y_pixel_start ) 00474 *y_pixel_start = int_y_min; 00475 if( int_y_max < *y_pixel_end ) 00476 *y_pixel_end = int_y_max; 00477 } 00478 00479 00480 /* ----------------------------- MNI Header ----------------------------------- 00481 @NAME : get_mapping 00482 @INPUT : volume 00483 origin 00484 x_axis 00485 y_axis 00486 x_translation 00487 y_translation 00488 x_scale 00489 y_scale 00490 @OUTPUT : pix_origin 00491 pix_x_axis 00492 pix_y_axis 00493 @RETURNS : 00494 @DESCRIPTION: Computes the mapping from pixels to voxels. 00495 @METHOD : 00496 @GLOBALS : 00497 @CALLS : 00498 @CREATED : 1993 David MacDonald 00499 @MODIFIED : 00500 ---------------------------------------------------------------------------- */ 00501 00502 public void get_mapping( 00503 Volume volume, 00504 Real origin[], 00505 Real x_axis[], 00506 Real y_axis[], 00507 Real x_translation, 00508 Real y_translation, 00509 Real x_scale, 00510 Real y_scale, 00511 Real pix_origin[], 00512 Real pix_x_axis[], 00513 Real pix_y_axis[] ) 00514 { 00515 int c, n_dims; 00516 Real separations[MAX_DIMENSIONS], x_len, y_len, comp; 00517 00518 n_dims = get_volume_n_dimensions( volume ); 00519 get_volume_separations( volume, separations ); 00520 00521 x_len = 0.0; 00522 y_len = 0.0; 00523 for_less( c, 0, n_dims ) 00524 { 00525 comp = x_axis[c] * separations[c]; 00526 x_len += comp * comp; 00527 comp = y_axis[c] * separations[c]; 00528 y_len += comp * comp; 00529 } 00530 00531 x_len = sqrt( x_len ); 00532 if( x_len == 0.0 ) 00533 x_len = 1.0; 00534 y_len = sqrt( y_len ); 00535 if( y_len == 0.0 ) 00536 y_len = 1.0; 00537 00538 x_scale *= x_len; 00539 y_scale *= y_len; 00540 00541 for_less( c, 0, n_dims ) 00542 { 00543 pix_x_axis[c] = x_axis[c] / x_scale; 00544 pix_y_axis[c] = y_axis[c] / y_scale; 00545 pix_origin[c] = origin[c] - pix_x_axis[c] * x_translation 00546 - pix_y_axis[c] * y_translation; 00547 } 00548 } 00549 00550 /* ----------------------------- MNI Header ----------------------------------- 00551 @NAME : map_pixel_to_voxel 00552 @INPUT : x_pixel 00553 y_pixel 00554 n_dims 00555 voxel_origin 00556 voxel_x_axis 00557 voxel_y_axis 00558 @OUTPUT : voxel 00559 @RETURNS : 00560 @DESCRIPTION: Converts a pixel to a voxel. 00561 @METHOD : 00562 @GLOBALS : 00563 @CALLS : 00564 @CREATED : 1993 David MacDonald 00565 @MODIFIED : 00566 ---------------------------------------------------------------------------- */ 00567 00568 private void map_pixel_to_voxel( 00569 Real x_pixel, 00570 Real y_pixel, 00571 int n_dims, 00572 Real voxel_origin[], 00573 Real voxel_x_axis[], 00574 Real voxel_y_axis[], 00575 Real voxel[] ) 00576 { 00577 int c; 00578 00579 for_less( c, 0, n_dims ) 00580 { 00581 voxel[c] = voxel_origin[c] + (Real) x_pixel * voxel_x_axis[c] + 00582 (Real) y_pixel * voxel_y_axis[c]; 00583 } 00584 } 00585 00586 /* ----------------------------- MNI Header ----------------------------------- 00587 @NAME : dot 00588 @INPUT : n 00589 v1 00590 v2 00591 @OUTPUT : 00592 @RETURNS : Dot product 00593 @DESCRIPTION: Computes the dot product of 2 n-dimensional vectors. 00594 @METHOD : 00595 @GLOBALS : 00596 @CALLS : 00597 @CREATED : 1993 David MacDonald 00598 @MODIFIED : 00599 ---------------------------------------------------------------------------- */ 00600 00601 private Real dot( 00602 int n, 00603 Real v1[], 00604 Real v2[] ) 00605 { 00606 int i; 00607 Real d; 00608 00609 d = 0.0; 00610 for_less( i, 0, n ) 00611 d += v1[i] * v2[i]; 00612 00613 return( d ); 00614 } 00615 00616 /* ----------------------------- MNI Header ----------------------------------- 00617 @NAME : get_two_axes_coordinates 00618 @INPUT : n 00619 vector 00620 x_axis 00621 y_axis 00622 @OUTPUT : x_pos 00623 y_pos 00624 @RETURNS : 00625 @DESCRIPTION: Computes the (x,y) coordinate given a vector and an x and y 00626 axis. Finds the values such that 00627 x_pos * x_axis + y_pos * y_axis = vector 00628 @METHOD : 00629 @GLOBALS : 00630 @CALLS : 00631 @CREATED : 1993 David MacDonald 00632 @MODIFIED : 00633 ---------------------------------------------------------------------------- */ 00634 00635 private void get_two_axes_coordinates( 00636 int n, 00637 Real vector[], 00638 Real x_axis[], 00639 Real y_axis[], 00640 Real *x_pos, 00641 Real *y_pos ) 00642 { 00643 Real x_dot_x, x_dot_v, x_dot_y, y_dot_v, y_dot_y, bottom; 00644 00645 x_dot_x = dot( n, x_axis, x_axis ); 00646 x_dot_v = dot( n, x_axis, vector ); 00647 x_dot_y = dot( n, x_axis, y_axis ); 00648 y_dot_y = dot( n, y_axis, y_axis ); 00649 y_dot_v = dot( n, y_axis, vector ); 00650 00651 bottom = x_dot_x * y_dot_y - x_dot_y * x_dot_y; 00652 00653 *x_pos = (x_dot_v * y_dot_y - x_dot_y * y_dot_v) / bottom; 00654 *y_pos = (y_dot_v * x_dot_x - x_dot_y * x_dot_v) / bottom; 00655 } 00656 00657 /* ----------------------------- MNI Header ----------------------------------- 00658 @NAME : map_voxel_to_pixel 00659 @INPUT : n 00660 voxel 00661 origin 00662 x_axis 00663 y_axis 00664 @OUTPUT : x_pixel 00665 y_pixel 00666 @RETURNS : 00667 @DESCRIPTION: Converts a voxel coordinate to a pixel coordinate. 00668 @METHOD : 00669 @GLOBALS : 00670 @CALLS : 00671 @CREATED : 1993 David MacDonald 00672 @MODIFIED : 00673 ---------------------------------------------------------------------------- */ 00674 00675 public void map_voxel_to_pixel( 00676 int n, 00677 Real voxel[], 00678 Real origin[], 00679 Real x_axis[], 00680 Real y_axis[], 00681 Real *x_pixel, 00682 Real *y_pixel ) 00683 { 00684 int c; 00685 Real vector[MAX_DIMENSIONS]; 00686 00687 for_less( c, 0, n ) 00688 vector[c] = voxel[c] - origin[c]; 00689 00690 get_two_axes_coordinates( n, vector, x_axis, y_axis, 00691 x_pixel, y_pixel ); 00692 } 00693 00694 /* ----------------------------- MNI Header ----------------------------------- 00695 @NAME : convert_slice_pixel_to_voxel 00696 @INPUT : volume 00697 x_pixel 00698 y_pixel 00699 slice_position 00700 x_axis_index 00701 y_axis_index 00702 x_translation 00703 y_translation 00704 x_scale 00705 y_scale 00706 @OUTPUT : pixel_slice_position 00707 @RETURNS : 00708 @DESCRIPTION: Converts pixel from the viewport to a voxel position. The 00709 centres of pixels correspond to .0, as opposed to .5. 00710 @METHOD : 00711 @GLOBALS : 00712 @CALLS : 00713 @CREATED : 1993 David MacDonald 00714 @MODIFIED : 00715 ---------------------------------------------------------------------------- */ 00716 00717 public BOOLEAN convert_slice_pixel_to_voxel( 00718 Volume volume, 00719 Real x_pixel, 00720 Real y_pixel, 00721 Real origin[], 00722 Real x_axis[], 00723 Real y_axis[], 00724 Real x_translation, 00725 Real y_translation, 00726 Real x_scale, 00727 Real y_scale, 00728 Real voxel[] ) 00729 { 00730 int n_dims; 00731 Real real_x_axis[MAX_DIMENSIONS], real_y_axis[MAX_DIMENSIONS]; 00732 Real real_origin[MAX_DIMENSIONS]; 00733 00734 get_mapping( volume, origin, x_axis, y_axis, 00735 x_translation, y_translation, x_scale, y_scale, 00736 real_origin, real_x_axis, real_y_axis ); 00737 00738 n_dims = get_volume_n_dimensions( volume ); 00739 00740 map_pixel_to_voxel( x_pixel, y_pixel, n_dims, real_origin, 00741 real_x_axis, real_y_axis, voxel ); 00742 00743 return( voxel_is_within_volume( volume, voxel ) ); 00744 } 00745 00746 /* ----------------------------- MNI Header ----------------------------------- 00747 @NAME : convert_voxel_to_slice_pixel 00748 @INPUT : volume 00749 voxel_position 00750 x_axis_index 00751 y_axis_index 00752 x_translation 00753 y_translation 00754 x_scale 00755 y_scale 00756 @OUTPUT : x_pixel 00757 y_pixel 00758 @RETURNS : 00759 @DESCRIPTION: Converts the voxel position (centres of voxels are at .0, 00760 not .5) to a pixel position in the viewport. To get the 00761 integer pixel position, take ROUND( x_pixel ). 00762 @METHOD : 00763 @GLOBALS : 00764 @CALLS : 00765 @CREATED : 1993 David MacDonald 00766 @MODIFIED : 00767 ---------------------------------------------------------------------------- */ 00768 00769 public void convert_voxel_to_slice_pixel( 00770 Volume volume, 00771 Real voxel[], 00772 Real origin[], 00773 Real x_axis[], 00774 Real y_axis[], 00775 Real x_translation, 00776 Real y_translation, 00777 Real x_scale, 00778 Real y_scale, 00779 Real *x_pixel, 00780 Real *y_pixel ) 00781 { 00782 int n_dims; 00783 Real real_x_axis[MAX_DIMENSIONS], real_y_axis[MAX_DIMENSIONS]; 00784 Real real_origin[MAX_DIMENSIONS]; 00785 00786 n_dims = get_volume_n_dimensions( volume ); 00787 00788 get_mapping( volume, origin, x_axis, y_axis, 00789 x_translation, y_translation, x_scale, y_scale, 00790 real_origin, real_x_axis, real_y_axis ); 00791 00792 map_voxel_to_pixel( n_dims, voxel, real_origin, 00793 real_x_axis, real_y_axis, x_pixel, y_pixel ); 00794 } 00795 00796 /* ----------------------------- MNI Header ----------------------------------- 00797 @NAME : resize_volume_slice 00798 @INPUT : old_x_viewport_size 00799 old_y_viewport_size 00800 old_used_x_viewport_size 00801 old_used_y_viewport_size 00802 new_x_viewport_size 00803 new_y_viewport_size 00804 @MODIFIED : x_translation 00805 y_translation 00806 x_scale 00807 y_scale 00808 @OUTPUT : used_x_viewport_size 00809 used_y_viewport_size 00810 @RETURNS : 00811 @DESCRIPTION: Computes the new transformation required after resizing the 00812 viewport containing a slice. 00813 @METHOD : 00814 @GLOBALS : 00815 @CALLS : 00816 @CREATED : 1993 David MacDonald 00817 @MODIFIED : 00818 ---------------------------------------------------------------------------- */ 00819 00820 public void resize_volume_slice( 00821 int old_x_viewport_size, 00822 int old_y_viewport_size, 00823 int old_used_x_viewport_size, 00824 int old_used_y_viewport_size, 00825 int new_x_viewport_size, 00826 int new_y_viewport_size, 00827 Real *x_translation, 00828 Real *y_translation, 00829 Real *x_scale, 00830 Real *y_scale, 00831 int *used_x_viewport_size, 00832 int *used_y_viewport_size ) 00833 { 00834 Real x_scale_factor, y_scale_factor, scale_factor; 00835 00836 if( old_used_x_viewport_size <= 0 ) 00837 old_used_x_viewport_size = 1; 00838 if( old_used_y_viewport_size <= 0 ) 00839 old_used_y_viewport_size = 1; 00840 00841 x_scale_factor = (Real) new_x_viewport_size / 00842 (Real) old_used_x_viewport_size; 00843 y_scale_factor = (Real) new_y_viewport_size / 00844 (Real) old_used_y_viewport_size; 00845 scale_factor = MIN( x_scale_factor, y_scale_factor ); 00846 00847 if( used_x_viewport_size != NULL ) 00848 { 00849 *used_x_viewport_size = ROUND( scale_factor * 00850 (Real) old_used_x_viewport_size); 00851 } 00852 if( used_y_viewport_size != NULL ) 00853 { 00854 *used_y_viewport_size = ROUND( scale_factor * 00855 (Real) old_used_y_viewport_size ); 00856 } 00857 00858 scale_slice_about_viewport_centre( scale_factor, 00859 old_x_viewport_size, old_y_viewport_size, 00860 x_translation, y_translation, 00861 x_scale, y_scale ); 00862 00863 *x_translation += (Real) (new_x_viewport_size - old_x_viewport_size) / 2.0; 00864 *y_translation += (Real) (new_y_viewport_size - old_y_viewport_size) / 2.0; 00865 } 00866 00867 /* ----------------------------- MNI Header ----------------------------------- 00868 @NAME : fit_volume_slice_to_viewport 00869 @INPUT : volume 00870 origin 00871 x_axis 00872 y_axis 00873 x_viewport_size 00874 y_viewport_size 00875 fraction_oversize 00876 @OUTPUT : x_translation 00877 y_translation 00878 x_scale 00879 y_scale 00880 used_x_viewport_size 00881 used_y_viewport_size 00882 @RETURNS : 00883 @DESCRIPTION: Computes the transformation required to fit a slice to the 00884 viewport. 00885 @METHOD : 00886 @GLOBALS : 00887 @CALLS : 00888 @CREATED : 1993 David MacDonald 00889 @MODIFIED : 00890 ---------------------------------------------------------------------------- */ 00891 00892 public void fit_volume_slice_to_viewport( 00893 Volume volume, 00894 Real origin[], 00895 Real x_axis[], 00896 Real y_axis[], 00897 int x_viewport_size, 00898 int y_viewport_size, 00899 Real fraction_oversize, 00900 Real *x_translation, 00901 Real *y_translation, 00902 Real *x_scale, 00903 Real *y_scale, 00904 int *used_x_viewport_size, 00905 int *used_y_viewport_size ) 00906 { 00907 Real x_min, x_max, y_min, y_max; 00908 00909 get_volume_mapping_range( volume, origin, x_axis, y_axis, 00910 0.0, 0.0, 1.0, 1.0, 00911 &x_min, &x_max, &y_min, &y_max ); 00912 00913 if( x_min == x_max || y_min == y_max ) 00914 { 00915 *x_translation = 0.0; 00916 *y_translation = 0.0; 00917 *x_scale = 0.0; 00918 *y_scale = 0.0; 00919 return; 00920 } 00921 00922 *x_scale = (Real) x_viewport_size / (x_max - x_min) / 00923 (1.0 + fraction_oversize); 00924 *y_scale = (Real) y_viewport_size / (y_max - y_min) / 00925 (1.0 + fraction_oversize); 00926 00927 if( *x_scale < *y_scale ) 00928 *y_scale = *x_scale; 00929 else 00930 *x_scale = *y_scale; 00931 00932 if( used_x_viewport_size != NULL ) 00933 { 00934 *used_x_viewport_size = (int) (*x_scale * (Real) (x_max - x_min) * 00935 (1.0 + fraction_oversize)); 00936 } 00937 00938 if( used_y_viewport_size != NULL ) 00939 { 00940 *used_y_viewport_size = (int) (*y_scale * (Real) (y_max - y_min) * 00941 (1.0 + fraction_oversize)); 00942 } 00943 00944 *x_translation = ((Real) x_viewport_size - *x_scale * (x_max - x_min))/2.0 - 00945 *x_scale * x_min; 00946 *y_translation = ((Real) y_viewport_size - *y_scale * (y_max - y_min))/2.0 - 00947 *y_scale * y_min; 00948 } 00949 00950 /* ----------------------------- MNI Header ----------------------------------- 00951 @NAME : scale_slice_about_viewport_centre 00952 @INPUT : scale_factor 00953 x_viewport_size 00954 y_viewport_size 00955 @OUTPUT : x_translation 00956 y_translation 00957 x_scale 00958 y_scale 00959 @RETURNS : 00960 @DESCRIPTION: Modifies the scale and translation to scale a slice about 00961 the centre pixel of a viewport. 00962 @METHOD : 00963 @GLOBALS : 00964 @CALLS : 00965 @CREATED : 1993 David MacDonald 00966 @MODIFIED : 00967 ---------------------------------------------------------------------------- */ 00968 00969 public void scale_slice_about_viewport_centre( 00970 Real scale_factor, 00971 int x_viewport_size, 00972 int y_viewport_size, 00973 Real *x_translation, 00974 Real *y_translation, 00975 Real *x_scale, 00976 Real *y_scale ) 00977 { 00978 Real x_centre, y_centre; 00979 00980 x_centre = (Real) x_viewport_size / 2.0; 00981 y_centre = (Real) y_viewport_size / 2.0; 00982 00983 *x_translation = x_centre - scale_factor * (x_centre - *x_translation); 00984 *y_translation = y_centre - scale_factor * (y_centre - *y_translation); 00985 *x_scale *= scale_factor; 00986 *y_scale *= scale_factor; 00987 } 00988 00989 /* ----------------------------- MNI Header ----------------------------------- 00990 @NAME : voxel_is_within_volume 00991 @INPUT : volume 00992 voxel_position 00993 @OUTPUT : 00994 @RETURNS : TRUE if voxel is within volume. 00995 @DESCRIPTION: Determines if a voxel position is within the volume. 00996 @CREATED : Mar 1993 David MacDonald 00997 @MODIFIED : 00998 ---------------------------------------------------------------------------- */ 00999 01000 public BOOLEAN voxel_is_within_volume( 01001 Volume volume, 01002 Real voxel_position[] ) 01003 { 01004 int i, sizes[MAX_DIMENSIONS]; 01005 BOOLEAN inside; 01006 01007 inside = TRUE; 01008 01009 get_volume_sizes( volume, sizes ); 01010 01011 for_less( i, 0, get_volume_n_dimensions(volume) ) 01012 { 01013 if( voxel_position[i] < -0.5 || 01014 voxel_position[i] >= (Real) sizes[i] - 0.5 ) 01015 { 01016 inside = FALSE; 01017 break; 01018 } 01019 } 01020 01021 return( inside ); 01022 } 01023 01024 /* ----------------------------- MNI Header ----------------------------------- 01025 @NAME : int_voxel_is_within_volume 01026 @INPUT : volume 01027 indices 01028 @OUTPUT : 01029 @RETURNS : TRUE if voxel within volume 01030 @DESCRIPTION: Determines if the voxel with integer coordinates is within the 01031 volume. 01032 @CREATED : Mar 1993 David MacDonald 01033 @MODIFIED : 01034 ---------------------------------------------------------------------------- */ 01035 01036 public BOOLEAN int_voxel_is_within_volume( 01037 Volume volume, 01038 int indices[] ) 01039 { 01040 int i, sizes[MAX_DIMENSIONS]; 01041 BOOLEAN inside; 01042 01043 inside = TRUE; 01044 01045 get_volume_sizes( volume, sizes ); 01046 01047 for_less( i, 0, get_volume_n_dimensions(volume) ) 01048 { 01049 if( indices[i] < 0 || indices[i] >= sizes[i] ) 01050 { 01051 inside = FALSE; 01052 break; 01053 } 01054 } 01055 01056 return( inside ); 01057 } 01058 01059 /* ----------------------------- MNI Header ----------------------------------- 01060 @NAME : convert_real_to_int_voxel 01061 @INPUT : n_dimensions 01062 voxel 01063 @OUTPUT : int_voxel 01064 @RETURNS : 01065 @DESCRIPTION: Converts real valued voxel positions to integer positions, by 01066 rounding. 01067 @METHOD : 01068 @GLOBALS : 01069 @CALLS : 01070 @CREATED : 1993 David MacDonald 01071 @MODIFIED : 01072 ---------------------------------------------------------------------------- */ 01073 01074 public void convert_real_to_int_voxel( 01075 int n_dimensions, 01076 Real voxel[], 01077 int int_voxel[] ) 01078 { 01079 int i; 01080 01081 for_less( i, 0, n_dimensions ) 01082 int_voxel[i] = ROUND( voxel[i] ); 01083 } 01084 01085 /* ----------------------------- MNI Header ----------------------------------- 01086 @NAME : convert_int_to_real_voxel 01087 @INPUT : n_dimensions 01088 int_voxel 01089 @OUTPUT : voxel 01090 @RETURNS : 01091 @DESCRIPTION: Converts int valued voxel positions to real positions. 01092 @METHOD : 01093 @GLOBALS : 01094 @CALLS : 01095 @CREATED : Apr. 16, 1996 David MacDonald 01096 @MODIFIED : 01097 ---------------------------------------------------------------------------- */ 01098 01099 public void convert_int_to_real_voxel( 01100 int n_dimensions, 01101 int int_voxel[], 01102 Real voxel[] ) 01103 { 01104 int i; 01105 01106 for_less( i, 0, n_dimensions ) 01107 voxel[i] = (Real) int_voxel[i]; 01108 } 01109 01110 /* ----------------------------- MNI Header ----------------------------------- 01111 @NAME : voxel_contains_range 01112 @INPUT : volume 01113 int_voxel 01114 target_value 01115 @OUTPUT : 01116 @RETURNS : TRUE if voxel contains this value 01117 @DESCRIPTION: Determines if the voxel contains this value, by assuming 01118 trilinear interpolation between the 8 corner values of this 01119 voxel. 01120 @CREATED : Mar 1993 David MacDonald 01121 @MODIFIED : 01122 ---------------------------------------------------------------------------- */ 01123 01124 public BOOLEAN voxel_contains_range( 01125 Volume volume, 01126 int int_voxel[], 01127 Real min_value, 01128 Real max_value ) 01129 { 01130 BOOLEAN less, greater; 01131 int i, n_values; 01132 Real value, values[1 << MAX_DIMENSIONS]; 01133 01134 less = FALSE; 01135 greater = FALSE; 01136 01137 get_volume_value_hyperslab( volume, 01138 int_voxel[0], 01139 int_voxel[1], 01140 int_voxel[2], 01141 int_voxel[3], 01142 int_voxel[4], 01143 2, 2, 2, 2, 2, values ); 01144 01145 n_values = 1 << get_volume_n_dimensions( volume ); 01146 01147 for_less( i, 0, n_values ) 01148 { 01149 value = values[i]; 01150 01151 if( value < min_value ) 01152 { 01153 if( greater ) 01154 return( TRUE ); 01155 less = TRUE; 01156 } 01157 else if( value > max_value ) 01158 { 01159 if( less ) 01160 return( TRUE ); 01161 greater = TRUE; 01162 } 01163 else 01164 return( TRUE ); 01165 } 01166 01167 return( FALSE ); 01168 } 01169 01170 /* ----------------------------- MNI Header ----------------------------------- 01171 @NAME : volumes_are_same_grid 01172 @INPUT : volume1 01173 volume2 01174 @OUTPUT : 01175 @RETURNS : TRUE or FALSE 01176 @DESCRIPTION: Tests if the volumes represent the same tessellation, i.e., 01177 the same number of voxels, and the same voxel to world transform. 01178 @METHOD : 01179 @GLOBALS : 01180 @CALLS : 01181 @CREATED : Aug. 1, 1995 David MacDonald 01182 @MODIFIED : 01183 ---------------------------------------------------------------------------- */ 01184 01185 public BOOLEAN volumes_are_same_grid( 01186 Volume volume1, 01187 Volume volume2 ) 01188 { 01189 int c, n_dims, i, j; 01190 int sizes1[MAX_DIMENSIONS], sizes2[MAX_DIMENSIONS]; 01191 General_transform *transform1, *transform2; 01192 Transform *lin_transform1, *lin_transform2; 01193 01194 n_dims = get_volume_n_dimensions( volume1 ); 01195 if( n_dims != get_volume_n_dimensions( volume2 ) ) 01196 return( FALSE ); 01197 01198 get_volume_sizes( volume1, sizes1 ); 01199 get_volume_sizes( volume2, sizes2 ); 01200 01201 for_less( c, 0, n_dims ) 01202 { 01203 if( sizes1[c] != sizes2[c] ) 01204 return( FALSE ); 01205 } 01206 01207 transform1 = get_voxel_to_world_transform( volume1 ); 01208 transform2 = get_voxel_to_world_transform( volume2 ); 01209 01210 if( get_transform_type(transform1) != get_transform_type(transform2) ) 01211 return( FALSE ); 01212 01213 /*--- for now, we will not support checking non-linear transforms for 01214 equality */ 01215 01216 if( get_transform_type(transform1) != LINEAR || 01217 get_transform_type(transform2) != LINEAR ) 01218 return( FALSE ); 01219 01220 lin_transform1 = get_linear_transform_ptr( transform1 ); 01221 lin_transform2 = get_linear_transform_ptr( transform2 ); 01222 01223 for_less( i, 0, 4 ) 01224 { 01225 for_less( j, 0, 4 ) 01226 { 01227 if( Transform_elem(*lin_transform1,i,j) != 01228 Transform_elem(*lin_transform2,i,j) ) 01229 { 01230 return( FALSE ); 01231 } 01232 } 01233 } 01234 01235 return( TRUE ); 01236 }