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

Linear interpolation for simulation. More...

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

Go to the source code of this file.

Functions

int tacInput2sim (TAC *itac, TAC *ttac, TAC *stac, TPCSTATUS *status)
 Modify input TAC based on tissue TAC, for use in simulations.
int tacVb (TAC *ttac, const int i, TAC *btac, double Vb, const int simVb, const int petVolume, TPCSTATUS *status)
 Correct TTACs for vascular blood, or simulate its effect.
int tacInterpolateToEqualLengthFrames (TAC *inp, double minfdur, double maxfdur, TAC *tac, TPCSTATUS *status)

Detailed Description

Linear interpolation for simulation.

Definition in file lisim.c.

Function Documentation

◆ tacInput2sim()

int tacInput2sim ( TAC * itac,
TAC * ttac,
TAC * stac,
TPCSTATUS * status )

Modify input TAC based on tissue TAC, for use in simulations.

Note
NOT FUNCTIONAL/TESTED YET!
See also
tacReadModelingInput, tacReadModelingData, tacInterpolateInto, tacInterpolateToEqualLengthFrames
Returns
enum tpcerror (TPCERROR_OK when successful).
Author
Vesa Oikonen
Parameters
itacPointer to original input TAC structure. Not modified. Data must not have missing values.
ttacPointer to tissue TAC structure. Not modified. Data must not have missing values.
stacPointer to the new input data, containing input data with new sample times.
statusPointer to status data; obligatory.

Definition at line 29 of file lisim.c.

38 {
39 if(status==NULL) return TPCERROR_FAIL;
40 int verbose=status->verbose;
41 if(verbose>0) printf("%s(itac, ttac, ...)\n", __func__);
42
43 /* Check the input */
44 if(itac==NULL || ttac==NULL || stac==NULL) {
45 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
46 return TPCERROR_FAIL;
47 }
48 if(itac->sampleNr<1 || ttac->sampleNr<1 || itac->tacNr<1) {
49 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
50 return TPCERROR_NO_DATA;
51 }
52
53 /* Delete any previous data */
54 tacFree(stac);
55
56 /* Make a working copy of input data */
57 TAC inp; tacInit(&inp);
58 status->error=tacDuplicate(itac, &inp); if(status->error!=TPCERROR_OK) return(status->error);
59 /* Make sure that input data contains also frame mid times */
60 if(inp.isframe) tacSetX(&inp, NULL);
61 /* Convert time units to TTAC units */
62 tacXUnitConvert(&inp, ttac->tunit, status);
63 if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
64
65 /* Sort the input data by sample times */
66 tacSortByTime(&inp, status); if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
67 /* Add zero sample if necessary */
68 tacAddZeroSample(&inp, status);
69 if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
70
71 /* Make a combined list of sample times based on input and tissue TACs */
72 unsigned int maxn=inp.sampleNr+ttac->sampleNr;
73 if(inp.isframe) maxn+=2*inp.sampleNr;
74 if(ttac->isframe) maxn+=ttac->sampleNr;
75 double x[maxn];
76
77 int n=0;
78 for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x[i])) {
79 int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x[i])<1.0E-20) break;
80 if(j==n) x[n++]=inp.x[i];
81 }
82 if(inp.isframe) {
83 for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x1[i])) {
84 int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x1[i])<1.0E-20) break;
85 if(j==n) x[n++]=inp.x1[i];
86 }
87 for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x2[i])) {
88 int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x2[i])<1.0E-20) break;
89 if(j==n) x[n++]=inp.x2[i];
90 }
91 }
92 if(ttac->isframe) {
93 for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
94 int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
95 if(j==n) x[n++]=ttac->x2[i];
96 }
97 for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
98 int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
99 if(j==n) x[n++]=ttac->x2[i];
100 }
101 } else {
102 for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x[i])) {
103 int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x[i])<1.0E-20) break;
104 if(j==n) x[n++]=ttac->x[i];
105 }
106 }
107 if(verbose>2) printf("n=%d\n", n);
108 statSortDouble(x, n, 0);
109
110 /* Put new sample times into output data */
111 status->error=tacDuplicate(itac, stac);
112 if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
113 status->error=tacAllocateMoreSamples(stac, n-stac->sampleNr);
114 if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
115 stac->sampleNr=n; stac->tunit=inp.tunit; stac->isframe=0;
116 for(int i=0; i<n; i++) stac->x[i]=x[i];
117
118 /* Interpolate input data into the new sample times */
119 for(int i=0; i<inp.tacNr; i++) {
120 status->error=liInterpolate(inp.x, inp.c[i].y, inp.sampleNr,
121 stac->x, stac->c[i].y, NULL, NULL, stac->sampleNr, 4, 1, verbose-10);
122 if(status->error!=TPCERROR_OK) break;
123 }
124 if(status->error!=TPCERROR_OK) {tacFree(&inp); tacFree(stac); return(status->error);}
125
126 tacFree(&inp);
127
128 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
129 return(TPCERROR_OK);
130}
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
statSortDouble
void statSortDouble(double *data, unsigned int n, int order)
Definition sort.c:99
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::x2
double * x2
Definition tpctac.h:101
TAC::tunit
unit tunit
Definition tpctac.h:109
TAC::x1
double * x1
Definition tpctac.h:99
TAC::tacNr
int tacNr
Definition tpctac.h:91
TPCSTATUS::verbose
int verbose
Verbose level, used by statusPrint() etc.
Definition tpcextensions.h:242
TPCSTATUS::error
tpcerror error
Error code.
Definition tpcextensions.h:245
tacFree
void tacFree(TAC *tac)
Definition tac.c:106
tacDuplicate
int tacDuplicate(TAC *tac1, TAC *tac2)
Make a duplicate of TAC structure.
Definition tac.c:356
tacInit
void tacInit(TAC *tac)
Definition tac.c:24
tacAllocateMoreSamples
int tacAllocateMoreSamples(TAC *tac, int addNr)
Allocate memory for more samples in TAC data.
Definition tac.c:435
tacSortByTime
int tacSortByTime(TAC *d, TPCSTATUS *status)
Definition tacorder.c:74
tacXUnitConvert
int tacXUnitConvert(TAC *tac, const int u, TPCSTATUS *status)
Definition tacunits.c:23
tacAddZeroSample
int tacAddZeroSample(TAC *d, TPCSTATUS *status)
Add an initial sample to TAC(s) with zero time and concentration.
Definition tacx.c:366
tacSetX
int tacSetX(TAC *d, TPCSTATUS *status)
Set TAC x values based on x1 and x2 values, or guess x1 and x2 values based on x values.
Definition tacx.c:653
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

