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

Functions for computing K1 pixel-by-pixel. More...

#include "libtpcmodext.h"

Go to the source code of this file.

Functions

int img_k1_using_ki (DFT *input, IMG *dyn_img, int frame_nr, IMG *ki_img, IMG *k1_img, IMG *k2k3_img, char *status, int verbose)
 

Detailed Description

Functions for computing K1 pixel-by-pixel.

Author
Vesa Oikonen
Todo
Needs to be rewritten, including processing of units.

Definition in file img_k1.c.

Function Documentation

◆ img_k1_using_ki()

int img_k1_using_ki ( DFT * input,
IMG * dyn_img,
int frame_nr,
IMG * ki_img,
IMG * k1_img,
IMG * k2k3_img,
char * status,
int verbose )

Computing pixel-by-pixel the K1 for irreversible PET tracers using previously determined Ki (K1*k3/(k2+k3) and bilinear regression.

Returns
Returns 0 if succesful, and >1 in case of an error.
Parameters
inputPointer to the TAC data to be used as model input. Sample times in minutes. Curve is interpolated to PET frame times, if necessary
dyn_imgPointer to dynamic PET image data. Image and input data must be in the same calibration units
frame_nrNr of frames that will be included in the fit [3-frame_nr]
ki_imgPointer to previously calculated Ki image
k1_imgPointer to initiated IMG structure where K1 values will be placed
k2k3_imgPointer to initiated IMG structure where (k2+k3) values will be placed; enter NULL, if not needed
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 16 of file img_k1.c.

37 {
38 int zi, yi, xi, fi, ret=0;
39 int nnls_n=2, nnls_m, m;
40 DFT tac;
41 clock_t fitStart, fitFinish;
42
43
44 if(verbose>0) {
45 printf("%s(input, dyn_img, %d, ki_img, k1_img, ", __func__, frame_nr);
46 if(k2k3_img==NULL) printf("NULL, "); else printf("k2k3_img, ");
47 if(status==NULL) printf("NULL)\n"); else printf("status)\n");
48 }
49 /* Initial check for the arguments */
50 if(status!=NULL) sprintf(status, "invalid data");
51 if(dyn_img->status!=IMG_STATUS_OCCUPIED || dyn_img->dimt<1) return(1);
52 if(ki_img->status!=IMG_STATUS_OCCUPIED || ki_img->dimt<1) return(1);
53 if(input==NULL || input->frameNr<1) return(1);
54 if(frame_nr>dyn_img->dimt) return(1);
55 if(k1_img==NULL) return(1);
56 if(frame_nr<3) {
57 if(status!=NULL) sprintf(status, "invalid fit time");
58 return(1);
59 }
60
61 /* Check that input contains samples until at least 80% of line fit range */
62 if(verbose>1) {
63 printf("input_last_sample_time := %g\n", input->x[input->frameNr-1]);
64 printf("k1_start_time := %g\n", dyn_img->start[0]/60.0);
65 printf("k1_end_time := %g\n", dyn_img->end[frame_nr-1]/60.0);
66 }
67 if(input->x[input->frameNr-1] <
68 (0.2*dyn_img->mid[0]+0.8*dyn_img->mid[frame_nr-1])/60.0)
69 {
70 sprintf(status, "too few input samples"); return(2);
71 }
72
73 fitStart=clock();
74
75 /* Allocate memory for interpolated and integrated input and
76 tissue pixel TACs */
77 dftInit(&tac); if(dftSetmem(&tac, frame_nr, 2)!=0) {
78 sprintf(status, "out of memory"); return(3);
79 }
80 strcpy(tac.voi[0].voiname, "input"); strcpy(tac.voi[1].voiname, "tissue");
81 tac.voiNr=2; tac.frameNr=frame_nr;
82 for(fi=0; fi<tac.frameNr; fi++) {
83 tac.x1[fi]=dyn_img->start[fi]/60.;
84 tac.x2[fi]=dyn_img->end[fi]/60.;
85 tac.x[fi]=dyn_img->mid[fi]/60.;
86 }
87
88 /* If input sample times are much different from PET frame times,
89 as they usually are if arterial plasma is used as input, then use
90 interpolate4pet(), otherwise with image-derived input use petintegral()
91 */
92 ret=0;
93 if(input->frameNr<dyn_img->dimt) ret=1;
94 for(fi=0; fi<tac.frameNr; fi++) {
95 if(verbose>8) printf(" %g vs %g\n", input->x1[fi], tac.x1[fi]);
96 if(input->x1[fi]>tac.x1[fi]+0.034 || input->x1[fi]<tac.x1[fi]-0.034) {
97 ret++; continue;}
98 if(input->x2[fi]>tac.x2[fi]+0.034 || input->x2[fi]<tac.x2[fi]-0.034) ret++;
99 }
100 if(ret>0) {
101 if(verbose>1) printf("using interpolate4pet() for input curve\n");
102 ret=interpolate4pet(
103 input->x, input->voi[0].y, input->frameNr,
104 tac.x1, tac.x2, tac.voi[0].y, tac.voi[0].y2, tac.voi[0].y3, tac.frameNr
105 );
106 } else {
107 if(verbose>1) printf("copying input curve and using petintegral()\n");
108 /* because times are the same, the values can be copied directly */
109 for(fi=0; fi<tac.frameNr; fi++) tac.voi[0].y[fi]=input->voi[0].y[fi];
110 ret=petintegral(tac.x1, tac.x2, tac.voi[0].y, tac.frameNr,
111 tac.voi[0].y2, tac.voi[0].y3);
112 }
113 if(ret) {
114 dftEmpty(&tac); sprintf(status, "cannot interpolate input data");
115 if(verbose>0) printf(" ret := %d\n", ret);
116 return(4);
117 }
118 if(verbose>3) dftPrint(&tac);
119
120 /*
121 * Allocate result images and fill the header info
122 */
123 /* K1 image */
124 imgEmpty(k1_img);
125 ret=imgAllocate(k1_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1);
126 if(ret) {
127 sprintf(status, "cannot allocate memory for K1 image");
128 imgEmpty(k1_img); dftEmpty(&tac); return(5);
129 }
130 ret=imgCopyhdr(dyn_img, k1_img);
131 if(ret) {
132 sprintf(status, "cannot copy header info for result image");
133 imgEmpty(k1_img); dftEmpty(&tac); return(5);
134 }
135 k1_img->unit=CUNIT_ML_PER_ML_PER_MIN;
136 k1_img->decayCorrection=IMG_DC_NONCORRECTED; k1_img->isWeight=0;
137 k1_img->start[0]=dyn_img->start[0]; k1_img->end[0]=dyn_img->end[frame_nr-1];
138 if(verbose>9) imgInfo(k1_img);
139
140 /* (k2+k3) image */
141 if(k2k3_img!=NULL) {
142 imgEmpty(k2k3_img);
143 ret=imgAllocate(k2k3_img, dyn_img->dimz, dyn_img->dimy, dyn_img->dimx, 1);
144 if(ret) {
145 sprintf(status, "cannot allocate memory for (k2+k3) image");
146 imgEmpty(k2k3_img); imgEmpty(k1_img); dftEmpty(&tac); return(6);
147 }
148 ret=imgCopyhdr(dyn_img, k2k3_img);
149 if(ret) {
150 sprintf(status, "cannot copy header info for result image");
151 imgEmpty(k2k3_img); imgEmpty(k1_img); dftEmpty(&tac); return(6);
152 }
153 k2k3_img->unit=CUNIT_PER_MIN;
154 k2k3_img->decayCorrection=IMG_DC_NONCORRECTED; k2k3_img->isWeight=0;
155 k2k3_img->start[0]=dyn_img->start[0];
156 k2k3_img->end[0]=dyn_img->end[frame_nr-1];
157 if(verbose>9) imgInfo(k2k3_img);
158 }
159
160 /*
161 * Allocate memory required by NNLS; C99 !!!
162 */
163 if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
164 nnls_m=frame_nr;
165 double *nnls_a[nnls_n], nnls_mat[nnls_n][nnls_m], nnls_b[nnls_m],
166 nnls_zz[nnls_m];
167 double nnls_x[nnls_n], nnls_wp[nnls_n], nnls_rnorm;
168 int nnls_index[nnls_n];
169 nnls_a[0]=nnls_mat[0];
170 nnls_a[1]=nnls_mat[1];
171
172 /*
173 * Compute pixel-by-pixel
174 */
175 if(verbose>1) printf("computing K1 pixel-by-pixel\n");
176 for(zi=0; zi<dyn_img->dimz; zi++) {
177 for(yi=0; yi<dyn_img->dimy; yi++) {
178 for(xi=0; xi<dyn_img->dimx; xi++) {
179 /* Initiate pixel output values */
180 k1_img->m[zi][yi][xi][0]=0.0;
181 if(k2k3_img!=NULL) k2k3_img->m[zi][yi][xi][0]=0.0;
182
183 /* Copy and integrate pixel curve */
184 for(m=0; m<nnls_m; m++) tac.voi[1].y[m]=dyn_img->m[zi][yi][xi][m];
185 ret=petintegral(tac.x1, tac.x2, tac.voi[1].y, tac.frameNr,
186 tac.voi[1].y2, NULL);
187 if(ret) continue;
188 /* if AUC at the end is <= 0, then do nothing */
189 if(tac.voi[1].y2[nnls_m-1]<=0.0) continue;
190
191 /* Fill the NNLS data matrix */
192 for(m=0; m<nnls_m; m++) {
193 nnls_a[0][m]=tac.voi[0].y2[m];
194 nnls_a[1][m]=ki_img->m[zi][yi][xi][0]*tac.voi[0].y3[m]-tac.voi[1].y2[m];
195 nnls_b[m]=tac.voi[1].y[m];
196 }
197
198 /* NNLS */
199 ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm,
200 nnls_wp, nnls_zz, nnls_index);
201 if(ret>1) { /* no solution is possible */
202 continue;
203 } else if(ret==1) { /* max iteration count exceeded */ }
204 k1_img->m[zi][yi][xi][0]=nnls_x[0];
205 if(k2k3_img!=NULL) k2k3_img->m[zi][yi][xi][0]=nnls_x[1];
206
207 } /* next column */
208 } /* next row */
209 } /* next plane */
210 dftEmpty(&tac);
211
212 fitFinish=clock(); //printf("CLOCKS_PER_SEC=%ld\n", CLOCKS_PER_SEC);
213 //printf("%ld - %ld\n", fitFinish, fitStart);
214 if(verbose>0 && fitFinish!=(clock_t)(-1)) printf("done in %g seconds.\n",
215 (double)(fitFinish - fitStart) / (double)CLOCKS_PER_SEC );
216
217 return(0);
218}
dftInit
void dftInit(DFT *data)
Definition dft.c:38
dftSetmem
int dftSetmem(DFT *data, int frameNr, int voiNr)
Definition dft.c:57
dftEmpty
void dftEmpty(DFT *data)
Definition dft.c:20
dftPrint
void dftPrint(DFT *data)
Definition dftio.c:538
imgInfo
void imgInfo(IMG *image)
Definition img.c:359
imgAllocate
int imgAllocate(IMG *image, int planes, int rows, int columns, int frames)
Definition img.c:194
imgCopyhdr
int imgCopyhdr(IMG *image1, IMG *image2)
Definition img.c:471
imgEmpty
void imgEmpty(IMG *image)
Definition img.c:121
petintegral
int petintegral(double *x1, double *x2, double *y, int nr, double *ie, double *iie)
Integrate PET TAC data to frame mid times.
Definition integr.c:771
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
IMG_STATUS_OCCUPIED
#define IMG_STATUS_OCCUPIED
Definition libtpcimgio.h:1479
IMG_DC_NONCORRECTED
#define IMG_DC_NONCORRECTED
Definition libtpcimgio.h:1533
nnls
int nnls(double **a, int m, int n, double *b, double *x, double *rnorm, double *w, double *zz, int *index)
Definition nnls.c:38
DFT::voi
Voi * voi
Definition libtpccurveio.h:108
DFT::x1
double * x1
Definition libtpccurveio.h:104
DFT::voiNr
int voiNr
Definition libtpccurveio.h:78
DFT::x2
double * x2
Definition libtpccurveio.h:106
DFT::frameNr
int frameNr
Definition libtpccurveio.h:76
DFT::x
double * x
Definition libtpccurveio.h:102
IMG::dimx
unsigned short int dimx
Definition libtpcimgio.h:1780
IMG::m
float **** m
Definition libtpcimgio.h:1799
IMG::decayCorrection
char decayCorrection
Definition libtpcimgio.h:1688
IMG::unit
char unit
Definition libtpcimgio.h:1673
IMG::status
char status
Definition libtpcimgio.h:1665
IMG::dimt
unsigned short int dimt
Definition libtpcimgio.h:1778
IMG::start
float * start
Definition libtpcimgio.h:1815
IMG::dimz
unsigned short int dimz
Definition libtpcimgio.h:1784
IMG::dimy
unsigned short int dimy
Definition libtpcimgio.h:1782
IMG::end
float * end
Definition libtpcimgio.h:1817
IMG::mid
float * mid
Definition libtpcimgio.h:1819
IMG::isWeight
char isWeight
Definition libtpcimgio.h:1825
Voi::y2
double * y2
Definition libtpccurveio.h:61
Voi::voiname
char voiname[MAX_REGIONSUBNAME_LEN+1]
Definition libtpccurveio.h:50
Voi::y
double * y
Definition libtpccurveio.h:59
Voi::y3
double * y3
Definition libtpccurveio.h:63
Generated on Fri Sep 19 2025 13:32:07 for TPCCLIB by doxygen 1.12.0