00001 #include "mrilib.h"
00002 
00003 /*--------------------------------------------------------------------------*/
00004 
00005 MRI_IMAGE * mri_fft2D( MRI_IMAGE * im , int mode )
00006 {
00007    MRI_IMAGE * cxim , * outim ;
00008    int nx,ny , nxup,nyup , ii,jj ;
00009    complex * cxar , * outar , * cpt , * opt ;
00010 
00011    if( im == NULL ) return NULL ;
00012 
00013    /* convert input to complex */
00014 
00015    cxim = mri_to_complex(im) ;
00016    cxar = MRI_COMPLEX_PTR(cxim) ;
00017 
00018    /* compute size of output */
00019 
00020    nx = cxim->nx ; nxup = csfft_nextup_one35(nx) ;
00021    ny = cxim->ny ; nyup = csfft_nextup_one35(ny) ;
00022 
00023    /* create output array */
00024 
00025    outim = mri_new( nxup , nyup , MRI_complex ) ;
00026    outar = MRI_COMPLEX_PTR(outim) ;
00027 
00028    /* copy input to output, zero padding along the way */
00029 
00030    opt = outar ;
00031    cpt = cxar  ;
00032    for( jj=0 ; jj < ny ; jj++ ){
00033       for( ii=0 ; ii < nx   ; ii++ ){opt->r=cpt->r; opt->i=cpt->i; opt++; cpt++;}
00034       for(      ; ii < nxup ; ii++ ){opt->r=opt->i=0.0; opt++;}
00035    }
00036    for( ; jj < nyup ; jj++ ){opt->r=opt->i=0.0; opt++;}
00037 
00038    mri_free(cxim) ;
00039 
00040    /* row FFTs */
00041 
00042    for( jj=0 ; jj < ny ; jj++ )
00043       csfft_cox( mode , nxup , outar+jj*nxup ) ;
00044 
00045    /* column FFTs */
00046 
00047    cxar = (complex *) malloc(sizeof(complex)*nyup) ;
00048 
00049    for( ii=0 ; ii < nxup ; ii++ ){
00050       for( jj=0 ; jj < nyup ; jj++ ) cxar[jj] = outar[ii+jj*nxup] ;
00051       csfft_cox( mode , nyup , cxar ) ;
00052       for( jj=0 ; jj < nyup ; jj++ ) outar[ii+jj*nxup] = cxar[jj] ;
00053    }
00054 
00055    free(cxar) ; return outim ;
00056 }
00057 
00058 /*--------------------------------------------------------------------------*/
00059 
00060 MRI_IMAGE * cx_scramble( MRI_IMAGE * ima , MRI_IMAGE * imb ,
00061                          float alpha     , float beta       )
00062 {
00063    int ii , npix ;
00064    double r1,r2 , t1,t2 , rr,tt ;
00065    complex * ar , * br , * cr ;
00066    double aa,aa1 , bb,bb1 ;
00067    MRI_IMAGE * imc ;
00068 
00069    if( ima == NULL        || ima->kind != MRI_complex ||
00070        imb == NULL        || imb->kind != MRI_complex ||
00071        ima->nx != imb->nx || ima->ny != imb->ny       ||
00072        alpha < 0.0        || alpha > 1.0              ||
00073        beta  < 0.0        || beta  > 1.0                ) return NULL ;
00074 
00075    npix = ima->nvox ;
00076    ar   = MRI_COMPLEX_PTR(ima) ;
00077    br   = MRI_COMPLEX_PTR(imb) ;
00078    imc  = mri_new_conforming( ima , MRI_complex ) ;
00079    cr   = MRI_COMPLEX_PTR(imc) ;
00080 
00081    aa   = alpha ; aa1 = 1.0 - aa ;
00082    bb   = beta  ; bb1 = 1.0 - bb ;
00083 
00084    for( ii=0 ; ii < npix ; ii++ ){
00085       r1 = CABS(ar[ii]) ; r2 = CABS(br[ii]) ; rr = pow(r1,aa)*pow(r2,aa1) ;
00086       t1 = CARG(ar[ii]) ; t2 = CARG(br[ii]) ; tt = t1-t2 ;
00087            if( tt < -PI ) t2 -= 2.0*PI ;
00088       else if( tt >  PI ) t2 += 2.0*PI ;
00089       tt = bb*t1 + bb1*t2 ;
00090       cr[ii].r = rr * cos(tt) ; cr[ii].i = rr * sin(tt) ;
00091    }
00092 
00093    return imc ;
00094 }
00095 
00096 /*--------------------------------------------------------------------------*/
00097 
00098 MRI_IMAGE * mri_scramble( MRI_IMAGE * ima , MRI_IMAGE * imb ,
00099                           float alpha     , float beta       )
00100 {
00101    MRI_IMAGE * cxa, * cxb, * cxc ;
00102    int nx,ny , nxup,nyup ;
00103 
00104    if( ima == NULL        || imb == NULL        ||
00105        ima->nx != imb->nx || ima->ny != imb->ny ||
00106        alpha < 0.0        || alpha > 1.0        ||
00107        beta  < 0.0        || beta  > 1.0          ) return NULL ;
00108 
00109    cxa = mri_fft2D( ima , -1 ) ;
00110    cxb = mri_fft2D( imb , -1 ) ;
00111    cxc = cx_scramble( cxa,cxb,alpha,beta ) ;
00112    mri_free(cxa) ; mri_free(cxb) ;
00113    cxa = mri_fft2D( cxc , 1 ) ;
00114    mri_free(cxc) ;
00115    cxb = mri_to_mri( ima->kind , cxa ) ;
00116    mri_free(cxa) ;
00117    return cxb ;
00118 }
00119 
00120 /*--------------------------------------------------------------------------*/
00121 
00122 int main( int argc , char * argv[] )
00123 {
00124    float alpha,beta ;
00125    MRI_IMAGE * ima , * imb , * imc , * imd ;
00126    int nx,ny,nxup,nyup ;
00127 
00128    if( argc < 6 ){
00129       printf("Usage: phace alpha beta im1 im2 fname\n") ; exit(0) ;
00130    }
00131 
00132    alpha = strtod(argv[1],NULL) ;
00133    beta  = strtod(argv[2],NULL) ;
00134    ima   = mri_read(argv[3]) ;
00135    imb   = mri_read(argv[4]) ;
00136    imc   = mri_scramble( ima,imb , alpha,beta )  ;
00137 
00138    nx   = ima->nx ; ny   = ima->ny ;
00139    nxup = imc->nx ; nyup = imc->ny ;
00140 
00141    if( nx < nxup || ny < nyup ){
00142       imd = mri_cut_2D( imc , 0,nx-1 , 0,ny-1 ) ;
00143       mri_free(imc) ;
00144       imc = imd ;
00145    }
00146 
00147    mri_free(ima) ; mri_free(imb) ;
00148 
00149    mri_write( argv[5] , imd ) ;
00150    exit(0) ;
00151 }