Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
retroicor.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <math.h>
00021
00022 #ifndef M_PI
00023 # define M_PI 3.14159265358979323846
00024 #endif
00025
00026 #include "retroicor.h"
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
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;
00052 int maxidx;
00053
00054
00055 if (cdata == NULL || startpoint >= numpoints || startpoint < 0 ||
00056 maxpoint == NULL || endpoint == NULL) {
00057
00058 return -1;
00059 }
00060
00061
00062 while (cdata[curidx] <= threshold) {
00063 curidx += 1;
00064
00065
00066 if (curidx >= numpoints) {
00067 return -1;
00068 }
00069 }
00070
00071
00072 maxidx = curidx;
00073 while (cdata[curidx] > threshold) {
00074
00075 if (cdata[curidx] > cdata[maxidx]) {
00076 maxidx = curidx;
00077 }
00078
00079 curidx += 1;
00080
00081
00082 if (curidx >= numpoints) {
00083 break;
00084 }
00085 }
00086
00087 *maxpoint = maxidx;
00088 *endpoint = curidx;
00089 return 0;
00090 }
00091
00092
00093
00094
00095
00096 MRI_IMAGE * RIC_ToCardiacPhase(const MRI_IMAGE * card, float threshold) {
00097
00098 int numSamps;
00099 MRI_IMAGE * cardphase;
00100 float * cpdata;
00101 float * cdata;
00102 double twoPI = 2 * M_PI;
00103 double twoPI_t2_t1;
00104 double phase;
00105 int lastpeakpt = 0;
00106 int t = 0;
00107 int t_1 = 0;
00108 int t_2;
00109
00110
00111 if (card == NULL || card->nx < 2 || card->kind != MRI_float) {
00112 return NULL;
00113 }
00114
00115
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
00122 while (_RIC_findNextCardiacPeak(cdata, numSamps, lastpeakpt, &t_2,
00123 &lastpeakpt, threshold) == 0) {
00124
00125
00126 twoPI_t2_t1 = twoPI / (t_2 - t_1);
00127 phase = 0.0;
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
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;
00150 MRI_IMAGE * respphase;
00151 float * rpdata;
00152 float * rdata;
00153 float * nrdata;
00154 int curr;
00155 int currdist;
00156 float maxr, minr;
00157 double hist[RIC_HISTSIZE];
00158 int curb;
00159 double binfact;
00160 double binshift;
00161 double leftsum, rightsum;
00162 double nisPI;
00163
00164
00165 if (resp == NULL || resp->nx < 2 || resp->kind != MRI_float || winsize<2) {
00166 return NULL;
00167 }
00168
00169
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
00183
00184
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
00195 for (curr = 0; curr < numSamps; curr += 1) {
00196 nrdata[curr] = rdata[curr] - minr;
00197 }
00198 maxr -= minr;
00199 minr = 0.0;
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 binfact = (RIC_HISTSIZE - 2 * RIC_HISTFUDGE) / maxr;
00211 binshift = 0.5 - RIC_HISTFUDGE;
00212 if (binfact <= 0.0 || binshift <= 0.0) {
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
00219 }
00220
00221
00222 for (curb = 1; curb < RIC_HISTSIZE; curb += 1) {
00223 hist[curb] += hist[curb - 1];
00224 }
00225
00226 nisPI = M_PI / numSamps;
00227 for (curb = 0; curb < RIC_HISTSIZE; curb += 1) {
00228 hist[curb] *= nisPI;
00229 }
00230
00231
00232
00233 for (curr = 0; curr < numSamps; curr += 1) {
00234
00235
00236
00237
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
00253 if ((rightsum - leftsum) < 0.0) {
00254 rpdata[curr] = - hist[(int)rint(nrdata[curr] * binfact
00255 - binshift)];
00256 } else {
00257 rpdata[curr] = hist[(int)rint(nrdata[curr] * binfact
00258 - binshift)];
00259 }
00260 }
00261
00262 free(nrdata);
00263 return respphase;
00264 }
00265
00266
00267
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
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 double * RIC_CalcVoxelMeans(const THD_3dim_dataset * dset, int ignore) {
00304
00305 double * avg;
00306 float scalefactor;
00307 int ival, nvals;
00308 int ivox, nvoxs;
00309
00310
00311 if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||
00312 ignore < 0 || ignore >= DSET_NVALS(dset)) {
00313
00314 return NULL;
00315 }
00316
00317
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
00327
00328
00329 for (ivox = 0; ivox < nvoxs; ivox += 1) {
00330 avg[ivox] = 0.0;
00331 }
00332
00333
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:
00348 free(avg);
00349 return NULL;
00350 }
00351 }
00352
00353
00354 nvals -= ignore;
00355 for (ivox = 0; ivox < nvoxs; ivox += 1) {
00356 avg[ivox] /= nvals;
00357 }
00358
00359 return avg;
00360 }
00361
00362
00363
00364 #define RIC_CALCCOEFFAB__DO_CALCNUMERDENOM(t) { \
00365 \
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 \
00375 for (m = 1; m <= M; m += 1) { \
00376 slicestart = 0; \
00377 \
00378 for (islice = 0; islice < nslices; islice += 1) { \
00379 sliceend = slicestart + nvoxpers; \
00380 \
00381 \
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 \
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 \
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;
00424 int ival, nvals;
00425 int ivox, nvoxs;
00426 double * newa, * newb;
00427 float * pdata;
00428 int m;
00429 double mp, cmp, smp;
00430 double sampd;
00431 double slicetime;
00432 int slicestart;
00433 int sliceend;
00434 double * denoma, * denomb;
00435 int idenom, ndenoms;
00436 int inumer, nnumers;
00437 int islice, nslices;
00438 int nvoxpers;
00439
00440
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
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
00493
00494 for (ival = ignore; ival < nvals; ival += 1) {
00495 scalefactor = DSET_BRICK_FACTOR(dset, ival);
00496
00497
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:
00509 free(newa); free(newb); free(denoma); free(denomb);
00510 return -1;
00511 }
00512 }
00513
00514
00515 idenom = 0;
00516 inumer = 0;
00517
00518 for (m = 1; m <= M; m += 1) {
00519 slicestart = 0;
00520
00521 for (islice = 0; islice < nslices; islice += 1) {
00522 sliceend = slicestart + nvoxpers;
00523
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
00540
00541 #define RIC_CORRECTDATASET__DO_CORRECT(t) { \
00542 \
00543 t * brickdset = DSET_BRICK_ARRAY(dset, ival); \
00544 if (brickdset == NULL) { \
00545 return -1; \
00546 } \
00547 \
00548 icoeff = 0; \
00549 \
00550 for (m = 1; m <= M; m += 1) { \
00551 slicestart = 0; \
00552 \
00553 for (islice = 0; islice < nslices; islice += 1) { \
00554 sliceend = slicestart + nvoxpers; \
00555 \
00556 \
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 \
00566 mp = m * pdata[(int)(slicetime / sampd)]; \
00567 cmp = cos(mp); smp = sin(mp); \
00568 \
00569 \
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;
00593 int ival, nvals;
00594 int ivox, nvoxs;
00595 float * pdata;
00596 int m;
00597 double mp, cmp, smp;
00598 double sampd;
00599 double slicetime;
00600 int slicestart;
00601 int sliceend;
00602 int icoeff, ncoeffs;
00603 int islice, nslices;
00604 int nvoxpers;
00605
00606
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
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
00633
00634 for (ival = ignore; ival < nvals; ival += 1) {
00635 scalefactor = DSET_BRICK_FACTOR(dset, ival);
00636
00637
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:
00649 return -1;
00650 }
00651 }
00652
00653 return 0;
00654 }