00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013 #ifdef HP
00014 # undef INLINE
00015 # pragma INLINE xxMRI_scaler
00016 # pragma INLINE xxMRI_rotfunc
00017 #endif
00018
00019 #define FINS(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00020 ? 0.0 : far[(i)+(j)*nx] )
00021
00022
00023
00024 static float sx_scale , sy_scale ;
00025
00026 INLINE void xxMRI_scaler( float xpr, float ypr, float *xx , float *yy )
00027 {
00028 *xx = sx_scale * xpr ;
00029 *yy = sy_scale * ypr ;
00030 return ;
00031 }
00032
00033
00034
00035 MRI_IMAGE *mri_warp( MRI_IMAGE *im , int nxnew , int nynew , int wtype ,
00036 void wfunc(float,float,float *,float *) )
00037 {
00038 switch( wtype ){
00039 case MRI_BILINEAR:
00040 return mri_warp_bilinear( im , nxnew , nynew , wfunc ) ;
00041
00042 case MRI_BICUBIC:
00043 return mri_warp_bicubic( im , nxnew , nynew , wfunc ) ;
00044
00045 default:
00046 fprintf( stderr , "mri_warp: illegal wtype %d\n" , wtype ) ;
00047 MRI_FATAL_ERROR ;
00048 }
00049 return NULL ;
00050 }
00051
00052
00053
00054 MRI_IMAGE *mri_resize_NN( MRI_IMAGE *im , int nxnew , int nynew )
00055 {
00056 int nx,ny , nnx,nny , ii,jj , pp,qq , bb ;
00057 float fx,fy ;
00058 MRI_IMAGE *nim ;
00059 char *nar , *ar ;
00060
00061 if( im == NULL ) return NULL ;
00062
00063 nx = im->nx ; ny = im->ny ;
00064 nnx = nxnew ; nny = nynew ;
00065 fx = ((float)nx) / (float)nnx ;
00066 fy = ((float)ny) / (float)nny ;
00067
00068 nim = mri_new( nnx , nny , im->kind ) ;
00069 nar = mri_data_pointer( nim ) ;
00070 ar = mri_data_pointer( im ) ;
00071 bb = im->pixel_size ;
00072
00073 for( jj=0 ; jj < nny ; jj++ ){
00074 qq = (int)( fy*jj ) ;
00075 for( ii=0 ; ii < nnx ; ii++ ){
00076 pp = (int)( fx*ii ) ;
00077 memcpy( nar + (ii+jj*nnx)*bb , ar + (pp+qq*nx)*bb , bb ) ;
00078 }
00079 }
00080
00081 MRI_COPY_AUX(nim,im) ;
00082 nim->dx *= fx ;
00083 nim->dy *= fy ;
00084 return nim ;
00085 }
00086
00087
00088
00089 MRI_IMAGE *mri_squareaspect( MRI_IMAGE *im )
00090 {
00091 int nx,ny , nxnew,nynew ;
00092 float dx,dy , rr ;
00093
00094 if( im == NULL ) return NULL ;
00095
00096 dx = fabs(im->dx) ; dy = fabs(im->dy) ;
00097 if( dx == 0.0 || dy == 0.0 ) return NULL ;
00098 rr = dy / dx ; if( rr == 1.0 ) return NULL ;
00099
00100 nx = im->nx ; ny = im->ny ;
00101
00102 if( rr < 1.0 ){
00103 nxnew = rint( nx/rr ) ; if( nxnew <= nx ) return NULL ;
00104 nynew = ny ;
00105 } else {
00106 nynew = rint( ny*rr ) ; if( nynew <= ny ) return NULL ;
00107 nxnew = nx ;
00108 }
00109
00110 return mri_resize_NN( im , nxnew , nynew ) ;
00111 }
00112
00113
00114
00115 MRI_IMAGE *mri_resize( MRI_IMAGE *im , int nxnew , int nynew )
00116 {
00117 int nx,ny , nnx,nny , wtype ;
00118
00119 nx = im->nx ; ny = im->ny ;
00120 nnx = nxnew ; nny = nynew ;
00121
00122 if( nnx <= 0 && nny <= 0 ){
00123 fprintf( stderr , "mri_resize: nxnew,nynew = %d %d\n",nxnew,nynew ) ;
00124 MRI_FATAL_ERROR ;
00125 }
00126
00127 sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ;
00128 sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ;
00129
00130 if( nnx <= 0 ){
00131 sx_scale = sy_scale ;
00132 nnx = sx_scale * nx ;
00133 } else if( nny <= 0 ){
00134 sy_scale = sx_scale ;
00135 nny = sy_scale * ny ;
00136 }
00137
00138
00139 #if 0
00140 wtype = MRI_BICUBIC ;
00141
00142 if( ( ((nnx>=nx) && (nnx%nx==0)) || ((nnx<nx) && (nx%nnx==0)) ) &&
00143 ( ((nny>=ny) && (nny%ny==0)) || ((nny<ny) && (ny%nny==0)) ) ){
00144
00145 wtype = MRI_BILINEAR ;
00146 } else {
00147 wtype = MRI_BICUBIC ;
00148 }
00149
00150 return mri_warp( im , nnx,nny , wtype , xxMRI_scaler ) ;
00151 #else
00152 return mri_warp_bicubic( im , nnx,nny , xxMRI_scaler ) ;
00153 #endif
00154 }
00155
00156
00157
00158 MRI_IMAGE *mri_warp_bicubic( MRI_IMAGE *im , int nxnew , int nynew ,
00159 void wf( float,float,float *,float *) )
00160 {
00161 MRI_IMAGE *imfl , *newImg ;
00162 float *far , *nar ;
00163 float xpr,ypr , xx,yy , fx,fy ;
00164 int ii,jj, nx,ny , ix,jy ;
00165 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00166 float bot,top,val ;
00167
00168 nx = im->nx ;
00169 ny = im->ny ;
00170
00171 nxnew = (nxnew > 0) ? nxnew : nx ;
00172 nynew = (nynew > 0) ? nynew : ny ;
00173
00174 switch( im->kind ){
00175
00176 case MRI_float:
00177 imfl = im ; break ;
00178
00179 default:
00180 imfl = mri_to_float(im) ; break ;
00181
00182 case MRI_short:{
00183 imfl = mri_to_float(im) ;
00184 newImg = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
00185 mri_free(imfl) ;
00186 imfl = mri_to_mri(MRI_short,newImg) ;
00187 mri_free(newImg) ; return imfl ;
00188 }
00189
00190 case MRI_byte:{
00191 imfl = mri_to_float(im) ;
00192 newImg = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
00193 mri_free(imfl) ;
00194 imfl = mri_to_mri(MRI_byte,newImg) ;
00195 mri_free(newImg) ; return imfl ;
00196 }
00197
00198 case MRI_rgb:{
00199 MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00200 MRI_IMAGE *rim,*gim,*bim ;
00201 rim = mri_warp_bicubic( IMARR_SUBIM(imar,0), nxnew,nynew, wf ) ;
00202 gim = mri_warp_bicubic( IMARR_SUBIM(imar,1), nxnew,nynew, wf ) ;
00203 bim = mri_warp_bicubic( IMARR_SUBIM(imar,2), nxnew,nynew, wf ) ;
00204 DESTROY_IMARR(imar) ;
00205 newImg = mri_3to_rgb( rim,gim,bim ) ;
00206 mri_free(rim); mri_free(gim); mri_free(bim); return newImg;
00207 }
00208
00209 }
00210
00211
00212
00213 far = mri_data_pointer( imfl ) ;
00214
00215 newImg = mri_new( nxnew , nynew , MRI_float ) ;
00216 nar = mri_data_pointer( newImg ) ;
00217
00218 bot = top = far[0] ;
00219 for( ii=1 ; ii < imfl->nvox ; ii++ ){
00220 if( far[ii] > top ) top = far[ii] ;
00221 else if( far[ii] < bot ) bot = far[ii] ;
00222 }
00223
00224
00225
00226 for( jj=0 ; jj < nynew ; jj++ ){
00227 ypr = jj ;
00228 for( ii=0 ; ii < nxnew ; ii++ ){
00229 xpr = ii ;
00230 wf( xpr,ypr , &xx,&yy ) ;
00231
00232 ix = floor( xx ) ; fx = xx - ix ;
00233 jy = floor( yy ) ; fy = yy - jy ;
00234
00235
00236
00237 #define P_M1(x) (-(x)*((x)-1)*((x)-2))
00238 #define P_00(x) (3*((x)+1)*((x)-1)*((x)-2))
00239 #define P_P1(x) (-3*(x)*((x)+1)*((x)-2))
00240 #define P_P2(x) ((x)*((x)+1)*((x)-1))
00241
00242 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;
00243 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
00244
00245 f_jm1 = wt_m1 * FINS(ix-1,jy-1)
00246 + wt_00 * FINS(ix ,jy-1)
00247 + wt_p1 * FINS(ix+1,jy-1)
00248 + wt_p2 * FINS(ix+2,jy-1) ;
00249
00250 f_j00 = wt_m1 * FINS(ix-1,jy)
00251 + wt_00 * FINS(ix ,jy)
00252 + wt_p1 * FINS(ix+1,jy)
00253 + wt_p2 * FINS(ix+2,jy) ;
00254
00255 f_jp1 = wt_m1 * FINS(ix-1,jy+1)
00256 + wt_00 * FINS(ix ,jy+1)
00257 + wt_p1 * FINS(ix+1,jy+1)
00258 + wt_p2 * FINS(ix+2,jy+1) ;
00259
00260 f_jp2 = wt_m1 * FINS(ix-1,jy+2)
00261 + wt_00 * FINS(ix ,jy+2)
00262 + wt_p1 * FINS(ix+1,jy+2)
00263 + wt_p2 * FINS(ix+2,jy+2) ;
00264
00265
00266
00267 val = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00268 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
00269
00270 if( val > top ) val = top ;
00271 else if( val < bot ) val = bot ;
00272
00273 nar[ii+jj*nxnew] = val ;
00274 }
00275 }
00276
00277
00278
00279 if( im != imfl ) mri_free(imfl) ;
00280 return newImg ;
00281 }
00282
00283
00284
00285 MRI_IMAGE *mri_warp_bilinear( MRI_IMAGE *im , int nxnew , int nynew ,
00286 void wf( float,float,float *,float *) )
00287 {
00288 MRI_IMAGE *imfl , *newImg ;
00289 float *far , *nar ;
00290 float xpr,ypr , xx,yy , fx,fx1,fy,fy1 , f00,f10,f01,f11 ;
00291 int ii,jj, nx,ny , ix,jy ;
00292
00293 nx = im->nx ;
00294 ny = im->ny ;
00295
00296 nxnew = (nxnew > 0) ? nxnew : nx ;
00297 nynew = (nynew > 0) ? nynew : ny ;
00298
00299 if( im->kind == MRI_float ){
00300 imfl = im ;
00301 } else {
00302 imfl = mri_to_float( im ) ;
00303 }
00304 far = mri_data_pointer( imfl ) ;
00305
00306 newImg = mri_new( nxnew , nynew , MRI_float ) ;
00307 nar = mri_data_pointer( newImg ) ;
00308
00309
00310
00311 for( jj=0 ; jj < nynew ; jj++ ){
00312 ypr = jj ;
00313 for( ii=0 ; ii < nxnew ; ii++ ){
00314 xpr = ii ;
00315 wf( xpr,ypr , &xx,&yy ) ;
00316
00317 ix = floor( xx ) ; fx = xx - ix ; fx1 = 1.0 - fx ;
00318 jy = floor( yy ) ; fy = yy - jy ; fy1 = 1.0 - fy ;
00319
00320 f00 = FINS(ix ,jy ) ;
00321 f01 = FINS(ix+1,jy ) ;
00322 f10 = FINS(ix ,jy+1) ;
00323 f11 = FINS(ix+1,jy+1) ;
00324
00325 nar[ii+jj*nxnew] = fx1 * ( fy1 * f00 + fy * f01 )
00326 +fx * ( fy1 * f10 + fy * f11 ) ;
00327
00328 }
00329 }
00330
00331
00332
00333 if( im != imfl ) mri_free(imfl) ;
00334 return newImg ;
00335 }
00336
00337
00338
00339 static float rot_dx , rot_dy , rot_cph , rot_sph ;
00340
00341 INLINE void xxMRI_rotfunc( float xpr , float ypr , float *xx , float *yy )
00342 {
00343 *xx = rot_cph * xpr + rot_sph * ypr + rot_dx ;
00344 *yy = -rot_sph * xpr + rot_cph * ypr + rot_dy ;
00345 return ;
00346 }
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356 MRI_IMAGE *mri_rotate( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
00357 {
00358 MRI_IMAGE *imwarp ;
00359 int nxnew , nynew ;
00360
00361 rot_cph = cos(phi) ;
00362 rot_sph = sin(phi) ;
00363
00364 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00365 -(0.5 * im->ny) * rot_sph ;
00366
00367 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00368 +(0.5 * im->ny) * (1.0-rot_cph) ;
00369
00370 if( scl <= 0.0 ){
00371 nxnew = nynew = 0 ;
00372 } else {
00373 nxnew = scl * im->nx + 0.49 ;
00374 nynew = scl * im->ny + 0.49 ;
00375 rot_cph /= scl ;
00376 rot_sph /= scl ;
00377 }
00378
00379 return mri_warp_bicubic( im , nxnew,nynew , xxMRI_rotfunc ) ;
00380 }
00381
00382 MRI_IMAGE *mri_rotate_bilinear( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
00383 {
00384 MRI_IMAGE *imwarp ;
00385 int nxnew , nynew ;
00386
00387 rot_cph = cos(phi) ;
00388 rot_sph = sin(phi) ;
00389
00390 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00391 -(0.5 * im->ny) * rot_sph ;
00392
00393 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00394 +(0.5 * im->ny) * (1.0-rot_cph) ;
00395
00396 if( scl <= 0.0 ){
00397 nxnew = nynew = 0 ;
00398 } else {
00399 nxnew = scl * im->nx + 0.49 ;
00400 nynew = scl * im->ny + 0.49 ;
00401 rot_cph /= scl ;
00402 rot_sph /= scl ;
00403 }
00404
00405 return mri_warp_bilinear( im , nxnew,nynew , xxMRI_rotfunc ) ;
00406 }
00407
00408 #undef WARP_POINT_ROUTINES
00409 #ifdef WARP_POINT_ROUTINES
00410
00411
00412
00413
00414
00415 #define INSIDE(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) ? 0.0 : \
00416 (kk==MRI_byte) ? (float) MRI_BYTE_PIXEL(im,i,j) : \
00417 (kk==MRI_short) ? (float) MRI_SHORT_PIXEL(im,i,j) : \
00418 (kk==MRI_int) ? (float) MRI_INT_PIXEL(im,i,j) : \
00419 (kk==MRI_float) ? (float) MRI_FLOAT_PIXEL(im,i,j) : \
00420 (kk==MRI_double) ? (float) MRI_DOUBLE_PIXEL(im,i,j) : \
00421 (kk==MRI_complex) ? (float) CABS(MRI_COMPLEX_PIXEL(im,i,j)) : 0.0 )
00422
00423 float mri_warp_bicubic_point( MRI_IMAGE *im , int ii , int jj ,
00424 void wf( float,float,float *,float *) )
00425 {
00426 float xx,yy , fx,fy , newPt ;
00427 int nx,ny , ix,jy , kk ;
00428 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00429
00430 nx = im->nx ;
00431 ny = im->ny ;
00432 kk = im->kind ;
00433
00434
00435
00436 wf( (float) ii , (float) jj , &xx,&yy ) ;
00437
00438 ix = floor( xx ) ; fx = xx - ix ;
00439 jy = floor( yy ) ; fy = yy - jy ;
00440
00441 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;
00442 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
00443
00444 f_jm1 = wt_m1 * INSIDE(ix-1,jy-1)
00445 + wt_00 * INSIDE(ix ,jy-1)
00446 + wt_p1 * INSIDE(ix+1,jy-1)
00447 + wt_p2 * INSIDE(ix+2,jy-1) ;
00448
00449 f_j00 = wt_m1 * INSIDE(ix-1,jy)
00450 + wt_00 * INSIDE(ix ,jy)
00451 + wt_p1 * INSIDE(ix+1,jy)
00452 + wt_p2 * INSIDE(ix+2,jy) ;
00453
00454 f_jp1 = wt_m1 * INSIDE(ix-1,jy+1)
00455 + wt_00 * INSIDE(ix ,jy+1)
00456 + wt_p1 * INSIDE(ix+1,jy+1)
00457 + wt_p2 * INSIDE(ix+2,jy+1) ;
00458
00459 f_jp2 = wt_m1 * INSIDE(ix-1,jy+2)
00460 + wt_00 * INSIDE(ix ,jy+2)
00461 + wt_p1 * INSIDE(ix+1,jy+2)
00462 + wt_p2 * INSIDE(ix+2,jy+2) ;
00463
00464
00465
00466 newPt = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00467 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
00468
00469 return newPt ;
00470 }
00471
00472
00473
00474
00475
00476 float mri_rotate_point( MRI_IMAGE *im, float aa, float bb, float phi, float scl ,
00477 int ix , int jy )
00478 {
00479 MRI_IMAGE *imwarp ;
00480 int nxnew , nynew ;
00481
00482 rot_cph = cos(phi) ;
00483 rot_sph = sin(phi) ;
00484
00485 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00486 -(0.5 * im->ny) * rot_sph ;
00487
00488 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00489 +(0.5 * im->ny) * (1.0-rot_cph) ;
00490
00491 if( scl > 0.0 ){
00492 rot_cph /= scl ;
00493 rot_sph /= scl ;
00494 }
00495
00496 return mri_warp_bicubic_point( im , ix,jy , xxMRI_rotfunc ) ;
00497 }
00498 #endif