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

BFM for single tissue compartment model. More...

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

Go to the source code of this file.

Functions

int bfm1TCM (TAC *input, TAC *tissue, int bfNr, const double k2min, const double k2max, const int distr, TAC *bf, TPCSTATUS *status)

Detailed Description

BFM for single tissue compartment model.

Definition in file bf_1tcm.c.

Function Documentation

◆ bfm1TCM()

int bfm1TCM ( TAC * input,
TAC * tissue,
int bfNr,
const double k2min,
const double k2max,
const int distr,
TAC * bf,
TPCSTATUS * status )

Calculate set of basis functions for generic radiowater model.

Todo
Add tests.
Returns
enum tpcerror (TPCERROR_OK when successful).
See also
bfmSRTM
Parameters
inputPointer to blood input TAC (not modified). If frame start and end times are available (->isframes!=0), then BTAC integral is calculated and used as input.
tissuePointer to PET TTAC data, which must contain the sample (frame) times.
bfNrNr of basis functions to calculate.
k2minMinimum of k2 (sec-1 or min-1, corresponding to TAC time units). If set to a negative value, |k2min| is used, but basis function with k2=0 is included.
k2maxMaximum of k2 (sec-1 or min-1, corresponding to TAC time units).
distrDistribution of k2 values: 0=log-based; 1=even.
bfPointer to output TAC, to be allocated and filled with basis functions in here.
statusPointer to status data; enter NULL if not needed.

Definition at line 27 of file bf_1tcm.c.

