TPCCLIB
Loading...
Searching...
No Matches
Functions
litac.c File Reference

Linear interpolation of TACs. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcli.h"
#include "tpctac.h"
#include "tpctacmod.h"

Go to the source code of this file.

Functions

int tacIntegrate (TAC *inp, TAC *out, TPCSTATUS *status)
 Integrate TACs from one TAC structure into a new TAC structure.
int tacIntegrateFE (TAC *inp, TAC *out, TAC *out2, TPCSTATUS *status)
 Integrate TACs from one TAC structure into a new TAC structure. Integrals are calculated at frame end times.
int tacInterpolate (TAC *inp, TAC *xinp, TAC *tac, TAC *itac, TAC *iitac, TPCSTATUS *status)
 Interpolate and/or integrate TACs from one TAC structure into a new TAC structure, using sample times from another TAC structure which is not changed.
int tacInterpolateInto (TAC *inp, TAC *tac, TAC *itac, TAC *iitac, TPCSTATUS *status)
 Add TACs from one TAC structure into another TAC structure, interpolating the input TACs and allocating space if necessary.
double tacAUC (TAC *tac, int ti, double t1, double t2, TPCSTATUS *status)
 Calculates TAC AUC from t1 to t2.

Detailed Description

Linear interpolation of TACs.

Definition in file litac.c.

Function Documentation

◆ tacAUC()

double tacAUC ( TAC * tac,
int ti,
double t1,
double t2,
TPCSTATUS * status )

Calculates TAC AUC from t1 to t2.

See also
tacInterpolate, liInterpolate, liIntegrate
Returns
AUC, or NaN in case of an error.
Parameters
tacPointer to TAC structure.
Precondition
Data must be sorted by increasing x. No gaps or overlap between frames allowed.
See also
tacSortByTime, tacSetXContiguous
Parameters
tiIndex [0..tacNr-1] of the TAC to calculate AUC from.
t1AUC start time.
t2AUC end time.
statusPointer to status data; enter NULL if not needed.

Definition at line 486 of file litac.c.

