#include "mrilib.h" /*** NOT 7D SAFE ***/ /****************************************************************************/ /** functions to "warp" 3D images -- not very efficient, but quite general **/ /****************************************************************************/ /*--------------------------------------------------------------------------*/ static int wtype = MRI_LINEAR ; void mri_warp3D_method( int mode ){ wtype = mode ; } /** set interpolation **/ /*--------------------------------------------------------------------------*/ static byte *womask = NULL ; /* 19 Nov 2004 */ void mri_warp3D_set_womask( MRI_IMAGE *wim ) { womask = (wim == NULL || wim->kind != MRI_byte) ? NULL : MRI_BYTE_PTR(wim) ; } #undef SKIP #define SKIP(i,j,k) (womask!=NULL && womask[(i)+(j)*nxnew+(k)*nxynew]==0) /*--------------------------------------------------------------------------*/ static int zout = 1 ; void mri_warp3D_zerout( int zzzz ){ zout = zzzz ; } /** outside = 0? **/ /*--------------------------------------------------------------------------*/ /* Macros for interpolation and access to data arrays. */ #undef FAR #define FAR(i,j,k) far[(i)+(j)*nx+(k)*nxy] /* input point */ #undef NAR #define NAR(i,j,k) nar[(i)+(j)*nxnew+(k)*nxynew] /* output point */ /** clip input mm to range [0..nn] (for indexing) **/ #undef CLIP #define CLIP(mm,nn) if(mm < 0)mm=0; else if(mm > nn)mm=nn /* define Lagrange cubic interpolation polynomials */ #undef P_M1 #undef P_00 #undef P_P1 #undef P_P2 #undef P_FACTOR #define P_M1(x) (-(x)*((x)-1)*((x)-2)) #define P_00(x) (3*((x)+1)*((x)-1)*((x)-2)) #define P_P1(x) (-3*(x)*((x)+1)*((x)-2)) #define P_P2(x) ((x)*((x)+1)*((x)-1)) #define P_FACTOR 4.62962963e-3 /* 1/216 = final scaling factor */ /*---------------------------------------------------------------------------*/ /*! Transform a 3D image geometrically, using cubic interpolation. See comments for mri_warp3D() for details about the parameters. -----------------------------------------------------------------------------*/ INLINE MRI_IMAGE *mri_warp3D_cubic( MRI_IMAGE *im, int nxnew, int nynew, int nznew, void wf(float,float,float,float *,float *,float *) ) { MRI_IMAGE *imfl , *newImg ; float *far , *nar ; float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ; int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ; float bot,top,val ; float nxh,nyh,nzh ; int ix_m1,ix_00,ix_p1,ix_p2 ; /* interpolation indices */ int jy_m1,jy_00,jy_p1,jy_p2 ; /* (input image) */ int kz_m1,kz_00,kz_p1,kz_p2 ; float wt_m1,wt_00,wt_p1,wt_p2 ; /* interpolation weights */ float f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, /* interpolants */ f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00, f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1, f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2, f_km1 , f_k00 , f_kp1 , f_kp2 ; ENTRY("mri_warp3D_cubic") ; /*--- sanity check inputs ---*/ if( im == NULL || wf == NULL ) RETURN(NULL) ; /*-- dimensional analysis --*/ nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */ ny = im->ny ; ny1 = ny-1 ; nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ; nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */ nynew = (nynew > 0) ? nynew : ny ; nznew = (nznew > 0) ? nznew : nz ; nxynew = nxnew*nynew ; /*----- Allow for different input image types, by breaking them into components, doing each one separately, and then reassembling. This is needed since we only warp float-valued images below. -----*/ switch( im->kind ){ case MRI_float: /* use input directly */ imfl = im ; break ; default:{ /* floatize input */ imfl = mri_to_float(im) ; newImg = mri_warp3D_cubic( imfl , nxnew,nynew,nznew , wf ) ; mri_free(imfl) ; imfl = mri_to_mri(im->kind,newImg) ; if( imfl != NULL ){ mri_free(newImg); newImg = imfl; } RETURN(newImg) ; } case MRI_rgb:{ /* break into 3 pieces */ MRI_IMARR *imar = mri_rgb_to_3float(im) ; MRI_IMAGE *rim,*gim,*bim ; rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; gim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; bim = mri_warp3D_cubic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,2) ) ; FREE_IMARR(imar) ; newImg = mri_3to_rgb( rim,gim,bim ) ; mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg); } case MRI_complex:{ /* break into 2 pieces */ MRI_IMARR *imar = mri_complex_to_pair(im) ; MRI_IMAGE *rim, *iim ; rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; iim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; FREE_IMARR(imar) ; newImg = mri_pair_to_complex( rim , iim ) ; mri_free(rim); mri_free(iim); RETURN(newImg); } } /* end of special cases of input datum */ /*----- at this point, imfl is in float format -----*/ far = MRI_FLOAT_PTR( imfl ) ; /* input image data */ newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */ nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */ bot = top = far[0] ; /* find input data range */ for( ii=1 ; ii < imfl->nvox ; ii++ ){ if( far[ii] > top ) top = far[ii] ; else if( far[ii] < bot ) bot = far[ii] ; } /*** loop over output points and warp to them ***/ nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ; for( kk=0 ; kk < nznew ; kk++ ){ zpr = kk ; for( jj=0 ; jj < nynew ; jj++ ){ ypr = jj ; for( ii=0 ; ii < nxnew ; ii++ ){ xpr = ii ; if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */ wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */ if( zout && (xx < -0.5 || xx > nxh || yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ NAR(ii,jj,kk) = 0.0 ; continue ; } ix = floor(xx) ; fx = xx - ix ; /* integer and */ jy = floor(yy) ; fy = yy - jy ; /* fractional coords */ kz = floor(zz) ; fz = zz - kz ; ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ; CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ; jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ; CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ; kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ; CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ; wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ; /* interpolation weights */ wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ; #undef XINT #define XINT(j,k) wt_m1*FAR(ix_m1,j,k)+wt_00*FAR(ix_00,j,k) \ +wt_p1*FAR(ix_p1,j,k)+wt_p2*FAR(ix_p2,j,k) /* interpolate to location ix+fx at each jy,kz level */ f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ; f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ; f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ; f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ; f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ; f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ; /* interpolate to jy+fy at each kz level */ wt_m1 = P_M1(fy) ; wt_00 = P_00(fy) ; wt_p1 = P_P1(fy) ; wt_p2 = P_P2(fy) ; f_km1 = wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 ; f_k00 = wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 ; f_kp1 = wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 ; f_kp2 = wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 ; /* interpolate to kz+fz to get output */ wt_m1 = P_M1(fz) ; wt_00 = P_00(fz) ; wt_p1 = P_P1(fz) ; wt_p2 = P_P2(fz) ; val = P_FACTOR * ( wt_m1 * f_km1 + wt_00 * f_k00 + wt_p1 * f_kp1 + wt_p2 * f_kp2 ) ; if( val > top ) val = top ; /* clip to input data range */ else if( val < bot ) val = bot ; NAR(ii,jj,kk) = val ; } } } /*** cleanup and return ***/ if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ RETURN(newImg); } /*---------------------------------------------------------------------------*/ /*! Transform a 3D image geometrically, using linear interpolation. See comments for mri_warp3D() for details about the parameters. -----------------------------------------------------------------------------*/ INLINE MRI_IMAGE *mri_warp3D_linear( MRI_IMAGE *im, int nxnew, int nynew, int nznew, void wf(float,float,float,float *,float *,float *) ) { MRI_IMAGE *imfl , *newImg ; float *far , *nar ; float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ; int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ; float val ; float nxh,nyh,nzh ; int ix_00,ix_p1 ; /* interpolation indices */ int jy_00,jy_p1 ; /* (input image) */ int kz_00,kz_p1 ; float wt_00,wt_p1 ; /* interpolation weights */ float f_j00_k00, f_jp1_k00, f_j00_kp1, f_jp1_kp1, f_k00, f_kp1 ; #if 0 int nzset ; float zzsum ; #endif ENTRY("mri_warp3D_linear") ; /*--- sanity check inputs ---*/ if( im == NULL || wf == NULL ) RETURN(NULL) ; /*-- dimensional analysis --*/ nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */ ny = im->ny ; ny1 = ny-1 ; nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ; #if 0 fprintf(stderr,"mri_warp3D_linear: nx=%d ny=%d nz=%d\n",nx,ny,nz) ; nzset = 0 ; zzsum = 0.0 ; #endif nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */ nynew = (nynew > 0) ? nynew : ny ; nznew = (nznew > 0) ? nznew : nz ; nxynew = nxnew*nynew ; /*----- allow for different input image types, by breaking them into components, doing each one separately, and then reassembling -----*/ switch( im->kind ){ case MRI_float: /* use input directly */ imfl = im ; break ; default:{ /* floatize-input */ imfl = mri_to_float(im) ; newImg = mri_warp3D_linear( imfl , nxnew,nynew,nznew , wf ) ; mri_free(imfl) ; imfl = mri_to_mri(im->kind,newImg) ; if( imfl != NULL ){ mri_free(newImg); newImg = imfl; } RETURN(newImg) ; } case MRI_rgb:{ MRI_IMARR *imar = mri_rgb_to_3float(im) ; MRI_IMAGE *rim,*gim,*bim ; rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; gim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; bim = mri_warp3D_linear( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,2) ) ; FREE_IMARR(imar) ; newImg = mri_3to_rgb( rim,gim,bim ) ; mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg); } case MRI_complex:{ MRI_IMARR *imar = mri_complex_to_pair(im) ; MRI_IMAGE *rim, *iim ; rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; iim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; FREE_IMARR(imar) ; newImg = mri_pair_to_complex( rim , iim ) ; mri_free(rim); mri_free(iim); RETURN(newImg); } } /* end of special cases of input datum */ /*----- at this point, imfl is in float format -----*/ far = MRI_FLOAT_PTR( imfl ) ; /* input image data */ newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */ nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */ /*** loop over output points and warp to them ***/ nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ; for( kk=0 ; kk < nznew ; kk++ ){ zpr = kk ; for( jj=0 ; jj < nynew ; jj++ ){ ypr = jj ; for( ii=0 ; ii < nxnew ; ii++ ){ xpr = ii ; if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */ wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */ if( zout && (xx < -0.5 || xx > nxh || yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ #if 0 nzset++ ; #endif NAR(ii,jj,kk) = 0.0 ; continue ; } ix = floor(xx) ; fx = xx - ix ; /* integer and */ jy = floor(yy) ; fy = yy - jy ; /* fractional coords */ kz = floor(zz) ; fz = zz - kz ; ix_00 = ix ; ix_p1 = ix+1 ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; jy_00 = jy ; jy_p1 = jy+1 ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; kz_00 = kz ; kz_p1 = kz+1 ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; wt_00 = 1.0-fx ; wt_p1 = fx ; #undef XINT #define XINT(j,k) wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k) /* interpolate to location ix+fx at each jy,kz level */ f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ; f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ; /* interpolate to jy+fy at each kz level */ wt_00 = 1.0-fy ; wt_p1 = fy ; f_k00 = wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 ; f_kp1 = wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 ; /* interpolate to kz+fz to get output */ val = (1.0-fz) * f_k00 + fz * f_kp1 ; NAR(ii,jj,kk) = val ; #if 0 zzsum += fabsf(val) ; #endif } } } /*** cleanup and return ***/ if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ #if 0 zzsum /= (nxnew*nynew*nznew) ; fprintf(stderr," nzset=%d zzsum=%g\n",nzset,zzsum) ; #endif RETURN(newImg); } /*---------------------------------------------------------------------------*/ /*! Transform a 3D image geometrically, using NN 'interpolation'. See comments for mri_warp3D() for details about the parameters. -----------------------------------------------------------------------------*/ INLINE MRI_IMAGE *mri_warp3D_NN( MRI_IMAGE *im, int nxnew, int nynew, int nznew, void wf(float,float,float,float *,float *,float *) ) { MRI_IMAGE *imfl , *newImg ; float *far , *nar ; float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ; int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ; float nxh,nyh,nzh ; ENTRY("mri_warp3D_NN") ; /*--- sanity check inputs ---*/ if( im == NULL || wf == NULL ) RETURN(NULL) ; /*-- dimensional analysis --*/ nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */ ny = im->ny ; ny1 = ny-1 ; nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ; nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */ nynew = (nynew > 0) ? nynew : ny ; nznew = (nznew > 0) ? nznew : nz ; nxynew = nxnew*nynew ; /*----- allow for different input image types, by breaking them into components, doing each one separately, and then reassembling -----*/ switch( im->kind ){ case MRI_float: /* use input directly */ imfl = im ; break ; default:{ /* floatize-input */ imfl = mri_to_float(im) ; newImg = mri_warp3D_NN( imfl , nxnew,nynew,nznew , wf ) ; mri_free(imfl) ; imfl = mri_to_mri(im->kind,newImg) ; if( imfl != NULL ){ mri_free(newImg); newImg = imfl; } RETURN(newImg) ; } case MRI_rgb:{ MRI_IMARR *imar = mri_rgb_to_3float(im) ; MRI_IMAGE *rim,*gim,*bim ; rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; gim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; bim = mri_warp3D_NN( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,2) ) ; FREE_IMARR(imar) ; newImg = mri_3to_rgb( rim,gim,bim ) ; mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg); } case MRI_complex:{ MRI_IMARR *imar = mri_complex_to_pair(im) ; MRI_IMAGE *rim, *iim ; rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; iim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; FREE_IMARR(imar) ; newImg = mri_pair_to_complex( rim , iim ) ; mri_free(rim); mri_free(iim); RETURN(newImg); } } /* end of special cases of input datum */ /*----- at this point, imfl is in float format -----*/ far = MRI_FLOAT_PTR( imfl ) ; /* input image data */ newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */ nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */ /*** loop over output points and warp to them ***/ nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ; for( kk=0 ; kk < nznew ; kk++ ){ zpr = kk ; for( jj=0 ; jj < nynew ; jj++ ){ ypr = jj ; for( ii=0 ; ii < nxnew ; ii++ ){ xpr = ii ; if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */ wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */ if( zout && (xx < -0.5 || xx > nxh || yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ NAR(ii,jj,kk) = 0.0 ; continue ; } ix = (int)rint(xx) ; jy = (int)rint(yy) ; kz = (int)rint(zz) ; CLIP(ix,nx1) ; CLIP(jy,ny1) ; CLIP(kz,nz1) ; NAR(ii,jj,kk) = FAR(ix,jy,kz) ; } } } /*** cleanup and return ***/ if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ RETURN(newImg); } /* define quintic interpolation polynomials (Lagrange) */ #undef Q_M2 #undef Q_M1 #undef Q_00 #undef Q_P1 #undef Q_P2 #undef Q_P3 #define Q_M2(x) (x*(x*x-1.0)*(2.0-x)*(x-3.0)*0.008333333) #define Q_M1(x) (x*(x*x-4.0)*(x-1.0)*(x-3.0)*0.041666667) #define Q_00(x) ((x*x-4.0)*(x*x-1.0)*(3.0-x)*0.083333333) #define Q_P1(x) (x*(x*x-4.0)*(x+1.0)*(x-3.0)*0.083333333) #define Q_P2(x) (x*(x*x-1.0)*(x+2.0)*(3.0-x)*0.041666667) #define Q_P3(x) (x*(x*x-1.0)*(x*x-4.0)*0.008333333) /*---------------------------------------------------------------------------*/ /*! Transform a 3D image geometrically, using quintic interpolation. See comments for mri_warp3D() for details about the parameters. - RWCox - 06 Aug 2003 -----------------------------------------------------------------------------*/ INLINE MRI_IMAGE *mri_warp3D_quintic( MRI_IMAGE *im, int nxnew, int nynew, int nznew, void wf(float,float,float,float *,float *,float *) ) { MRI_IMAGE *imfl , *newImg ; float *far , *nar ; float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ; int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ; float bot,top,val ; float nxh,nyh,nzh ; int ix_m2,ix_m1,ix_00,ix_p1,ix_p2,ix_p3 ; /* interpolation indices */ int jy_m2,jy_m1,jy_00,jy_p1,jy_p2,jy_p3 ; /* (input image) */ int kz_m2,kz_m1,kz_00,kz_p1,kz_p2,kz_p3 ; float wt_m2,wt_m1,wt_00,wt_p1,wt_p2,wt_p3 ; /* interpolation weights */ float f_jm2_km2, f_jm1_km2, f_j00_km2, f_jp1_km2, f_jp2_km2, f_jp3_km2, f_jm2_km1, f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, f_jp3_km1, f_jm2_k00, f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00, f_jp3_k00, f_jm2_kp1, f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1, f_jp3_kp1, f_jm2_kp2, f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2, f_jp3_kp2, f_jm2_kp3, f_jm1_kp3, f_j00_kp3, f_jp1_kp3, f_jp2_kp3, f_jp3_kp3, f_km2 , f_km1 , f_k00 , f_kp1 , f_kp2 , f_kp3 ; ENTRY("mri_warp3D_quintic") ; /*--- sanity check inputs ---*/ if( im == NULL || wf == NULL ) RETURN(NULL) ; /*-- dimensional analysis --*/ nx = im->nx ; nx1 = nx-1 ; /* input image dimensions */ ny = im->ny ; ny1 = ny-1 ; nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ; nxnew = (nxnew > 0) ? nxnew : nx ; /* output image dimensions */ nynew = (nynew > 0) ? nynew : ny ; nznew = (nznew > 0) ? nznew : nz ; nxynew = nxnew*nynew ; /*----- Allow for different input image types, by breaking them into components, doing each one separately, and then reassembling. This is needed since we only warp float-valued images below. -----*/ switch( im->kind ){ case MRI_float: /* use input directly */ imfl = im ; break ; default:{ /* floatize input */ imfl = mri_to_float(im) ; newImg = mri_warp3D_quintic( imfl , nxnew,nynew,nznew , wf ) ; mri_free(imfl) ; imfl = mri_to_mri(im->kind,newImg) ; if( imfl != NULL ){ mri_free(newImg); newImg = imfl; } RETURN(newImg) ; } case MRI_rgb:{ /* break into 3 pieces */ MRI_IMARR *imar = mri_rgb_to_3float(im) ; MRI_IMAGE *rim,*gim,*bim ; rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; gim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; bim = mri_warp3D_quintic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,2) ) ; FREE_IMARR(imar) ; newImg = mri_3to_rgb( rim,gim,bim ) ; mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg); } case MRI_complex:{ /* break into 2 pieces */ MRI_IMARR *imar = mri_complex_to_pair(im) ; MRI_IMAGE *rim, *iim ; rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,0) ) ; iim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ; mri_free( IMARR_SUBIM(imar,1) ) ; FREE_IMARR(imar) ; newImg = mri_pair_to_complex( rim , iim ) ; mri_free(rim); mri_free(iim); RETURN(newImg); } } /* end of special cases of input datum */ /*----- at this point, imfl is in float format -----*/ far = MRI_FLOAT_PTR( imfl ) ; /* input image data */ newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ; /* make output image */ nar = MRI_FLOAT_PTR( newImg ) ; /* output image data */ bot = top = far[0] ; /* find input data range */ for( ii=1 ; ii < imfl->nvox ; ii++ ){ if( far[ii] > top ) top = far[ii] ; else if( far[ii] < bot ) bot = far[ii] ; } /*** loop over output points and warp to them ***/ nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ; for( kk=0 ; kk < nznew ; kk++ ){ zpr = kk ; for( jj=0 ; jj < nynew ; jj++ ){ ypr = jj ; for( ii=0 ; ii < nxnew ; ii++ ){ xpr = ii ; if( SKIP(ii,jj,kk) ) continue ; /* 19 Nov 2004 */ wf( xpr,ypr,zpr , &xx,&yy,&zz ) ; /* get xx,yy,zz in original image */ if( zout && (xx < -0.5 || xx > nxh || yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){ NAR(ii,jj,kk) = 0.0 ; continue ; } ix = floor(xx) ; fx = xx - ix ; /* integer and */ jy = floor(yy) ; fy = yy - jy ; /* fractional coords */ kz = floor(zz) ; fz = zz - kz ; /* in input image */ /* compute indexes from which to interpolate (-2,-1,0,+1,+2,+3), but clipped to lie inside input image volume */ ix_m1 = ix-1 ; ix_00 = ix ; ix_p1 = ix+1 ; ix_p2 = ix+2 ; CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ; ix_m2 = ix-2 ; ix_p3 = ix+3 ; CLIP(ix_m2,nx1) ; CLIP(ix_p3,nx1) ; jy_m1 = jy-1 ; jy_00 = jy ; jy_p1 = jy+1 ; jy_p2 = jy+2 ; CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ; jy_m2 = jy-2 ; jy_p3 = jy+3 ; CLIP(jy_m2,ny1) ; CLIP(jy_p3,ny1) ; kz_m1 = kz-1 ; kz_00 = kz ; kz_p1 = kz+1 ; kz_p2 = kz+2 ; CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ; kz_m2 = kz-2 ; kz_p3 = kz+3 ; CLIP(kz_m2,nz1) ; CLIP(kz_p3,nz1) ; wt_m1 = Q_M1(fx) ; wt_00 = Q_00(fx) ; /* interpolation weights */ wt_p1 = Q_P1(fx) ; wt_p2 = Q_P2(fx) ; /* in x-direction */ wt_m2 = Q_M2(fx) ; wt_p3 = Q_P3(fx) ; #undef XINT #define XINT(j,k) wt_m2*FAR(ix_m2,j,k)+wt_m1*FAR(ix_m1,j,k) \ +wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k) \ +wt_p2*FAR(ix_p2,j,k)+wt_p3*FAR(ix_p3,j,k) /* interpolate to location ix+fx at each jy,kz level */ f_jm2_km2 = XINT(jy_m2,kz_m2) ; f_jm1_km2 = XINT(jy_m1,kz_m2) ; f_j00_km2 = XINT(jy_00,kz_m2) ; f_jp1_km2 = XINT(jy_p1,kz_m2) ; f_jp2_km2 = XINT(jy_p2,kz_m2) ; f_jp3_km2 = XINT(jy_p3,kz_m2) ; f_jm2_km1 = XINT(jy_m2,kz_m1) ; f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ; f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ; f_jp3_km1 = XINT(jy_p3,kz_m1) ; f_jm2_k00 = XINT(jy_m2,kz_00) ; f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ; f_jp3_k00 = XINT(jy_p3,kz_00) ; f_jm2_kp1 = XINT(jy_m2,kz_p1) ; f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ; f_jp3_kp1 = XINT(jy_p3,kz_p1) ; f_jm2_kp2 = XINT(jy_m2,kz_p2) ; f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ; f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ; f_jp3_kp2 = XINT(jy_p3,kz_p2) ; f_jm2_kp3 = XINT(jy_m2,kz_p3) ; f_jm1_kp3 = XINT(jy_m1,kz_p3) ; f_j00_kp3 = XINT(jy_00,kz_p3) ; f_jp1_kp3 = XINT(jy_p1,kz_p3) ; f_jp2_kp3 = XINT(jy_p2,kz_p3) ; f_jp3_kp3 = XINT(jy_p3,kz_p3) ; /* interpolate to jy+fy at each kz level */ wt_m1 = Q_M1(fy) ; wt_00 = Q_00(fy) ; wt_p1 = Q_P1(fy) ; wt_p2 = Q_P2(fy) ; wt_m2 = Q_M2(fy) ; wt_p3 = Q_P3(fy) ; f_km2 = wt_m2 * f_jm2_km2 + wt_m1 * f_jm1_km2 + wt_00 * f_j00_km2 + wt_p1 * f_jp1_km2 + wt_p2 * f_jp2_km2 + wt_p3 * f_jp3_km2 ; f_km1 = wt_m2 * f_jm2_km1 + wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1 + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 + wt_p3 * f_jp3_km1 ; f_k00 = wt_m2 * f_jm2_k00 + wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 + wt_p3 * f_jp3_k00 ; f_kp1 = wt_m2 * f_jm2_kp1 + wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 + wt_p3 * f_jp3_kp1 ; f_kp2 = wt_m2 * f_jm2_kp2 + wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2 + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 + wt_p3 * f_jp3_kp2 ; f_kp3 = wt_m2 * f_jm2_kp3 + wt_m1 * f_jm1_kp3 + wt_00 * f_j00_kp3 + wt_p1 * f_jp1_kp3 + wt_p2 * f_jp2_kp3 + wt_p3 * f_jp3_kp3 ; /* interpolate to kz+fz to get output */ wt_m1 = Q_M1(fz) ; wt_00 = Q_00(fz) ; wt_p1 = Q_P1(fz) ; wt_p2 = Q_P2(fz) ; wt_m2 = Q_M2(fz) ; wt_p3 = Q_P3(fz) ; val = wt_m2 * f_km2 + wt_m1 * f_km1 + wt_00 * f_k00 + wt_p1 * f_kp1 + wt_p2 * f_kp2 + wt_p3 * f_kp3 ; if( val > top ) val = top ; /* clip to input data range */ else if( val < bot ) val = bot ; NAR(ii,jj,kk) = val ; } } } /*** cleanup and return ***/ if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ RETURN(newImg); } /*--------------------------------------------------------------------------*/ /*! Generic warping function. Calls the specific one for the desired interpolation method, which was set by mri_warp3D_method(). - im = input image (3D). - nxnew,nynew,nznew = size of output image grid. - wf( inew,jnew,knew , iold,jold,kold ) is a function that takes as input 3 indices in [0..nxnew-1,0..nynew-1,0..nznew-1] (the output image index space) and returns indices in the input image space. - indices in and out of wf() are floats. - by default, values outside the original image volume are treated as zero; this mode can be turned off using mri_warp3D_zerout(). ----------------------------------------------------------------------------*/ INLINE MRI_IMAGE *mri_warp3D( MRI_IMAGE *im, int nxnew, int nynew, int nznew, void wf(float,float,float,float *,float *,float *) ) { switch( wtype ){ default: case MRI_CUBIC: return mri_warp3D_cubic ( im , nxnew,nynew,nznew , wf ) ; case MRI_LINEAR: return mri_warp3D_linear ( im , nxnew,nynew,nznew , wf ) ; case MRI_NN: return mri_warp3D_NN ( im , nxnew,nynew,nznew , wf ) ; case MRI_QUINTIC: /* 06 Aug 2003 */ return mri_warp3D_quintic( im , nxnew,nynew,nznew , wf ) ; } return NULL ; /* unreachable */ } /*---------------------------------------------------------------------------*/ static float sx_scale, sy_scale, sz_scale ; /* global scaler data */ /*! Internal transform function for resize "warp". */ static INLINE void w3dMRI_scaler( float a,float b,float c, float *x,float *y,float *z ) { *x = sx_scale * a; *y = sy_scale * b; *z = sz_scale * c; } /*--------------------------------------------------------------------------*/ /*! Special image warp just for resizing. ----------------------------------------------------------------------------*/ MRI_IMAGE *mri_warp3D_resize( MRI_IMAGE *im , int nxnew, int nynew, int nznew ) { int nx,ny,nz , nnx,nny,nnz ; if( im == NULL ) return NULL ; nx = im->nx ; ny = im->ny ; nz = im->nz ; nnx = nxnew ; nny = nynew ; nnz = nznew ; if( nnx <= 0 && nny <= 0 && nnz <= 0 ) return NULL ; sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ; /* global variables */ sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ; sz_scale = (nnz > 0) ? ((float)nz)/nnz : 0.0 ; if( nnx <= 0 ){ sx_scale = MAX(sy_scale,sz_scale) ; nnx = (int)rint(sx_scale*nx) ; } if( nny <= 0 ){ sy_scale = MAX(sx_scale,sz_scale) ; nny = (int)rint(sy_scale*ny) ; } if( nnz <= 0 ){ sz_scale = MAX(sx_scale,sy_scale) ; nnz = (int)rint(sz_scale*nz) ; } return mri_warp3D( im , nnx,nny,nnz , w3dMRI_scaler ) ; } /*---------------------------------------------------------------------------*/ static float a11_aff,a12_aff,a13_aff,a14_aff , a21_aff,a22_aff,a23_aff,a24_aff , a31_aff,a32_aff,a33_aff,a34_aff ; /*! Internal transform function for affine "warp". */ static INLINE void w3dMRI_affine( float a,float b,float c, float *x,float *y,float *z ) { *x = a11_aff*a + a12_aff*b + a13_aff*c + a14_aff ; *y = a21_aff*a + a22_aff*b + a23_aff*c + a24_aff ; *z = a31_aff*a + a32_aff*b + a33_aff*c + a34_aff ; } /*--------------------------------------------------------------------------*/ /*! Special image warp for affine transformation of indices. ----------------------------------------------------------------------------*/ MRI_IMAGE * mri_warp3D_affine( MRI_IMAGE *im , THD_vecmat aff ) { if( im == NULL ) return NULL ; a11_aff = aff.mm.mat[0][0] ; a12_aff = aff.mm.mat[0][1] ; a13_aff = aff.mm.mat[0][2] ; a14_aff = aff.vv.xyz[0] ; a21_aff = aff.mm.mat[1][0] ; a22_aff = aff.mm.mat[1][1] ; a23_aff = aff.mm.mat[1][2] ; a24_aff = aff.vv.xyz[1] ; a31_aff = aff.mm.mat[2][0] ; a32_aff = aff.mm.mat[2][1] ; a33_aff = aff.mm.mat[2][2] ; a34_aff = aff.vv.xyz[2] ; return mri_warp3D( im , 0,0,0 , w3dMRI_scaler ) ; } /*===========================================================================*/ /*===========================================================================*/ /******** Functions for warping a dataset using coordinate transforms, *******/ /******** rather than the index transformation of the functions above. *******/ /*===========================================================================*/ /*===========================================================================*/ static THD_vecmat ijk_to_dicom_in, ijk_to_dicom_out, dicom_to_ijk_in ; static void (*warp_out_to_in)(float ,float ,float , float *,float *,float * ) ; /*----------------------------------------------------------------------------*/ /*! This is the function passed to mri_warp3D() for the ijk_out to ijk_in transformation. It does the ijk_out to dicom_out transform here, then calls the coordinate transform function for dicom_out to dicom_in, then does the dicom_in to ijk_in transform here. ------------------------------------------------------------------------------*/ static INLINE void w3d_warp_func( float xout, float yout, float zout, float *xin, float *yin, float *zin ) { THD_fvec3 xxx,yyy ; float xi,yi,zi ; LOAD_FVEC3(xxx,xout,yout,zout) ; yyy = VECMAT_VEC( ijk_to_dicom_out , xxx ) ; /* ijk_out to dicom_out */ warp_out_to_in( yyy.xyz[0],yyy.xyz[1],yyy.xyz[2], /* dicom_out to dicom_in */ &(xxx.xyz[0]) , &(xxx.xyz[1]), &(xxx.xyz[2]) ) ; yyy = VECMAT_VEC( dicom_to_ijk_in , xxx ) ; /* dicom_in to ijk_in */ *xin = yyy.xyz[0] ; *yin = yyy.xyz[1] ; *zin = yyy.xyz[2] ; } /*--------------------------------------------------------------------------*/ /*! - Find the 8 corners of the input dataset (voxel edges, not centers). - Warp each one using the provided wfunc(). - Return the min and max (x,y,z) coordinates of these warped points. ----------------------------------------------------------------------------*/ static void warp_corners( THD_3dim_dataset *inset, void wfunc(float,float,float,float *,float *,float *), float *xb , float *xt , float *yb , float *yt , float *zb , float *zt ) { THD_dataxes *daxes = inset->daxes ; THD_fvec3 corner , wcorn ; float nx0 = -0.5 , ny0 = -0.5 , nz0 = -0.5 ; float nx1 = daxes->nxx-0.5, ny1 = daxes->nyy-0.5, nz1 = daxes->nzz-0.5 ; float xx,yy,zz , xbot,ybot,zbot , xtop,ytop,ztop ; LOAD_FVEC3( corner , nx0,ny0,nz0 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = xtop = xx ; ybot = ytop = yy ; zbot = ztop = zz ; LOAD_FVEC3( corner , nx1,ny0,nz0 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx0,ny1,nz0 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx1,ny1,nz0 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx0,ny0,nz1 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx1,ny0,nz1 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx0,ny1,nz1 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); LOAD_FVEC3( corner , nx1,ny1,nz1 ) ; wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ; wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ; xbot = MIN(xx,xbot); xtop = MAX(xx,xtop); ybot = MIN(yy,ybot); ytop = MAX(yy,ytop); zbot = MIN(zz,zbot); ztop = MAX(zz,ztop); *xb = xbot; *xt = xtop; *yb = ybot; *yt = ytop; *zb = zbot; *zt = ztop; return; } /*--------------------------------------------------------------------------*/ /*! Geometrically transform a 3D dataset, producing a new dataset. - w_in2out transforms DICOM coords from input grid to output grid (only needed if newggg != NULL and flag == WARP3D_NEWGRID; otherwise, w_in2out can be NULL) - w_out2in is the inverse of w_in2out (cannot be NULL) - newggg is a pointer to something that determines the new grid on which the output dataset will be computed - NULL means compute on the inset grid - otherwise, flag determines how newggg is interpreted: - WARP3D_NEWGRID => newggg is a "float *" to a new grid size - WARP3D_NEWDSET => newggg is a "THD_3dim_dataset *" to a dataset whose grid will be used for the output grid - prefix = new dataset prefix (if NULL, then "warped") - zpad = number of planes to zeropad on each face of inset (>= 0) - flag = see above - Interpolation method can be set using mri_warp3D_method(). - At end, input dataset is unloaded from memory. - If input is 3D+time, the output dataset won't have slice-wise time offsets, even if the input did. ----------------------------------------------------------------------------*/ INLINE THD_3dim_dataset * THD_warp3D( THD_3dim_dataset *inset , void w_in2out(float,float,float,float *,float *,float *), void w_out2in(float,float,float,float *,float *,float *), void *newggg , char *prefix , int zpad , int flag ) { THD_3dim_dataset *outset , *qset , *newgset=NULL ; int nxin,nyin,nzin,nvals , ival ; int nxout,nyout,nzout ; float xbot,xtop , ybot,ytop , zbot,ztop ; int use_newgrid , use_oldgrid ; float ddd_newgrid , fac ; MRI_IMAGE *inim , *outim , *wim ; ENTRY("THD_warp3D") ; /*-- decide which grid the output will be computed on --*/ use_newgrid = ((flag & WARP3D_GRIDMASK) == WARP3D_NEWGRID) && (newggg != NULL) ; if( use_newgrid ){ float *gg = (float *)newggg ; if( thd_floatscan(1,gg) == 0 && *gg > 0.0 ) ddd_newgrid = *gg ; else use_newgrid = 0 ; } else if( newggg != NULL && (flag & WARP3D_GRIDMASK) == WARP3D_NEWDSET ){ newgset = (THD_3dim_dataset *) newggg ; if( !ISVALID_DSET(newgset) ) newgset = NULL ; } use_oldgrid = !use_newgrid && (newgset == NULL) ; /*-- see if inputs are valid --*/ if( !ISVALID_DSET(inset) || w_out2in == NULL || (w_in2out == NULL && use_newgrid) ){ fprintf(stderr,"** ERROR: THD_warp3D has bad inputs!\n") ; RETURN(NULL); } /*-- zeropad and replace input dataset, if desired --*/ if( zpad > 0 ){ qset = THD_zeropad( inset , zpad,zpad,zpad,zpad,zpad,zpad , "Quetzalcoatl" , ZPAD_PURGE ) ; DSET_unload(inset) ; if( qset == NULL ){ fprintf(stderr,"** ERROR: THD_warp3D can't zeropad!\n"); RETURN(NULL); } } else { qset = inset ; } /*-- read input data from disk, if not already present --*/ DSET_load(qset) ; if( !DSET_LOADED(qset) ){ fprintf(stderr,"** ERROR: THD_warp3D can't load input dataset!\n") ; if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ; RETURN(NULL); } /*-- compute mapping from input dataset (i,j,k) to DICOM coords --*/ { THD_vecmat ijk_to_xyz , xyz_to_dicom ; LOAD_DIAG_MAT( ijk_to_xyz.mm , qset->daxes->xxdel, qset->daxes->yydel, qset->daxes->zzdel ); LOAD_FVEC3 ( ijk_to_xyz.vv , qset->daxes->xxorg, qset->daxes->yyorg, qset->daxes->zzorg ); xyz_to_dicom.mm = qset->daxes->to_dicomm ; LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ; ijk_to_dicom_in = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ; dicom_to_ijk_in = INV_VECMAT( ijk_to_dicom_in ) ; } /*-- make empty output dataset --*/ nxin = DSET_NX(qset) ; nyin = DSET_NY(qset) ; nzin = DSET_NZ(qset) ; nvals = DSET_NVALS(qset) ; if( nxin*nyin*nzin <= 1 ){ fprintf(stderr,"** ERROR: THD_warp3D has nxin=%d nyin=%d nzin=%d!\n", nxin,nyin,nzin ) ; if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ; RETURN(NULL) ; } if( use_oldgrid ){ /*-- output is on same grid as input --*/ outset = EDIT_empty_copy( qset ) ; } else if( newgset != NULL ){ /*-- output on same grid as newgset --*/ outset = EDIT_empty_copy( newgset ) ; EDIT_dset_items( outset , ADN_nvals , nvals , ADN_type , qset->type , ADN_view_type , qset->view_type , ADN_func_type , qset->func_type , ADN_malloc_type , DATABLOCK_MEM_MALLOC , ADN_none ) ; if( DSET_NUM_TIMES(qset) > 1 ) /* and some time structure? */ EDIT_dset_items( outset , ADN_ntt , nvals , ADN_tunits , DSET_TIMEUNITS(qset) , ADN_ttorg , DSET_TIMEORIGIN(qset) , ADN_ttdel , DSET_TR(qset) , ADN_ttdur , DSET_TIMEDURATION(qset) , ADN_none ) ; THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ; /* 08 Jun 2004 */ INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ; } else { /*-- output is on new grid --*/ float xmid,ymid,zmid ; THD_ivec3 nxyz , orixyz ; THD_fvec3 dxyz , orgxyz ; /* compute DICOM coordinates of warped corners */ warp_corners( qset, w_in2out, &xbot,&xtop, &ybot,&ytop, &zbot,&ztop ) ; nxout = (int)( (xtop-xbot)/ddd_newgrid+0.999 ); if( nxout < 1 ) nxout = 1; nyout = (int)( (ytop-ybot)/ddd_newgrid+0.999 ); if( nyout < 1 ) nyout = 1; nzout = (int)( (ztop-zbot)/ddd_newgrid+0.999 ); if( nzout < 1 ) nzout = 1; xmid = 0.5*(xbot+xtop); ymid = 0.5*(ybot+ytop); zmid = 0.5*(zbot+ztop); xbot = xmid-0.5*(nxout-1)*ddd_newgrid; xtop = xbot+(nxout-1)*ddd_newgrid; ybot = ymid-0.5*(nyout-1)*ddd_newgrid; ytop = ybot+(nyout-1)*ddd_newgrid; zbot = zmid-0.5*(nzout-1)*ddd_newgrid; ztop = zbot+(nzout-1)*ddd_newgrid; #if 0 if( verb ) fprintf(stderr,"++ Transformed grid:\n" "++ xbot = %10.4g xtop = %10.4g nx = %d\n" "++ ybot = %10.4g ytop = %10.4g ny = %d\n" "++ zbot = %10.4g ztop = %10.4g nz = %d\n" , xbot,xtop,nxout , ybot,ytop,nyout , zbot,ztop,nzout ) ; #endif if( nxout*nyout*nzout <= 1 ){ fprintf(stderr,"** ERROR: THD_warp3D has nxout=%d nyout=%d nzout=%d!\n", nxout,nyout,nzout ) ; if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ; RETURN(NULL) ; } nxyz.ijk[0] = nxout ; dxyz.xyz[0] = ddd_newgrid ; /* setup axes */ nxyz.ijk[1] = nyout ; dxyz.xyz[1] = ddd_newgrid ; nxyz.ijk[2] = nzout ; dxyz.xyz[2] = ddd_newgrid ; orixyz.ijk[0] = ORI_R2L_TYPE ; orgxyz.xyz[0] = xbot ; orixyz.ijk[1] = ORI_A2P_TYPE ; orgxyz.xyz[1] = ybot ; orixyz.ijk[2] = ORI_I2S_TYPE ; orgxyz.xyz[2] = zbot ; /** create dataset and mangle it into the desired shape **/ outset = EDIT_empty_copy( NULL ) ; /* void and formless dataset */ EDIT_dset_items( outset , /* give it some structure! */ ADN_nxyz , nxyz , ADN_xyzdel , dxyz , ADN_xyzorg , orgxyz , ADN_xyzorient , orixyz , ADN_malloc_type , DATABLOCK_MEM_MALLOC , ADN_nvals , nvals , ADN_type , qset->type , ADN_view_type , qset->view_type , ADN_func_type , qset->func_type , ADN_none ) ; if( DSET_NUM_TIMES(qset) > 1 ) /* and some time structure? */ EDIT_dset_items( outset , ADN_ntt , nvals , ADN_tunits , DSET_TIMEUNITS(qset) , ADN_ttorg , DSET_TIMEORIGIN(qset) , ADN_ttdel , DSET_TR(qset) , ADN_ttdur , DSET_TIMEDURATION(qset) , ADN_none ) ; THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ; /* 08 Jun 2004 */ INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ; } /*-- end of warping to new grid --*/ nxout = DSET_NX(outset) ; nyout = DSET_NY(outset) ; nzout = DSET_NZ(outset) ; /*-- compute mapping from output dataset (i,j,k) to DICOM coords --*/ { THD_vecmat ijk_to_xyz , xyz_to_dicom ; LOAD_DIAG_MAT( ijk_to_xyz.mm, outset->daxes->xxdel, outset->daxes->yydel, outset->daxes->zzdel ); LOAD_FVEC3 ( ijk_to_xyz.vv, outset->daxes->xxorg, outset->daxes->yyorg, outset->daxes->zzorg ); xyz_to_dicom.mm = outset->daxes->to_dicomm ; LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ; ijk_to_dicom_out = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ; } /*-- add prefix to new dataset --*/ if( !THD_filename_ok(prefix) ) prefix = "warped" ; EDIT_dset_items( outset , ADN_prefix,prefix , ADN_none ) ; /*-- loop over bricks and warp them --*/ warp_out_to_in = w_out2in ; /* for use in w3d_warp_func(), supra */ for( ival=0 ; ival < nvals ; ival++ ){ inim = DSET_BRICK(qset,ival) ; fac = DSET_BRICK_FACTOR(qset,ival) ; if( fac > 0.0 && fac != 1.0 ) wim = mri_scale_to_float( fac , inim ) ; else wim = inim ; outim = mri_warp3D( wim , nxout,nyout,nzout , w3d_warp_func ) ; if( outim == NULL ){ fprintf(stderr,"** ERROR: THD_warp3D fails at ival=%d\n",ival); DSET_delete(outset); if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ; RETURN(NULL) ; } if( wim != inim ){ /* if we scaled input */ float gtop , fimfac ; mri_free(wim) ; /* 20 Oct 2003: rescale if input was an integer type */ if( outim->kind == MRI_float && MRI_IS_INT_TYPE(inim->kind) ){ fimfac = fac ; wim = mri_scalize( outim , inim->kind , &fimfac ) ; EDIT_BRICK_FACTOR( outset , ival , fimfac ) ; mri_free(outim) ; outim = wim ; } else { /* input not an integer: */ EDIT_BRICK_FACTOR( outset , ival , 0.0 ) ; /* output=unscaled float */ } } EDIT_substitute_brick( outset, ival,outim->kind,mri_data_pointer(outim) ); DSET_unload_one( qset , ival ) ; } /*-- done!!! --*/ if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ; THD_load_statistics( outset ) ; RETURN( outset ) ; } /*--------------------------------------------------------------------------*/ static THD_vecmat aff_out2in , aff_in2out ; /* data for afo2i() and afi2o() */ /*! Internal transform for THD_warp3D_affine(). ----------------------------*/ static INLINE void afo2i( float a,float b,float c, float *x,float *y,float *z ) { THD_fvec3 xxx , yyy ; LOAD_FVEC3( xxx , a,b,c ) ; yyy = VECMAT_VEC( aff_out2in , xxx ) ; UNLOAD_FVEC3( yyy , *x , *y , *z ) ; } /*! Internal transform for THD_warp3D_affine(). ----------------------------*/ static INLINE void afi2o( float a,float b,float c, float *x,float *y,float *z ) { THD_fvec3 xxx , yyy ; LOAD_FVEC3( xxx , a,b,c ) ; yyy = VECMAT_VEC( aff_in2out , xxx ) ; UNLOAD_FVEC3( yyy , *x , *y , *z ) ; } /*--------------------------------------------------------------------------*/ /*! A version of THD_warp3D() that uses an affine map, specified in the out2in parameter. The main advantage over THD_warp3D() is that the transformation functions internal -- afo2i() and afi2o() -- which makes it simpler on the caller and also somewhat faster. ----------------------------------------------------------------------------*/ THD_3dim_dataset * THD_warp3D_affine( THD_3dim_dataset *inset , THD_vecmat out2in , void *newggg , char *prefix , int zpad , int flag ) { aff_out2in = out2in ; aff_in2out = INV_VECMAT(aff_out2in) ; return THD_warp3D( inset , afi2o,afo2i , newggg,prefix,zpad,flag ) ; } /*--------------------------------------------------------------------------*/ /*! Internal transform functions for TT <-> MNI coords. Both are RAI, as per AFNI internal logic. 11 Mar 2004 - RW Cox [Jury Duty Day] cf. http://www.mrc-cbu.cam.ac.uk/Imaging/Common/mnispace.shtml ----------------------------------------------------------------------------*/ static INLINE void w3d_mni2tta( float mx , float my , float mz , float *tx, float *ty, float *tz ) { *tx = 0.99 * mx ; if( mz > 0.0 ){ *ty = 0.9688 * my + 0.0460 * mz ; *tz = -0.0485 * my + 0.9189 * mz ; } else { *ty = 0.9688 * my + 0.0420 * mz ; *tz = -0.0485 * my + 0.8390 * mz ; } } /*------- Should be the inverse of the above! -----------------------------*/ static INLINE void w3d_tta2mni( float tx , float ty , float tz , float *mx, float *my, float *mz ) { *mx = 1.01010 * tx ; *my = 1.02962 * ty - 0.05154 * tz ; *mz = 0.05434 * ty + 1.08554 * tz ; if( *mz < 0.0 ) *mz *= 1.09523 ; } /*--------------------------------------------------------------------------*/ THD_3dim_dataset * THD_warp3D_tta2mni( THD_3dim_dataset *inset , void *newggg , char *prefix , int zpad , int flag ) { return THD_warp3D( inset , w3d_tta2mni,w3d_mni2tta , NULL,prefix,0,0 ) ; } /*------- Should be the inverse of the above! -----------------------------*/ THD_3dim_dataset * THD_warp3D_mni2tta( THD_3dim_dataset *inset , void *newggg , char *prefix , int zpad , int flag ) { return THD_warp3D( inset , w3d_mni2tta,w3d_tta2mni , newggg,prefix,zpad,flag ) ; }