Referenced by bfm1TCM().

◆ tacInterpolateToEqualLengthFrames()

int tacInterpolateToEqualLengthFrames ( TAC * inp,
double minfdur,
double maxfdur,
TAC * tac,
TPCSTATUS * status )

Interpolate TAC data into data with equal frame lengths, starting from zero.

By default, the shortest frame length or sampling interval in the input data is used.

Remarks
This is needed for TAC convolution.
See also
tacInterpolate, liInterpolateForPET, tacGetSampleInterval, simIsSteadyInterval
Returns
enum tpcerror (TPCERROR_OK when successful).
Parameters
inpPointer to source TAC struct. Missing values are not allowed.
Precondition
Data must be sorted by increasing x.
Parameters
minfdurMinimum frame length; if NaN or <=0 then not applied, otherwise frame lengths are not allowed to be lower.
maxfdurMaximum frame length; if NaN or <=0 then not applied, otherwise frame lengths are not allowed to be higher.
tacPointer to target TAC struct into which the interpolated TACs are written.
Precondition
Struct must be initiated.
Parameters
statusPointer to status data; enter NULL if not needed.

Definition at line 214 of file lisim.c.

229 {
230 int verbose=0; if(status!=NULL) verbose=status->verbose;
231 if(verbose>0) {printf("%s(tac1, %g, %g, ...)\n", __func__, minfdur, maxfdur); fflush(stdout);}
232 if(tac==NULL || inp==NULL || inp->sampleNr<1 || inp->tacNr<1) {
233 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
234 return TPCERROR_FAIL;
235 }
236 if(minfdur>0.0 && maxfdur>0.0 && minfdur>maxfdur) {
237 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
238 return TPCERROR_INVALID_X;
239 }
240 if(!tacIsX(inp)) {
241 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
242 return TPCERROR_INVALID_X;
243 }
244
245 /* Get the minimum sample interval (that is higher than 0) */
246 double freq=nan("");
247 int ret=tacGetSampleInterval(inp, 1.0E-08, &freq, NULL);
248 //printf("dataset based freq: %g\n", freq);
249 if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
250 if(minfdur>0.0 && freq<minfdur) freq=minfdur;
251 if(maxfdur>0.0 && freq>maxfdur) freq=maxfdur;
252 //printf("refined freq: %g\n", freq);
253
254 /* Get the x range */
255 double xmax=nan("");
256 if(tacXRange(inp, NULL, &xmax)!=0) {
257 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
258 return TPCERROR_INVALID_XRANGE;
259 }
260 //printf("xmax: %g\n", xmax);
261 if(freq>0.5*xmax) {
262 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
263 return TPCERROR_INVALID_XRANGE;
264 }
265
266
267 /* Calculate the number of required interpolated samples */
268 int nreq=ceil(xmax/freq); //printf("required_n := %d\n", nreq);
269 /* Check that it is not totally stupid */
270 if(nreq<1 || nreq>2E+06) {
271 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_BIG);
272 return TPCERROR_TOO_BIG;
273 }
274
275 /* Setup the output TAC */
276 if(tacAllocate(tac, nreq, inp->tacNr)!=TPCERROR_OK) {
277 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
278 return TPCERROR_OUT_OF_MEMORY;
279 }
280 tac->tacNr=inp->tacNr; tac->sampleNr=nreq;
281 /* Copy header information */
282 tacCopyHdr(inp, tac);
283 for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
284
285 /* Set X (sample time) information */
286 tac->isframe=1;
287 for(int i=0; i<nreq; i++) tac->x1[i]=(double)i*freq;
288 for(int i=0; i<nreq; i++) tac->x2[i]=(double)(i+1)*freq;
289 for(int i=0; i<nreq; i++) tac->x[i]=0.5*(tac->x1[i]+tac->x2[i]);
290
291 /* Calculate mean value during each interval */
292 if(verbose>2) {printf(" interpolating...\n"); fflush(stdout);}
293 for(int i=0; i<inp->tacNr; i++) {
294 if(liInterpolateForPET(inp->x, inp->c[i].y, inp->sampleNr, tac->x1, tac->x2,
295 tac->c[i].y, NULL, NULL, tac->sampleNr, 2, 1, 0)!=0) {
296 tacFree(tac);
297 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
298 return TPCERROR_INVALID_XRANGE;
299 }
300 }
301 if(verbose>2) {printf(" ... done.\n"); fflush(stdout);}
302
303 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
304 return(TPCERROR_OK);
305}
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
tacGetSampleInterval
int tacGetSampleInterval(TAC *d, double ilimit, double *minfdur, double *maxfdur)
Get the shortest and longest sampling intervals or frame lengths in TAC structure.
Definition tacx.c:832
tacIsX
int tacIsX(TAC *d)
Verify if TAC structure contains reasonable x values (times).
Definition tacx.c:226
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_TOO_BIG
@ TPCERROR_TOO_BIG
File is too big.
Definition tpcextensions.h:209
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_INVALID_X
@ TPCERROR_INVALID_X
Invalid sample time.
Definition tpcextensions.h:213

