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 }

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 }