#include "mrilib.h" static int verb = 1 ; #undef FREEIF #define FREEIF(x) if((x)!=NULL)free((void *)(x)) THD_3dim_dataset * THD_deghoster( THD_3dim_dataset *inset , THD_3dim_dataset *filset , int pe , int fe , int se ) ; void orfilt_vector( int nvec , float *vec ) ; static int orfilt_len = 11 ; /*----------------------------------------------------------------------------*/ int main( int argc , char *argv[] ) { char *prefix = "Deghost" ; int iarg ; int fe=1 , pe=2 , se=3 , nvals ; THD_3dim_dataset *inset=NULL , *outset , *filset=NULL ; if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ printf( "Usage: 3dDeghost [options] dataset\n" "\n" "* This program tries do remove N/2 (AKA Nyquist) ghosts from an EPI\n" " magnitude time series dataset.\n" "* If you apply it to some other kind of dataset (e.g., spiral), weird\n" " things will probably transpire.\n" "* The input EPI dataset should NOT be filtered, masked, cropped,\n" " registered, or pre-processed in any way!\n" "* This program will not work well if the input EPI dataset is heavily\n" " 'shaded' -- that is, its intensity varies dramatically inside the brain.\n" "* The output dataset is always stored in float format.\n" "* Only the Amitabha Buddha knows if this program is actually useful.\n" "\n" "========\n" "OPTIONS:\n" "========\n" " -input dataset = Another way to specify the input dataset\n" " -prefix pp = Use 'pp' for prefix of output dataset\n" " -FPS abc = Define the Frequency, Phase, and Slice\n" " directions in the dataset based on the\n" " axis orientations inside the dataset header\n" " (e.g., see the output of 3dinfo). The 'abc'\n" " code is a permutaton of the digits '123'.\n" " * The first digit 'a' specifies which dataset\n" " axis/index is the Frequency encoding direction.\n" " * The second digit 'b' specifies which dataset\n" " direction is the Phase encoding direction.\n" " * The third digit 'c' specifies which dataset\n" " direction is the Slice encoding direction.\n" " -->>** The default value for 'abc' is '123'; that is,\n" " the dataset is ordered so that the first index\n" " (x-axis) is frequency, the second index is phase,\n" " and the third index is slice. In most cases,\n" " this is how the reconstruction software will\n" " store the images. Only in unusual cases should\n" " you need the '-FPS' option!\n" " -filt N = Length of time series filter to apply when\n" " estimating ghosting parameters. Set N to 0 or 1\n" " to turn this feature off; otherwise, N should be an\n" " odd positive integer from 3 to 19 [default N=%d].\n" " * Longer filter lengths ARE allowed, but will be slow\n" " (cases with N <= 19 are hand coded for speed).\n" " * Datasets with fewer than 4 time points will not\n" " be filtered. For longer datasets, if the filter\n" " length is too big, it will be shortened ruthlessly.\n" "=======\n" "METHOD:\n" "=======\n" "Would you believe me if I said magic? Would you accept secret algorithms\n" "known only to the Olmecs? How about something so ad hoc that it cannot\n" "be described without embarrassment and shame?\n" "\n" "-- Feb 2014 - Zhark the Phantasmal\n" , orfilt_len ) ; PRINT_COMPILE_DATE ; exit(0) ; } mainENTRY("3dDeghost main"); machdep(); AFNI_logger("3dDeghost",argc,argv); PRINT_VERSION("3dDeghost") ; /*-- scan command line --*/ iarg = 1 ; while( iarg < argc && argv[iarg][0] == '-' ){ /*---*/ if( strcasecmp(argv[iarg],"-quiet") == 0 ){ verb = 0 ; iarg++ ; continue ; } if( strcasecmp(argv[iarg],"-verb") == 0 ){ verb++ ; iarg++ ; continue ; } /*---*/ if( strcasecmp(argv[iarg],"-filt") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ; orfilt_len = (int)strtod(argv[iarg],NULL) ; if( orfilt_len > 1 && orfilt_len%2 == 0 ){ orfilt_len++ ; INFO_message("-filt %d has been adjusted to %d (must be odd)" , orfilt_len-1 , orfilt_len) ; } if( orfilt_len > 19 ) WARNING_message("-filt %d is over the recommended limit of 19",orfilt_len) ; iarg++ ; continue ; } /*---*/ if( strcasecmp(argv[iarg],"-prefix") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ; prefix = argv[iarg] ; if( !THD_filename_ok(prefix) ) ERROR_exit("Illegal value after -prefix!\n"); iarg++ ; continue ; } /*---*/ if( strcasecmp(argv[iarg],"-input") == 0 || strcasecmp(argv[iarg],"-inset") == 0 ){ if( ++iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ; if( inset != NULL ) ERROR_exit("You can't give the input dataset twice!") ; inset = THD_open_dataset( argv[iarg] ) ; CHECK_OPEN_ERROR(inset,argv[iarg]) ; DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ; iarg++ ; continue ; } /*---*/ if( strcasecmp(argv[iarg],"-FPS") == 0 ){ /* stolen from 3dAllineate.c */ char *fps ; if( ++iarg >= argc ) ERROR_exit("Need argument after option '%s'",argv[iarg-1]) ; fps = argv[iarg] ; if( strlen(fps) < 3 ) ERROR_exit("Code '%s' after '%s' is too short", fps , argv[iarg-1] ) ; switch( fps[0] ){ default: ERROR_exit("Illegal '%s' F code '%c' :-(" , argv[iarg-1],fps[0] ); case 'i': case 'I': case 'x': case 'X': case '1': fe = 1; break; case 'j': case 'J': case 'y': case 'Y': case '2': fe = 2; break; case 'k': case 'K': case 'z': case 'Z': case '3': fe = 3; break; } switch( fps[1] ){ default: ERROR_exit("Illegal '%s' P code '%c' :-(" , argv[iarg-1],fps[1] ); case 'i': case 'I': case 'x': case 'X': case '1': pe = 1; break; case 'j': case 'J': case 'y': case 'Y': case '2': pe = 2; break; case 'k': case 'K': case 'z': case 'Z': case '3': pe = 3; break; } switch( fps[2] ){ default: ERROR_exit("Illegal '%s' S code '%c' :-(" , argv[iarg-1],fps[2] ); case 'i': case 'I': case 'x': case 'X': case '1': se = 1; break; case 'j': case 'J': case 'y': case 'Y': case '2': se = 2; break; case 'k': case 'K': case 'z': case 'Z': case '3': se = 3; break; } if( fe+pe+se != 6 ) ERROR_exit("Code '%s' after '%s' is nonsensical", fps , argv[iarg-1] ) ; iarg++ ; continue ; } /*---*/ ERROR_exit("Unknown option: %s\n",argv[iarg]); } if( inset == NULL && iarg >= argc ) ERROR_exit("No dataset name on command line?\n"); /*-- read input if needed --*/ if( inset == NULL ){ inset = THD_open_dataset( argv[iarg] ) ; CHECK_OPEN_ERROR(inset,argv[iarg]) ; DSET_load( inset ) ; CHECK_LOAD_ERROR(inset) ; } /*-- filter input? --*/ nvals = DSET_NVALS(inset) ; if( orfilt_len > nvals/2 ){ orfilt_len = nvals/2 ; if( orfilt_len%2 == 0 ) orfilt_len++ ; } if( orfilt_len > 1 && nvals > 1 ){ MRI_vectim *invect ; int ii ; if( verb ) INFO_message("Filtering input dataset: filter length=%d",orfilt_len) ; invect = THD_dset_to_vectim(inset,NULL,0) ; THD_vectim_applyfunc( invect , orfilt_vector ) ; filset = EDIT_empty_copy( inset ) ; for( ii=0 ; ii < nvals ; ii++ ) EDIT_substitute_brick( filset , ii , MRI_float , NULL ) ; THD_vectim_to_dset( invect , filset ) ; VECTIM_destroy(invect) ; } else { if( verb ) INFO_message("Time series filtering is turned off") ; } /***** outsource the work *****/ outset = THD_deghoster( inset , (filset!=NULL)?filset:inset , pe,fe,se ) ; if( outset == NULL ) ERROR_exit("THD_deghoster fails :-(((") ; if( filset != NULL ) DSET_delete(filset) ; EDIT_dset_items( outset , ADN_prefix,prefix , ADN_none ) ; tross_Copy_History( inset , outset ) ; tross_Make_History( "3dDeghost" , argc,argv , outset ) ; DSET_write(outset) ; WROTE_DSET(outset) ; exit(0) ; } /*============================================================================*/ /*** Stuff for the time series filtering ***/ #undef DTYPE #define DTYPE float #include "cs_qsort_small.h" #undef OFILT #define OFILT(j) 0.2f*(ar[j-1]+ar[j+1]+3.0f*ar[j]) /* odd lengths */ #undef EFILT #define EFILT(j) 0.2f*(ar[j-2]+ar[j+1])+0.3f*(ar[j-1]+ar[j]) /* even lengths */ static float orfilt( int n , float *ar ) { int nby2 ; switch( n ){ /* fast cases */ case 1: return ar[0] ; case 2: return 0.5f*(ar[0]+ar[1]) ; case 3: qsort3_float(ar) ; return OFILT(1) ; case 5: qsort5_float(ar) ; return OFILT(2) ; case 7: qsort7_float(ar) ; return OFILT(3) ; case 9: qsort9_float(ar) ; return OFILT(4) ; case 11: qsort11_float(ar) ; return OFILT(5) ; case 13: qsort13_float(ar) ; return OFILT(6) ; case 15: qsort15_float(ar) ; return OFILT(7) ; case 17: qsort17_float(ar) ; return OFILT(8) ; case 19: qsort19_float(ar) ; return OFILT(9) ; case 21: qsort21_float(ar) ; return OFILT(10) ; case 25: qsort25_float(ar) ; return OFILT(12) ; case 27: qsort27_float(ar) ; return OFILT(13) ; case 4: qsort4_float(ar) ; return EFILT(2) ; case 6: qsort6_float(ar) ; return EFILT(3) ; case 8: qsort8_float(ar) ; return EFILT(4) ; case 10: qsort10_float(ar) ; return EFILT(5) ; case 12: qsort12_float(ar) ; return EFILT(6) ; case 14: qsort14_float(ar) ; return EFILT(7) ; case 16: qsort16_float(ar) ; return EFILT(8) ; case 18: qsort18_float(ar) ; return EFILT(9) ; case 20: qsort20_float(ar) ; return EFILT(10) ; } /* general case for n not in above list -- will be slower */ qsort_float(n,ar) ; nby2 = n/2 ; return (n%2==0) ? EFILT(nby2) : OFILT(nby2) ; } #undef OFILT #undef EFILT /*----------------------------------------------------------------------------*/ void orfilt_vector( int nvec , float *vec ) { static float *ar=NULL ; static int nar=0 ; static float *vv=NULL ; static int nvv=0 ; int ii , ibot,itop , nby2=orfilt_len/2 , nv1=nvec-1,nii ; if( orfilt_len == 0 ){ if( ar != NULL ){ free(ar); nar = 0; ar = NULL; } if( vv != NULL ){ free(vv); nvv = 0; vv = NULL; } return ; } if( orfilt_len == 1 ) return ; if( ar == NULL || nar < orfilt_len ){ nar = orfilt_len ; ar = (float *)realloc(ar,sizeof(float)*nar) ; } if( vv == NULL || nvv < nvec ){ nvv = nvec ; vv = (float *)realloc(vv,sizeof(float)*nvv) ; } for( ii=0 ; ii < nvec ; ii++ ){ ibot = ii-nby2 ; if( ibot < 0 ) ibot = 0 ; itop = ii+nby2 ; if( itop > nv1 ) itop = nv1 ; nii = itop - ibot + 1 ; memcpy( ar , vec+ibot , sizeof(float)*nii ) ; vv[ii] = orfilt( nii , ar ) ; } memcpy( vec , vv , sizeof(float)*nvec ) ; return ; } /*============================================================================*/ /***------------------------------------------------------------------------***/ /* 'vec' variables are in smask = brain voxels whose N/2 location is outside */ static int nvec ; static float *bvec=NULL, *gvec=NULL, *xvec=NULL, *yvec=NULL, *ctvec=NULL, *stvec=NULL ; static byte *smask=NULL ; /* 'vim' variables cover the whole slice */ static int nvim ; static float *bvim=NULL, *gvim=NULL, *xvim=NULL, *yvim=NULL, *ctvim=NULL, *stvim=NULL ; static float noise_estimate=0.0f ; /* initial estimate of noise level */ /* these parameters define theta(x,y) */ static int ntheta_par=2 ; static double theta_par[2] , d_par ; void compute_theta( int npt , float *xpt , float *ypt , float *cpt , float *spt ) { float th0=theta_par[0] , th1=theta_par[1] , thth ; int ii ; for( ii=0 ; ii < npt ; ii++ ){ thth = th1*(xpt[ii]-th0) ; cpt[ii] = fabsf(cosf(thth)) ; spt[ii] = fabsf(sinf(thth)) ; } return ; } void compute_thvec(void) /* simplest model */ { compute_theta( nvec,xvec,yvec , ctvec,stvec ) ; } compute_thvim(void) { compute_theta( nvim,xvim,yvim , ctvim,stvim ) ; } /***------------------------------------------------------------------------***/ /*** Given Iy, Iyn, theta (in form of cos and sin), and noise level d, find M that is the solution to F(M) = Iy*cos(theta)^2 / sqrt( M^2 * cos(theta)^2 + d^2 ) + Iyn*sin(theta)^2 / sqrt( M^2 * sin(theta)^2 + d^2 ) - 1 = 0 If d=0, M is solved for analytically (mzero). If d > 0, then 3 steps of Newton's method are used -- unless d is too big, in which case M=0 is the return value. *//*------------------------------------------------------------------------***/ #undef COMPUTE_ff_df #define COMPUTE_ff_df(mm) \ { mq = (mm)*(mm); mcd = sqrt(mq*ctq+dq) ; msd = sqrt(mq*stq+dq) ; \ ff = iy*ctq/mcd + iyn*stq/msd - 1.0f ; \ df = -iy*(mm)*ctq*ctq/(mcd*mcd*mcd) - iyn*(mm)*stq*stq/(msd*msd*msd) ; } float find_mhat( float iy , float iyn , float ct , float st , float d ) { float mzero, ff,df, ctq=ct*ct, stq=st*st, mq, dq=d*d , mval , msd,mcd ; mzero = iy*ct + iyn*st ; /* analytic solution for d=0 */ if( d <= 0.0f ) return mzero ; /* easiest case */ if( d >= iy*ctq + iyn*stq ) return 0.0f ; /* unlikely */ COMPUTE_ff_df(mzero) ; if( ff >= 0.0f || df >= 0.0f ) return mzero ; /* should not happen */ mval = mzero - ff/df ; /* Newton step #1 */ if( mval <= 0.0 ) return 0.0 ; /* should not happen */ COMPUTE_ff_df(mval) ; if( df >= 0.0f ) return mval ; /* should not happen */ mval = mval - ff/df ; /* Newton step #2 */ if( mval <= 0.0 ) return 0.0 ; /* should not happen */ COMPUTE_ff_df(mval) ; if( df >= 0.0f ) return mval ; /* should not happen */ mval = mval - ff/df ; /* Newton step #3 */ if( mval <= 0.0 ) return 0.0 ; /* should not happen */ return mval ; } /***------------------------------------------------------------------------***/ /*** Given Iy, Iyn, theta (in form of cos and sin), and noise level d, find My and Myn that are the solution to Iy^2 = My^2 * cos(theta)^2 + Myn^2 * sin(theta)^2 + d^2 Iyn^2 = My^2 * sin(theta)^2 + Myn^2 * cos(theta)^2 + d^2 *//*------------------------------------------------------------------------***/ float_pair find_mpair( float iy , float iyn , float ct , float st , float d ) { float ctq=ct*ct , stq=st*st , dq=d*d , iyt,iynt , my,myn , den ; float_pair result={0.0f,0.0f} ; den = ctq*ctq - stq*stq ; if( den <= 0.0f ) return result ; /* bad inputs */ iyt = iy *iy - dq ; if( iyt < 0.0f ) iyt = 0.0f ; iynt = iyn*iyn - dq ; if( iynt < 0.0f ) iynt = 0.0f ; my = ( iyt*ctq - iynt*stq) / den ; if( my < 0.0f ) my = 0.0f ; myn = (-iyt*stq + iynt*ctq) / den ; if( myn < 0.0f ) myn = 0.0f ; result.a = sqrtf(my) ; result.b = sqrtf(myn) ; return result ; } /***------------------------------------------------------------------------***/ /* this function is the target for powell_newuoa_con() */ double theta_func( int npar , double *thpar ) { int ii ; double sum=0.0 ; float mhat,e1,e2 ; theta_par[0] = thpar[0] ; theta_par[1] = thpar[1] ; d_par = thpar[2] ; compute_thvec() ; /** INFO_message("========== theta_func(%g,%g,%g) ==========",thpar[0],thpar[1],thpar[2]) ; **/ for( ii=0 ; ii < nvec ; ii++ ){ mhat = find_mhat( bvec[ii],gvec[ii] , ctvec[ii],stvec[ii] , d_par ) ; e1 = bvec[ii]-sqrt(mhat*mhat*ctvec[ii]*ctvec[ii]+d_par*d_par) ; e2 = gvec[ii]-sqrt(mhat*mhat*stvec[ii]*stvec[ii]+d_par*d_par) ; sum += e1*e1 + e2*e2 ; /** ININFO_message(" bvec=%g gvec=%g ctvec=%g stvec=%g ==> mhat=%g e1=%g e2=%g", bvec[ii],gvec[ii] , ctvec[ii],stvec[ii] , mhat,e1,e2) ; **/ } /** penalty for nonzero parameters **/ #if 1 sum += 0.022*nvec*noise_estimate*noise_estimate * ( fabs(thpar[0]) + 99.9*fabs(thpar[1]) ) ; #endif /** ININFO_message("RESULT = %g",sum) ; **/ return sum ; } /***------------------------------------------------------------------------***/ void optimize_theta(void) { double thpar[3] , thbot[3] , thtop[3] ; int nite ; thpar[0] = 0.00 ; thbot[0] = -9.99 ; thtop[0] = 9.99 ; thpar[1] = 0.000; /* initial values */ thbot[1] = -0.002; /* lower limits */ thtop[1] = 0.099; /* upper limits */ thpar[2] = 0.111*noise_estimate ; thbot[2] = 0.001*noise_estimate ; thtop[2] = 3.333*noise_estimate ; nite = powell_newuoa_con( 3 , thpar,thbot,thtop , 299 , 0.111 , 0.0001 , 999 , theta_func ) ; if( fabs(thpar[1]) < 0.00111 ) thpar[0] = thpar[1] = 0.0 ; theta_par[0] = thpar[0] ; theta_par[1] = thpar[1] ; d_par = thpar[2] ; return ; } /***------------------------------------------------------------------------***/ #define CLFRAC 0.4f byte * DEG_automask_image( MRI_IMAGE *im ) { byte *bmask , *cmask ; float *iar , cval ; int nx=im->nx , ny=im->ny , nz=im->nz , nvox=im->nvox , ii ; THD_automask_set_clipfrac(CLFRAC) ; bmask = mri_automask_image(im) ; cmask = (byte *)malloc(sizeof(byte)*nx*ny*nz) ; memcpy(cmask,bmask,sizeof(byte)*nx*ny*nz) ; iar = MRI_FLOAT_PTR(im) ; cval = THD_cliplevel(im,CLFRAC) ; THD_mask_dilate(nx,ny,nz,cmask,3,2) ; /* embiggen the mask */ THD_mask_dilate(nx,ny,nz,cmask,3,2) ; THD_mask_dilate(nx,ny,nz,cmask,3,2) ; THD_mask_dilate(nx,ny,nz,cmask,3,2) ; for( ii=0 ; ii < nvox ; ii++ ) /* enlittle it now (cromulently) */ if( !bmask[ii] && cmask[ii] && iar[ii] < cval ) cmask[ii] = 0 ; free(bmask) ; return cmask ; } /***------------------------------------------------------------------------***/ #undef FREEUP #define FREEUP \ do{ mri_free(medim); FREEIF(bmask); FREEIF(amask); \ FREEIF(smask); FREEIF(bvec ); FREEIF(gvec ); FREEIF(xvec ); \ FREEIF(yvec ); FREEIF(ctvec); FREEIF(stvec); FREEIF(bvim ); \ FREEIF(gvim ); FREEIF(xvim ); FREEIF(yvim ); FREEIF(ctvim); \ FREEIF(stvim); orfilt_vector(0,NULL); \ FREEIF(xzero_t); FREEIF(thet1_t); FREEIF(dparr_t); \ } while(0) THD_3dim_dataset * THD_deghoster( THD_3dim_dataset *inset , THD_3dim_dataset *filset, int pe , int fe , int se ) { MRI_IMAGE *medim=NULL , *tim=NULL , *oim=NULL ; float cval, *mar=NULL , *tar=NULL , *oar=NULL ; float *xzero_t=NULL , *thet1_t=NULL , *dparr_t=NULL , t1med,t1bmv; byte *bmask=NULL , *amask=NULL , sm ; int nvox , nx,ny,nz , dp=0,df=0,ds=0 , np=0,nf=0,ns=0,np2,nf2 ; int pp,ff,ss,nfp , ii , ppg , nsm,ism , vv,nv , iim , sskip ; THD_3dim_dataset *outset=NULL ; float iy,iyn ; float_pair mp ; /* create brain mask (bmask) */ medim = THD_median_brick(inset) ; bmask = DEG_automask_image(medim) ; /* brain mask (we hope) */ nx = medim->nx ; ny = medim->ny ; nz = medim->nz ; nvox = medim->nvox ; nv = DSET_NVALS(inset) ; /* estimate noise level from data outside the mask (crudely) */ mar = MRI_FLOAT_PTR(medim) ; cval = THD_cliplevel(medim,CLFRAC) ; for( noise_estimate=0.0f,iim=ii=0 ; ii < nvox ; ii++ ){ if( !bmask[ii] && mar[ii] < cval ){ noise_estimate += mar[ii] ; iim++ ; } } if( iim < 9 ){ FREEUP; return NULL; } /* should not happen */ noise_estimate /= iim ; /* initial estimate of noise level */ if( verb > 1 ) INFO_message("Global crude noise_estimate = %g",noise_estimate) ; /* chop out all sub-threshold voxels (amask) */ amask = (byte *)malloc(sizeof(byte)*nvox) ; /* clipped brain mask */ memcpy(amask,bmask,sizeof(byte)*nvox) ; for( ii=0 ; ii < nvox ; ii++ ) if( amask[ii] && mar[ii] < cval ) amask[ii] = 0 ; /* setting up slice coordinates f,p,s */ if( pe == 1 ){ dp = 1 ; np = nx ; } else if( pe == 2 ){ dp = nx ; np = ny ; } else if( pe == 3 ){ dp = nx*ny ; np = ns ; } if( fe == 1 ){ df = 1 ; nf = nx ; } else if( fe == 2 ){ df = nx ; nf = ny ; } else if( fe == 3 ){ df = nx*ny ; nf = nz ; } if( se == 1 ){ ds = 1 ; ns = nx ; } else if( se == 2 ){ ds = nx ; ns = ny ; } else if( se == 3 ){ ds = nx*ny ; ns = nz ; } #undef IJK #define IJK(f,p,s) ((f)*df+(p)*dp+(s)*ds) nvim = nfp = nf * np ; np2 = np / 2 ; nf2 = nf / 2 ; smask = (byte * )malloc(sizeof(byte) *nfp) ; bvec = (float *)malloc(sizeof(float)*nfp) ; gvec = (float *)malloc(sizeof(float)*nfp) ; xvec = (float *)malloc(sizeof(float)*nfp) ; yvec = (float *)malloc(sizeof(float)*nfp) ; ctvec = (float *)malloc(sizeof(float)*nfp) ; stvec = (float *)malloc(sizeof(float)*nfp) ; bvim = (float *)malloc(sizeof(float)*nvim) ; gvim = (float *)malloc(sizeof(float)*nvim) ; xvim = (float *)malloc(sizeof(float)*nvim) ; yvim = (float *)malloc(sizeof(float)*nvim) ; ctvim = (float *)malloc(sizeof(float)*nvim) ; stvim = (float *)malloc(sizeof(float)*nvim) ; xzero_t = (float *)malloc(sizeof(float)*nv) ; thet1_t = (float *)malloc(sizeof(float)*nv) ; dparr_t = (float *)malloc(sizeof(float)*nv) ; /* copy input to output (will be ghost edited later) */ outset = EDIT_empty_copy(inset) ; for( vv=0 ; vv < nv ; vv++ ){ oim = THD_extract_float_brick(vv,inset) ; oar = MRI_FLOAT_PTR(oim) ; EDIT_BRICK_FACTOR( outset , vv , 0.0f ) ; EDIT_substitute_brick( outset , vv , MRI_float , oar ) ; mri_clear_and_free(oim) ; } /* loop over slices */ for( ss=0 ; ss < ns ; ss++ ){ /* make copy of brain mask in this slice, then edit it down to voxels in the brain whose N/2 point is outside the brain (smask) */ for( iim=nsm=pp=0 ; pp < np ; pp++ ){ if( pp >= np2 ) ppg = pp-np2 ; else ppg = pp+np2 ; for( ff=0 ; ff < nf ; ff++,iim++ ){ smask[iim] = sm = amask[IJK(ff,pp,ss)] && !bmask[IJK(ff,ppg,ss)] ; xvim[iim] = ff-nf2 ; yvim[iim] = pp-np2 ; if( sm ){ xvec[nsm] = xvim[iim]; yvec[nsm] = yvim[iim]; nsm++; } } } if( nsm < nfp/20 ){ /* skip this slice */ if( verb ) INFO_message("deghost: skipping slice #%d -- too few points in smask",ss) ; continue ; } nvec = nsm ; if( verb ) INFO_message("deghost: processing slice #%d",ss) ; /* smask is now the mask of brain voxels whose Nyquist ghost locations are NOT in the brain mask */ /* loop over time points, estimate the ghost correction parameters */ for( vv=0 ; vv < nv ; vv++ ){ tim = THD_extract_float_brick(vv,filset) ; tar = MRI_FLOAT_PTR(tim) ; /* extract the vector of image values in smask, and the vector of image values at the ghost locations */ for( iim=ism=pp=0 ; pp < np ; pp++ ){ if( pp >= np2 ) ppg = pp-np2 ; else ppg = pp+np2 ; for( ff=0 ; ff < nf ; ff++,iim++ ){ bvim[iim] = tar[IJK(ff,pp,ss)] ; gvim[iim] = tar[IJK(ff,ppg,ss)] ; if( smask[iim] ){ bvec[ism] = bvim[iim]; gvec[ism++] = gvim[iim]; } } } /* fit the theta parameters from the smask region and save them */ optimize_theta() ; xzero_t[vv] = theta_par[0] ; thet1_t[vv] = theta_par[1] ; dparr_t[vv] = d_par ; mri_free(tim) ; tim = NULL ; } /* now check the slice parameters for reasonability */ sskip = 0 ; if( nv > 4 ){ orfilt_len = 3 ; orfilt_vector(nv,xzero_t) ; orfilt_vector(nv,thet1_t) ; orfilt_vector(nv,dparr_t) ; qmedmadbmv_float(nv,thet1_t,&t1med,NULL,&t1bmv) ; if( verb ) ININFO_message(" slice #%d -- median(theta1)=%g stdev=%g ratio=%g", ss,t1med,t1bmv,(t1bmv>0.0f)?t1med/t1bmv:0.0f) ; if( t1med == 0.0f || fabsf(t1med) <= 0.111f*t1bmv ){ sskip = 1 ; ININFO_message(" skipping slice #%d -- theta1 too small",ss) ; } } if( sskip ) continue ; /* skip processing this slice */ /* loop over time points, estimate the un-ghosted image */ for( vv=0 ; vv < nv ; vv++ ){ tim = THD_extract_float_brick(vv,inset) ; /* input data for slice */ tar = MRI_FLOAT_PTR(tim) ; oar = DSET_ARRAY(outset,vv) ; /* output data for volume */ /* compute theta at each voxel */ if( thet1_t[vv] == 0.0f ) continue ; /* ghost amplitude is 0 ==> skip this time point */ if( verb > 1 ) ININFO_message(" slice=%d index=%d theta = %g %g %g", ss,vv,xzero_t[vv],thet1_t[vv],dparr_t[vv]) ; theta_par[0] = xzero_t[vv]; theta_par[1] = thet1_t[vv]; d_par = dparr_t[vv]; compute_thvim() ; /* compute output values at each voxel: (a) inside the smask == voxel in brain, N/2 ghost isn't (b) not in the smask but in the bmask == voxel && N/2 ghost are in brain (c) otherwise == voxel is unimportant effluvium */ for( iim=pp=0 ; pp < np ; pp++ ){ if( pp >= np2 ) ppg = pp-np2 ; else ppg = pp+np2 ; for( ff=0 ; ff < nf ; ff++,iim++ ){ iy = tar[IJK(ff,pp,ss)] ; iyn = tar[IJK(ff,ppg,ss)] ; if( smask[iim] ){ oar[IJK(ff,pp,ss)] = find_mhat( iy,iyn , ctvim[iim],stvim[iim] , d_par ) ; oar[IJK(ff,ppg,ss)] = 0.0f ; } else if( bmask[IJK(ff,pp,ss)] && bmask[IJK(ff,ppg,ss)] ){ if( ppg > pp ){ mp = find_mpair( iy,iyn , ctvim[iim],stvim[iim] , d_par ) ; oar[IJK(ff,pp,ss)] = mp.a ; oar[IJK(ff,ppg,ss)] = mp.b ; } } else { /* nada: output is already a copy of input */ } } } } } /* end of loop over slices */ FREEUP ; return outset ; }