Doxygen Source Code Documentation
3dDespike.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | TFAC 0.1 |
#define | ITFAC 10.0 |
Functions | |
INLINE float | mytanh (float x) |
int | main (int argc, char *argv[]) |
Define Documentation
|
Definition at line 9 of file 3dDespike.c. Referenced by main(). |
|
Definition at line 8 of file 3dDespike.c. Referenced by main(). |
Function Documentation
|
compute the overall minimum and maximum voxel values for a dataset Definition at line 24 of file 3dDespike.c. References ADN_brick_fac, ADN_datum_all, ADN_func_type, ADN_none, ADN_prefix, AFNI_logger(), argc, BYTEIZE, cl1_solve(), DSET_ARRAY, DSET_BRICK_BYTES, DSET_BRICK_TYPE, DSET_BRIKNAME, DSET_datum_constant, DSET_delete, DSET_HEADNAME, DSET_index_to_ix, DSET_index_to_jy, DSET_index_to_kz, DSET_load, DSET_LOADED, DSET_NUM_TIMES, DSET_NVOX, DSET_NX, DSET_NY, DSET_NZ, DSET_unload_one, DSET_write, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), far, fit, free, FUNC_FIM_TYPE, ISVALID_DSET, ITFAC, machdep(), mainENTRY, malloc, mytanh(), PRINT_VERSION, qg(), qmed_float(), ref, strtod(), TFAC, THD_automask(), THD_countmask(), THD_filename_ok(), THD_is_file(), THD_mask_dilate(), THD_open_dataset(), tross_Copy_History(), tross_Make_History(), and var.
00025 { 00026 THD_3dim_dataset *dset , *oset=NULL , *tset=NULL ; 00027 int nvals , iv , nxyz , ii,jj,kk , iarg , kz,kzold ; 00028 float cut1=2.5,cut2=4.0 , sq2p , fq , val , fsig ; 00029 MRI_IMAGE *flim ; 00030 float *far, *dar , *var , *fitar ; 00031 char *prefix="despike" , *tprefix=NULL ; 00032 00033 int corder=-1 , nref , ignore=0 , polort=2 , nuse , nomask=0 ; 00034 int nspike, nbig, nproc ; 00035 float **ref ; 00036 float *fit , *ssp , snew , c21,ic21 , pspike,pbig , cls ; 00037 short *sar , *qar ; 00038 byte *tar , *mask=NULL ; 00039 float *zar , *yar ; 00040 int datum ; 00041 00042 /*----- Read command line -----*/ 00043 00044 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){ 00045 printf("Usage: 3dDespike [options] dataset\n" 00046 "Removes 'spikes' from the 3D+time input dataset and writes\n" 00047 "a new dataset with the spike values replaced by something\n" 00048 "more pleasing.\n" 00049 "\n" 00050 "Method:\n" 00051 " * L1 fit a smooth-ish curve to each voxel time series\n" 00052 " [see -corder option for description of the curve].\n" 00053 " * Compute the MAD of the difference between the curve and\n" 00054 " the data time series (the residuals).\n" 00055 " * Estimate the standard deviation 'sigma' of the residuals\n" 00056 " as sqrt(PI/2)*MAD.\n" 00057 " * For each voxel value, define s = (value-curve)/sigma.\n" 00058 " * Values with s > c1 are replaced with a value that yields\n" 00059 " a modified s' = c1+(c2-c1)*tanh((s-c1)/(c2-c1)).\n" 00060 " * c1 is the threshold value of s for a 'spike' [default c1=2.5].\n" 00061 " * c2 is the upper range of the allowed deviation from the curve:\n" 00062 " s=[c1..infinity) is mapped to s'=[c1..c2) [default c2=4].\n" 00063 "\n" 00064 "Options:\n" 00065 " -ignore I = Ignore the first I points in the time series:\n" 00066 " these values will just be copied to the\n" 00067 " output dataset [default I=0].\n" 00068 " -corder L = Set the curve fit order to L:\n" 00069 " the curve that is fit to voxel data v(t) is\n" 00070 "\n" 00071 " k=L [ (2*PI*k*t) (2*PI*k*t) ]\n" 00072 " f(t) = a+b*t+c*t*t + SUM [ d * sin(--------) + e * cos(--------) ]\n" 00073 " k=1 [ k ( T ) k ( T ) ]\n" 00074 "\n" 00075 " where T = duration of time series;\n" 00076 " the a,b,c,d,e parameters are chosen to minimize\n" 00077 " the sum over t of |v(t)-f(t)| (L1 regression);\n" 00078 " this type of fitting is is insensitive to large\n" 00079 " spikes in the data. The default value of L is\n" 00080 " NT/30, where NT = number of time points.\n" 00081 "\n" 00082 " -cut c1 c2 = Alter default values for the spike cut values\n" 00083 " [default c1=2.5, c2=4.0].\n" 00084 " -prefix pp = Save de-spiked dataset with prefix 'pp'\n" 00085 " [default pp='despike']\n" 00086 " -ssave ttt = Save 'spikiness' measure s for each voxel into a\n" 00087 " 3D+time dataset with prefix 'ttt' [default=no save]\n" 00088 " -nomask = Process all voxels\n" 00089 " [default=use a mask of high-intensity voxels, ]\n" 00090 " [as created via '3dAutomask -dilate 4 dataset'].\n" 00091 "\n" 00092 "Caveats:\n" 00093 "* Despiking may interfere with image registration, since head\n" 00094 " movement may produce 'spikes' at the edge of the brain, and\n" 00095 " this information would be used in the registration process.\n" 00096 " This possibility has not been explored.\n" 00097 "* Check your data visually before and after despiking and\n" 00098 " registration!\n" 00099 " [Hint: open 2 AFNI controllers, and turn Time Lock on.]\n" 00100 ) ; 00101 exit(0) ; 00102 } 00103 00104 /** AFNI package setup and logging **/ 00105 00106 mainENTRY("3dDespike main"); machdep(); AFNI_logger("3dDespike",argc,argv); 00107 PRINT_VERSION("3dDespike") ; 00108 00109 /** parse options **/ 00110 00111 iarg = 1 ; 00112 while( iarg < argc && argv[iarg][0] == '-' ){ 00113 00114 /** don't use masking **/ 00115 00116 if( strcmp(argv[iarg],"-nomask") == 0 ){ 00117 nomask = 1 ; 00118 iarg++ ; continue ; 00119 } 00120 00121 /** output dataset prefix **/ 00122 00123 if( strcmp(argv[iarg],"-prefix") == 0 ){ 00124 prefix = argv[++iarg] ; 00125 if( !THD_filename_ok(prefix) ){ 00126 fprintf(stderr,"** ERROR: -prefix is not good!\n"); exit(1); 00127 } 00128 iarg++ ; continue ; 00129 } 00130 00131 /** ratio dataset prefix **/ 00132 00133 if( strcmp(argv[iarg],"-ssave") == 0 ){ 00134 tprefix = argv[++iarg] ; 00135 if( !THD_filename_ok(tprefix) ){ 00136 fprintf(stderr,"** ERROR: -ssave prefix is not good!\n"); exit(1); 00137 } 00138 iarg++ ; continue ; 00139 } 00140 00141 /** trigonometric polynomial order **/ 00142 00143 if( strcmp(argv[iarg],"-corder") == 0 ){ 00144 corder = strtol( argv[++iarg] , NULL , 10 ) ; 00145 if( corder < 0 ){ 00146 fprintf(stderr,"** Illegal value of -corder!\n"); exit(1); 00147 } 00148 iarg++ ; continue ; 00149 } 00150 00151 /** how much to ignore at start **/ 00152 00153 if( strcmp(argv[iarg],"-ignore") == 0 ){ 00154 ignore = strtol( argv[++iarg] , NULL , 10 ) ; 00155 if( ignore < 0 ){ 00156 fprintf(stderr,"** Illegal value of -ignore!\n"); exit(1); 00157 } 00158 iarg++ ; continue ; 00159 } 00160 00161 /** thresholds for s ratio **/ 00162 00163 if( strcmp(argv[iarg],"-cut") == 0 ){ 00164 cut1 = strtod( argv[++iarg] , NULL ) ; 00165 cut2 = strtod( argv[++iarg] , NULL ) ; 00166 if( cut1 < 1.0 || cut2 < cut1+0.5 ){ 00167 fprintf(stderr,"** Illegal values after -cut!\n"); exit(1); 00168 } 00169 iarg++ ; continue ; 00170 } 00171 00172 fprintf(stderr,"** Unknown option: %s\n",argv[iarg]) ; exit(1) ; 00173 } 00174 00175 c21 = cut2-cut1 ; ic21 = 1.0/c21 ; 00176 00177 /*----- read input dataset -----*/ 00178 00179 if( iarg >= argc ){ 00180 fprintf(stderr,"** No input dataset!?\n"); exit(1); 00181 } 00182 00183 dset = THD_open_dataset( argv[iarg] ) ; 00184 if( !ISVALID_DSET(dset) ){ 00185 fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]); exit(1); 00186 } 00187 datum = DSET_BRICK_TYPE(dset,0) ; 00188 if( (datum != MRI_short && datum != MRI_float) || !DSET_datum_constant(dset) ){ 00189 fprintf(stderr,"** Can't process non-short, non-float dataset!\n") ; exit(1) ; 00190 } 00191 fprintf(stderr,"Input data type = %s\n",MRI_TYPE_name[datum]) ; 00192 nvals = DSET_NUM_TIMES(dset) ; nuse = nvals - ignore ; 00193 if( nuse < 15 ){ 00194 fprintf(stderr,"** Can't use dataset with < 15 time points per voxel!\n") ; 00195 exit(1) ; 00196 } 00197 fprintf(stderr,"++ ignoring first %d time points, using last %d\n",ignore,nuse); 00198 if( corder > 0 && 4*corder+2 > nuse ){ 00199 fprintf(stderr,"** -corder %d is too big for NT=%d!\n",corder,nvals) ; 00200 exit(1) ; 00201 } else if( corder < 0 ){ 00202 corder = rint(nuse/30.0) ; 00203 fprintf(stderr,"++ using %d time points => -corder %d\n",nuse,corder) ; 00204 } else { 00205 fprintf(stderr,"++ -corder %d set from command line\n",corder) ; 00206 } 00207 nxyz = DSET_NVOX(dset) ; 00208 fprintf(stderr,"++ Loading dataset %s\n",argv[iarg]) ; 00209 DSET_load(dset) ; 00210 if( !DSET_LOADED(dset) ){ 00211 fprintf(stderr,"** Can't read dataset from disk!\n") ; exit(1) ; 00212 } 00213 00214 /*-- create automask --*/ 00215 00216 if( !nomask ){ 00217 mask = THD_automask( dset ) ; 00218 for( ii=0 ; ii < 4 ; ii++ ) 00219 THD_mask_dilate( DSET_NX(dset), DSET_NY(dset), DSET_NZ(dset), mask, 3 ) ; 00220 ii = THD_countmask( DSET_NVOX(dset) , mask ) ; 00221 fprintf(stderr,"++ %d voxels in the mask [out of %d in dataset]\n",ii,DSET_NVOX(dset)) ; 00222 } else { 00223 fprintf(stderr,"++ processing all %d voxels in dataset\n",DSET_NVOX(dset)) ; 00224 } 00225 00226 /*-- create empty despiked dataset --*/ 00227 00228 oset = EDIT_empty_copy( dset ) ; 00229 EDIT_dset_items( oset , 00230 ADN_prefix , prefix , 00231 ADN_brick_fac , NULL , 00232 ADN_datum_all , datum , 00233 ADN_none ) ; 00234 00235 if( THD_is_file(DSET_HEADNAME(oset)) ){ 00236 fprintf(stderr,"** ERROR: output dataset already exists: %s\n",DSET_HEADNAME(oset)); 00237 exit(1); 00238 } 00239 00240 tross_Copy_History( oset , dset ) ; 00241 tross_Make_History( "3dDespike" , argc , argv , oset ) ; 00242 00243 /* create bricks (will be filled with zeros) */ 00244 00245 for( iv=0 ; iv < nvals ; iv++ ) 00246 EDIT_substitute_brick( oset , iv , datum , NULL ) ; 00247 00248 /* copy the ignored bricks */ 00249 00250 switch( datum ){ 00251 case MRI_short: 00252 for( iv=0 ; iv < ignore ; iv++ ){ 00253 sar = DSET_ARRAY(oset,iv) ; 00254 qar = DSET_ARRAY(dset,iv) ; 00255 memcpy( sar , qar , DSET_BRICK_BYTES(dset,iv) ) ; 00256 DSET_unload_one(dset,iv) ; 00257 } 00258 break ; 00259 case MRI_float: 00260 for( iv=0 ; iv < ignore ; iv++ ){ 00261 zar = DSET_ARRAY(oset,iv) ; 00262 yar = DSET_ARRAY(dset,iv) ; 00263 memcpy( zar , yar , DSET_BRICK_BYTES(dset,iv) ) ; 00264 DSET_unload_one(dset,iv) ; 00265 } 00266 break ; 00267 } 00268 00269 /*-- setup to save a threshold statistic dataset, if desired --*/ 00270 00271 if( tprefix != NULL ){ 00272 float *fac ; 00273 tset = EDIT_empty_copy( dset ) ; 00274 fac = (float *) malloc( sizeof(float) * nvals ) ; 00275 for( ii=0 ; ii < nvals ; ii++ ) fac[ii] = TFAC ; 00276 EDIT_dset_items( tset , 00277 ADN_prefix , tprefix , 00278 ADN_brick_fac , fac , 00279 ADN_datum_all , MRI_byte , 00280 ADN_func_type , FUNC_FIM_TYPE , 00281 ADN_none ) ; 00282 free(fac) ; 00283 00284 tross_Copy_History( tset , dset ) ; 00285 tross_Make_History( "3dDespike" , argc , argv , tset ) ; 00286 00287 if( THD_is_file(DSET_HEADNAME(tset)) ){ 00288 fprintf(stderr,"** ERROR: -ssave dataset already exists!\n"); exit(1); 00289 } 00290 00291 tross_Copy_History( tset , dset ) ; 00292 tross_Make_History( "3dDespike" , argc , argv , tset ) ; 00293 00294 for( iv=0 ; iv < nvals ; iv++ ) 00295 EDIT_substitute_brick( tset , iv , MRI_byte , NULL ) ; 00296 } 00297 00298 /*-- setup to find spikes --*/ 00299 00300 sq2p = sqrt(0.5*PI) ; 00301 var = (float *) malloc( sizeof(float) * nvals ) ; 00302 far = (float *) malloc( sizeof(float) * nvals ) ; 00303 dar = (float *) malloc( sizeof(float) * nvals ) ; 00304 fitar = (float *) malloc( sizeof(float) * nvals ) ; 00305 ssp = (float *) malloc( sizeof(float) * nvals ) ; 00306 00307 /* make ref functions */ 00308 00309 nref = 2*corder+3 ; 00310 ref = (float **) malloc( sizeof(float *) * nref ) ; 00311 for( jj=0 ; jj < nref ; jj++ ) 00312 ref[jj] = (float *) malloc( sizeof(float) * nuse ) ; 00313 00314 fit = (float *) malloc( sizeof(float) * nref ) ; 00315 00316 /* r(t) = 1 */ 00317 00318 for( iv=0 ; iv < nuse ; iv++ ) ref[0][iv] = 1.0 ; 00319 jj = 1 ; 00320 00321 /* r(t) = t - tmid */ 00322 00323 { float tm = 0.5 * (nuse-1.0) ; float fac = 2.0 / nuse ; 00324 for( iv=0 ; iv < nuse ; iv++ ) ref[1][iv] = (iv-tm)*fac ; 00325 jj = 2 ; 00326 00327 /* r(t) = (t-tmid)**jj */ 00328 00329 for( ; jj <= polort ; jj++ ) 00330 for( iv=0 ; iv < nuse ; iv++ ) 00331 ref[jj][iv] = pow( (iv-tm)*fac , (double)jj ) ; 00332 } 00333 00334 for( kk=1 ; kk <= corder ; kk++ ){ 00335 fq = (2.0*PI*kk)/nuse ; 00336 00337 /* r(t) = sin(2*PI*k*t/N) */ 00338 00339 for( iv=0 ; iv < nuse ; iv++ ) 00340 ref[jj][iv] = sin(fq*iv) ; 00341 jj++ ; 00342 00343 /* r(t) = cos(2*PI*k*t/N) */ 00344 00345 for( iv=0 ; iv < nuse ; iv++ ) 00346 ref[jj][iv] = cos(fq*iv) ; 00347 jj++ ; 00348 } 00349 00350 /*--- loop over voxels and do work ---*/ 00351 00352 fprintf(stderr,"++ edit thresholds: %.1f .. %.1f standard deviations\n",cut1,cut2) ; 00353 fprintf(stderr,"++ [ %.4f%% .. %.4f%% of normal distribution]\n", 00354 200.0*qg(cut1) , 200.0*qg(cut2) ) ; 00355 fprintf(stderr,"++ %d slices to process\n",DSET_NZ(dset)) ; 00356 kzold = -1 ; 00357 nspike = 0 ; nbig = 0 ; nproc = 0 ; 00358 for( ii=0 ; ii < nxyz ; ii++ ){ /* ii = voxel index */ 00359 00360 if( mask != NULL && mask[ii] == 0 ) continue ; /* skip this voxel */ 00361 00362 kz = DSET_index_to_kz(dset,ii) ; /* starting a new slice */ 00363 if( kz != kzold ){ 00364 fprintf(stderr, "++ start slice %2d",kz ) ; 00365 if( nproc > 0 ){ 00366 pspike = (100.0*nspike)/nproc ; 00367 pbig = (100.0*nbig )/nproc ; 00368 fprintf(stderr,"; done %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]", 00369 nproc,nspike,pspike,nbig,pbig ) ; 00370 } 00371 fprintf(stderr,"\n") ; 00372 kzold = kz ; 00373 } 00374 00375 /*** extract ii-th time series into far[] ***/ 00376 00377 switch( datum ){ 00378 case MRI_short: 00379 for( iv=0 ; iv < nuse ; iv++ ){ 00380 qar = DSET_ARRAY(dset,iv+ignore) ; /* skip ignored data */ 00381 far[iv] = (float) qar[ii] ; 00382 } 00383 break ; 00384 case MRI_float: 00385 for( iv=0 ; iv < nuse ; iv++ ){ 00386 zar = DSET_ARRAY(dset,iv+ignore) ; 00387 far[iv] = zar[ii] ; 00388 } 00389 break ; 00390 } 00391 memcpy(dar,far,sizeof(float)*nuse) ; /* copy time series into dar[] */ 00392 00393 /*** solve for L1 fit ***/ 00394 00395 cls = cl1_solve( nuse , nref , far , ref , fit,0 ); 00396 00397 if( cls < 0.0 ){ /* fit failed! */ 00398 #if 0 00399 fprintf(stderr,"curve fit fails at voxel %d %d %d\n", 00400 DSET_index_to_ix(dset,ii) , 00401 DSET_index_to_jy(dset,ii) , 00402 DSET_index_to_kz(dset,ii) ) ; 00403 #endif 00404 continue ; /* skip this voxel */ 00405 } 00406 00407 for( iv=0 ; iv < nuse ; iv++ ){ /* detrend */ 00408 val = fit[0] 00409 + fit[1]*ref[1][iv] /* quadratic part of curve fit */ 00410 + fit[2]*ref[2][iv] ; 00411 for( jj=3 ; jj < nref ; jj++ ) /* rest of curve fit */ 00412 val += fit[jj] * ref[jj][iv] ; 00413 00414 fitar[iv] = val ; /* save curve value */ 00415 var[iv] = dar[iv]-val ; /* remove fitted value = resid */ 00416 far[iv] = fabs(var[iv]) ; /* abs value of resid */ 00417 } 00418 00419 /*** compute estimate standard deviation of detrended data ***/ 00420 00421 fsig = sq2p * qmed_float(nuse,far) ; /* also mangles far array */ 00422 00423 /*** process time series for spikes, editing data in dar[] ***/ 00424 00425 if( fsig > 0.0 ){ /* data wasn't fit perfectly */ 00426 00427 /* find spikiness for each point in time */ 00428 00429 fq = 1.0 / fsig ; 00430 for( iv=0 ; iv < nuse ; iv++ ){ 00431 ssp[iv] = fq * var[iv] ; /* spikiness s = how many sigma out */ 00432 } 00433 00434 /* save spikiness in -ssave datset */ 00435 00436 if( tset != NULL ){ 00437 for( iv=0 ; iv < nuse ; iv++ ){ 00438 tar = DSET_ARRAY(tset,iv+ignore) ; 00439 snew = ITFAC*fabs(ssp[iv]) ; /* scale for byte storage */ 00440 tar[ii] = BYTEIZE(snew) ; /* cf. mrilib.h */ 00441 } 00442 } 00443 00444 /* process values of |s| > cut1, editing dar[] */ 00445 00446 for( iv=0 ; iv < nuse ; iv++ ){ 00447 if( ssp[iv] > cut1 ){ 00448 snew = cut1 + c21*mytanh((ssp[iv]-cut1)*ic21) ; /* edit s down */ 00449 dar[iv] = fitar[iv] + snew*fsig ; 00450 nspike++ ; if( ssp[iv] > cut2 ) nbig++ ; /* count edits */ 00451 } else if( ssp[iv] < -cut1 ){ 00452 snew = -cut1 + c21*mytanh((ssp[iv]+cut1)*ic21) ; /* edit s up */ 00453 dar[iv] = fitar[iv] + snew*fsig ; 00454 nspike++ ; if( ssp[iv] < -cut2 ) nbig++ ; 00455 } 00456 } 00457 nproc += nuse ; /* number data points processed */ 00458 00459 } /* end of processing time series when fsig is positive */ 00460 00461 /* put dar[] time series (possibly edited above) into output bricks */ 00462 00463 switch( datum ){ 00464 case MRI_short: 00465 for( iv=0 ; iv < nuse ; iv++ ){ 00466 sar = DSET_ARRAY(oset,iv+ignore) ; /* output brick */ 00467 sar[ii] = (short) dar[iv] ; /* original or mutated data */ 00468 } 00469 break ; 00470 case MRI_float: 00471 for( iv=0 ; iv < nuse ; iv++ ){ 00472 zar = DSET_ARRAY(oset,iv+ignore) ; /* output brick */ 00473 zar[ii] = dar[iv] ; /* original or mutated data */ 00474 } 00475 break ; 00476 } 00477 00478 } /* end of loop over voxels #ii */ 00479 00480 /*--- finish up ---*/ 00481 00482 DSET_delete(dset) ; /* delete input dataset */ 00483 00484 if( nproc > 0 ){ 00485 pspike = (100.0*nspike)/nproc ; 00486 pbig = (100.0*nbig )/nproc ; 00487 fprintf(stderr,"++ FINAL: %d data points, %d edits [%.3f%%], %d big edits [%.3f%%]\n", 00488 nproc,nspike,pspike,nbig,pbig ) ; 00489 } else { 00490 fprintf(stderr,"++ FINAL: no good voxels found to process!!!!\n") ; 00491 } 00492 00493 /* write results */ 00494 00495 DSET_write(oset) ; 00496 fprintf(stderr,"++ Wrote output dataset %s\n",DSET_BRIKNAME(oset)) ; 00497 DSET_delete(oset) ; 00498 00499 if( tset != NULL ){ 00500 DSET_write(tset) ; 00501 fprintf(stderr,"++ Wrote -ssave dataset %s\n",DSET_BRIKNAME(tset)) ; 00502 } 00503 00504 exit(0) ; 00505 } |
|
Definition at line 13 of file 3dDespike.c. Referenced by main().
00014 { 00015 register float ex , exi ; 00016 if( x > 7.0 ) return 1.0 ; /* 03 Sep: check for stupid inputs */ 00017 else if( x < -7.0 ) return -1.0 ; 00018 ex = exp(x) ; exi = 1.0/ex ; 00019 return (ex-exi)/(ex+exi) ; 00020 } |