{              \
    t * brickdset = DSET_ARRAY(dset, ival);             \
    if (brickdset == NULL) {                            \
        free(avg);                                      \
        return NULL;                                    \
    }                                                   \
    if (scalefactor == 0.0) {                           \
        for (ivox = 0; ivox < nvoxs; ivox += 1) {       \
            avg[ivox] += brickdset[ivox];               \
        }                                               \
    } else {                                            \
        for (ivox = 0; ivox < nvoxs; ivox += 1) {       \
            avg[ivox] += brickdset[ivox] * scalefactor; \
        }                                               \
    }                                                   \
}

00049                                               {
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 }

00421                                        {
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 }

00303                                                                        {
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 }

00590                                           {
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 }

00096                                                                         {
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 }

00147                                                                  {
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 }