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

Functions for calculation of basis functions for PET modelling. More...

#include "libtpcmodext.h"

Go to the source code of this file.

Functions

int bf_srtm (double *t, double *cr, int n, int bfNr, double t3min, double t3max, DFT *bf)
 
int bfRadiowater (DFT *input, DFT *tissue, DFT *bf, int bfNr, double k2min, double k2max, char *status, int verbose)
 
int bfIrr2TCM (DFT *input, DFT *tissue, DFT *bf, int bfNr, double thetamin, double thetamax, char *status, int verbose)
 

Detailed Description

Functions for calculation of basis functions for PET modelling.

Author
Vesa Oikonen

Definition in file bf_model.c.

Function Documentation

◆ bf_srtm()

int bf_srtm ( double * t,
double * cr,
int n,
int bfNr,
double t3min,
double t3max,
DFT * bf )

Calculates set of basis functions for SRTM.

Returns
Returns 0 if successful, otherwise non-zero.
See also
bfRadiowater, bfIrr2TCM
Parameters
tPET frame mid times.
crNon-decay corrected Cr(t).
nNr of PET frames.
bfNrNr of basis functions to calculate.
t3mintheta3 min.
t3maxtheta3 max.
bfdata for basis functions is allocated and filled here.

Definition at line 16 of file bf_model.c.

