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

Functions for simulating image of vessel (tube). More...

#include "libtpcidi.h"

Go to the source code of this file.

Functions

int idiSimulateTubeVol (VOL *vol, int zi, double cx, double cy, double r, double FWHM, double cbkg, double cblo)
int idiSimulateTubeImg (IMG *img, int zi, double cx, double cy, double r, double *cbkg, double *cblo)
int idiSimulateTubeImgPlane (int simmet, IMG *img, int zi, double cx, double cy, double r, double FWHM, double *cbkg, double *cblo)

Detailed Description

Functions for simulating image of vessel (tube).

Author
Vesa Oikonen

Definition in file vessel.c.

Function Documentation

◆ idiSimulateTubeImg()

int idiSimulateTubeImg ( IMG * img,
int zi,
double cx,
double cy,
double r,
double * cbkg,
double * cblo )

Simulate dynamic image surrounding a circular object extending across image planes (vessel).

Postcondition
Apply 2D Gaussian smoothing.
See also
imgGaussianFIRFilter, idiSimulateTubeVol, idiSimulateTubeImgPlane, imgSimulateSphere
Returns
Returns 0 if successful, or <> 0 in case of an error.
Parameters
imgPointer to allocated dynamic image; volume must contain pixel sizes and dimensions, and the same time frames as the TAC data.
ziPlane index [0..dimz-1].
cxX distance of circle centre (mm) from the upper left corner of the image.
cyY distance of circle centre (mm) from the upper left corner of the image.
rRadius of vessel (mm).
cbkgArray of background activities; size must be the same as img frame number.
cbloArray of vessel activities; can be higher or lower than background; size must be the same as img frame number.

Definition at line 81 of file vessel.c.

98 {
99 double dx1, dy1, dx2, dy2, d1, d2, d3, d4, r2;
100 int xi, yi, fi, n;
101
102 if(img->status!=IMG_STATUS_OCCUPIED) return(1);
103 if(zi<0 || zi>=img->dimz) return(2);
104 if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
105
106 /* Simulate the circle */
107 r2=r*r;
108 for(yi=0; yi<img->dimy; yi++) {
109 dy1=dy2=(0.5+(double)yi)*img->sizey - cy;
110 dy1-=0.25*img->sizey; dy2+=0.25*img->sizey;
111 for(xi=0; xi<img->dimx; xi++) {
112 dx1=dx2=(0.5+(double)xi)*img->sizex - cx;
113 dx1-=0.25*img->sizex; dx2+=0.25*img->sizex;
114 d1=dx1*dx1+dy1*dy1; d2=dx1*dx1+dy2*dy2;
115 d3=dx2*dx2+dy1*dy1; d4=dx2*dx2+dy2*dy2;
116 n=0; if(d1<=r2) n++; if(d2<=r2) n++; if(d3<=r2) n++; if(d4<=r2) n++;
117 for(fi=0; fi<img->dimt; fi++)
118 img->m[zi][yi][xi][fi]=0.25*((double)n*cblo[fi]+(double)(4-n)*cbkg[fi]);
119 }
120 }
121
122 return(0);
123}
IMG_STATUS_OCCUPIED
#define IMG_STATUS_OCCUPIED
Definition libtpcimgio.h:1479
IMG::sizex
float sizex
Definition libtpcimgio.h:1714
IMG::dimx
unsigned short int dimx
Definition libtpcimgio.h:1780
IMG::m
float **** m
Definition libtpcimgio.h:1799
IMG::status
char status
Definition libtpcimgio.h:1665
IMG::dimt
unsigned short int dimt
Definition libtpcimgio.h:1778
IMG::sizey
float sizey
Definition libtpcimgio.h:1716
IMG::dimz
unsigned short int dimz
Definition libtpcimgio.h:1784
IMG::dimy
unsigned short int dimy
Definition libtpcimgio.h:1782

◆ idiSimulateTubeImgPlane()

int idiSimulateTubeImgPlane ( int simmet,
IMG * img,
int zi,
double cx,
double cy,
double r,
double FWHM,
double * cbkg,
double * cblo )

Simulate image surrounding a circular object extending across image planes (bar or vessel). This function applies currently two methods to simulate the spill-over and spill-in effects:

Method 1. is based on Germano et al. JNM 1992; 33: 613-620. This equation can not be used to estimate recovery coefficient or the true activity concentration inside the vessel, but only to fit the radius of the vessel (see Germano), which then can be used to calculate the RC.

Method 2. simulates just the vessel without PVE; add 2D Gaussian smoothing later.

To simulate the circular vessel correctly in 2D image matrix use the equations in Brix et al. Nuklearmedizin 2002;41:184-190 instead; however, that would require numerical solution to double integrals, which may be either slow or error-prone for fitting purposes.

See also
idiSimulateTubeVol, idiSimulateTubeImg, imgGaussianFIRFilter, imgCircleMask
Returns
Returns 0 if successful, or <> 0 in case of an error.
Parameters
simmetSimulation method: 0=Germano, -1 or 1=No PVE.
imgPointer to allocated dynamic image; volume must contain pixel sizes and dimensions, and the same time frames as the TAC data.
ziPlane index [0..dimz-1].
cxX distance of circle centre (mm) from the upper left corner of the image.
cyY distance of circle centre (mm) from the upper left corner of the image.
rRadius of vessel (mm).
FWHMFWHM (mm).
cbkgArray of background activities; size must be the same as img frame number.
cbloArray of vessel activities; can be higher or lower than background; size must be the same as img frame number.

Definition at line 143 of file vessel.c.