500 {
501 int verbose=0; if(status!=NULL) verbose=status->verbose;
502 if(verbose>0) printf("%s(tac, %d, %g, %g)\n", __func__, ti, t1, t2);
503 /* Check the function input */
504 if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
505 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
506 return(nan(""));
507 }
508 if(ti<0 || ti>=tac->tacNr) {
509 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
510 return(nan(""));
511 }
512 if(t1>=t2) {
513 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
514 return(nan(""));
515 }
516
517 /* Get and check the x range */
518 double xmin=nan(""), xmax=nan("");
519 if(tacXRange(tac, &xmin, &xmax)!=0 || t1>=xmax || t2<=xmin) {
520 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
521 return(nan(""));
522 }
523
524 /* Make temporary x and yi array */
525 double tx[2]; tx[0]=t1; tx[1]=t2;
526 double tyi[2];
527
528 /* If TAC contains just frame mid times, then simple dot-to-dot integration */
529 if(tac->isframe==0) {
530 if(liInterpolate(tac->x, tac->c[ti].y, tac->sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
531 {
532 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
533 return(nan(""));
534 }
535 /* That's it then */
536 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
537 return(tyi[1]-tyi[0]);
538 }
539
540 /* We have frame start and end times */
541 /* Convert the data into stepwise dot-to-dot data */
542 TAC stac; tacInit(&stac);
543 {
544 int ret=tacFramesToSteps(tac, &stac, NULL);
545 if(ret!=TPCERROR_OK) {
546 tacFree(&stac);
547 statusSet(status, __func__, __FILE__, __LINE__, ret);
548 return(nan(""));
549 }
550 }
551 /* Then simple dot-to-dot integration */
552 if(liInterpolate(stac.x, stac.c[ti].y, stac.sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
553 {
554 tacFree(&stac);
555 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
556 return(nan(""));
557 }
558 tacFree(&stac);
559
560 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
561 return(tyi[1]-tyi[0]);
562}
liInterpolate
int liInterpolate(double *x, double *y, const int nr, double *newx, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
Linear interpolation and/or integration with trapezoidal method.
Definition interpolate.c:139
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
TACC::y
double * y
Definition tpctac.h:75
TAC
Definition tpctac.h:87
TAC::x
double * x
Definition tpctac.h:97
TAC::sampleNr
int sampleNr
Definition tpctac.h:89
TAC::isframe
int isframe
Definition tpctac.h:95
TAC::c
TACC * c
Definition tpctac.h:117
TAC::tacNr
int tacNr
Definition tpctac.h:91
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
tacFree
void tacFree(TAC *tac)
Definition tac.c:106
tacInit
void tacInit(TAC *tac)
Definition tac.c:24
tacFramesToSteps
int tacFramesToSteps(TAC *inp, TAC *out, TPCSTATUS *status)
Transform TAC with frames into TAC with frames represented with stepwise changing dot-to-dot data.
Definition tacx.c:942
tacXRange
int tacXRange(TAC *d, double *xmin, double *xmax)
Get the range of x values (times) in TAC structure.
Definition tacx.c:124
TPCERROR_INVALID_XRANGE
@ TPCERROR_INVALID_XRANGE
Invalid sample time range.
Definition tpcextensions.h:214
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203

◆ tacIntegrate()

int tacIntegrate ( TAC * inp,
TAC * out,
TPCSTATUS * status )

Integrate TACs from one TAC structure into a new TAC structure.

PET frame lengths are taken into account, if available.

See also
liIntegrate, liIntegratePET, tacInterpolate, tacInterpolateInto
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
inpPointer to source TAC structure.
Precondition
Data must be sorted by increasing x.
Parameters
outPointer to target TAC structure into which the integrated TAC(s) are written. Structure must be initiated; any previous contents are deleted.
statusPointer to status data; enter NULL if not needed.

Definition at line 27 of file litac.c.

36 {
37 int verbose=0; if(status!=NULL) verbose=status->verbose;
38 if(verbose>0) printf("%s()\n", __func__);
39 if(inp==NULL || out==NULL) {
40 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
41 return TPCERROR_FAIL;
42 }
43 tacFree(out);
44 int ret=TPCERROR_OK;
45 ret=tacDuplicate(inp, out);
46 if(ret!=TPCERROR_OK) {
47 statusSet(status, __func__, __FILE__, __LINE__, ret);
48 return(ret);
49 }
50 if(inp->isframe==0) {
51 for(int i=0; i<inp->tacNr; i++) {
52 ret=liIntegrate(inp->x, inp->c[i].y, inp->sampleNr, out->c[i].y, 3, 0);
53 if(ret!=TPCERROR_OK) break;
54 }
55 } else {
56 for(int i=0; i<inp->tacNr; i++) {
57 ret=liIntegratePET(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, NULL, 0);
58 if(ret!=TPCERROR_OK) break;
59 }
60 }
61 if(ret!=TPCERROR_OK) tacFree(out);
62 statusSet(status, __func__, __FILE__, __LINE__, ret);
63 return(ret);
64}
liIntegrate
int liIntegrate(double *x, double *y, const int nr, double *yi, const int se, const int verbose)
Linear integration of TAC with trapezoidal method.
Definition integrate.c:33
liIntegratePET
int liIntegratePET(double *x1, double *x2, double *y, int nr, double *ie, double *iie, const int verbose)
Calculate PET TAC AUC from start to each time frame, as averages during each frame.
Definition integrate.c:120
TAC::x2
double * x2
Definition tpctac.h:101
TAC::x1
double * x1
Definition tpctac.h:99
tacDuplicate
int tacDuplicate(TAC *tac1, TAC *tac2)
Make a duplicate of TAC structure.
Definition tac.c:356

◆ tacIntegrateFE()

int tacIntegrateFE ( TAC * inp,
TAC * out,
TAC * out2,
TPCSTATUS * status )

Integrate TACs from one TAC structure into a new TAC structure. Integrals are calculated at frame end times.

TAC must contain frame start and end times.

See also
tacIntegrate, liIntegrateFE, tacIsXContiguous, tacSetXContiguous, tacMultipleSamples
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
inpPointer to source TAC structure.
Precondition
Data must be sorted by increasing x. Frames must not overlap.
Parameters
outPointer to target TAC structure into which the integrated TAC(s) are written. Structure must be initiated; any previous contents are deleted.
out2Pointer to target TAC structure into which the 2nd integrals are written. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
statusPointer to status data; enter NULL if not needed.

Definition at line 75 of file litac.c.

88 {
89 int verbose=0; if(status!=NULL) verbose=status->verbose;
90 if(verbose>0) printf("%s()\n", __func__);
91 if(inp==NULL || out==NULL) {
92 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
93 return TPCERROR_FAIL;
94 }
95 if(inp->isframe==0) {
96 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
97 return(TPCERROR_INVALID_X);
98 }
99
100 /* Make space */
101 tacFree(out); if(out2!=NULL) tacFree(out2);
102 int ret=TPCERROR_OK;
103 ret=tacDuplicate(inp, out);
104 if(ret!=TPCERROR_OK) {
105 statusSet(status, __func__, __FILE__, __LINE__, ret);
106 return(ret);
107 }
108 if(out2!=NULL) {
109 ret=tacDuplicate(inp, out2);
110 if(ret!=TPCERROR_OK) {
111 tacFree(out);
112 statusSet(status, __func__, __FILE__, __LINE__, ret);
113 return(ret);
114 }
115 }
116
117 for(int i=0; i<inp->tacNr; i++) {
118 if(out2==NULL)
119 ret=liIntegrateFE(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, NULL, 0);
120 else
121 ret=liIntegrateFE(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, out2->c[i].y, 0);
122 if(ret!=TPCERROR_OK) break;
123 }
124 if(ret!=TPCERROR_OK) {tacFree(out); if(out2!=NULL) tacFree(out2);}
125
126 statusSet(status, __func__, __FILE__, __LINE__, ret);
127 return(ret);
128}
liIntegrateFE
int liIntegrateFE(double *x1, double *x2, double *y, int nr, double *ie, double *iie, const int verbose)
Linear integration of PET TAC to frame end times.
Definition integrate.c:219
TPCERROR_INVALID_X
@ TPCERROR_INVALID_X
Invalid sample time.
Definition tpcextensions.h:213

◆ tacInterpolate()

int tacInterpolate ( TAC * inp,
TAC * xinp,
TAC * tac,
TAC * itac,
TAC * iitac,
TPCSTATUS * status )

Interpolate and/or integrate TACs from one TAC structure into a new TAC structure, using sample times from another TAC structure which is not changed.

PET frame lengths of output TACs are taken into account in interpolation, if available. Input frame lengths are taken into account if the framing is same as with output TACs, otherwise frame middle times are used.

See also
tacInterpolateInto, tacAUC, tacFramesToSteps, tacXCopy, liInterpolate
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
inpPointer to source TAC structure; make sure that x (time) and y (concentration) units are the same as in output TAC structure, because input time range is verified to cover the output data range. Detailed verification of suitability of TAC for interpolation must be done elsewhere.
Precondition
Data must be sorted by increasing x.
Parameters
xinpPointer to TAC structure which contains the x values (sample times) for the interpolation and integration. Data in this TAC structure is not modified.
Precondition
Data must be sorted by increasing x.
Parameters
tacPointer to target TAC structure into which the interpolated input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
itacPointer to target TAC structure into which the integrated input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
iitacPointer to target TAC structure into which 2nd integrals of input TAC(s) are added. Structure must be initiated; any previous contents are deleted. Enter NULL if not needed.
statusPointer to status data; enter NULL if not needed.

Definition at line 141 of file litac.c.

163 {
164 int verbose=0; if(status!=NULL) verbose=status->verbose;
165 if(verbose>0) printf("%s()\n", __func__);
166 if(inp==NULL || xinp==NULL || xinp==tac || xinp==itac || xinp==iitac) {
167 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
168 return TPCERROR_FAIL;
169 }
170 if(tac==NULL && itac==NULL && iitac==NULL) { // there's nothing to do
171 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
172 return TPCERROR_OK;
173 }
174 /* Check the function input */
175 if(inp->sampleNr<1 || inp->tacNr<1 || xinp->sampleNr<1) {
176 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
177 return TPCERROR_NO_DATA;
178 }
179
180 /* Check that input does not have any missing values */
181 if(tacNaNs(inp)>0 || tacXNaNs(xinp)>0) {
182 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
183 return TPCERROR_MISSING_DATA;
184 }
185
186
187 /* If input data have similar x values as requested already, then
188 copy x1 and x2 into input if available */
189 if(xinp->isframe && tacXMatch(inp, xinp, verbose-2) && inp->sampleNr<=xinp->sampleNr) {
190 for(int i=0; i<inp->sampleNr; i++) {
191 inp->x[i]=xinp->x[i]; inp->x1[i]=xinp->x1[i]; inp->x2[i]=xinp->x2[i];
192 }
193 inp->isframe=xinp->isframe;
194 }
195
196
197 /* Check that there is no need for excess extrapolation */
198 {
199 double start1, end1, start2, end2, range2;
200 if(tacXRange(inp, &start1, &end1) || tacXRange(xinp, &start2, &end2)) {
201 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
202 return TPCERROR_INVALID_X;
203 }
204 range2=end2-start2;
205 if(verbose>2) printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
206 if(start1>start2+0.2*range2 || end1<end2-0.3*range2) {
207 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
208 return TPCERROR_INVALID_XRANGE;
209 }
210 }
211
212 /* Delete any previous contents in output structs */
213 if(tac!=NULL) tacFree(tac);
214 if(itac!=NULL) tacFree(itac);
215 if(iitac!=NULL) tacFree(iitac);
216
217 /* Allocate memory for output data */
218 int ret=TPCERROR_OK;
219 if(tac!=NULL) ret=tacAllocate(tac, xinp->sampleNr, inp->tacNr);
220 if(ret==TPCERROR_OK && itac!=NULL)
221 ret=tacAllocate(itac, xinp->sampleNr, inp->tacNr);
222 if(ret==TPCERROR_OK && iitac!=NULL)
223 ret=tacAllocate(iitac, xinp->sampleNr, inp->tacNr);
224 if(ret!=TPCERROR_OK) {
225 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
226 return TPCERROR_OUT_OF_MEMORY;
227 }
228 if(tac!=NULL) {tac->tacNr=inp->tacNr; tac->sampleNr=xinp->sampleNr;}
229 if(itac!=NULL) {itac->tacNr=inp->tacNr; itac->sampleNr=xinp->sampleNr;}
230 if(iitac!=NULL) {iitac->tacNr=inp->tacNr; iitac->sampleNr=xinp->sampleNr;}
231
232
233 /* Copy header information */
234 if(tac!=NULL) {
235 tacCopyHdr(inp, tac);
236 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
237 }
238 if(itac!=NULL) {
239 tacCopyHdr(inp, itac);
240 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, itac->c+i);
241 }
242 if(iitac!=NULL) {
243 tacCopyHdr(inp, iitac);
244 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+i);
245 }
246 /* Copy X (time) information */
247 if(tac!=NULL) {
248 tac->isframe=xinp->isframe; tacXCopy(xinp, tac, 0, tac->sampleNr-1);
249 }
250 if(itac!=NULL) {
251 itac->isframe=xinp->isframe; tacXCopy(xinp, itac, 0, itac->sampleNr-1);
252 }
253 if(iitac!=NULL) {
254 iitac->isframe=xinp->isframe; tacXCopy(xinp, iitac, 0, iitac->sampleNr-1);
255 }
256
257 /* Check if TACs have the same frame times already */
258 if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=xinp->sampleNr) {
259 if(verbose>2) printf(" same frame times\n");
260 int ret=0, ri, fi;
261 double *i1, *i2;
262 /* copy the values directly and integrate in place, if requested */
263 for(ri=0; ri<inp->tacNr && ret==0; ri++) {
264 if(tac!=NULL)
265 for(fi=0; fi<tac->sampleNr; fi++)
266 tac->c[ri].y[fi]=inp->c[ri].y[fi];
267 /* If integrals not requested then that's it for this TAC */
268 if(itac==NULL && iitac==NULL) continue;
269 /* otherwise set pointers for integrals and then integrate */
270 if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;
271 if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
272 if(inp->isframe)
273 ret=liIntegratePET(inp->x1, inp->x2, inp->c[ri].y,
274 xinp->sampleNr, i1, i2, verbose-2);
275 else
276 ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
277 xinp->x, NULL, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
278 }
279 if(ret) {
280 if(tac!=NULL) tacFree(tac);
281 if(itac!=NULL) tacFree(itac);
282 if(iitac!=NULL) tacFree(iitac);
283 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
284 return TPCERROR_INVALID_VALUE;
285 }
286 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
287 return(TPCERROR_OK); // that's it then
288 }
289
290 /* Interpolate and/or integrate inp data to tac sample times,
291 taking into account tac frame lengths if available */
292 if(verbose>2) printf(" different frame times\n");
293 ret=0;
294 double *i0, *i1, *i2;
295 for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
296 /* set pointers for output y arrays */
297 if(tac!=NULL) i0=tac->c[ri].y; else i0=NULL;
298 if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;
299 if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
300 if(xinp->isframe)
301 ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
302 xinp->x1, xinp->x2, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
303 else
304 ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
305 xinp->x, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
306 }
307 if(ret) {
308 if(tac!=NULL) tacFree(tac);
309 if(itac!=NULL) tacFree(itac);
310 if(iitac!=NULL) tacFree(iitac);
311 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
312 return TPCERROR_INVALID_VALUE;
313 }
314
315 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
316 return(TPCERROR_OK);
317}
liInterpolateForPET
int liInterpolateForPET(double *x, double *y, const int nr, double *newx1, double *newx2, double *newy, double *newyi, double *newyii, const int newnr, const int se, const int ee, const int verbose)
Linear TAC interpolation and/or integration to PET frames.
Definition interpolate.c:317
tacAllocate
int tacAllocate(TAC *tac, int sampleNr, int tacNr)
Definition tac.c:130
tacCopyTacchdr
int tacCopyTacchdr(TACC *d1, TACC *d2)
Definition tac.c:282
tacCopyHdr
int tacCopyHdr(TAC *tac1, TAC *tac2)
Copy TAC header data from tac1 to tac2.
Definition tac.c:310
tacNaNs
int tacNaNs(TAC *tac)
Definition tacnan.c:71
tacXNaNs
int tacXNaNs(TAC *tac)
Definition tacnan.c:23
tacXMatch
int tacXMatch(TAC *d1, TAC *d2, const int verbose)
Check whether sample (frame) times are the same (or very close to) in two TAC structures.
Definition tacx.c:249
tacXCopy
int tacXCopy(TAC *tac1, TAC *tac2, int i1, int i2)
Definition tacx.c:24
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_MISSING_DATA
@ TPCERROR_MISSING_DATA
File contains missing values.
Definition tpcextensions.h:207

