00001 /***********************************************************************
00002   plug_retroicor.c: AFNI plugin 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   Fred Tam
00011   Sunnybrook & Women's College Health Sciences Centre
00012   April 15, 2002
00013 
00014   Copyright (C) 2002 Sunnybrook & Women's College Health Sciences Centre,
00015   Toronto, ON, Canada. As part of the AFNI software package, this software
00016   is distributed under the GNU General Public License, Version 2.
00017   See the file README.copyright for details.
00018 ************************************************************************/
00019 
00020 /*****************************************************************************
00021    Major portions of this software are copyrighted by the Medical College
00022    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00023    License, Version 2.  See the file README.Copyright for details.
00024 ******************************************************************************/
00025 
00026 #include "afni.h"
00027 
00028 #ifndef ALLOW_PLUGINS
00029 #  error "Plugins not properly set up -- see machdep.h"
00030 #endif
00031 
00032 #include "retroicor.h"
00033 
00034 static char * PRIC_main( PLUGIN_interface * ) ;
00035 
00036 static char helpstring[] =
00037   " Purpose: Perform Retrospective Image Correction for physiological\n"
00038   "   motion effects, using a slightly modified version of the RETROICOR\n"
00039   "   algorithm described in:\n"
00040   " \n"
00041   "     Glover, G. H., Li, T., & Ress, D. (2000). Image-based method for\n"
00042   "   retrospective correction of physiological motion effects in fMRI:\n"
00043   "   RETROICOR. Magnetic Resonance in Medicine, 44, 162-167.\n"
00044   " \n"
00045   " Inputs:\n"
00046   " \n"
00047   "   Datasets: Input  = 3D+time dataset to process\n"
00048   "             Output = Prefix for new, corrected dataset\n"
00049   "             Ignore = Number of initial timepoints to ignore in input\n"
00050   "                      (These points will be passed through uncorrected)\n"
00051   " \n"
00052   "   Cardiac: Input     = 1D cardiac waveform data for correction\n"
00053   "            Output    = Filename for 1D cardiac phase output\n"
00054   "                        (Leave blank if you don't want this output)\n"
00055   "            Threshold = Threshold for detection of R-wave peaks in input\n"
00056   "                        (Make sure it's above the background noise level;\n"
00057   "                        Try 3/4 or 4/5 times range plus minimum)\n"
00058   " \n"
00059   "   Resp: Input  = 1D respiratory waveform data for correction\n"
00060   "         Output = Filename for 1D resp phase output\n"
00061   "                  (Leave blank if you don't want this output)\n"
00062 /*-- removed winsize ui
00063   "         Window = Window size for input point estimate of slope\n"
00064   "                  (Try 1/3 or 1/2 times respiratory sampling rate in Hz)\n"
00065 removed winsize ui --*/
00066   " \n"
00067   "   Params: Order = The order of the correction\n"
00068   "           (2 is typical; higher-order terms yield little improvement\n"
00069   "            according to Glover et al.)\n"
00070   " \n"
00071   "   ** The input and output datasets and at least one input for correction\n"
00072   "      (cardiac and/or resp) are required.\n"
00073   " \n"
00074   " NOTES\n"
00075   " -----\n"
00076   " \n"
00077   " The durations of the physiological inputs are assumed to equal the\n"
00078   " duration of the dataset. Any constant sampling rate may be used, but\n"
00079   " 40 Hz seems to be acceptable. This plugin's cardiac peak detection\n"
00080   " algorithm is rather simplistic, so you might try using the scanner's\n"
00081   " cardiac gating output (transform it to a spike wave if necessary).\n"
00082   " \n"
00083   " This plugin uses slice timing information embedded in the dataset to\n"
00084   " estimate the proper cardiac/respiratory phase for each slice. It makes\n"
00085   " sense to run this plugin before any program that may destroy the slice\n"
00086   " timings (e.g. 3dvolreg for motion correction).\n"
00087   " \n"
00088   " Author -- Fred Tam, April 2002"
00089 ;
00090 
00091 #define PRIC_C_DEF_THRESHOLD 10
00092 #define PRIC_R_DEF_WINSIZE 20
00093 #define PRIC_M_DEF_ORDER 2
00094 #define PRIC_I_DEF_IGNORE 0
00095 
00096 /***********************************************************************
00097    Set up the interface to the user
00098 ************************************************************************/
00099 
00100 PLUGIN_interface * PLUGIN_init( int ncall )
00101 {
00102    PLUGIN_interface * plint ;
00103 
00104    if( ncall > 0 ) return NULL ;  /* only one interface */
00105 
00106    /*-- set titles and call point --*/
00107 
00108    plint = PLUTO_new_interface( "RETROICOR" ,
00109                                 "Physio Correction of a 3D+time Dataset" ,
00110                                 helpstring ,
00111                                 PLUGIN_CALL_VIA_MENU , PRIC_main  ) ;
00112 
00113    PLUTO_add_hint( plint , "Physio Correction of a 3D+time Dataset" ) ;
00114 
00115    PLUTO_set_sequence( plint , "A:newdset:retroicor" ) ;
00116 
00117    /*-- first line of input: Datasets --*/
00118 
00119    PLUTO_add_option( plint , "Datasets" , "Datasets" , TRUE ) ;
00120 
00121    PLUTO_add_dataset( plint , "Input" ,
00122                                     ANAT_ALL_MASK , FUNC_FIM_MASK ,
00123                                     DIMEN_4D_MASK | BRICK_ALLREAL_MASK ) ;
00124    PLUTO_add_hint( plint , "Choose 3D+time input" ) ;
00125 
00126    PLUTO_add_string( plint , "Output" , 0, NULL , 19 ) ;
00127    PLUTO_add_hint( plint , "Prefix for corrected 3D+time output" ) ;
00128 
00129    PLUTO_add_number( plint , "Ignore" , 0, 100, 0, PRIC_I_DEF_IGNORE, TRUE ) ;
00130    PLUTO_add_hint( plint , "Number of initial input timepoints to ignore" ) ;
00131 
00132    /*-- second line of input: Cardiac --*/
00133 
00134    PLUTO_add_option( plint , "Cardiac", "Cardiac" , FALSE ) ;
00135 
00136    PLUTO_add_timeseries( plint , "Input" );
00137    PLUTO_add_hint( plint , "Choose 1D cardiac waveform input");
00138 
00139    PLUTO_add_string( plint , "Output" , 0, NULL , 19 ) ;
00140    PLUTO_add_hint( plint , "Filename for 1D cardiac phase output (optional)" );
00141 
00142    PLUTO_add_number( plint , "Threshold" , -1280, 1270, 1,
00143                      PRIC_C_DEF_THRESHOLD, TRUE ) ;
00144    PLUTO_add_hint( plint , "Threshold for input R-wave peak detection" ) ;
00145 
00146    /*-- third line of input: Resp --*/
00147 
00148    PLUTO_add_option( plint , "Resp", "Resp" , FALSE ) ;
00149 
00150    PLUTO_add_timeseries( plint , "Input" );
00151    PLUTO_add_hint( plint , "Choose 1D resp waveform input");
00152 
00153    PLUTO_add_string( plint , "Output" , 0, NULL , 19 ) ;
00154    PLUTO_add_hint( plint , "Filename for 1D resp phase output (optional)" ) ;
00155    /*-- removed winsize ui
00156    PLUTO_add_number( plint , "Window" , 2, 100, 0, PRIC_R_DEF_WINSIZE, TRUE ) ;
00157    PLUTO_add_hint( plint , "Window size for input point estimate of slope" ) ;
00158    removed winsize ui --*/
00159    /*-- fourth line of input: Params --*/
00160 
00161    PLUTO_add_option( plint , "Params", "Params" , FALSE ) ;
00162 
00163    PLUTO_add_number( plint , "Order" , 1, 5, 0, PRIC_M_DEF_ORDER, FALSE ) ;
00164    PLUTO_add_hint( plint , "Order of correction" ) ;
00165 
00166    return plint ;
00167 }
00168 
00169 /***************************************************************************
00170   Main routine for this plugin (will be called from AFNI).
00171 ****************************************************************************/
00172 
00173 static char * PRIC_main( PLUGIN_interface * plint )
00174 {
00175    char * tag , * new_prefix , * cphase1d = NULL, * rphase1d = NULL;
00176    MCW_idcode * idc ;
00177    THD_3dim_dataset * dset , * new_dset;
00178    double * avg = NULL;
00179    double * ca , * cb, * ra, * rb;
00180    MRI_IMAGE * card = NULL, * resp = NULL;
00181    MRI_IMAGE * cardphase = NULL, * respphase = NULL;
00182    float threshold, ignore_input, M_input;
00183    /*-- removed winsize ui
00184    , winsize_input;
00185    removed winsize ui --*/
00186    int ignore;
00187    int M = PRIC_M_DEF_ORDER;
00188    int winsize;
00189    float tr;
00190    int ival, nvals;
00191    FILE * fp;
00192    float * cpdata, * rpdata;
00193    char * histstring;
00194 
00195    /*--------------------------------------------------------------------*/
00196    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00197 
00198    if( plint == NULL )
00199       return "*********************\n"
00200              "PRIC_main: NULL input\n"
00201              "*********************"    ;
00202 
00203    /*-- first line of input: Datasets (required) --*/
00204 
00205    tag = PLUTO_get_optiontag(plint) ;
00206    if( tag==NULL || strcmp(tag,"Datasets") != 0 )
00207       return "*******************************\n"
00208              "PRIC_main: bad Input option tag\n"
00209              "*******************************"    ;
00210 
00211    idc  = PLUTO_get_idcode(plint) ;
00212    dset = PLUTO_find_dset(idc) ;
00213    if( dset == NULL )
00214       return "****************************\n"
00215              "PRIC_main: bad input dataset\n"
00216              "****************************"   ;
00217 
00218    new_prefix = PLUTO_get_string(plint) ;
00219    if( ! PLUTO_prefix_ok(new_prefix) )
00220        return "*********************\n"
00221               "PRIC_main: bad prefix\n"
00222               "*********************"   ;
00223 
00224    ignore_input = PLUTO_get_number(plint) ;
00225    ignore = rint(ignore_input);
00226    if( ignore_input == BAD_NUMBER || ignore < 0 )
00227       return "***************************\n"
00228              "PRIC_main: bad ignore input\n"
00229              "***************************"   ;
00230 
00231    /*-- second line of input: Cardiac --*/
00232 
00233    tag = PLUTO_get_optiontag(plint) ;
00234 
00235    if( tag != NULL && strcmp(tag, "Cardiac") == 0 ) {
00236 
00237        card = PLUTO_get_timeseries(plint) ;
00238        if( card == NULL )
00239            return "****************************\n"
00240                   "PRIC_main: bad cardiac input\n"
00241                   "****************************"   ;
00242 
00243        /* It's okay if this is zero length--it means no output */
00244        cphase1d = PLUTO_get_string(plint) ;
00245        if( cphase1d == NULL ||
00246            (strlen(cphase1d) > 0 && ! THD_filename_ok(cphase1d)) )
00247            return "*****************************\n"
00248                   "PRIC_main: bad CPhase 1D name\n"
00249                   "*****************************"   ;
00250 
00251        threshold = PLUTO_get_number(plint) ;
00252        if( threshold == BAD_NUMBER )
00253            return "******************************\n"
00254                   "PRIC_main: bad threshold input\n"
00255                   "******************************"   ;
00256 
00257        tag = PLUTO_get_optiontag(plint) ;
00258    }
00259 
00260    /*-- third line of input: Resp --*/
00261 
00262    /* By this point we already got the next option tag */
00263 
00264    if( tag != NULL && strcmp(tag, "Resp") == 0 ) {
00265 
00266        resp = PLUTO_get_timeseries(plint) ;
00267        if( resp == NULL )
00268            return "**************************\n"
00269                   "PRIC_main: bad resp input\n"
00270                   "*************************"   ;
00271 
00272        /* It's okay if this is zero length--it means no output */
00273        rphase1d = PLUTO_get_string(plint) ;
00274        if( rphase1d == NULL ||
00275            (strlen(rphase1d) > 0 && ! THD_filename_ok(rphase1d)) )
00276            return "*****************************\n"
00277                   "PRIC_main: bad RPhase 1D name\n"
00278                   "*****************************"   ;
00279 
00280        /*-- removed winsize ui
00281        winsize_input = PLUTO_get_number(plint) ;
00282        winsize = rint(winsize_input);
00283        if( winsize_input == BAD_NUMBER || winsize < 2 )
00284            return "********************************\n"
00285                   "PRIC_main: bad window size input\n"
00286                   "********************************"   ;
00287        removed winsize ui --*/
00288 
00289        tag = PLUTO_get_optiontag(plint) ;
00290    }
00291 
00292    /* Check that at least one of Cardiac and Resp were selected */
00293    if (card == NULL && resp == NULL)
00294        return "****************************************************\n"
00295               "PRIC_main: at least one of Cardiac and Resp required\n"
00296               "****************************************************"   ;
00297 
00298    /*-- fourth line of input: Params --*/
00299 
00300    /* By this point we already got the next option tag */
00301 
00302    if( tag != NULL && strcmp(tag, "Params") == 0 ) {
00303        M_input = PLUTO_get_number(plint) ;
00304        M = rint(M_input);
00305        if( M_input == BAD_NUMBER || M < 1)
00306            return "**************************\n"
00307                   "PRIC_main: bad Order input\n"
00308                   "**************************"   ;
00309    }
00310 
00311    /*------------------------------------------------------*/
00312    /*---------- At this point, the inputs are OK ----------*/
00313 
00314    PLUTO_popup_meter(plint);
00315 
00316    /*-- copy the image data for editing in place --*/
00317 
00318    new_dset = PLUTO_copy_dset( dset , new_prefix );
00319    if( new_dset == NULL )
00320        return "********************************\n"
00321               "PRIC_main: error copying dataset\n"
00322               "********************************"   ;
00323    tross_Copy_History(dset, new_dset); /* Copy and add to new_dset history */
00324    histstring = PLUTO_commandstring(plint);
00325    tross_Append_History(new_dset, histstring);
00326    free(histstring);
00327    DSET_unload( dset ) ;  /* We won't need the old dataset anymore */
00328 
00329    PLUTO_set_meter(plint, 10);
00330 
00331    /*-- calculate cardiac correction coefficients if requested --*/
00332 
00333    if (card != NULL) {
00334        /*-- convert cardiac waveform to phase --*/
00335        cardphase = RIC_ToCardiacPhase(card, threshold) ;
00336        if (cardphase == NULL) {
00337            THD_delete_3dim_dataset( new_dset , False ) ;
00338            return "******************************************\n"
00339                   "PRIC_main: error transforming cardiac data\n"
00340                   "******************************************"   ;
00341        }
00342        PLUTO_set_meter(plint, 20);
00343 
00344        /*-- calculate dataset voxel means --*/
00345        avg = RIC_CalcVoxelMeans(new_dset, ignore);
00346        if (avg == NULL) {
00347            THD_delete_3dim_dataset( new_dset , False ) ;
00348            mri_free(cardphase);
00349            return "************************************************\n"
00350                   "PRIC_main: error calculating dataset voxel means\n"
00351                   "************************************************"   ;
00352        }
00353        PLUTO_set_meter(plint, 33);
00354 
00355        /*-- calculate coefficients for each voxel --*/
00356        if (RIC_CalcCoeffAB(new_dset, cardphase, avg, &ca, &cb, M, ignore)
00357            != 0) {
00358 
00359            THD_delete_3dim_dataset( new_dset , False ) ;
00360            mri_free(cardphase);
00361            return "******************************************************\n"
00362                   "PRIC_main: error calculating cardiac a, b coefficients\n"
00363                   "******************************************************"   ;
00364        }
00365    }
00366    PLUTO_set_meter(plint, 30);
00367 
00368    /*-- calculate respiratory correction coefficients if requested --*/
00369 
00370    if (resp != NULL) {
00371        /*-- Set winsize to 1/2 sampling rate of resp in Hz --*/
00372        tr = new_dset->taxis->ttdel;
00373        switch (new_dset->taxis->units_type) {
00374        case UNITS_MSEC_TYPE: tr /= 1000; break;
00375        case UNITS_SEC_TYPE:  break;
00376        case UNITS_HZ_TYPE:   tr = 1 / tr; break;
00377        default:
00378            THD_delete_3dim_dataset( new_dset , False ) ;
00379            return "*****************************************\n"
00380                   "PRIC_main: bad time units type in dataset\n"
00381                   "*****************************************"   ;
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            return "***************************************\n"
00390                   "PRIC_main: error transforming resp data\n"
00391                   "***************************************"   ;
00392        }
00393        PLUTO_set_meter(plint, 40);
00394 
00395        /*-- calculate dataset voxel means if not already done --*/
00396        if (avg == NULL) {
00397            avg = RIC_CalcVoxelMeans(new_dset, ignore);
00398            if (avg == NULL) {
00399                THD_delete_3dim_dataset( new_dset , False ) ;
00400                mri_free(respphase);
00401                return "**************************************************\n"
00402                       "PRIC_main: error calculating dataset voxel means 2\n"
00403                       "**************************************************"   ;
00404            }
00405        }
00406        PLUTO_set_meter(plint, 43);
00407 
00408        /*-- calculate coefficients for each voxel --*/
00409        if (RIC_CalcCoeffAB(new_dset, respphase, avg, &ra, &rb, M, ignore)
00410            != 0) {
00411 
00412            THD_delete_3dim_dataset( new_dset , False ) ;
00413            mri_free(respphase);
00414            return "***************************************************\n"
00415                   "PRIC_main: error calculating resp a, b coefficients\n"
00416                   "***************************************************"   ;
00417        }
00418    }
00419    PLUTO_set_meter(plint, 50);
00420 
00421    /*-- do cardiac correction if requested --*/
00422 
00423    if (card != NULL) {
00424        /*-- correct the image data --*/
00425        if (RIC_CorrectDataset(new_dset, cardphase, ca, cb, M, ignore) != 0) {
00426            THD_delete_3dim_dataset( new_dset , False ) ;
00427            mri_free(cardphase);
00428            free(ca); free(cb);
00429            return "********************************************\n"
00430                   "PRIC_main: error applying cardiac correction\n"
00431                   "********************************************"   ;
00432        }
00433        PLUTO_set_meter(plint, 60);
00434 
00435        /*-- if requested, write phase data to file and pass to AFNI --*/
00436        if ( THD_filename_ok(cphase1d) ) {
00437            /* Write the file */
00438            fp = fopen(cphase1d, "w");
00439            nvals = cardphase->nx;
00440            cpdata = MRI_FLOAT_PTR(cardphase);
00441            for (ival = 0; ival < nvals; ival += 1) {
00442                fprintf(fp, "%f\n", cpdata[ival]);
00443            }
00444            fclose(fp);
00445 
00446            /* Pass to afni */
00447            PLUTO_register_timeseries( cphase1d , cardphase ) ;
00448        }
00449 
00450        mri_free(cardphase);
00451        free(ca); free(cb); free(avg); avg = NULL;
00452    }
00453    PLUTO_set_meter(plint, 70);
00454 
00455    /*-- do resp correction if requested --*/
00456 
00457    if (resp != NULL) {
00458        /*-- correct the image data --*/
00459        if (RIC_CorrectDataset(new_dset, respphase, ra, rb, M, ignore) != 0) {
00460            THD_delete_3dim_dataset( new_dset , False ) ;
00461            mri_free(respphase);
00462            free(ra); free(rb);
00463            return "*****************************************\n"
00464                   "PRIC_main: error applying resp correction\n"
00465                   "*****************************************"   ;
00466        }
00467        PLUTO_set_meter(plint, 80);
00468 
00469        /*-- if requested, write phase data to file and pass to AFNI --*/
00470        if ( THD_filename_ok(rphase1d) ) {
00471            /* Write the file */
00472            fp = fopen(rphase1d, "w");
00473            nvals = respphase->nx;
00474            rpdata = MRI_FLOAT_PTR(respphase);
00475            for (ival = 0; ival < nvals; ival += 1) {
00476                fprintf(fp, "%f\n", rpdata[ival]);
00477            }
00478            fclose(fp);
00479 
00480            /* Pass to afni */
00481            PLUTO_register_timeseries( rphase1d , respphase ) ;
00482        }
00483 
00484        mri_free(respphase);
00485        free(ra); free(rb); if (avg != NULL) free(avg);
00486    }
00487    PLUTO_set_meter(plint, 90);
00488 
00489    /*-- give new dataset to AFNI --*/
00490 
00491    if ( PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ) {
00492        THD_delete_3dim_dataset( new_dset , False ) ;
00493        return "*********************************************\n"
00494               "PRIC_main: failure to add new dataset to AFNI\n"
00495               "*********************************************" ;
00496    }
00497    PLUTO_set_meter(plint, 100);
00498 
00499    /*-- done successfully!!! --*/
00500 
00501    return NULL ;
00502 }