BICPL: interpolate.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/interpolate.c,v 1.11 2000/02/06 15:30:55 stever Exp $"; 00020 #endif 00021 00022 /* ----------------------------- MNI Header ----------------------------------- 00023 @NAME : interpolate_volume_to_slice 00024 @INPUT : volume1 00025 n_dims1 00026 origin1 00027 x_axis1 00028 y_axis1 00029 volume2 00030 n_dims2 00031 origin2 00032 x_axis2 00033 y_axis2 00034 degrees_continuity 00035 cmode_colour_map 00036 rgb_colour_map 00037 empty_colour 00038 @OUTPUT : pixels 00039 @RETURNS : 00040 @DESCRIPTION: Interpolates the volume to the given slice, using either 00041 nearest neighbour, linear or cubic interpolation. 00042 @METHOD : 00043 @GLOBALS : 00044 @CALLS : 00045 @CREATED : 1993 David MacDonald 00046 @MODIFIED : 00047 ---------------------------------------------------------------------------- */ 00048 00049 /* ARGSUSED */ 00050 00051 public void interpolate_volume_to_slice( 00052 Volume volume1, 00053 int n_dims1, 00054 Real origin1[], 00055 Real x_axis1[], 00056 Real y_axis1[], 00057 Volume volume2, 00058 int n_dims2, 00059 Real origin2[], 00060 Real x_axis2[], 00061 Real y_axis2[], 00062 int x_pixel_start, 00063 int x_pixel_end, 00064 int y_pixel_start, 00065 int y_pixel_end, 00066 int degrees_continuity, 00067 unsigned short **cmode_colour_map, 00068 Colour **rgb_colour_map, 00069 Colour empty_colour, 00070 pixels_struct *pixels ) 00071 { 00072 int dim, x, y; 00073 int int_voxel_value1, int_voxel_value2; 00074 Real outside_value1, outside_value2; 00075 Real start_voxel1[MAX_DIMENSIONS], voxel1[MAX_DIMENSIONS]; 00076 Real start_voxel2[MAX_DIMENSIONS], voxel2[MAX_DIMENSIONS]; 00077 Real value1, voxel_value1, value2, voxel_value2; 00078 Real min_voxel1, max_voxel1, min_voxel2, max_voxel2; 00079 unsigned short *cmode_ptr; 00080 Colour *rgb_ptr; 00081 BOOLEAN inside1, inside2; 00082 Pixel_types pixel_type; 00083 #ifdef REPORT_PROGRESS 00084 progress_struct progress; 00085 #endif 00086 00087 if( is_an_rgb_volume( volume1 ) || 00088 (volume2 != NULL && is_an_rgb_volume( volume2 )) ) 00089 degrees_continuity = -1; 00090 00091 get_volume_voxel_range( volume1, &min_voxel1, &max_voxel1 ); 00092 00093 pixel_type = pixels->pixel_type; 00094 00095 for_less( dim, 0, n_dims1 ) 00096 { 00097 start_voxel1[dim] = origin1[dim] + (Real) x_pixel_start * x_axis1[dim] + 00098 (Real) y_pixel_start * y_axis1[dim]; 00099 } 00100 00101 outside_value1 = 0.0; 00102 00103 if( volume2 != NULL ) 00104 { 00105 get_volume_voxel_range( volume2, &min_voxel2, &max_voxel2 ); 00106 for_less( dim, 0, n_dims2 ) 00107 { 00108 start_voxel2[dim] = origin2[dim] + 00109 (Real) x_pixel_start * x_axis2[dim] + 00110 (Real) y_pixel_start * y_axis2[dim]; 00111 } 00112 00113 outside_value2 = 0.0; 00114 } 00115 00116 #ifdef REPORT_PROGRESS 00117 initialize_progress_report( &progress, FALSE, 00118 y_pixel_end - y_pixel_start + 1, 00119 "Creating Slice" ); 00120 #endif 00121 00122 for_inclusive( y, y_pixel_start, y_pixel_end ) 00123 { 00124 for_less( dim, 0, n_dims1 ) 00125 voxel1[dim] = start_voxel1[dim]; 00126 00127 if( volume2 != NULL ) 00128 { 00129 for_less( dim, 0, n_dims2 ) 00130 voxel2[dim] = start_voxel2[dim]; 00131 } 00132 00133 if( pixel_type == RGB_PIXEL ) 00134 rgb_ptr = &PIXEL_RGB_COLOUR(*pixels,x_pixel_start,y); 00135 else 00136 cmode_ptr = &PIXEL_COLOUR_INDEX_16(*pixels,x_pixel_start,y); 00137 00138 for_inclusive( x, x_pixel_start, x_pixel_end ) 00139 { 00140 inside1 = voxel_is_within_volume( volume1, voxel1 ); 00141 00142 if( inside1 ) 00143 { 00144 (void) evaluate_volume( volume1, voxel1, NULL, 00145 degrees_continuity, 00146 FALSE, outside_value1, 00147 &value1, NULL, NULL ); 00148 00149 voxel_value1 = convert_value_to_voxel( volume1, value1 ); 00150 00151 if( voxel_value1 < min_voxel1 ) 00152 voxel_value1 = min_voxel1; 00153 else if( voxel_value1 > max_voxel1 ) 00154 voxel_value1 = max_voxel1; 00155 } 00156 00157 for_less( dim, 0, n_dims1 ) 00158 voxel1[dim] += x_axis1[dim]; 00159 00160 if( volume2 != NULL ) 00161 { 00162 inside2 = voxel_is_within_volume( volume2, voxel2 ); 00163 00164 if( inside2 ) 00165 { 00166 (void) evaluate_volume( volume2, voxel2, NULL, 00167 degrees_continuity, 00168 FALSE, outside_value2, 00169 &value2, NULL, NULL ); 00170 00171 voxel_value2 = convert_value_to_voxel( volume2, value2 ); 00172 00173 if( voxel_value2 < min_voxel2 ) 00174 voxel_value2 = min_voxel2; 00175 else if( voxel_value2 > max_voxel2 ) 00176 voxel_value2 = max_voxel2; 00177 } 00178 00179 for_less( dim, 0, n_dims2 ) 00180 voxel2[dim] += x_axis2[dim]; 00181 00182 if( pixel_type == RGB_PIXEL ) 00183 { 00184 if( inside1 && inside2 ) 00185 { 00186 int_voxel_value1 = ROUND( voxel_value1 ); 00187 int_voxel_value2 = ROUND( voxel_value2 ); 00188 *rgb_ptr = rgb_colour_map[int_voxel_value1] 00189 [int_voxel_value2]; 00190 } 00191 else 00192 *rgb_ptr = empty_colour; 00193 00194 ++rgb_ptr; 00195 } 00196 else 00197 { 00198 if( inside1 && inside2 ) 00199 { 00200 int_voxel_value1 = ROUND( voxel_value1 ); 00201 int_voxel_value2 = ROUND( voxel_value2 ); 00202 *cmode_ptr = cmode_colour_map[int_voxel_value1] 00203 [int_voxel_value2]; 00204 } 00205 else 00206 *cmode_ptr = (unsigned short) empty_colour; 00207 00208 ++cmode_ptr; 00209 } 00210 } 00211 else 00212 { 00213 if( pixel_type == RGB_PIXEL ) 00214 { 00215 if( inside1 ) 00216 { 00217 if( rgb_colour_map == NULL ) 00218 *rgb_ptr = (Colour) voxel_value1; 00219 else 00220 { 00221 int_voxel_value1 = ROUND( voxel_value1 ); 00222 *rgb_ptr = rgb_colour_map[0][int_voxel_value1]; 00223 } 00224 } 00225 else 00226 *rgb_ptr = empty_colour; 00227 00228 ++rgb_ptr; 00229 } 00230 else 00231 { 00232 if( inside1 ) 00233 { 00234 int_voxel_value1 = ROUND( voxel_value1 ); 00235 if( cmode_colour_map == NULL ) 00236 *cmode_ptr = (unsigned short) int_voxel_value1; 00237 else 00238 { 00239 *cmode_ptr = cmode_colour_map[0][int_voxel_value1]; 00240 } 00241 } 00242 else 00243 *cmode_ptr = (unsigned short) empty_colour; 00244 00245 ++cmode_ptr; 00246 } 00247 } 00248 } 00249 00250 for_less( dim, 0, n_dims1 ) 00251 start_voxel1[dim] += y_axis1[dim]; 00252 00253 if( volume2 != NULL ) 00254 { 00255 for_less( dim, 0, n_dims2 ) 00256 start_voxel2[dim] += y_axis2[dim]; 00257 } 00258 00259 #ifdef REPORT_PROGRESS 00260 update_progress_report( &progress, y-y_pixel_start+1 ); 00261 #endif 00262 } 00263 00264 #ifdef REPORT_PROGRESS 00265 terminate_progress_report( &progress ); 00266 #endif 00267 }