31 {
32 int bi, fi, ret;
33 double a, b, c;
34
35 // printf("\n%s(*t, *cr, %d, %d, %g, %g, *bf)\n", __func__, n, bfNr, t3min, t3max);
36
37 /* Check the parameters */
38 if(t==NULL || cr==NULL || n<2 || bfNr<1 || t3min<1.0E-10 || t3min>=t3max) return(1);
39 if(bf==NULL || bf->voiNr>0) return(1);
40
41 /* Allocate meory for basis functions */
42 ret=dftSetmem(bf, n, bfNr); if(ret) return(2);
43
44 /* Copy and set information fields */
45 bf->voiNr=bfNr; bf->frameNr=n;
46 bf->_type=DFT_FORMAT_STANDARD;
47 for(bi=0; bi<bf->voiNr; bi++) {
48 snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
49 strcpy(bf->voi[bi].hemisphere, ".");
50 strcpy(bf->voi[bi].place, ".");
51 strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
52 }
53 for(fi=0; fi<bf->frameNr; fi++) bf->x[fi]=t[fi];
54
55 /* Compute theta3 values to size fields */
56 a=log10(t3min); b=log10(t3max); c=(b-a)/(double)(bfNr-1);
57 for(bi=0; bi<bf->voiNr; bi++) {
58 bf->voi[bi].size=pow(10.0, (double)bi*c+a);
59 }
60
61 /* Calculate the functions */
62 for(bi=0; bi<bf->voiNr; bi++) {
63 a=bf->voi[bi].size;
64 ret=simC1(t, cr, n, 1.0, a, bf->voi[bi].y);
65 if(ret) return(4);
66 }
67
68 return(0);
69}
dftSetmem
int dftSetmem(DFT *data, int frameNr, int voiNr)
Definition dft.c:57
DFT_FORMAT_STANDARD
#define DFT_FORMAT_STANDARD
Definition libtpccurveio.h:405
simC1
int simC1(double *t, double *ca, int nr, double k1, double k2, double *ct)
Definition simulate.c:1317
DFT::_type
int _type
Definition libtpccurveio.h:125
DFT::voi
Voi * voi
Definition libtpccurveio.h:108
DFT::voiNr
int voiNr
Definition libtpccurveio.h:78
DFT::frameNr
int frameNr
Definition libtpccurveio.h:76
DFT::x
double * x
Definition libtpccurveio.h:102
Voi::size
double size
Definition libtpccurveio.h:56
Voi::voiname
char voiname[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:50
Voi::y
double * y
Definition libtpccurveio.h:59
Voi::name
char name[MAX_REGIONNAME_LEN+1]
Definition libtpccurveio.h:48
Voi::hemisphere
char hemisphere[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:52
Voi::place
char place[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:54

◆ bfIrr2TCM()

int bfIrr2TCM ( DFT * input,
DFT * tissue,
DFT * bf,
int bfNr,
double thetamin,
double thetamax,
char * status,
int verbose )

Calculates set of basis functions for irreversible 2TCM.

Returns
Returns 0 if successful, otherwise non-zero.
See also
bf_srtm, bfRadiowater
Parameters
inputArterial PTAC (not modified).
tissuePET TTAC (not modified, just to get frame times).
bfPlace for basis functions (initiated DFT struct, allocated and filled here).
bfNrNr of basis functions to calculate.
thetaminMinimum of theta=k2+k3 (sec-1 or min-1, corresponding to TAC time units).
thetamaxMaximum of theta=k2+k3 (sec-1 or min-1, corresponding to TAC time units).
statusPointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 217 of file bf_model.c.

235 {
236 int bi, fi, ret;
237 double a, b, c;
238
239 if(verbose>0)
240 printf("\n%s(*inp, *tis, *bf, %d, %g, %g, status, %d)\n", __func__,
241 bfNr, thetamin, thetamax, verbose);
242
243 /* Check the parameters */
244 if(input==NULL || tissue==NULL || bf==NULL) {
245 if(status!=NULL) strcpy(status, "program error");
246 else if(verbose>0) fprintf(stderr, "invalid function parameters\n");
247 return(1);
248 }
249 if(input->frameNr<3 || input->voiNr<1) {
250 if(status!=NULL) strcpy(status, "no input data");
251 else if(verbose>0) fprintf(stderr, "invalid input data\n");
252 return(2);
253 }
254 if(tissue->frameNr<1) {
255 if(status!=NULL) strcpy(status, "no pet data");
256 else if(verbose>0) fprintf(stderr, "invalid PET data\n");
257 return(3);
258 }
259 if(input->timeunit!=tissue->timeunit) {
260 if(status!=NULL) strcpy(status, "invalid time units");
261 else if(verbose>0) fprintf(stderr, "invalid time units\n");
262 return(4);
263 }
264 if(bfNr<2) {
265 if(status!=NULL) strcpy(status, "invalid nr of basis functions");
266 else if(verbose>0) fprintf(stderr, "invalid number of basis functions\n");
267 return(5);
268 }
269 if(thetamin<0.0) thetamin=0.0;
270 if(thetamin>=thetamax) {
271 if(status!=NULL) strcpy(status, "invalid theta range");
272 else if(verbose>0) fprintf(stderr, "invalid theta range\n");
273 return(6);
274 }
275 if(verbose>1) {
276 printf("input timerange: %g - %g\n", input->x[0], input->x[input->frameNr-1]);
277 printf("tissue timerange: %g - %g\n", tissue->x[0], tissue->x[tissue->frameNr-1]);
278 }
279
280 /* Allocate memory for basis functions */
281 if(verbose>1) printf("allocating memory for basis functions\n");
282 ret=dftSetmem(bf, tissue->frameNr, bfNr);
283 if(ret) {
284 if(status!=NULL) strcpy(status, "out of memory");
285 else if(verbose>0) fprintf(stderr, "out of memory\n");
286 return(10);
287 }
288
289 /* Copy and set information fields */
290 bf->voiNr=bfNr; bf->frameNr=tissue->frameNr;
291 bf->_type=tissue->_type;
292 dftCopymainhdr2(tissue, bf, 1);
293 for(bi=0; bi<bf->voiNr; bi++) {
294 snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
295 strcpy(bf->voi[bi].hemisphere, ".");
296 strcpy(bf->voi[bi].place, ".");
297 strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
298 }
299 for(fi=0; fi<bf->frameNr; fi++) {
300 bf->x[fi]=tissue->x[fi];
301 bf->x1[fi]=tissue->x1[fi];
302 bf->x2[fi]=tissue->x2[fi];
303 }
304
305 /* Compute the range of theta values to size fields */
306 if(verbose>1) printf("computing theta values\n");
307 a=thetamin; b=thetamax; c=(b-a)/(double)(bfNr-1);
308 if(verbose>20) printf("a=%g b=%g, c=%g\n", a, b, c);
309 for(bi=0; bi<bf->voiNr; bi++) bf->voi[bi].size=(double)bi*c+a;
310 if(verbose>2) {
311 printf("final BF theta range: %g - %g\n", bf->voi[0].size, bf->voi[bf->voiNr-1].size);
312 printf("theta step size: %g\n", c);
313 }
314
315 /* Allocate memory for simulated TAC */
316 double *sim;
317 sim=(double*)malloc(input->frameNr*sizeof(double));
318 if(sim==NULL) {
319 if(status!=NULL) strcpy(status, "out of memory");
320 else if(verbose>0) fprintf(stderr, "out of memory\n");
321 dftEmpty(bf); return(11);
322 }
323
324 /* Calculate the basis functions at input time points */
325 if(verbose>1) printf("computing basis functions at input sample times\n");
326 for(bi=0; bi<bf->voiNr; bi++) {
327 a=bf->voi[bi].size;
328 ret=simC1(input->x, input->voi[0].y, input->frameNr, 1.0, a, sim);
329 if(ret) {
330 if(status!=NULL) strcpy(status, "simulation problem");
331 else if(verbose>0) fprintf(stderr, "simulation problem\n");
332 free(sim); dftEmpty(bf); return(20);
333 }
334 if(verbose>100) {
335 printf("\ntheta := %g\n", a);
336 printf("simulated TAC:\n");
337 for(fi=0; fi<input->frameNr; fi++)
338 printf(" %12.6f %12.3f\n", input->x[fi], sim[fi]);
339 }
340 /* interpolate to PET time frames */
341 if(tissue->timetype==DFT_TIME_STARTEND)
342 ret=interpolate4pet(input->x, sim, input->frameNr, tissue->x1, tissue->x2,
343 bf->voi[bi].y, NULL, NULL, bf->frameNr);
344 else
345 ret=interpolate(input->x, sim, input->frameNr, tissue->x,
346 bf->voi[bi].y, NULL, NULL, bf->frameNr);
347 if(ret) {
348 if(status!=NULL) strcpy(status, "simulation problem");
349 else if(verbose>0) fprintf(stderr, "simulation problem\n");
350 free(sim); dftEmpty(bf); return(20);
351 }
352
353 } // next basis function
354
355 free(sim);
356 if(verbose>1) printf("%s() done.\n\n", __func__);
357 if(status!=NULL) strcpy(status, "ok");
358 return(0);
359}
dftCopymainhdr2
int dftCopymainhdr2(DFT *dft1, DFT *dft2, int ow)
Definition dft.c:587
dftEmpty
void dftEmpty(DFT *data)
Definition dft.c:20
interpolate
int interpolate(double *x, double *y, int nr, double *newx, double *newy, double *newyi, double *newyii, int newnr)
Linear interpolation and integration.
Definition integr.c:28
interpolate4pet
int interpolate4pet(double *x, double *y, int nr, double *newx1, double *newx2, double *newy, double *newyi, double *newyii, int newnr)
Interpolate and integrate TAC to PET frames.
Definition integr.c:510
DFT_TIME_STARTEND
#define DFT_TIME_STARTEND
Definition libtpccurveio.h:438
DFT::timetype
int timetype
Definition libtpccurveio.h:99
DFT::timeunit
int timeunit
Definition libtpccurveio.h:85
DFT::x1
double * x1
Definition libtpccurveio.h:104
DFT::x2
double * x2
Definition libtpccurveio.h:106

◆ bfRadiowater()

int bfRadiowater ( DFT * input,
DFT * tissue,
DFT * bf,
int bfNr,
double k2min,
double k2max,
char * status,
int verbose )

Calculates set of basis functions for generic radiowater model.

Returns
Returns 0 if successful, otherwise non-zero.
See also
bf_srtm, bfIrr2TCM
Parameters
inputArterial blood input TAC (not modified).
tissuePET TACs (not modified, just to get frame times).
bfPlace for basis functions (initiated DFT struct, allocated and filled here).
bfNrNr of basis functions to calculate.
k2minMinimum of k2 (sec-1 or min-1, corresponding to TAC time units).
k2maxMaximum of k2 (sec-1 or min-1, corresponding to TAC time units).
statusPointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed.
verboseVerbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 77 of file bf_model.c.

95 {
96 if(verbose>0)
97 printf("\n%s(*inp, *tis, *bf, %d, %g, %g, status, %d)\n", __func__,
98 bfNr, k2min, k2max, verbose);
99
100 /* Check the parameters */
101 if(input==NULL || tissue==NULL || bf==NULL) {
102 if(status!=NULL) strcpy(status, "program error");
103 return 1;
104 }
105 if(input->frameNr<3 || input->voiNr<1) {
106 if(status!=NULL) strcpy(status, "no input data");
107 return 2;
108 }
109 if(tissue->frameNr<1) {
110 if(status!=NULL) strcpy(status, "no pet data");
111 return 3;
112 }
113 if(input->timeunit!=tissue->timeunit) {
114 if(status!=NULL) strcpy(status, "invalid time units");
115 return 4;
116 }
117 if(bfNr<2) {
118 if(status!=NULL) strcpy(status, "invalid nr of basis functions");
119 return 5;
120 }
121 if(k2min<1.0E-10) k2min=1.0E-10; // range calculation does not work otherwise
122 if(k2min>=k2max || k2min<0.0) {
123 if(status!=NULL) strcpy(status, "invalid k2 range");
124 return 6;
125 }
126 if(verbose>1) {
127 printf("input timerange: %g - %g\n", input->x[0], input->x[input->frameNr-1]);
128 printf("tissue timerange: %g - %g\n", tissue->x[0], tissue->x[tissue->frameNr-1]);
129 }
130
131 /* Allocate memory for basis functions */
132 if(verbose>1) printf("allocating memory for basis functions\n");
133 if(dftSetmem(bf, tissue->frameNr, bfNr)) {
134 if(status!=NULL) strcpy(status, "out of memory");
135 return 10;
136 }
137
138 /* Copy and set information fields */
139 bf->voiNr=bfNr; bf->frameNr=tissue->frameNr;
140 bf->_type=tissue->_type;
141 dftCopymainhdr2(tissue, bf, 1);
142 for(int bi=0; bi<bf->voiNr; bi++) {
143 snprintf(bf->voi[bi].voiname, 6, "B%5.5d", bi+1);
144 strcpy(bf->voi[bi].hemisphere, ".");
145 strcpy(bf->voi[bi].place, ".");
146 strcpy(bf->voi[bi].name, bf->voi[bi].voiname);
147 }
148 for(int fi=0; fi<bf->frameNr; fi++) {
149 bf->x[fi]=tissue->x[fi];
150 bf->x1[fi]=tissue->x1[fi];
151 bf->x2[fi]=tissue->x2[fi];
152 }
153
154 /* Compute the range of k2 values to size fields */
155 if(verbose>1) printf("computing k2 values\n");
156 double a, b, c;
157 a=log10(k2min); b=log10(k2max); c=(b-a)/(double)(bfNr-1);
158 if(verbose>20) printf("a=%g b=%g, c=%g\n", a, b, c);
159 for(int bi=0; bi<bf->voiNr; bi++) {
160 bf->voi[bi].size=pow(10.0, (double)bi*c+a);
161 }
162 if(verbose>2) {
163 printf("final BF k2 range: %g - %g\n", bf->voi[0].size, bf->voi[bf->voiNr-1].size);
164 }
165
166 /* Allocate memory for simulated TAC */
167 double *sim;
168 sim=(double*)malloc(input->frameNr*sizeof(double));
169 if(sim==NULL) {
170 if(status!=NULL) strcpy(status, "out of memory");
171 dftEmpty(bf); return 11;
172 }
173
174 /* Calculate the basis functions at input time points */
175 if(verbose>1) printf("computing basis functions at input sample times\n");
176 for(int bi=0; bi<bf->voiNr; bi++) {
177 a=bf->voi[bi].size;
178 if(simC1(input->x, input->voi[0].y, input->frameNr, 1.0, a, sim)) {
179 if(status!=NULL) strcpy(status, "simulation problem");
180 free(sim); dftEmpty(bf);
181 return(20);
182 }
183 if(verbose>100) {
184 printf("\nk2 := %g\n", a);
185 printf("simulated TAC:\n");
186 for(int fi=0; fi<input->frameNr; fi++)
187 printf(" %12.6f %12.3f\n", input->x[fi], sim[fi]);
188 }
189 /* interpolate to PET time frames */
190 int ret=0;
191 if(tissue->timetype==DFT_TIME_STARTEND)
192 ret=interpolate4pet(input->x, sim, input->frameNr, tissue->x1, tissue->x2,
193 bf->voi[bi].y, NULL, NULL, bf->frameNr);
194 else
195 ret=interpolate(input->x, sim, input->frameNr, tissue->x,
196 bf->voi[bi].y, NULL, NULL, bf->frameNr);
197 if(ret) {
198 if(status!=NULL) strcpy(status, "simulation problem");
199 free(sim); dftEmpty(bf);
200 return(20);
201 }
202
203 } // next basis function
204
205 free(sim);
206 if(verbose>1) printf("%s() done.\n\n", __func__);
207 if(status!=NULL) strcpy(status, "ok");
208 return(0);
209}
Generated on Fri Sep 19 2025 13:32:07 for TPCCLIB by doxygen 1.12.0