164 {
165 int xi, yi, fi;
166
167 if(img->status!=IMG_STATUS_OCCUPIED) return(1);
168 if(zi<0 || zi>=img->dimz) return(2);
169 if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
170 if(cbkg==NULL || cblo==NULL) return(4);
171
172 double d, s=FWHM/2.354820; // 2*sqrt(2*ln(2))
173
174 if(simmet==0) { // Germano
175 double v, w1, w2, w3, w4, dy, dx;
176 for(yi=0; yi<img->dimy; yi++) {
177 dy=(0.5+(double)yi)*img->sizey - cy;
178 for(xi=0; xi<img->dimx; xi++) {
179 dx=(0.5+(double)xi)*img->sizex - cx;
180 d=hypot(dx-0.25*img->sizex, dy-0.25*img->sizey);
181 w1=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
182 d=hypot(dx-0.25*img->sizex, dy+0.25*img->sizey);
183 w2=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
184 d=hypot(dx+0.25*img->sizex, dy-0.25*img->sizey);
185 w3=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
186 d=hypot(dx+0.25*img->sizex, dy+0.25*img->sizey);
187 w4=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
188 for(fi=0; fi<img->dimt; fi++) {
189 v=cblo[fi]-cbkg[fi];
190 img->m[zi][yi][xi][fi]= 0.25*v*(w1+w2+w3+w4) + cbkg[fi];
191 }
192 }
193 }
194 } else { // Vessel without PVE
195 double dx1, dy1, dx2, dy2, d1, d2, d3, d4, r2;
196 int n;
197 /* Simulate the circle */
198 r2=r*r;
199 for(yi=0; yi<img->dimy; yi++) {
200 dy1=dy2=(0.5+(double)yi)*img->sizey - cy;
201 dy1-=0.25*img->sizey; dy2+=0.25*img->sizey;
202 for(xi=0; xi<img->dimx; xi++) {
203 dx1=dx2=(0.5+(double)xi)*img->sizex - cx;
204 dx1-=0.25*img->sizex; dx2+=0.25*img->sizex;
205 d1=dx1*dx1+dy1*dy1; d2=dx1*dx1+dy2*dy2;
206 d3=dx2*dx2+dy1*dy1; d4=dx2*dx2+dy2*dy2;
207 n=0; if(d1<=r2) n++; if(d2<=r2) n++; if(d3<=r2) n++; if(d4<=r2) n++;
208 for(fi=0; fi<img->dimt; fi++)
209 img->m[zi][yi][xi][fi]=0.25*((double)n*cblo[fi]+(double)(4-n)*cbkg[fi]);
210 }
211 }
212 }
213
214 return(0);
215}

◆ idiSimulateTubeVol()

int idiSimulateTubeVol ( VOL * vol,
int zi,
double cx,
double cy,
double r,
double FWHM,
double cbkg,
double cblo )

Simulate image volume surrounding a circular object extending across image planes (bar or vessel). Based on Germano et al. JNM 1992; 33: 613-620.

This equation can not be used to estimate recovery coefficient or the true activity concentration inside the vessel, but only to fit the radius of the vessel (see Germano), which then can be used to calculate the RC.

To simulate the circular vessel correctly in 2D image matrix use the equations in Brix et al. Nuklearmedizin 2002;41:184-190 instead; however, that would require numerical solution to double integrals, which may be either slow or error-prone for fitting purposes.

In this version, the activity is calculated as an average of four samples inside the pixel.

See also
idiSimulateTubeImg, idiSimulateTubeImgPlane, imgGaussianFIRFilter, imgCircleMask
Returns
Returns 0 if successful, or <> 0 in case of an error.
Parameters
volPointer to allocated image volume; volume must contain pixel sizes and dimensions
ziPlane index [0..dimz-1].
cxX distance of circle centre (mm) from the upper left corner of the image.
cyY distance of circle centre (mm) from the upper left corner of the image.
rRadius of vessel (mm).
FWHMFWHM (mm).
cbkgBackground activity.
cbloVessel activity; can be higher or lower than background.

Definition at line 26 of file vessel.c.

43 {
44 double s, dx, dy, d, v;
45 int xi, yi;
46
47 if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
48 if(zi<0 || zi>=vol->dimz) return(2);
49 if(vol->sizey<=0.0 || vol->sizex<=0.0) return(3);
50
51 s=FWHM/2.354820; // 2*sqrt(2*ln(2))
52 v=cblo-cbkg;
53 for(yi=0; yi<vol->dimy; yi++) {
54 dy=(0.5+(double)yi)*vol->sizey - cy;
55 for(xi=0; xi<vol->dimx; xi++) {
56 dx=(0.5+(double)xi)*vol->sizex - cx;
57 vol->v[zi][yi][xi]=0.0;
58 d=hypot(dx-0.25*vol->sizex, dy-0.25*vol->sizey);
59 vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
60 d=hypot(dx-0.25*vol->sizex, dy+0.25*vol->sizey);
61 vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
62 d=hypot(dx+0.25*vol->sizex, dy-0.25*vol->sizey);
63 vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
64 d=hypot(dx+0.25*vol->sizex, dy+0.25*vol->sizey);
65 vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
66 vol->v[zi][yi][xi]*=0.25;
67 vol->v[zi][yi][xi]+= cbkg;
68 }
69 }
70
71 return(0);
72}
VOL::status
char status
Definition libtpcimgio.h:1906
VOL::dimx
unsigned short int dimx
Definition libtpcimgio.h:1918
VOL::dimy
unsigned short int dimy
Definition libtpcimgio.h:1920
VOL::dimz
unsigned short int dimz
Definition libtpcimgio.h:1922
VOL::sizex
float sizex
Definition libtpcimgio.h:1912
VOL::v
float *** v
Definition libtpcimgio.h:1934
VOL::sizey
float sizey
Definition libtpcimgio.h:1914
Generated on for TPCCLIB by doxygen 1.15.0