◆ tacInterpolateInto()

int tacInterpolateInto ( TAC * inp,
TAC * tac,
TAC * itac,
TAC * iitac,
TPCSTATUS * status )

Add TACs from one TAC structure into another TAC structure, interpolating the input TACs and allocating space if necessary.

PET frame lengths of output TACs are taken into account in interpolation, if available. Input frame lengths are taken into account if the framing is same as with output TACs, otherwise frame middle times are used.

See also
tacInterpolate, tacAUC, tacInterpolateToEqualLengthFrames, liInterpolateForPET
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
inpPointer to source TAC struct; make sure that x (time) and y (concentration) units are the same as in output TAC struct, because input time range is verified to cover the output data range. Detailed verification of suitability of TAC for interpolation must be done elsewhere.
Precondition
Data must be sorted by increasing x.
Parameters
tacPointer to target TAC struct into which the input TAC(s) are added. Struct must contain at least one TAC, and the x or x1 and x2 values.
Precondition
Data must be sorted by increasing x.
Parameters
itacPointer to target TAC struct into which integrated input TAC(s) are added; enter NULL if not needed.
iitacPointer to target TAC struct into which 2nd integrals of input TAC(s) are added; enter NULL if not needed.
statusPointer to status data; enter NULL if not needed.

Definition at line 330 of file litac.c.

