1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** NOT 7D SAFE ***/
10
11 /*** functions to "warp" images -- not very efficient, but quite general ***/
12
13 #ifdef HP
14 # undef INLINE
15 # pragma INLINE xxMRI_scaler
16 # pragma INLINE xxMRI_rotfunc
17 #endif
18
19 #define FINS(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
20 ? 0.0 : far[(i)+(j)*nx] )
21
22 /**************************************************************************/
23
24 static float sx_scale , sy_scale ; /* global scaler data */
25
26 /* inline func with static vars should be static 24 Sep 2015 [rickr] */
xxMRI_scaler(float xpr,float ypr,float * xx,float * yy)27 static INLINE void xxMRI_scaler( float xpr, float ypr, float *xx , float *yy )
28 {
29 *xx = sx_scale * xpr ;
30 *yy = sy_scale * ypr ;
31 return ;
32 }
33
34 /**************************************************************************/
35
mri_warp(MRI_IMAGE * im,int nxnew,int nynew,int wtype,void wfunc (float,float,float *,float *))36 MRI_IMAGE *mri_warp( MRI_IMAGE *im , int nxnew , int nynew , int wtype ,
37 void wfunc(float,float,float *,float *) )
38 {
39 switch( wtype ){
40 case MRI_BILINEAR:
41 return mri_warp_bilinear( im , nxnew , nynew , wfunc ) ;
42
43 case MRI_BICUBIC:
44 return mri_warp_bicubic( im , nxnew , nynew , wfunc ) ;
45
46 default:
47 fprintf( stderr , "mri_warp: illegal wtype %d\n" , wtype ) ;
48 MRI_FATAL_ERROR ;
49 }
50 return NULL ;
51 }
52
53 /****************************************************************************/
54
mri_resize_NN(MRI_IMAGE * im,int nxnew,int nynew)55 MRI_IMAGE *mri_resize_NN( MRI_IMAGE *im , int nxnew , int nynew ) /* 08 Jun 2004 */
56 {
57 int nx,ny , nnx,nny , ii,jj , pp,qq , bb ;
58 float fx,fy ;
59 MRI_IMAGE *nim ;
60 char *nar , *ar ;
61
62 if( im == NULL ) return NULL ;
63
64 nx = im->nx ; ny = im->ny ;
65 nnx = nxnew ; nny = nynew ;
66 fx = ((float)nx) / (float)nnx ;
67 fy = ((float)ny) / (float)nny ;
68
69 nim = mri_new( nnx , nny , im->kind ) ;
70 nar = mri_data_pointer( nim ) ;
71 ar = mri_data_pointer( im ) ;
72 bb = im->pixel_size ;
73
74 for( jj=0 ; jj < nny ; jj++ ){
75 qq = (int)( fy*jj ) ;
76 for( ii=0 ; ii < nnx ; ii++ ){
77 pp = (int)( fx*ii ) ;
78 memcpy( nar + (ii+jj*nnx)*bb , ar + (pp+qq*nx)*bb , bb ) ;
79 }
80 }
81
82 MRI_COPY_AUX(nim,im) ;
83 nim->dx *= fx ;
84 nim->dy *= fy ;
85 return nim ;
86 }
87
88 /****************************************************************************/
89
mri_squareaspect(MRI_IMAGE * im)90 MRI_IMAGE *mri_squareaspect( MRI_IMAGE *im ) /* 08 Jun 2004 */
91 {
92 int nx,ny , nxnew,nynew ;
93 float dx,dy , rr ;
94
95 if( im == NULL ) return NULL ;
96
97 dx = fabs(im->dx) ; dy = fabs(im->dy) ;
98 if( dx == 0.0 || dy == 0.0 ) return NULL ;
99 rr = dy / dx ; if( rr == 1.0 ) return NULL ;
100
101 nx = im->nx ; ny = im->ny ;
102
103 if( rr < 1.0 ){
104 nxnew = rint( nx/rr ) ; if( nxnew <= nx ) return NULL ;
105 nynew = ny ;
106 } else {
107 nynew = rint( ny*rr ) ; if( nynew <= ny ) return NULL ;
108 nxnew = nx ;
109 }
110
111 return mri_resize_NN( im , nxnew , nynew ) ;
112 }
113
114 /****************************************************************************/
115
mri_resize(MRI_IMAGE * im,int nxnew,int nynew)116 MRI_IMAGE *mri_resize( MRI_IMAGE *im , int nxnew , int nynew )
117 {
118 int nx,ny , nnx,nny , wtype ;
119
120 nx = im->nx ; ny = im->ny ;
121 nnx = nxnew ; nny = nynew ;
122
123 if( nnx <= 0 && nny <= 0 ){
124 fprintf( stderr , "mri_resize: nxnew,nynew = %d %d\n",nxnew,nynew ) ;
125 MRI_FATAL_ERROR ;
126 }
127
128 sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ;
129 sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ;
130
131 if( nnx <= 0 ){
132 sx_scale = sy_scale ;
133 nnx = sx_scale * nx ;
134 } else if( nny <= 0 ){
135 sy_scale = sx_scale ;
136 nny = sy_scale * ny ;
137 }
138
139
140 #if 0
141 wtype = MRI_BICUBIC ;
142
143 if( ( ((nnx>=nx) && (nnx%nx==0)) || ((nnx<nx) && (nx%nnx==0)) ) &&
144 ( ((nny>=ny) && (nny%ny==0)) || ((nny<ny) && (ny%nny==0)) ) ){
145
146 wtype = MRI_BILINEAR ;
147 } else {
148 wtype = MRI_BICUBIC ;
149 }
150
151 return mri_warp( im , nnx,nny , wtype , xxMRI_scaler ) ;
152 #else
153 return mri_warp_bicubic( im , nnx,nny , xxMRI_scaler ) ;
154 #endif
155 }
156
157 /**************************************************************************/
158
mri_warp_bicubic(MRI_IMAGE * im,int nxnew,int nynew,void wf (float,float,float *,float *))159 MRI_IMAGE *mri_warp_bicubic( MRI_IMAGE *im , int nxnew , int nynew ,
160 void wf( float,float,float *,float *) )
161 {
162 MRI_IMAGE *imfl , *newImg ;
163 float *far , *nar ;
164 float xpr,ypr , xx,yy , fx,fy ;
165 int ii,jj, nx,ny , ix,jy ;
166 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
167 float bot,top,val ; /* 29 Mar 2003 */
168
169 nx = im->nx ; /* input image dimensions, for convenience */
170 ny = im->ny ;
171
172 nxnew = (nxnew > 0) ? nxnew : nx ; /* default output image sizes */
173 nynew = (nynew > 0) ? nynew : ny ;
174
175 switch( im->kind ){ /* 29 Mar 2003: allow for different input types */
176 /* by doing components 1 at a time */
177 case MRI_float:
178 imfl = im ; break ;
179
180 default:
181 imfl = mri_to_float(im) ; break ;
182
183 case MRI_short:{
184 imfl = mri_to_float(im) ;
185 newImg = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
186 mri_free(imfl) ;
187 imfl = mri_to_mri(MRI_short,newImg) ;
188 mri_free(newImg) ; return imfl ;
189 }
190
191 case MRI_byte:{
192 imfl = mri_to_float(im) ;
193 newImg = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
194 mri_free(imfl) ;
195 imfl = mri_to_mri(MRI_byte,newImg) ;
196 mri_free(newImg) ; return imfl ;
197 }
198
199 case MRI_rgb:{
200 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
201 MRI_IMAGE *rim,*gim,*bim ;
202 rim = mri_warp_bicubic( IMARR_SUBIM(imar,0), nxnew,nynew, wf ) ;
203 gim = mri_warp_bicubic( IMARR_SUBIM(imar,1), nxnew,nynew, wf ) ;
204 bim = mri_warp_bicubic( IMARR_SUBIM(imar,2), nxnew,nynew, wf ) ;
205 DESTROY_IMARR(imar) ;
206 newImg = mri_3to_rgb( rim,gim,bim ) ;
207 mri_free(rim); mri_free(gim); mri_free(bim); return newImg;
208 }
209
210 }
211
212 /* at this point, imfl is in MRI_float format */
213
214 far = mri_data_pointer( imfl ) ; /* easy access to float data */
215
216 newImg = mri_new( nxnew , nynew , MRI_float ) ; /* output image */
217 nar = mri_data_pointer( newImg ) ; /* output image data */
218
219 bot = top = far[0] ; /* 29 Mar 2003: */
220 for( ii=1 ; ii < imfl->nvox ; ii++ ){ /* clip output data range */
221 if( far[ii] > top ) top = far[ii] ;
222 else if( far[ii] < bot ) bot = far[ii] ;
223 }
224
225 /*** loop over output points and warp to them ***/
226
227 for( jj=0 ; jj < nynew ; jj++ ){
228 ypr = jj ;
229 for( ii=0 ; ii < nxnew ; ii++ ){
230 xpr = ii ;
231 wf( xpr,ypr , &xx,&yy ) ; /* get xx,yy in original image */
232
233 ix = floor( xx ) ; fx = xx - ix ; /* integer and */
234 jy = floor( yy ) ; fy = yy - jy ; /* fractional coords */
235
236 /* define cubic interpolation polynomials and data from original grid */
237
238 #define P_M1(x) (-(x)*((x)-1)*((x)-2))
239 #define P_00(x) (3*((x)+1)*((x)-1)*((x)-2))
240 #define P_P1(x) (-3*(x)*((x)+1)*((x)-2))
241 #define P_P2(x) ((x)*((x)+1)*((x)-1))
242
243 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ; /* weights for interpolating */
244 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ; /* in the x-direction */
245
246 f_jm1 = wt_m1 * FINS(ix-1,jy-1) /* interpolate to ix + fx */
247 + wt_00 * FINS(ix ,jy-1) /* at levels jy-1 .. jy+2 */
248 + wt_p1 * FINS(ix+1,jy-1)
249 + wt_p2 * FINS(ix+2,jy-1) ;
250
251 f_j00 = wt_m1 * FINS(ix-1,jy)
252 + wt_00 * FINS(ix ,jy)
253 + wt_p1 * FINS(ix+1,jy)
254 + wt_p2 * FINS(ix+2,jy) ;
255
256 f_jp1 = wt_m1 * FINS(ix-1,jy+1)
257 + wt_00 * FINS(ix ,jy+1)
258 + wt_p1 * FINS(ix+1,jy+1)
259 + wt_p2 * FINS(ix+2,jy+1) ;
260
261 f_jp2 = wt_m1 * FINS(ix-1,jy+2)
262 + wt_00 * FINS(ix ,jy+2)
263 + wt_p1 * FINS(ix+1,jy+2)
264 + wt_p2 * FINS(ix+2,jy+2) ;
265
266 /* interpolate between y-levels to jy+fy */
267
268 val = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
269 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
270
271 if( val > top ) val = top ; /* 29 Mar 2003 */
272 else if( val < bot ) val = bot ;
273
274 nar[ii+jj*nxnew] = val ;
275 }
276 }
277
278 /*** cleanup and return ***/
279
280 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
281 return newImg ;
282 }
283
284 /*************************************************************************/
285
mri_warp_bilinear(MRI_IMAGE * im,int nxnew,int nynew,void wf (float,float,float *,float *))286 MRI_IMAGE *mri_warp_bilinear( MRI_IMAGE *im , int nxnew , int nynew ,
287 void wf( float,float,float *,float *) )
288 {
289 MRI_IMAGE *imfl , *newImg ;
290 float *far , *nar ;
291 float xpr,ypr , xx,yy , fx,fx1,fy,fy1 , f00,f10,f01,f11 ;
292 int ii,jj, nx,ny , ix,jy ;
293
294 nx = im->nx ; /* dimensions, for convenience */
295 ny = im->ny ;
296
297 nxnew = (nxnew > 0) ? nxnew : nx ; /* default output image sizes */
298 nynew = (nynew > 0) ? nynew : ny ;
299
300 if( im->kind == MRI_float ){ /* convert input to float, if needed */
301 imfl = im ;
302 } else {
303 imfl = mri_to_float( im ) ;
304 }
305 far = mri_data_pointer( imfl ) ; /* easy access to float data */
306
307 newImg = mri_new( nxnew , nynew , MRI_float ) ; /* output image */
308 nar = mri_data_pointer( newImg ) ; /* output image data */
309
310 /*** loop over output points and warp to them ***/
311
312 for( jj=0 ; jj < nynew ; jj++ ){
313 ypr = jj ;
314 for( ii=0 ; ii < nxnew ; ii++ ){
315 xpr = ii ;
316 wf( xpr,ypr , &xx,&yy ) ; /* get xx,yy in original image */
317
318 ix = floor( xx ) ; fx = xx - ix ; fx1 = 1.0 - fx ;
319 jy = floor( yy ) ; fy = yy - jy ; fy1 = 1.0 - fy ;
320
321 f00 = FINS(ix ,jy ) ;
322 f01 = FINS(ix+1,jy ) ;
323 f10 = FINS(ix ,jy+1) ;
324 f11 = FINS(ix+1,jy+1) ;
325
326 nar[ii+jj*nxnew] = fx1 * ( fy1 * f00 + fy * f01 )
327 +fx * ( fy1 * f10 + fy * f11 ) ;
328
329 }
330 }
331
332 /*** cleanup and return ***/
333
334 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
335 return newImg ;
336 }
337
338 /**********************************************************************/
339
340 static float rot_dx , rot_dy , rot_cph , rot_sph ; /* global rotfunc data */
341
xxMRI_rotfunc(float xpr,float ypr,float * xx,float * yy)342 static INLINE void xxMRI_rotfunc(float xpr , float ypr , float *xx , float *yy)
343 {
344 *xx = rot_cph * xpr + rot_sph * ypr + rot_dx ;
345 *yy = -rot_sph * xpr + rot_cph * ypr + rot_dy ;
346 return ;
347 }
348
349 /**--------------------------------------------------------------------------------
350 Rotate and shift an image, using bicubic interpolation:
351 aa = shift in x
352 bb = shift in y
353 phi = angle in radians
354 scl = size scale factor (0.0 --> leave the same size)
355 -----------------------------------------------------------------------------------**/
356
mri_rotate(MRI_IMAGE * im,float aa,float bb,float phi,float scl)357 MRI_IMAGE *mri_rotate( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
358 {
359 MRI_IMAGE *imwarp ;
360 int nxnew , nynew ;
361
362 rot_cph = cos(phi) ;
363 rot_sph = sin(phi) ;
364
365 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
366 -(0.5 * im->ny) * rot_sph ;
367
368 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
369 +(0.5 * im->ny) * (1.0-rot_cph) ;
370
371 if( scl <= 0.0 ){
372 nxnew = nynew = 0 ;
373 } else {
374 nxnew = scl * im->nx + 0.49 ;
375 nynew = scl * im->ny + 0.49 ;
376 rot_cph /= scl ;
377 rot_sph /= scl ;
378 }
379
380 return mri_warp_bicubic( im , nxnew,nynew , xxMRI_rotfunc ) ;
381 }
382
mri_rotate_bilinear(MRI_IMAGE * im,float aa,float bb,float phi,float scl)383 MRI_IMAGE *mri_rotate_bilinear( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
384 {
385 MRI_IMAGE *imwarp ;
386 int nxnew , nynew ;
387
388 rot_cph = cos(phi) ;
389 rot_sph = sin(phi) ;
390
391 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
392 -(0.5 * im->ny) * rot_sph ;
393
394 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
395 +(0.5 * im->ny) * (1.0-rot_cph) ;
396
397 if( scl <= 0.0 ){
398 nxnew = nynew = 0 ;
399 } else {
400 nxnew = scl * im->nx + 0.49 ;
401 nynew = scl * im->ny + 0.49 ;
402 rot_cph /= scl ;
403 rot_sph /= scl ;
404 }
405
406 return mri_warp_bilinear( im , nxnew,nynew , xxMRI_rotfunc ) ;
407 }
408
409 #undef WARP_POINT_ROUTINES
410 #ifdef WARP_POINT_ROUTINES
411
412 /*--------------------------------------------------------------------
413 Return one point from a warped image
414 ----------------------------------------------------------------------*/
415
416 #define INSIDE(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) ? 0.0 : \
417 (kk==MRI_byte) ? (float) MRI_BYTE_PIXEL(im,i,j) : \
418 (kk==MRI_short) ? (float) MRI_SHORT_PIXEL(im,i,j) : \
419 (kk==MRI_int) ? (float) MRI_INT_PIXEL(im,i,j) : \
420 (kk==MRI_float) ? (float) MRI_FLOAT_PIXEL(im,i,j) : \
421 (kk==MRI_double) ? (float) MRI_DOUBLE_PIXEL(im,i,j) : \
422 (kk==MRI_complex) ? (float) CABS(MRI_COMPLEX_PIXEL(im,i,j)) : 0.0 )
423
mri_warp_bicubic_point(MRI_IMAGE * im,int ii,int jj,void wf (float,float,float *,float *))424 float mri_warp_bicubic_point( MRI_IMAGE *im , int ii , int jj ,
425 void wf( float,float,float *,float *) )
426 {
427 float xx,yy , fx,fy , newPt ;
428 int nx,ny , ix,jy , kk ;
429 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
430
431 nx = im->nx ; /* input image dimensions, for convenience */
432 ny = im->ny ;
433 kk = im->kind ;
434
435 /*** warp to output point ***/
436
437 wf( (float) ii , (float) jj , &xx,&yy ) ; /* get xx,yy in original image */
438
439 ix = floor( xx ) ; fx = xx - ix ; /* integer and */
440 jy = floor( yy ) ; fy = yy - jy ; /* fractional coords */
441
442 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ; /* weights for interpolating */
443 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ; /* in the x-direction */
444
445 f_jm1 = wt_m1 * INSIDE(ix-1,jy-1) /* interpolate to ix + fx */
446 + wt_00 * INSIDE(ix ,jy-1) /* at levels jy-1 .. jy+2 */
447 + wt_p1 * INSIDE(ix+1,jy-1)
448 + wt_p2 * INSIDE(ix+2,jy-1) ;
449
450 f_j00 = wt_m1 * INSIDE(ix-1,jy)
451 + wt_00 * INSIDE(ix ,jy)
452 + wt_p1 * INSIDE(ix+1,jy)
453 + wt_p2 * INSIDE(ix+2,jy) ;
454
455 f_jp1 = wt_m1 * INSIDE(ix-1,jy+1)
456 + wt_00 * INSIDE(ix ,jy+1)
457 + wt_p1 * INSIDE(ix+1,jy+1)
458 + wt_p2 * INSIDE(ix+2,jy+1) ;
459
460 f_jp2 = wt_m1 * INSIDE(ix-1,jy+2)
461 + wt_00 * INSIDE(ix ,jy+2)
462 + wt_p1 * INSIDE(ix+1,jy+2)
463 + wt_p2 * INSIDE(ix+2,jy+2) ;
464
465 /* interpolate between y-levels to jy+fy */
466
467 newPt = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
468 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
469
470 return newPt ;
471 }
472
473 /*--------------------------------------------------------------------
474 Return one point from a rotated image
475 ----------------------------------------------------------------------*/
476
mri_rotate_point(MRI_IMAGE * im,float aa,float bb,float phi,float scl,int ix,int jy)477 float mri_rotate_point( MRI_IMAGE *im, float aa, float bb, float phi, float scl ,
478 int ix , int jy )
479 {
480 MRI_IMAGE *imwarp ;
481 int nxnew , nynew ;
482
483 rot_cph = cos(phi) ;
484 rot_sph = sin(phi) ;
485
486 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
487 -(0.5 * im->ny) * rot_sph ;
488
489 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
490 +(0.5 * im->ny) * (1.0-rot_cph) ;
491
492 if( scl > 0.0 ){
493 rot_cph /= scl ;
494 rot_sph /= scl ;
495 }
496
497 return mri_warp_bicubic_point( im , ix,jy , xxMRI_rotfunc ) ;
498 }
499 #endif /* WARP_POINT_ROUTINES */
500