00001 /***********************************************************************
00002   retroicor.c: Implementation of 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 #include <math.h>
00021 
00022 #ifndef M_PI
00023 # define M_PI           3.14159265358979323846  /* pi */
00024 #endif
00025 
00026 #include "retroicor.h"
00027 
00028 /* Find the next local maximum value in the float data sequence that is
00029    above the given threshold, starting at the given point.
00030 
00031    PRE:
00032    cdata is an array of float data;
00033    numpoints is the total number of datapoints in the cdata array;
00034    startpoint is the index of the first point to consider (zero indexed);
00035    maxpoint and endpoint are valid non-NULL pointers to integer values;
00036    threshold is the value which data points must exceed to be considered;
00037    0 <= startpoint < numpoints;
00038 
00039    POST:
00040    return value is -1 on error or if there was no data with value > threshold,
00041    else return value is 0 and *maxpoint is the index of the point with maximum
00042    value in the next peak after and including startpoint and *endpoint is the
00043    index of the first point after the peak, where a peak is a contiguous
00044    series of datapoints with value > threshold;
00045 */
00046 int _RIC_findNextCardiacPeak(const float * cdata,
00047                              int numpoints, int startpoint,
00048                              int * maxpoint, int * endpoint,
00049                              float threshold) {
00050 
00051     int curidx = startpoint;  /* Index of current point under examination */
00052     int maxidx;               /* Index of maximum point found so far */
00053 
00054     /* Quick check of arguments */
00055     if (cdata == NULL || startpoint >= numpoints || startpoint < 0 ||
00056         maxpoint == NULL || endpoint == NULL) {
00057 
00058         return -1;
00059     }
00060 
00061     /* Skip over data <= threshold, if any */
00062     while (cdata[curidx] <= threshold) {
00063         curidx += 1;
00064 
00065         /* Check if we ran out of data without finding anything > threshold */
00066         if (curidx >= numpoints) {
00067             return -1;
00068         }
00069     }
00070 
00071     /* Find the max point in this region > threshold */
00072     maxidx = curidx;
00073     while (cdata[curidx] > threshold) {
00074         /* Check if this is a new max */
00075         if (cdata[curidx] > cdata[maxidx]) {
00076             maxidx = curidx;
00077         }
00078 
00079         curidx += 1;
00080 
00081         /* Stop if we ran out of data */
00082         if (curidx >= numpoints) {
00083             break;
00084         }
00085     }
00086 
00087     *maxpoint = maxidx;
00088     *endpoint = curidx;
00089     return 0;
00090 }
00091 
00092 /* Implementation Note: Some of the calculations in the loops can be
00093    factored out, but I choose not to for clarity. This is all O(N) in the
00094    number of timepoints anyways (I think).
00095 */
00096 MRI_IMAGE * RIC_ToCardiacPhase(const MRI_IMAGE * card, float threshold) {
00097 
00098     int numSamps;            /* Number of samples in vector */
00099     MRI_IMAGE * cardphase;   /* The cardiac phase vector to return */
00100     float * cpdata;          /* Pointer to cardiac phase vector data */
00101     float * cdata;           /* Pointer to cardiac vector data */
00102     double twoPI = 2 * M_PI; /* 2 * PI  (intermediate value for calculation) */
00103     double twoPI_t2_t1;      /* 2 * PI / (t_2 - t_1)  (intermediate value) */
00104     double phase;            /* card phase: 2 * PI / (t_2 - t_1) * (t - t_1) */
00105     int lastpeakpt = 0;      /* The last cdata timepoint searched for a peak */
00106     int t = 0;               /* The current time */
00107     int t_1 = 0;             /* The previous cardiac peak time */
00108     int t_2;                 /* The next cardiac peak time */
00109 
00110     /* Quick check of arguments */
00111     if (card == NULL || card->nx < 2 || card->kind != MRI_float) {
00112         return NULL;
00113     }
00114 
00115     /* Initialize */
00116     numSamps = card->nx;
00117     cardphase = mri_new(numSamps, 1, MRI_float);
00118     cpdata = MRI_FLOAT_PTR(cardphase);
00119     cdata = MRI_FLOAT_PTR(card);
00120 
00121     /* Iterate over the cardiac peaks, assuming data start is peak */
00122     while (_RIC_findNextCardiacPeak(cdata, numSamps, lastpeakpt, &t_2,
00123                                     &lastpeakpt, threshold) == 0) {
00124 
00125         /* Fill in the cardiac phase values between peaks */
00126         twoPI_t2_t1 = twoPI / (t_2 - t_1);
00127         phase = 0.0;  /* Since we always have t == t_1 at this point) */
00128         for ( ; t < t_2; t += 1) {
00129             cpdata[t] = phase;
00130             phase += twoPI_t2_t1;
00131         }
00132 
00133         t_1 = t_2;
00134     }
00135 
00136     /* Fill in any remaining phase values, assuming end of data is peak */
00137     twoPI_t2_t1 = twoPI / (numSamps - t_1);
00138     phase = 0.0;
00139     for ( ; t < numSamps; t += 1) {
00140         cpdata[t] = phase;
00141         phase += twoPI_t2_t1;
00142     }
00143 
00144     return cardphase;
00145 }
00146 
00147 MRI_IMAGE * RIC_ToRespPhase(const MRI_IMAGE * resp, int winsize) {
00148 
00149     int numSamps;           /* Number of samples in input vector */
00150     MRI_IMAGE * respphase;    /* The resp phase vector to return */
00151     float * rpdata;           /* Pointer to resp phase vector data */
00152     float * rdata;            /* Pointer to resp vector data */
00153     float * nrdata;           /* Pointer to normalized resp data */
00154     int curr;                 /* Current resp data vector item */
00155     int currdist;             /* Distance from current resp data vector item */
00156     float maxr, minr;         /* Max and min resp data vector values */
00157     double hist[RIC_HISTSIZE];/* Histogram of resp data values */
00158     int curb;                 /* Current histogram bin */
00159     double binfact;           /* Multiply value by this to get bin number */
00160     double binshift;          /* Subtract from bin # to centre bins on range */
00161     double leftsum, rightsum; /* Sum of resp values left, right of point */
00162     double nisPI;             /* Intermediate value: M_PI / numInSamps */
00163 
00164     /* Quick check of arguments */
00165     if (resp == NULL || resp->nx < 2 || resp->kind != MRI_float || winsize<2) {
00166         return NULL;
00167     }
00168 
00169     /* Initialize */
00170     numSamps = resp->nx;
00171     nrdata = malloc(sizeof(float) * numSamps);
00172     if (nrdata == NULL) {
00173         return NULL;
00174     }
00175     respphase = mri_new(numSamps, 1, MRI_float);
00176     rpdata = MRI_FLOAT_PTR(respphase);
00177     rdata = MRI_FLOAT_PTR(resp);
00178     for (curb = 0; curb < RIC_HISTSIZE; curb += 1) {
00179         hist[curb] = 0.0;
00180     }
00181 
00182     /* Collect stats */
00183 
00184     /* Find max and min (for normalization and histogram range) */
00185     maxr = minr = rdata[0];
00186     for (curr = 1; curr < numSamps; curr += 1) {
00187         if (rdata[curr] > maxr) {
00188             maxr = rdata[curr];
00189         } else if (rdata[curr] < minr) {
00190             minr = rdata[curr];
00191         }
00192     }
00193 
00194     /* Normalize the data (just subtract minr to make everything >= 0) */
00195     for (curr = 0; curr < numSamps; curr += 1) {
00196         nrdata[curr] = rdata[curr] - minr;
00197     }
00198     maxr -= minr;
00199     minr = 0.0;
00200 
00201     /***** We'll use nrdata from now on instead of rdata *****/
00202 
00203     /* Create histogram (the following makes sense only if all data >= 0) */
00204     /* Actually this is broken; rounding means the probability of a value
00205        falling into the first and last buckets is less than it should be
00206        --but this is good enough and takes few brain cells to program */
00207     /* Added a fudge factor so that the first and last histogram bins are
00208        almost as wide as the others but not so wide that rounding the bin
00209        number could cause it to refer to a nonexistant bin */
00210     binfact = (RIC_HISTSIZE - 2 * RIC_HISTFUDGE) / maxr;
00211     binshift = 0.5 - RIC_HISTFUDGE;
00212     if (binfact <= 0.0 || binshift <= 0.0) {  /* Check just in case! */
00213         free(nrdata);
00214         return NULL;
00215     }
00216     for (curr = 0; curr < numSamps; curr += 1) {
00217         hist[(int)rint(nrdata[curr] * binfact - binshift)] += 1.0;
00218         /* a.k.a. nrdata[curr] / maxr * RIC_HISTSIZE */
00219     }
00220 
00221     /* Turn the histogram into a vector of percentile values */
00222     for (curb = 1; curb < RIC_HISTSIZE; curb += 1) {
00223         hist[curb] += hist[curb - 1];
00224     }
00225     /* Actually, make that PI * percentile, since we use that later */
00226     nisPI = M_PI / numSamps;
00227     for (curb = 0; curb < RIC_HISTSIZE; curb += 1) {
00228         hist[curb] *= nisPI;
00229     }
00230 
00231     /* Iterate over the resp waveform samples, converting each to phase */
00232 
00233     for (curr = 0; curr < numSamps; curr += 1) {
00234         /* Estimate of dR/dt: Slope of N-point means on either side of point */
00235         /* Further Simplification: We don't actually need the mean to find
00236            the slope, just the sums on either side of the point */
00237         /* There are better ways to do this but YES WE HAVE NO BANANAS */
00238         rightsum = leftsum = 0;
00239         for (currdist = 0; currdist < winsize; currdist += 1) {
00240             if (curr + currdist < numSamps) {
00241                 rightsum += nrdata[curr + currdist];
00242             } else {
00243                 rightsum += nrdata[curr];
00244             }
00245             if (curr - currdist >= 0) {
00246                 leftsum += nrdata[curr - currdist];
00247             } else {
00248                 leftsum += nrdata[curr];
00249             }
00250         }
00251 
00252         /* Hooray!! Now we get the result! */
00253         if ((rightsum - leftsum) < 0.0) {   /* Negative slope */
00254             rpdata[curr] = - hist[(int)rint(nrdata[curr] * binfact
00255                                             - binshift)];
00256         } else {                            /* Positive slope */
00257             rpdata[curr] =   hist[(int)rint(nrdata[curr] * binfact
00258                                             - binshift)];
00259         }
00260     }
00261 
00262     free(nrdata);
00263     return respphase;
00264 }
00265 
00266 /* I know the "== 0.0" is unsafe, but the AFNI docs say that a zero scale
00267    factor means the data is not scaled. Compare with epsilon? */
00268 #define RIC_CALCVOXELMEANS__DO_VOXSUM(t) {              \
00269     t * brickdset = DSET_ARRAY(dset, ival);             \
00270     if (brickdset == NULL) {                            \
00271         free(avg);                                      \
00272         return NULL;                                    \
00273     }                                                   \
00274     if (scalefactor == 0.0) {                           \
00275         for (ivox = 0; ivox < nvoxs; ivox += 1) {       \
00276             avg[ivox] += brickdset[ivox];               \
00277         }                                               \
00278     } else {                                            \
00279         for (ivox = 0; ivox < nvoxs; ivox += 1) {       \
00280             avg[ivox] += brickdset[ivox] * scalefactor; \
00281         }                                               \
00282     }                                                   \
00283 }
00284 
00285 /* Calculates the average value for each voxel in dset across all timepoints.
00286 
00287    PRE:
00288    dset is a valid 3D+T dataset with at least 1 timepoint;
00289    ignore is the number of timepoints to ignore at the beginning of dset;
00290    0 <= ignore < number of timepoints in dset;
00291 
00292    POST:
00293    return value is NULL on error,
00294    else, return value is an array of floats with the same dimensions as the
00295    subbricks of dset, containing the voxel value averages;
00296 
00297    Note: The caller is responsible for free()ing the returned block of memory
00298    when done.
00299 
00300    Note2: The complex datatype is not supported, and any such bricks will
00301    result in an error (NULL return value).
00302 */
00303 double * RIC_CalcVoxelMeans(const THD_3dim_dataset * dset, int ignore) {
00304 
00305     double * avg;       /* The voxel averages to be returned */
00306     float scalefactor; /* Current dset brick scaling factor */
00307     int ival, nvals;   /* Current, number of dset timepoints */
00308     int ivox, nvoxs;   /* Current, number of dset brick voxels */
00309 
00310     /* Quick check of arguments */
00311     if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||
00312         ignore < 0 || ignore >= DSET_NVALS(dset)) {
00313 
00314         return NULL;
00315     }
00316 
00317     /* Initialize */
00318     DSET_load(dset);
00319     nvals = DSET_NVALS(dset);
00320     nvoxs = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;
00321     avg = malloc(sizeof(double) * nvoxs);
00322     if (avg == NULL) {
00323         return NULL;
00324     }
00325 
00326     /* Calculate the voxel averages; treat matrices as 1-D arrays */
00327 
00328     /* Zero the voxel sums */
00329     for (ivox = 0; ivox < nvoxs; ivox += 1) {
00330         avg[ivox] = 0.0;
00331     }
00332 
00333     /* Sum each voxel across time (and hope there are not too many points) */
00334     for (ival = ignore; ival < nvals; ival += 1) {
00335         scalefactor = DSET_BRICK_FACTOR(dset, ival);
00336 
00337         switch (DSET_BRICK_TYPE(dset, ival)) {
00338         case MRI_short:
00339             RIC_CALCVOXELMEANS__DO_VOXSUM(short);
00340             break;
00341         case MRI_byte:
00342             RIC_CALCVOXELMEANS__DO_VOXSUM(byte);
00343             break;
00344         case MRI_float:
00345             RIC_CALCVOXELMEANS__DO_VOXSUM(float);
00346             break;
00347         default: /* Unsupported datatype */
00348             free(avg);
00349             return NULL;
00350         }
00351     }
00352 
00353     /* Divide by number of timepoints to get average */
00354     nvals -= ignore;  /* We do not average over the ignored timepoints */
00355     for (ivox = 0; ivox < nvoxs; ivox += 1) {
00356         avg[ivox] /= nvals;
00357     }
00358 
00359     return avg;
00360 }
00361 
00362 /* I know the "== 0.0" is unsafe, but the AFNI docs say that a zero scale
00363    factor means the data is not scaled. Compare with epsilon? */
00364 #define RIC_CALCCOEFFAB__DO_CALCNUMERDENOM(t) {                               \
00365     /* Load and check the dataset brick */                                    \
00366     t * brickdset = DSET_ARRAY(dset, ival);                                   \
00367     if (brickdset == NULL) {                                                  \
00368         free(newa); free(newb); free(denoma); free(denomb);                   \
00369         return -1;                                                            \
00370     }                                                                         \
00371                                                                               \
00372     idenom = 0;                                                               \
00373     inumer = 0;                                                               \
00374     /* m loop outside voxel loop <- coeff arrays bigger than one volume */    \
00375     for (m = 1; m <= M; m += 1) {                                             \
00376         slicestart = 0;                                                       \
00377         /* Iterate over slices within this TR */                              \
00378         for (islice = 0; islice < nslices; islice += 1) {                     \
00379             sliceend = slicestart + nvoxpers;                                 \
00380                                                                               \
00381             /* Get the slice time in milliseconds */                          \
00382             slicetime = THD_timeof_slice(ival, islice, dset);                 \
00383             switch (dset->taxis->units_type) {                                \
00384             case UNITS_MSEC_TYPE: break;                                      \
00385             case UNITS_SEC_TYPE:  slicetime *= 1000; break;                   \
00386             case UNITS_HZ_TYPE:   slicetime = 1000 / slicetime; break;        \
00387             default: free(newa); free(newb); free(denoma); free(denomb);      \
00388                 return -1;                                                    \
00389             }                                                                 \
00390                                                                               \
00391             /* Find phase at slice time and add up the denominators */        \
00392             mp = m * pdata[(int)(slicetime / sampd)];                         \
00393             cmp = cos(mp); smp = sin(mp);                                     \
00394             denoma[idenom] += cmp * cmp; denomb[idenom] += smp * smp;         \
00395                                                                               \
00396             /* Add up the numerators, using scale factor if necessary */      \
00397             if (scalefactor == 0.0) {                                         \
00398                 for (ivox = slicestart; ivox < sliceend; ivox += 1) {         \
00399                     newa[inumer] += (brickdset[ivox] - avg[ivox]) * cmp;      \
00400                     newb[inumer] += (brickdset[ivox] - avg[ivox]) * smp;      \
00401                     inumer += 1;                                              \
00402                 }                                                             \
00403             } else {                                                          \
00404                 for (ivox = slicestart; ivox < sliceend; ivox += 1) {         \
00405                     newa[inumer] +=                                           \
00406                         ((brickdset[ivox]*scalefactor) - avg[ivox]) * cmp;    \
00407                     newb[inumer] +=                                           \
00408                         ((brickdset[ivox]*scalefactor) - avg[ivox]) * smp;    \
00409                     inumer += 1;                                              \
00410                 }                                                             \
00411             }                                                                 \
00412                                                                               \
00413             slicestart = sliceend;                                            \
00414             idenom += 1;                                                      \
00415         }                                                                     \
00416     }                                                                         \
00417 }
00418 
00419 int RIC_CalcCoeffAB(THD_3dim_dataset * dset, const MRI_IMAGE * phase,
00420                     const double * avg, double ** a, double ** b,
00421                     int M, int ignore) {
00422 
00423     float scalefactor;         /* Scaling factor for current dset brick */
00424     int ival, nvals;           /* Current, number of dset timepoints */
00425     int ivox, nvoxs;           /* Current, number of voxels in dset subbrick */
00426     double * newa, * newb;     /* Coefficients a and b, to be calculated */
00427     float * pdata;             /* Phase data */
00428     int m;                     /* Current order number (1, M) */
00429     double mp, cmp, smp;       /* m * current phase; cos(mp); sin(mp) */
00430     double sampd;              /* Intersample time in ms for phase */
00431     double slicetime;          /* Time of current slice in milliseconds */
00432     int slicestart;            /* Voxel index of 1st voxel in current slice */
00433     int sliceend;              /* End voxel index of current slice */
00434     double * denoma, * denomb; /* Coeff denominators for each m and slice */
00435     int idenom, ndenoms;       /* Current, number of coeff denominators */
00436     int inumer, nnumers;       /* Current, number of coeff numerators */
00437     int islice, nslices;       /* Current, number of slices in a subbrick */
00438     int nvoxpers;              /* The number of voxels per slice */
00439 
00440     /* Quick check of arguments */
00441     if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||
00442         !ISVALID_TIMEAXIS(dset->taxis) ||
00443         phase == NULL || phase->nx < 1 || phase->ny != 1 ||
00444         avg == NULL || a == NULL || b == NULL || M < 1 ||
00445         ignore < 0 || ignore >= DSET_NVALS(dset)) {
00446 
00447         return -1;
00448     }
00449 
00450     /* Initialize */
00451     DSET_load(dset);
00452     nvals = DSET_NVALS(dset);
00453     nvoxpers = dset->daxes->nxx * dset->daxes->nyy;
00454     nslices = dset->daxes->nzz;
00455     nvoxs = nvoxpers * nslices;
00456     ndenoms = nslices * M;
00457     nnumers = nvoxs * M;
00458     sampd = dset->taxis->ttdel * nvals / phase->nx;
00459     switch (dset->taxis->units_type) {
00460     case UNITS_MSEC_TYPE: break;
00461     case UNITS_SEC_TYPE:  sampd *= 1000; break;
00462     case UNITS_HZ_TYPE:   sampd = 1000 / sampd; break;
00463     default: return -1;
00464     }
00465     pdata = MRI_FLOAT_PTR(phase);
00466     newa = malloc(sizeof(double) * nnumers);
00467     if (newa == NULL) {
00468         return -1;
00469     }
00470     newb = malloc(sizeof(double) * nnumers);
00471     if (newb == NULL) {
00472         free(newa);
00473         return -1;
00474     }
00475     for (inumer = 0; inumer < nnumers; inumer += 1) {
00476         newa[inumer] = 0.0; newb[inumer] = 0.0;
00477     }
00478     denoma = malloc(sizeof(double) * ndenoms);
00479     if (denoma == NULL) {
00480         free(newa); free(newb);
00481         return -1;
00482     }
00483     denomb = malloc(sizeof(double) * ndenoms);
00484     if (denomb == NULL) {
00485         free(newa); free(newb); free(denoma);
00486         return -1;
00487     }
00488     for (idenom = 0; idenom < ndenoms; idenom += 1) {
00489         denoma[idenom] = 0.0; denomb[idenom] = 0.0;
00490     }
00491 
00492     /* Calculate the coeff numerators and denominators */
00493     /* Iterate over dset timepoints (TRs) => BRIK read from disk once */
00494     for (ival = ignore; ival < nvals; ival += 1) {
00495         scalefactor = DSET_BRICK_FACTOR(dset, ival);
00496 
00497         /* Calculate coeff numerators and denominators over all data */
00498         switch (DSET_BRICK_TYPE(dset, ival)) {
00499         case MRI_short:
00500             RIC_CALCCOEFFAB__DO_CALCNUMERDENOM(short);
00501             break;
00502         case MRI_byte:
00503             RIC_CALCCOEFFAB__DO_CALCNUMERDENOM(byte);
00504             break;
00505         case MRI_float:
00506             RIC_CALCCOEFFAB__DO_CALCNUMERDENOM(float);
00507             break;
00508         default: /* Unsupported datatype */
00509             free(newa); free(newb); free(denoma); free(denomb);
00510             return -1;
00511         }
00512     }
00513 
00514     /* Divide the coeff numerators by their denominators */
00515     idenom = 0;
00516     inumer = 0;
00517     /* m loop outside voxel loop <- coeff arrays bigger than one volume */
00518     for (m = 1; m <= M; m += 1) {
00519         slicestart = 0;
00520         /* Iterate over slices */
00521         for (islice = 0; islice < nslices; islice += 1) {
00522             sliceend = slicestart + nvoxpers;
00523             /* Divide numerator by denominator */
00524             for (ivox = slicestart; ivox < sliceend; ivox += 1) {
00525                 newa[inumer] /= denoma[idenom];
00526                 newb[inumer] /= denomb[idenom];
00527                 inumer += 1;
00528             }
00529             slicestart = sliceend;
00530             idenom += 1;
00531         }
00532     }
00533 
00534     *a = newa; *b = newb;
00535     free(denoma); free(denomb);
00536     return 0;
00537 }
00538 
00539 /* I know the "== 0.0" is unsafe, but the AFNI docs say that a zero scale
00540    factor means the data is not scaled. Compare with epsilon? */
00541 #define RIC_CORRECTDATASET__DO_CORRECT(t) {                                   \
00542     /* Load and check the dataset brick */                                    \
00543     t * brickdset = DSET_BRICK_ARRAY(dset, ival);                             \
00544     if (brickdset == NULL) {                                                  \
00545         return -1;                                                            \
00546     }                                                                         \
00547                                                                               \
00548     icoeff = 0;                                                               \
00549     /* m loop outside voxel loop <- coeff arrays bigger than one volume */    \
00550     for (m = 1; m <= M; m += 1) {                                             \
00551         slicestart = 0;                                                       \
00552         /* Iterate over slices within this TR */                              \
00553         for (islice = 0; islice < nslices; islice += 1) {                     \
00554             sliceend = slicestart + nvoxpers;                                 \
00555                                                                               \
00556             /* Get the slice time in milliseconds */                          \
00557             slicetime = THD_timeof_slice(ival, islice, dset);                 \
00558             switch (dset->taxis->units_type) {                                \
00559             case UNITS_MSEC_TYPE: break;                                      \
00560             case UNITS_SEC_TYPE:  slicetime *= 1000; break;                   \
00561             case UNITS_HZ_TYPE:   slicetime = 1000 / slicetime; break;        \
00562             default: return -1;                                               \
00563             }                                                                 \
00564                                                                               \
00565             /* Find phase at slice time and calculate intermediate value */   \
00566             mp = m * pdata[(int)(slicetime / sampd)];                         \
00567             cmp = cos(mp); smp = sin(mp);                                     \
00568                                                                               \
00569             /* Apply the corrections to the voxels */                         \
00570             if (scalefactor == 0.0) {                                         \
00571                 for (ivox = slicestart; ivox < sliceend; ivox += 1) {         \
00572                     brickdset[ivox] -= a[icoeff] * cmp + b[icoeff] * smp;     \
00573                     icoeff += 1;                                              \
00574                 }                                                             \
00575             } else {                                                          \
00576                 for (ivox = slicestart; ivox < sliceend; ivox += 1) {         \
00577                     brickdset[ivox] -= (a[icoeff] * cmp + b[icoeff] * smp)    \
00578                         / scalefactor;                                        \
00579                     icoeff += 1;                                              \
00580                 }                                                             \
00581             }                                                                 \
00582                                                                               \
00583             slicestart = sliceend;                                            \
00584         }                                                                     \
00585     }                                                                         \
00586 }
00587 
00588 int RIC_CorrectDataset(THD_3dim_dataset * dset, const MRI_IMAGE * phase,
00589                        const double * a, const double * b,
00590                        int M, int ignore) {
00591 
00592     float scalefactor;         /* Scaling factor for current dset brick */
00593     int ival, nvals;           /* Current, number of dset timepoints */
00594     int ivox, nvoxs;           /* Current, number of voxels in dset subbrick */
00595     float * pdata;             /* Phase data */
00596     int m;                     /* Current order number (1, M) */
00597     double mp, cmp, smp;       /* m * phase; cos(mp); sin(mp) */
00598     double sampd;              /* Intersample time in ms for phase */
00599     double slicetime;          /* Time of current slice in milliseconds */
00600     int slicestart;            /* Voxel index of 1st voxel in current slice */
00601     int sliceend;              /* End voxel index of current slice */
00602     int icoeff, ncoeffs;       /* Current, number of a,b coefficients */
00603     int islice, nslices;       /* Current, number of slices in a subbrick */
00604     int nvoxpers;              /* The number of voxels per slice */
00605 
00606     /* Quick check of arguments */
00607     if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||
00608         !ISVALID_TIMEAXIS(dset->taxis) ||
00609         phase == NULL || phase->nx < 1 || phase->ny != 1 ||
00610         a == NULL || b == NULL || M < 1 ||
00611         ignore < 0 || ignore >= DSET_NVALS(dset)) {
00612 
00613         return -1;
00614     }
00615 
00616     /* Initialize */
00617     DSET_load(dset);
00618     nvals = DSET_NVALS(dset);
00619     nvoxpers = dset->daxes->nxx * dset->daxes->nyy;
00620     nslices = dset->daxes->nzz;
00621     nvoxs = nvoxpers * nslices;
00622     ncoeffs = nvoxs * M;
00623     sampd = dset->taxis->ttdel * nvals / phase->nx;
00624     switch (dset->taxis->units_type) {
00625     case UNITS_MSEC_TYPE: break;
00626     case UNITS_SEC_TYPE:  sampd *= 1000; break;
00627     case UNITS_HZ_TYPE:   sampd = 1000 / sampd; break;
00628     default: return -1;
00629     }
00630     pdata = MRI_FLOAT_PTR(phase);
00631 
00632     /* Correct each subbrick in dset */
00633     /* Iterate over dset timepoints (TRs) => BRIK read from disk once */
00634     for (ival = ignore; ival < nvals; ival += 1) {
00635         scalefactor = DSET_BRICK_FACTOR(dset, ival);
00636 
00637         /* Apply each order of correction to each voxel of the subbrick */
00638         switch (DSET_BRICK_TYPE(dset, ival)) {
00639         case MRI_short:
00640             RIC_CORRECTDATASET__DO_CORRECT(short);
00641             break;
00642         case MRI_byte:
00643             RIC_CORRECTDATASET__DO_CORRECT(byte);
00644             break;
00645         case MRI_float:
00646             RIC_CORRECTDATASET__DO_CORRECT(float);
00647             break;
00648         default: /* Unsupported datatype */
00649             return -1;
00650         }
00651     }
00652 
00653     return 0;
00654 }