00001 /*****************************************************************************
00002   3dretroicor.c: AFNI program to perform Retrospective Image
00003   Correction for physiological motion effects, using a slightly
00004   modified version of the RETROICOR algorithm described in:
00005 
00006     Glover, G. H., Li, T., & Ress, D. (2000). Image-based method for
00007   retrospective correction of physiological motion effects in fMRI:
00008   RETROICOR. Magnetic Resonance in Medicine, 44, 162-167.
00009 
00010   Adapted from plug_retroicor.c.
00011 
00012   Fred Tam
00013   Sunnybrook & Women's College Health Sciences Centre
00014   August 1, 2002
00015 
00016   Copyright (C) 2002 Sunnybrook & Women's College Health Sciences Centre,
00017   Toronto, ON, Canada. As part of the AFNI software package, this software
00018   is distributed under the GNU General Public License, Version 2.
00019   See the file README.copyright for details.
00020 *****************************************************************************/
00021 
00022 /*****************************************************************************
00023    Major portions of this software are copyrighted by the Medical College
00024    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00025    License, Version 2.  See the file README.Copyright for details.
00026 ******************************************************************************/
00027 
00028 #include "mrilib.h"
00029 
00030 #include "retroicor.h"
00031 
00032 #define TRIC_C_DEF_THRESHOLD 1
00033 #define TRIC_R_DEF_WINSIZE 20
00034 #define TRIC_M_DEF_ORDER 2
00035 #define TRIC_I_DEF_IGNORE 0
00036 #define TRIC_O_DEF_NEWPREFIX "retroicor"
00037 
00038 /*****************************************************************************
00039   Print out a help message to stdout
00040 ******************************************************************************/
00041 
00042 void TRIC_printhelp(void)
00043 {
00044     printf("Usage: 3dretroicor [options] dataset\n"
00045            "\n"
00046            "Performs Retrospective Image Correction for physiological\n"
00047            "motion effects, using a slightly modified version of the\n"
00048            "RETROICOR algorithm described in:\n"
00049            "\n"
00050            "  Glover, G. H., Li, T., & Ress, D. (2000). Image-based method\n"
00051            "for retrospective correction of physiological motion effects in\n"
00052            "fMRI: RETROICOR. Magnetic Resonance in Medicine, 44, 162-167.\n"
00053            "\n"
00054            "Options (defaults in []'s):\n"
00055            "\n"
00056            " -ignore    = The number of initial timepoints to ignore in the\n"
00057            "              input (These points will be passed through\n"
00058            "              uncorrected) [0]\n"
00059            " -prefix    = Prefix for new, corrected dataset [retroicor]\n"
00060            "\n"
00061            " -card      = 1D cardiac data file for cardiac correction\n"
00062            " -cardphase = Filename for 1D cardiac phase output\n"
00063            " -threshold = Threshold for detection of R-wave peaks in input\n"
00064            "              (Make sure it's above the background noise level;\n"
00065            "              Try 3/4 or 4/5 times range plus minimum) [1]\n"
00066            "\n"
00067            " -resp      = 1D respiratory waveform data for correction\n"
00068            " -respphase = Filename for 1D resp phase output\n"
00069            /*-- removed winsize ui
00070            " -window    = Window size for input point estimate of slope\n"
00071            "              (Try 1/3 or 1/2 times respiratory sampling rate in\n"
00072            "              Hz) [50]\n"
00073            removed winsize ui --*/
00074            "\n"
00075            " -order     = The order of the correction (2 is typical;\n"
00076            "              higher-order terms yield little improvement\n"
00077            "              according to Glover et al.) [2]\n"
00078            "\n"
00079            " -help      = Display this message and stop (must be first arg)\n"
00080            "\n"
00081            "Dataset: 3D+time dataset to process\n"
00082            "\n"
00083            "** The input dataset and at least one of -card and -resp are\n"
00084            "    required.\n"
00085            "\n"
00086            "NOTES\n"
00087            "-----\n"
00088            "\n"
00089            "The durations of the physiological inputs are assumed to equal\n"
00090            "the duration of the dataset. Any constant sampling rate may be\n"
00091            "used, but 40 Hz seems to be acceptable. This program's cardiac\n"
00092            "peak detection algorithm is rather simplistic, so you might try\n"
00093            "using the scanner's cardiac gating output (transform it to a\n"
00094            "spike wave if necessary).\n"
00095            "\n"
00096            "This program uses slice timing information embedded in the\n"
00097            "dataset to estimate the proper cardiac/respiratory phase for\n"
00098            "each slice. It makes sense to run this program before any\n"
00099            "program that may destroy the slice timings (e.g. 3dvolreg for\n"
00100            "motion correction).\n"
00101            "\n"
00102            "Author -- Fred Tam, August 2002\n"
00103            "\n"
00104            MASTER_SHORTHELP_STRING
00105            );
00106 }
00107 
00108 /*****************************************************************************
00109   Monolithic Mainline
00110 ******************************************************************************/
00111 
00112 int main(int argc, char * argv[])
00113 {
00114    char cphase1d[256] = "\0", rphase1d[256] = "\0";
00115    char new_prefix[THD_MAX_PREFIX] = TRIC_O_DEF_NEWPREFIX;
00116    MCW_idcode * idc ;
00117    THD_3dim_dataset * dset , * new_dset;
00118    double * avg = NULL;
00119    double * ca , * cb, * ra, * rb;
00120    MRI_IMAGE * card = NULL, * resp = NULL;
00121    MRI_IMAGE * cardphase = NULL, * respphase = NULL;
00122    float threshold = TRIC_C_DEF_THRESHOLD;
00123    int ignore = TRIC_I_DEF_IGNORE;
00124    int M = TRIC_M_DEF_ORDER;
00125    int winsize = TRIC_R_DEF_WINSIZE;
00126    float tr;
00127    int ival, nvals;
00128    FILE * fp;
00129    float * cpdata, * rpdata;
00130    int argi = 1;
00131 
00132    /*-------------------------------------------------------------*/
00133    /*----- Check arguments to see if they are reasonable-ish -----*/
00134 
00135    if (argc < 2 || strcmp(argv[1], "-help") == 0) {
00136        TRIC_printhelp();
00137        exit(-1);
00138    }
00139 
00140    mainENTRY("3dretroicor main"); PRINT_VERSION("3dretroicor") ;
00141    machdep();
00142    AFNI_logger("3dretroicor", argc, argv);
00143 
00144    /* Iterate over commandline args */
00145    while (argi < argc && argv[argi][0] == '-') {
00146 
00147        /*-- ignore --*/
00148 
00149        if (strcmp(argv[argi], "-ignore") == 0) {
00150            if (argi + 1 >= argc) {
00151                fprintf(stderr, "*** -ignore needs an argument!\n");
00152                exit(1);
00153            }
00154            argi += 1;
00155            ignore = atoi(argv[argi]);
00156            if (ignore < 0) {
00157                fprintf(stderr, "*** %i is not a valid number to ignore!\n",
00158                        ignore);
00159                exit(1);
00160            }
00161            argi += 1;
00162            continue;  /* Skip the rest of the loop and go to the next arg */
00163        }
00164 
00165        /*-- prefix --*/
00166 
00167        if (strcmp(argv[argi], "-prefix") == 0) {
00168            if (argi + 1 >= argc) {
00169                fprintf(stderr, "*** -prefix needs an argument!\n");
00170                exit(1);
00171            }
00172            argi += 1;
00173            MCW_strncpy(new_prefix, argv[argi], THD_MAX_PREFIX);
00174            if (! THD_filename_ok(new_prefix)) {
00175                fprintf(stderr, "*** %s is not a valid prefix!\n", new_prefix);
00176                exit(1);
00177            }
00178            argi += 1;
00179            continue;  /* Skip the rest of the loop and go to the next arg */
00180        }
00181 
00182        /*-- card --*/
00183 
00184        if (strcmp(argv[argi], "-card") == 0) {
00185            if (argi + 1 >= argc) {
00186                fprintf(stderr, "*** -card needs an argument!\n");
00187                exit(1);
00188            }
00189            argi += 1;
00190            card = mri_read_1D(argv[argi]);
00191            if (card == NULL) {
00192                fprintf(stderr, "*** Can't read -card %s\n", argv[argi]);
00193                exit(1);
00194            }
00195            argi += 1;
00196            continue;  /* Skip the rest of the loop and go to the next arg */
00197        }
00198 
00199        /*-- cardphase --*/
00200 
00201        if (strcmp(argv[argi], "-cardphase") == 0) {
00202            if (argi + 1 >= argc) {
00203                fprintf(stderr, "*** -cardphase needs an argument!\n");
00204                exit(1);
00205            }
00206            argi += 1;
00207            MCW_strncpy(cphase1d, argv[argi], 256);
00208            if (! THD_filename_ok(cphase1d)) {
00209                fprintf(stderr, "*** Bad name argument for -cardphase\n");
00210                exit(1);
00211            }
00212            argi += 1;
00213            continue;  /* Skip the rest of the loop and go to the next arg */
00214        }
00215 
00216        /*-- threshold --*/
00217 
00218        if (strcmp(argv[argi], "-threshold") == 0) {
00219            if (argi + 1 >= argc) {
00220                fprintf(stderr, "*** -threshold needs an argument!\n");
00221                exit(1);
00222            }
00223            argi += 1;
00224            threshold = atof(argv[argi]);
00225            argi += 1;
00226            continue;  /* Skip the rest of the loop and go to the next arg */
00227        }
00228 
00229        /*-- resp --*/
00230 
00231        if (strcmp(argv[argi], "-resp") == 0) {
00232            if (argi + 1 >= argc) {
00233                fprintf(stderr, "*** -resp needs an argument!\n");
00234                exit(1);
00235            }
00236            argi += 1;
00237            resp = mri_read_1D(argv[argi]);
00238            if (resp == NULL) {
00239                fprintf(stderr, "*** Can't read -resp %s\n", argv[argi]);
00240                exit(1);
00241            }
00242            argi += 1;
00243            continue;  /* Skip the rest of the loop and go to the next arg */
00244        }
00245 
00246        /*-- respphase --*/
00247 
00248        if (strcmp(argv[argi], "-respphase") == 0) {
00249            if (argi + 1 >= argc) {
00250                fprintf(stderr, "*** -respphase needs an argument!\n");
00251                exit(1);
00252            }
00253            argi += 1;
00254            MCW_strncpy(rphase1d, argv[argi], 256);
00255            if (! THD_filename_ok(rphase1d)) {
00256                fprintf(stderr, "*** Bad name argument for -respphase\n");
00257                exit(1);
00258            }
00259            argi += 1;
00260            continue;  /* Skip the rest of the loop and go to the next arg */
00261        }
00262 
00263        /*-- window --*/
00264 
00265        /*-- removed winsize ui --*/
00266 #if 0
00267        if (strcmp(argv[argi], "-window") == 0) {
00268            if (argi + 1 >= argc) {
00269                fprintf(stderr, "*** -window needs an argument!\n");
00270                exit(1);
00271            }
00272            argi += 1;
00273            winsize = atoi(argv[argi]);
00274            argi += 1;
00275            continue;  /* Skip the rest of the loop and go to the next arg */
00276        }
00277 #endif
00278        /*-- removed winsize ui --*/
00279 
00280        /*-- order --*/
00281 
00282        if (strcmp(argv[argi], "-order") == 0) {
00283            if (argi + 1 >= argc) {
00284                fprintf(stderr, "*** -order needs an argument!\n");
00285                exit(1);
00286            }
00287            argi += 1;
00288            M = atoi(argv[argi]);
00289            if (M < 1) {
00290                fprintf(stderr, "*** %i is not a valid order number\n", M);
00291                exit(1);
00292            }
00293            argi += 1;
00294            continue;  /* Skip the rest of the loop and go to the next arg */
00295        }
00296    } /* End iterating over commandline args */
00297 
00298    /* Check that at least one of Cardiac and Resp were selected */
00299    if (card == NULL && resp == NULL) {
00300        fprintf(stderr, "*** Need at least one correction (-card, -resp)\n");
00301        exit(1);
00302    }
00303 
00304    /*-- Open and check input dataset --*/
00305 
00306    if (argi >= argc) {
00307        fprintf(stderr, "*** No input dataset!?\n");
00308        exit(1);
00309    } else if (argi < argc - 1) {
00310        fprintf(stderr, "*** Too many input datasets?!\n");
00311        exit(1);
00312    }
00313 
00314    dset = THD_open_dataset(argv[argi]);
00315    if (! ISVALID_3DIM_DATASET(dset)) {
00316        fprintf(stderr, "*** Can't open dataset %s\n", argv[argi]);
00317        exit(1);
00318    }
00319    if (DSET_NUM_TIMES(dset) < 2) {
00320        fprintf(stderr, "*** Input dataset is not 3D+time!\n");
00321        exit(1);
00322    }
00323 
00324    /*----------------------------------------------------------*/
00325    /*---------- At this point, the inputs are OK-ish ----------*/
00326 
00327    /*-- copy the image data for editing in place --*/
00328 
00329    new_dset = EDIT_full_copy( dset , new_prefix );
00330    if( new_dset == NULL ) {
00331        fprintf(stderr, "*** Error copying dataset\n");
00332        exit(1);
00333    }
00334    tross_Copy_History(dset, new_dset); /* Copy and add to new_dset history */
00335    tross_Make_History("3dretroicor", argc, argv, new_dset);
00336    DSET_unload( dset ) ;  /* We won't need the old dataset anymore */
00337 
00338    /*-- calculate cardiac correction coefficients if requested --*/
00339 
00340    if (card != NULL) {
00341        /*-- convert cardiac waveform to phase --*/
00342        cardphase = RIC_ToCardiacPhase(card, threshold) ;
00343        if (cardphase == NULL) {
00344            THD_delete_3dim_dataset( new_dset , False ) ;
00345            fprintf(stderr, "*** Error transforming cardiac data\n");
00346            exit(2);
00347        }
00348 
00349        /*-- calculate dataset voxel means --*/
00350        avg = RIC_CalcVoxelMeans(new_dset, ignore);
00351        if (avg == NULL) {
00352            THD_delete_3dim_dataset( new_dset , False ) ;
00353            mri_free(cardphase);
00354            fprintf(stderr, "*** Error calculating dataset voxel means\n");
00355            exit(2);
00356        }
00357 
00358        /*-- calculate coefficients for each voxel --*/
00359        if (RIC_CalcCoeffAB(new_dset, cardphase, avg, &ca, &cb, M, ignore)
00360            != 0) {
00361 
00362            THD_delete_3dim_dataset( new_dset , False ) ;
00363            mri_free(cardphase);
00364            fprintf(stderr, "*** Error calculating cardiac a b coefficients\n");
00365            exit(2);
00366        }
00367    }
00368 
00369    /*-- calculate respiratory correction coefficients if requested --*/
00370 
00371    if (resp != NULL) {
00372        /*-- Set winsize to 1/2 sampling rate of resp in Hz --*/
00373        tr = new_dset->taxis->ttdel;
00374        switch (new_dset->taxis->units_type) {
00375        case UNITS_MSEC_TYPE: tr /= 1000; break;
00376        case UNITS_SEC_TYPE:  break;
00377        case UNITS_HZ_TYPE:   tr = 1 / tr; break;
00378        default:
00379            THD_delete_3dim_dataset( new_dset , False ) ;
00380            fprintf(stderr, "*** Bad time units type in dataset\n");
00381            exit(2);
00382        }
00383        winsize = ceil(resp->nx / (tr * DSET_NVALS(new_dset)) / 2.0);
00384 
00385        /*-- convert respiratory waveform to phase --*/
00386        respphase = RIC_ToRespPhase(resp, winsize) ;
00387        if (respphase == NULL) {
00388            THD_delete_3dim_dataset( new_dset , False ) ;
00389            fprintf(stderr, "*** Error transforming resp data\n");
00390            exit(2);
00391        }
00392 
00393        /*-- calculate dataset voxel means if not already done --*/
00394        if (avg == NULL) {
00395            avg = RIC_CalcVoxelMeans(new_dset, ignore);
00396            if (avg == NULL) {
00397                THD_delete_3dim_dataset( new_dset , False ) ;
00398                mri_free(respphase);
00399                fprintf(stderr, "*** Error calculating dataset voxel means2\n");
00400                exit(2);
00401            }
00402        }
00403 
00404        /*-- calculate coefficients for each voxel --*/
00405        if (RIC_CalcCoeffAB(new_dset, respphase, avg, &ra, &rb, M, ignore)
00406            != 0) {
00407 
00408            THD_delete_3dim_dataset( new_dset , False ) ;
00409            mri_free(respphase);
00410            fprintf(stderr, "*** Error calculating resp a, b coefficients\n");
00411            exit(2);
00412        }
00413    }
00414 
00415    /*-- do cardiac correction if requested --*/
00416 
00417    if (card != NULL) {
00418        /*-- correct the image data --*/
00419        if (RIC_CorrectDataset(new_dset, cardphase, ca, cb, M, ignore) != 0) {
00420            THD_delete_3dim_dataset( new_dset , False ) ;
00421            mri_free(cardphase);
00422            free(ca); free(cb);
00423            fprintf(stderr, "*** Error applying cardiac correction\n");
00424            exit(3);
00425        }
00426 
00427        /*-- if requested, write phase data to file and pass to AFNI --*/
00428        if ( THD_filename_ok(cphase1d) ) {
00429            /* Write the file */
00430            fp = fopen(cphase1d, "w");
00431            nvals = cardphase->nx;
00432            cpdata = MRI_FLOAT_PTR(cardphase);
00433            for (ival = 0; ival < nvals; ival += 1) {
00434                fprintf(fp, "%f\n", cpdata[ival]);
00435            }
00436            fclose(fp);
00437        }
00438 
00439        mri_free(cardphase);
00440        free(ca); free(cb); free(avg); avg = NULL;
00441    }
00442 
00443    /*-- do resp correction if requested --*/
00444 
00445    if (resp != NULL) {
00446        /*-- correct the image data --*/
00447        if (RIC_CorrectDataset(new_dset, respphase, ra, rb, M, ignore) != 0) {
00448            THD_delete_3dim_dataset( new_dset , False ) ;
00449            mri_free(respphase);
00450            free(ra); free(rb);
00451            fprintf(stderr, "*** Error applying resp correction\n");
00452            exit(3);
00453        }
00454 
00455        /*-- if requested, write phase data to file and pass to AFNI --*/
00456        if ( THD_filename_ok(rphase1d) ) {
00457            /* Write the file */
00458            fp = fopen(rphase1d, "w");
00459            nvals = respphase->nx;
00460            rpdata = MRI_FLOAT_PTR(respphase);
00461            for (ival = 0; ival < nvals; ival += 1) {
00462                fprintf(fp, "%f\n", rpdata[ival]);
00463            }
00464            fclose(fp);
00465        }
00466 
00467        mri_free(respphase);
00468        free(ra); free(rb); if (avg != NULL) free(avg);
00469    }
00470 
00471    /*-- write out new dataset --*/
00472 
00473    DSET_write(new_dset);
00474    fprintf(stderr,"++ output dataset: %s\n",DSET_BRIKNAME(new_dset)) ;
00475 
00476    /*-- done successfully!!! --*/
00477 
00478    exit(0);
00479 }