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

Iterative deconvolution method for image PVC. More...

Go to the source code of this file.

Functions

int imgPVCRRL (IMG *img1, double FWHMxy, double FWHMz, int reblur, double betaFactor, double alpha, int maxIterNr, double termLimit, IMG *img2, int verbose)
int imgPVCRVC (IMG *img1, double FWHMxy, double FWHMz, int reblur, double betaFactor, double alpha, int maxIterNr, double termLimit, IMG *img2, int verbose)

Detailed Description

Iterative deconvolution method for image PVC.

Copyright
(c) Turku PET Centre
Author
Vesa Oikonen

Definition in file imgidpvc.c.

Function Documentation

◆ imgPVCRRL()

int imgPVCRRL ( IMG * img1,
double FWHMxy,
double FWHMz,
int reblur,
double betaFactor,
double alpha,
int maxIterNr,
double termLimit,
IMG * img2,
int verbose )

Iterative deconvolution for partial volume correction applying Richardson-Lucy method.

This function processes image frames independently.

References:

  • Tohka J, Reilhac A. Deconvolution-based partial volume correction in Raclopride-PET and Monte Carlo comparison to MR-based method. Neuroimage 2008;39(4):1570-1584. https://doi.org/10.1016/j.neuroimage.2007.10.038
  • Boussion N et al. Incorporation of wavelet-based denoising in iterative deconvolution for partial volume correction in whole-body PET imaging. Eur J Nucl Med Mol Imaging 2009;36(7):1064-1075. https://doi.org/10.1007/s00259-009-1065-5
See also
imgPVCRVC, imgGaussianFIRFilter, imgRead
Returns
Returns 0 if ok.
Parameters
img1Pointer to the original 3D or 4D image; not changed.
FWHMxyFWHM in x,y dimensions in pixels.
FWHMzFWHM in z dimension in pixels.
reblurReblurring in each iteration.
betaFactorUpper limit for deconvolved image as a factor of image maximum value.
alphaConverging rate parameter (0
maxIterNrMaximum number of iterations; enter 0 to use the default.
termLimitLimit for termination (0<termLimit<1); more iterations with smaller limit.
img2Pointer to the corrected image.
Precondition
IMG structure must be initiated.
See also
imgInit, imgEmpty
Parameters
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 87 of file imgidpvc.c.

110 {
111 if(verbose>0) {
112 printf("%s(i, %g, %g, %d, %g, %g, %d, %g, ...)\n",
113 __func__, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit);
114 fflush(stdout);
115 }
116 if(img1==NULL || img2==NULL) return(1);
117 if(img1->dimz<1 || img1->dimy<8 || img1->dimx<8 || img1->dimt<1) return(2);
118 if(!(FWHMxy>0.0)) return(3);
119 if(img1->dimz<2) FWHMz=0.0;
120 if(!(betaFactor>0.0)) betaFactor=1.4;
121 if(maxIterNr<1) maxIterNr=25;
122 if(!(alpha>0.0)) alpha=1.6;
123 if(!(termLimit>0.0) || termLimit>=1) termLimit=0.005;
124
125 /* Convert FWHM into Gaussian SD */
126 float SDxy=FWHMxy/2.355;
127 float SDz=FWHMz/2.355;
128
129 int x, dimx=img1->dimx;
130 int y, dimy=img1->dimy;
131 int z, dimz=img1->dimz;
132 int t, dimt=img1->dimt;
133
134 /* Allocate output image and set headers */
135 {
136 int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, dimt, img1);
137 if(ret) return(10+ret);
138 }
139
140 /* Allocate image frame, for 'true' and correction matrices */
141 IMG tkm, cm; imgInit(&tkm); imgInit(&cm);
142 {
143 int ret=imgAllocateWithHeader(&tkm, dimz, dimy, dimx, 1, img1);
144 if(!ret) ret=imgAllocateWithHeader(&cm, dimz, dimy, dimx, 1, img1);
145 if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(5);}
146 }
147
148 /*
149 * One image frame at a time
150 */
151 for(t=0; t<dimt; t++) {
152 if(verbose>2 && dimt>1) {printf(" frame %d\n", 1+t); fflush(stdout);}
153
154 /* First, copy the original measured contents to "true" image matrix */
155 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
156 tkm.m[z][y][x][0]=img1->m[z][y][x][t];
157
158 /* Get image frame minimum and maximum for setting beta */
159 float omin, omax;
160 if(imgMinMax(&tkm, &omin, &omax)) {imgEmpty(&tkm); imgEmpty(&cm); return(6);}
161 if(!(omax>omin) || !(omax>0.0)) continue;
162 double beta=betaFactor*omax;
163 if(verbose>4) {
164 printf(" image_min := %g\n image_max := %g\n", omin, omax);
165 printf(" beta := %g\n", beta);
166 fflush(stdout);
167 }
168
169 /* Iterations */
170 int ret=0, iterNr=0;
171 for(iterNr=0; iterNr<maxIterNr; iterNr++) {
172 if(verbose>3) {printf(" iteration %d\n", 1+iterNr); fflush(stdout);}
173 /* Convolution of current estimate of true image with PSF (Gaussian filter) */
174 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
175 cm.m[z][y][x][0]=tkm.m[z][y][x][0];
176 ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
177 if(ret) {ret+=100; break;}
178 /* Divide the measured image matrix with it */
179 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
180 float f=img1->m[z][y][x][t]/cm.m[z][y][x][0];
181 if(!isfinite(f)) {
182 if(fabs(img1->m[z][y][x][t])<1.0E-06) f=0.0; else f=1.0;
183 }
184 cm.m[z][y][x][0]=f;
185 }
186 /* Re-blurring the correction matrix with PSF (Gaussian filter) */
187 if(reblur) {
188 ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
189 if(ret) {ret+=200; break;}
190 }
191 /* Multiply the previous estimate of true image matrix; */
192 /* and calculate factor to check for termination */
193 double dsum=0.0, psum=0.0;
194 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
195 double prev=tkm.m[z][y][x][0];
196 tkm.m[z][y][x][0]*=cm.m[z][y][x][0];
197 if(tkm.m[z][y][x][0]<=0.0 && img1->m[z][y][x][t]>0.0)
198 tkm.m[z][y][x][0]=0.01*img1->m[z][y][x][t];
199 if(tkm.m[z][y][x][0]>beta) tkm.m[z][y][x][0]=beta;
200 double d=tkm.m[z][y][x][0]-prev;
201 dsum+=d*d; psum+=prev*prev;
202 }
203 double dssqr=sqrt(dsum)/sqrt(psum); if(verbose>4) printf(" dssqr := %g\n", dssqr);
204 if(dssqr<termLimit) break;
205 } // next iteration
206 if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(ret);}
207 if(verbose>3) {printf(" iterNr := %d\n", (iterNr>=maxIterNr?maxIterNr:1+iterNr)); fflush(stdout);}
208
209 /* Copy the estimate of true matrix to output image */
210 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
211 img2->m[z][y][x][t]=tkm.m[z][y][x][0];
212
213 } // next frame
214
215 imgEmpty(&tkm); imgEmpty(&cm);
216 return(0);
217}
imgAllocateWithHeader
int imgAllocateWithHeader(IMG *image, int planes, int rows, int columns, int frames, IMG *image_from)
Definition img.c:279
imgEmpty
void imgEmpty(IMG *image)
Definition img.c:121
imgInit
void imgInit(IMG *image)
Definition img.c:60
imgGaussianFIRFilter
int imgGaussianFIRFilter(IMG *img, float xsd, float ysd, float zsd, double tolerance, int verbose)
Definition imgfilter.c:18
imgMinMax
int imgMinMax(IMG *img, float *minvalue, float *maxvalue)
Definition imgminmax.c:154
IMG::dimx
unsigned short int dimx
Definition libtpcimgio.h:1780
IMG::m
float **** m
Definition libtpcimgio.h:1799
IMG::dimt
unsigned short int dimt
Definition libtpcimgio.h:1778
IMG::dimz
unsigned short int dimz
Definition libtpcimgio.h:1784
IMG::dimy
unsigned short int dimy
Definition libtpcimgio.h:1782