46 {
47 int verbose=0; if(status!=NULL) verbose=status->verbose;
48 if(verbose>0) printf("%s(itac, ttac, %d, %g, %g, %d, bf)\n", __func__, bfNr, k2min, k2max, distr);
49 if(input==NULL || tissue==NULL || bf==NULL) {
50 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
51 return TPCERROR_FAIL;
52 }
53 if(input->sampleNr<1 || input->tacNr<1 || tissue->sampleNr<1) {
54 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
55 return TPCERROR_NO_DATA;
56 }
57 if(input->sampleNr<3 || bfNr<2) {
58 if(verbose>1) printf("invalid sample or BF number\n");
59 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
60 return TPCERROR_TOO_FEW;
61 }
62 if(fabs(k2min)<1.0E-10 || fabs(k2min)>=k2max) { // range cannot be computed otherwise
63 if(verbose>1) printf("invalid k2 range\n");
64 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
65 return TPCERROR_INVALID_VALUE;
66 }
67
68 int ret;
69
70 /* Check units */
71 if(status==NULL || !status->forgiving) {
72 if(!unitIsTime(input->tunit) || input->tunit!=tissue->tunit) {
73 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_UNIT);
74 return TPCERROR_UNKNOWN_UNIT;
75 }
76 }
77
78 /* Allocate memory for basis functions */
79 if(verbose>1) printf("allocating memory for basis functions\n");
80 tacFree(bf);
81 ret=tacAllocate(bf, tissue->sampleNr, bfNr);
82 if(ret!=TPCERROR_OK) {
83 statusSet(status, __func__, __FILE__, __LINE__, ret);
84 return ret;
85 }
86 /* Copy and set information fields */
87 bf->tacNr=bfNr; bf->sampleNr=tissue->sampleNr;
88 bf->format=TAC_FORMAT_TSV_UK;
89 bf->isframe=tissue->isframe; tacXCopy(tissue, bf, 0, tissue->sampleNr-1);
90 bf->tunit=tissue->tunit; bf->cunit=tissue->cunit;
91 for(int bi=0; bi<bf->tacNr; bi++) sprintf(bf->c[bi].name, "B%5.5d", bi+1);
92 /* Compute the range of k2 values to size fields */
93 if(verbose>1) printf("computing k2 values\n");
94 int zero=0; if(k2min<0.0) {zero=1; bf->c[0].size=0.0;}
95 if(distr==0) { // printf(" log-based distribution\n");
96 double a, b, c;
97 a=log10(fabs(k2min)); b=log10(k2max); c=(b-a)/(double)(bfNr-1-zero);
98 for(int bi=zero; bi<bf->tacNr; bi++) {
99 bf->c[bi].size=pow(10.0, (double)(bi-zero)*c+a);
100 }
101 } else { // printf(" even distribution\n");
102 double c=(k2max-fabs(k2min))/(double)(bfNr-1-zero);
103 for(int bi=zero; bi<bf->tacNr; bi++) {
104 bf->c[bi].size=fabs(k2min)+(double)(bi-zero)*c;
105 }
106 }
107
108
109#if(0) // The replacement works better with image-derived input with long frames
110
111 /* Allocate memory for simulated TAC */
112 double *sim;
113 sim=(double*)malloc(input->sampleNr*sizeof(double));
114 if(sim==NULL) {
115 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
116 tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
117 }
118
119 /* Simulate the basis functions at input time points, and interpolate to tissue sample times */
120 double *cbi=NULL;
121 if(input->isframe!=0) {
122 // Time frames are available:
123 // Integrate BTAC to frame middle times and use it as input */
124 if(verbose>1) printf("using BTAC integral as simulation input\n");
125 cbi=(double*)malloc(input->sampleNr*sizeof(double));
126 if(cbi==NULL) {
127 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
128 tacFree(bf); free(sim); return TPCERROR_OUT_OF_MEMORY;
129 }
130 ret=liIntegratePET(input->x1, input->x2, input->c[0].y, input->sampleNr, cbi, NULL, verbose-10);
131 if(ret) {
132 if(verbose>1) printf("liIntegratePET() = %d\n", ret);
133 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
134 tacFree(bf); free(sim); free(cbi); return TPCERROR_INVALID_VALUE;
135 }
136 }
137 if(verbose>1) printf("computing basis functions\n");
138 ret=0;
139 for(int bi=0; bi<bf->tacNr && !ret; bi++) {
140 if(input->isframe==0) { // time frames not available
141 ret=simC1(input->x, input->c[0].y, input->sampleNr, 1.0, bf->c[bi].size, sim);
142 } else { // time frames are available
143 ret=simC1_i(input->x, cbi, input->sampleNr, 1.0, bf->c[bi].size, sim);
144 }
145 if(ret) break;
146 if(verbose>100) {
147 printf("\nk2 := %g\n", bf->c[bi].size);
148 printf("simulated TAC:\n");
149 for(int i=0; i<input->sampleNr; i++) printf(" %12.6f %12.3f\n", input->x[i], sim[i]);
150 }
151 /* interpolate to PET time frames */
152 if(tissue->isframe)
153 ret=liInterpolateForPET(input->x, sim, input->sampleNr,
154 tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
155 else
156 ret=liInterpolate(input->x, sim, input->sampleNr, tissue->x,
157 bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
158 if(ret) break;
159 } // next basis function
160 free(sim); free(cbi);
161 if(ret) {
162 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
163 return TPCERROR_INVALID_VALUE;
164 }
165
166#else
167
168 /* Interpolate the input function to times based on both input and tissue data */
169 /* When input TAC has long frames, the last simulated tissue frame will be biased without this. */
170 TAC inp; tacInit(&inp); TPCSTATUS status2; statusInit(&status2);
171 tacInput2sim(input, tissue, &inp, &status2);
172 statusFree(&status2);
173
174 /* Allocate memory for simulated TAC */
175 double *sim;
176 sim=(double*)malloc(inp.sampleNr*sizeof(double));
177 if(sim==NULL) {
178 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
179 tacFree(bf); return TPCERROR_OUT_OF_MEMORY;
180 }
181
182 if(verbose>1) printf("computing basis functions\n");
183 ret=0;
184 for(int bi=0; bi<bf->tacNr && !ret; bi++) {
185 ret=simC1(inp.x, inp.c[0].y, inp.sampleNr, 1.0, bf->c[bi].size, sim);
186 if(ret) break;
187 if(verbose>100) {
188 printf("\nk2 := %g\n", bf->c[bi].size);
189 printf("simulated TAC:\n");
190 for(int i=0; i<inp.sampleNr; i++) printf(" %12.6f %12.3f\n", inp.x[i], sim[i]);
191 }
192 /* interpolate to PET time frames */
193 if(tissue->isframe)
194 ret=liInterpolateForPET(inp.x, sim, inp.sampleNr,
195 tissue->x1, tissue->x2, bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
196 else
197 ret=liInterpolate(inp.x, sim, inp.sampleNr, tissue->x,
198 bf->c[bi].y, NULL, NULL, bf->sampleNr, 4, 1, verbose-8);
199 if(ret) break;
200 } // next basis function
201 free(sim);
202 if(ret) {
203 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
204 return TPCERROR_INVALID_VALUE;
205 }
206
207 tacFree(&inp);
208#endif
209
210
211 statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
212 return(TPCERROR_OK);
213}
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
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
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
tacInput2sim
int tacInput2sim(TAC *itac, TAC *ttac, TAC *stac, TPCSTATUS *status)
Modify input TAC based on tissue TAC, for use in simulations.
Definition lisim.c:29
simC1
int simC1(double *t, double *ca, const int nr, const double k1, const double k2, double *ct)
Definition sim1cm.c:93
simC1_i
int simC1_i(double *t, double *cai, const int nr, const double k1, const double k2, double *ct)
Definition sim1cm.c:164
statusInit
void statusInit(TPCSTATUS *s)
Definition statusmsg.c:104
statusSet
void statusSet(TPCSTATUS *s, const char *func, const char *srcfile, int srcline, tpcerror error)
Definition statusmsg.c:142
statusFree
void statusFree(TPCSTATUS *s)
Definition statusmsg.c:126
TACC::name
char name[MAX_TACNAME_LEN+1]
Definition tpctac.h:81
TACC::y
double * y
Definition tpctac.h:75
TACC::size
double size
Definition tpctac.h:71
TAC
Definition tpctac.h:87
TAC::x
double * x
Definition tpctac.h:97
TAC::cunit
unit cunit
Definition tpctac.h:105
TAC::format
tacformat format
Definition tpctac.h:93
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::forgiving
int forgiving
Force level, 0 for strict tests for data units etc.
Definition tpcextensions.h:244
tacFree
void tacFree(TAC *tac)
Definition tac.c:106
tacAllocate
int tacAllocate(TAC *tac, int sampleNr, int tacNr)
Definition tac.c:130
tacInit
void tacInit(TAC *tac)
Definition tac.c:24
tacXCopy
int tacXCopy(TAC *tac1, TAC *tac2, int i1, int i2)
Definition tacx.c:24
TPCERROR_INVALID_VALUE
@ TPCERROR_INVALID_VALUE
Invalid value.
Definition tpcextensions.h:220
TPCERROR_UNKNOWN_UNIT
@ TPCERROR_UNKNOWN_UNIT
Unknown data unit.
Definition tpcextensions.h:226
TPCERROR_FAIL
@ TPCERROR_FAIL
General error.
Definition tpcextensions.h:192
TPCERROR_OUT_OF_MEMORY
@ TPCERROR_OUT_OF_MEMORY
Cannot allocate memory.
Definition tpcextensions.h:193
TPCERROR_OK
@ TPCERROR_OK
No error.
Definition tpcextensions.h:191
TPCERROR_NO_DATA
@ TPCERROR_NO_DATA
File contains no data.
Definition tpcextensions.h:203
TPCERROR_TOO_FEW
@ TPCERROR_TOO_FEW
File contains too few samples.
Definition tpcextensions.h:208
unitIsTime
int unitIsTime(int u)
Definition units.c:359
TAC_FORMAT_TSV_UK
@ TAC_FORMAT_TSV_UK
UK TSV (point as decimal separator).
Definition tpctac.h:37
Generated on for TPCCLIB by doxygen 1.15.0