◆ tacVb()

int tacVb ( TAC * ttac,
const int i,
TAC * btac,
double Vb,
const int simVb,
const int petVolume,
TPCSTATUS * status )

Correct TTACs for vascular blood, or simulate its effect.

See also
tacDelay, tacReadModelingInput, tacInterpolate, imgVb
Returns
enum tpcerror (TPCERROR_OK when successful).
Parameters
ttacPointer to TAC data to process.
iIndex of TAC to be processed; enter <0 to process all.
btacPointer to BTAC data to subtract or add; must contain exactly same sample times as TTAC and y values must be in the same units.
VbVb fraction [0,1].
simVbSwitch to either subtract vascular volume (0) or to simulate it (1).
petVolumeSwitch to model vascular volume as either 0 : Cpet = (1-Vb)*Ct + Vb*Cb, or 1 : Cpet = Ct + Vb*Cb
statusPointer to status data; obligatory.

Definition at line 138 of file lisim.c.

156 {
157 int verbose=0; if(status!=NULL) verbose=status->verbose;
158 if(verbose>0) printf("%s(ttac, %d, btac, %g, %d, %d)\n", __func__, i, Vb, simVb, petVolume);
159 if(ttac==NULL || btac==NULL) {
160 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
161 return TPCERROR_FAIL;
162 }
163 if(ttac->tacNr<1 || btac->tacNr<1 || ttac->sampleNr<1) {
164 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
165 return TPCERROR_NO_DATA;
166 }
167 if(ttac->sampleNr>btac->sampleNr) {
168 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
169 return TPCERROR_INCOMPATIBLE_DATA;
170 }
171 if(!(Vb>=0.0) || !(Vb<=1.0)) {
172 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
173 return TPCERROR_INVALID_VALUE;
174 }
175 if(Vb==1.0 && petVolume==0 && simVb==0) { // combination would lead to divide by zero
176 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
177 return TPCERROR_INVALID_VALUE;
178 }
179 if(i>=ttac->tacNr) {
180 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
181 return TPCERROR_INVALID_VALUE;
182 }
183
184
185 if(simVb==0) {
186 if(verbose>1) printf("subtracting blood\n");
187 for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
188 for(int fi=0; fi<ttac->sampleNr; fi++) {
189 ttac->c[ri].y[fi]-=Vb*btac->c[0].y[fi];
190 if(petVolume==0) ttac->c[ri].y[fi]/=(1.0-Vb);
191 }
192 }
193 } else {
194 if(verbose>1) printf("adding blood\n");
195 for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
196 for(int fi=0; fi<ttac->sampleNr; fi++) {
197 if(petVolume==0) ttac->c[ri].y[fi]*=(1.0-Vb);
198 ttac->c[ri].y[fi]+=Vb*btac->c[0].y[fi];
199 }
200 }
201 }
202
203 return(TPCERROR_OK);
204}
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_INCOMPATIBLE_DATA
@ TPCERROR_INCOMPATIBLE_DATA
Incompatible data.
Definition tpcextensions.h:216
Generated on for TPCCLIB by doxygen 1.15.0