◆ imgPVCRVC()

int imgPVCRVC ( IMG * img1,
double FWHMxy,
double FWHMz,
int reblur,
double betaFactor,
double alpha,
int maxIterNr,
double termLimit,
IMG * img2,
int verbose )

Iterative deconvolution for partial volume correction applying reblurred van Cittert method.

This function processes image frames independently.

References:

  • Jansson PA (ed.): Deconvolution of images and spectra, 2nd ed., Academic Press, 1996. ISBN: 0-12-380222-9.
  • Tohka J, Reilhac A. A Monte Carlo study of deconvolution algorithms for partial volume correction in quantitative PET. NSSMIC 2006. https://doi.org/10.1109/NSSMIC.2006.353719
See also
imgGaussianFIRFilter, imgRead
Returns
Returns 0 if ok.
Parameters
img1Pointer to the original 3D or 4D image; not changed.
FWHMxyFWHM in x,y dimensions in pixels.
FWHMzFWHM in z dimension in pixels.
reblurReblurring in each iteration.
betaFactorUpper limit for deconvolved image as a factor of image maximum value.
alphaConverging rate parameter (0
maxIterNrMaximum number of iterations; enter 0 to use the default.
termLimitLimit for termination (0<termLimit<1); more iterations with smaller limit.
img2Pointer to the corrected image.
Precondition
IMG structure must be initiated.
See also
imgInit, imgEmpty
Parameters
verboseVerbose level; if zero, then only warnings are printed into stderr.

Definition at line 234 of file imgidpvc.c.

257 {
258 if(verbose>0) {
259 printf("%s(i, %g, %g, %d, %g, %g, %d, %g, ...)\n",
260 __func__, FWHMxy, FWHMz, reblur, betaFactor, alpha, maxIterNr, termLimit);
261 fflush(stdout);
262 }
263 if(img1==NULL || img2==NULL) return(1);
264 if(img1->dimz<1 || img1->dimy<8 || img1->dimx<8 || img1->dimt<1) return(2);
265 if(!(FWHMxy>0.0)) return(3);
266 if(img1->dimz<2) FWHMz=0.0;
267 if(!(betaFactor>0.0)) betaFactor=1.4;
268 if(maxIterNr<1) maxIterNr=25;
269 if(!(alpha>0.0)) alpha=1.6;
270 if(!(termLimit>0.0) || termLimit>=1) termLimit=0.005;
271
272 /* Convert FWHM into Gaussian SD */
273 float SDxy=FWHMxy/2.355;
274 float SDz=FWHMz/2.355;
275
276 int x, dimx=img1->dimx;
277 int y, dimy=img1->dimy;
278 int z, dimz=img1->dimz;
279 int t, dimt=img1->dimt;
280
281 /* Allocate output image and set headers */
282 {
283 int ret=imgAllocateWithHeader(img2, dimz, dimy, dimx, dimt, img1);
284 if(ret) return(10+ret);
285 }
286
287 /* Allocate image frame, for 'true' and correction matrices */
288 IMG tkm, cm; imgInit(&tkm); imgInit(&cm);
289 {
290 int ret=imgAllocateWithHeader(&tkm, dimz, dimy, dimx, 1, img1);
291 if(!ret) ret=imgAllocateWithHeader(&cm, dimz, dimy, dimx, 1, img1);
292 if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(5);}
293 }
294
295 /*
296 * One image frame at a time
297 */
298 for(t=0; t<dimt; t++) {
299 if(verbose>2 && dimt>1) {printf(" frame %d\n", 1+t); fflush(stdout);}
300
301 /* First, copy the original measured contents to "true" image matrix */
302 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
303 tkm.m[z][y][x][0]=img1->m[z][y][x][t];
304
305 /* Get image frame minimum and maximum for setting beta */
306 float omin, omax;
307 if(imgMinMax(&tkm, &omin, &omax)) {imgEmpty(&tkm); imgEmpty(&cm); return(6);}
308 if(!(omax>omin) || !(omax>0.0)) continue;
309 double beta=betaFactor*omax;
310 if(verbose>4) {
311 printf(" image_min := %g\n image_max := %g\n", omin, omax);
312 printf(" beta := %g\n", beta);
313 fflush(stdout);
314 }
315
316 /* Iterations */
317 int ret=0, iterNr=0;
318 for(iterNr=0; iterNr<maxIterNr; iterNr++) {
319 if(verbose>3) {printf(" iteration %d\n", 1+iterNr); fflush(stdout);}
320 /* Convolution of current estimate of true image with PSF (Gaussian filter) */
321 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
322 cm.m[z][y][x][0]=tkm.m[z][y][x][0];
323 ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
324 if(ret) {ret+=100; break;}
325 /* Subtract from measured image matrix */
326 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
327 cm.m[z][y][x][0]=img1->m[z][y][x][t]-cm.m[z][y][x][0];
328 /* Re-blurring the difference with PSF (Gaussian filter) */
329 if(reblur) {
330 ret=imgGaussianFIRFilter(&cm, SDxy, SDxy, SDz, 1.0E-03, verbose-5);
331 if(ret) {ret+=200; break;}
332 }
333 /* Add to the previous estimate of true image matrix, applying relaxation function; */
334 /* and calculate factor to check for termination */
335 double dsum=0.0, psum=0.0;
336 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++) {
337 double r=alpha*(1.0 - (2.0/beta)*fabs(tkm.m[z][y][x][0] - beta/2.0));
338 double d=r*cm.m[z][y][x][0]; dsum+=d*d; psum+=tkm.m[z][y][x][0]*tkm.m[z][y][x][0];
339 tkm.m[z][y][x][0]+=r*cm.m[z][y][x][0];
340 }
341 double dssqr=sqrt(dsum)/sqrt(psum); if(verbose>4) printf(" dssqr := %g\n", dssqr);
342 if(dssqr<termLimit) break;
343 } // next iteration
344 if(ret) {imgEmpty(&tkm); imgEmpty(&cm); return(ret);}
345 if(verbose>3) {printf(" iterNr := %d\n", (iterNr>=maxIterNr?maxIterNr:1+iterNr)); fflush(stdout);}
346
347 /* Copy the estimate of true matrix to output image */
348 for(z=0; z<dimz; z++) for(y=0; y<dimy; y++) for(x=0; x<dimx; x++)
349 img2->m[z][y][x][t]=tkm.m[z][y][x][0];
350
351 } // next frame
352
353 imgEmpty(&tkm); imgEmpty(&cm);
354 return(0);
355}
Generated on for TPCCLIB by doxygen 1.15.0