348 {
349 int verbose=0; if(status!=NULL) verbose=status->verbose;
350 if(verbose>0) printf("%s()\n", __func__);
351 if(tac==NULL || inp==NULL) {
352 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
353 return TPCERROR_FAIL;
354 }
355
356 /* Check the function input */
357 if(tac->sampleNr<1 || inp->sampleNr<1 || inp->tacNr<1) {
358 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
359 return TPCERROR_NO_DATA;
360 }
361 if(itac!=NULL && itac->sampleNr!=tac->sampleNr) {
362 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
363 return TPCERROR_INVALID_X;
364 }
365 if(iitac!=NULL && iitac->sampleNr!=tac->sampleNr) {
366 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
367 return TPCERROR_INVALID_X;
368 }
369
370
371 /* Check that input does not have any missing values */
372 /* Check that there are no missing sample times in output TACs */
373 if(tacNaNs(inp)>0 || tacXNaNs(tac)>0) {
374 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
375 return TPCERROR_MISSING_DATA;
376 }
377
378 /* Check that there is no need for excess extrapolation */
379 {
380 double start1, end1, start2, end2, range2;
381 if(tacXRange(inp, &start1, &end1) || tacXRange(tac, &start2, &end2)) {
382 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
383 return TPCERROR_INVALID_X;
384 }
385 range2=end2-start2;
386 if(verbose>2)
387 printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
388 if(start1>(start2+0.2*range2) || end1<(end2-0.3*range2)) {
389 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
390 return TPCERROR_INVALID_XRANGE;
391 }
392 }
393
394 /* Allocate memory for interpolated data */
395 if(tacAllocateMore(tac, inp->tacNr)) {
396 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
397 return TPCERROR_OUT_OF_MEMORY;
398 }
399 /* Copy header information */
400 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+tac->tacNr+i);
401 /* Same for the integrals */
402 if(itac!=NULL) {
403 if(tacAllocateMore(itac, inp->tacNr)) {
404 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
405 return TPCERROR_OUT_OF_MEMORY;
406 }
407 for(int i=0; i<inp->tacNr; i++)
408 tacCopyTacchdr(inp->c+i, itac->c+itac->tacNr+i);
409 }
410 if(iitac!=NULL) {
411 if(tacAllocateMore(iitac, inp->tacNr)) {
412 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
413 return TPCERROR_OUT_OF_MEMORY;
414 }
415 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+iitac->tacNr+i);
416 }
417
418 /* Check if TACs have the same frame times already */
419 if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=tac->sampleNr) {
420 if(verbose>2) printf("same frame times\n");
421 int ret=0, ri, fi;
422 double *i1, *i2;
423 /* copy the values directly and integrate in place, if requested */
424 for(ri=0; ri<inp->tacNr && ret==0; ri++) {
425 for(fi=0; fi<tac->sampleNr; fi++)
426 tac->c[tac->tacNr+ri].y[fi]=inp->c[ri].y[fi];
427 /* If integrals not requested then that's it for this TAC */
428 if(itac==NULL && iitac==NULL) continue;
429 /* otherwise set pointers for integrals and then integrate */
430 if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;
431 if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
432 if(tac->isframe)
433 ret=liIntegratePET(tac->x1, tac->x2, tac->c[tac->tacNr+ri].y,
434 tac->sampleNr, i1, i2, verbose-2);
435 else
436 ret=liInterpolate(tac->x, tac->c[tac->tacNr+ri].y, tac->sampleNr,
437 tac->x, NULL, i1, i2, tac->sampleNr, 4, 1, verbose-2);
438 }
439 if(ret) {
440 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
441 return TPCERROR_INVALID_VALUE;
442 }
443 tac->tacNr+=inp->tacNr;
444 if(itac!=NULL) itac->tacNr+=inp->tacNr;
445 if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
446 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
447 return(TPCERROR_OK); // that's it then
448 }
449
450 /* Interpolate and, if requested, integrate inp data to tac sample times,
451 taking into account tac frame lengths if available */
452 int ret=0;
453 double *i1, *i2;
454 for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
455 /* set pointers for integrals */
456 if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;
457 if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
458 if(tac->isframe)
459 ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
460 tac->x1, tac->x2, tac->c[tac->tacNr+ri].y,
461 i1, i2, tac->sampleNr, 4, 1, verbose-2);
462 else
463 ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
464 tac->x, tac->c[tac->tacNr+ri].y,
465 i1, i2, tac->sampleNr, 4, 1, verbose-2);
466 }
467 if(ret) {
468 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
469 return TPCERROR_INVALID_VALUE;
470 }
471 tac->tacNr+=inp->tacNr;
472 if(itac!=NULL) itac->tacNr+=inp->tacNr;
473 if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
474
475 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
476 return(TPCERROR_OK);
477}
tacAllocateMore
int tacAllocateMore(TAC *tac, int tacNr)
Definition tac.c:178

Referenced by tacReadModelingData(), and tacReadModelingInput().

Generated on for TPCCLIB by doxygen 1.15.0