libtpcidi.h File Reference

Header file for libtpcidi. More...

#include "tpcclibConfig.h"
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpcsvg.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"

Go to the source code of this file.

Functions
int	imgCircleMask (IMG *img, int zi, double cx, double cy, double r, double mv, double *smv, int verbose)
int	imgRingMask (IMG *img, int zi, double cx, double cy, double r1, double r2, double mv, double *smv, int verbose)
int	imgSimulateRing (IMG *img, int fi, int zi, double cx, double cy, double r1, double r2, double vr, double vi, double vo, int verbose)
int	simMyocDiameterCurve (DFT *dft, double t1, double t2, double hbr, double maxdiam, double mindiam)
int	imgSimulateSphere (IMG *img, int fi, double cx, double cy, double cz, double r1, double r2, double vr, double vi, double vo, int verbose)
int	heartRecoverySpilloverCorrectionFactors (double R, double d, double s, double Vb, double *FMM, double *FMB, double *FBM, double *FBB)
int	imgMaskPixelTACs (IMG *img, IMG *mask, double thrs, DFT *dft, int verbose)
int	imgGetConcWeightedPeakPos (IMG *img, float thrs, IMG_PIXEL *pos, int verbose)
double	rcPeakPET (double FWHM, double R)
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

Header file for libtpcidi.

Author

Vesa Oikonen

Definition in file libtpcidi.h.

Function Documentation

heartRecoverySpilloverCorrectionFactors()

int heartRecoverySpilloverCorrectionFactors ( double R, double d, double s, double Vb, double * FMM, double * FMB, double * FBM, double * FBB )

extern

Calculate recovery and spillover correction coefficients.

Based on Henze E, Huang S-C, Ratib O, Hoffman E, Phelps ME, Schelbert HR. Measurements of regional tissue and blood-pool radiotracer concentrations from serial tomographic images of the heart. J Nucl Med. 1983;24:987-996.

See also

rcPeakPET, simMyocDiameterCurve, imgSimulateSphere

Returns

Returns 0 if successful, or <> 0 in case of an error.

Parameters

R	radius of cavity and circular ROI (mm).
d	thickness of myocardium (mm).
s	spatial resolution (mm); s = FWHM/(2*SQRT(2*LN(2))).
Vb	Vascular volume fraction in myocardium; Henze et al assumed 0.1.
FMM	Pointer to resulting correction coefficient FMM, see the article.
FMB	Pointer to resulting correction coefficient FMB, see the article.
FBM	Pointer to resulting correction coefficient FBM, see the article.
FBB	Pointer to resulting correction coefficient FBB, see the article.

Definition at line 37 of file heartcorr.c.

  {
  double hp;
 
  if(R<=0.0 || d<=0.0 || s<=0.0) return(1);
  if(FMM==NULL || FMB==NULL || FBM==NULL || FBB==NULL) return(2);
  if(Vb<0.0 || Vb>=1.0) return(3);
  hp = exp(-R*R/(2.0*s*s));
  *FMB = hp - exp(-(R+d)*(R+d)/(2.0*s*s));
  *FMM = erf( (d/2.0) / (M_SQRT2*s) );
  *FBB = 1.0 - hp;
  *FBM = Vb + 0.5*(1.0 - *FMM);  //0.6 - 0.5 * *FMM;
  return(0);
}

idiSimulateTubeImg()

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

extern

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

img	Pointer to allocated dynamic image; volume must contain pixel sizes and dimensions, and the same time frames as the TAC data.
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r	Radius of vessel (mm).
cbkg	Array of background activities; size must be the same as img frame number.
cblo	Array 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.

  {
  double dx1, dy1, dx2, dy2, d1, d2, d3, d4, r2;
  int xi, yi, fi, n;
 
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(zi<0 || zi>=img->dimz) return(2);
  if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
 
  /* Simulate the circle */
  r2=r*r;
  for(yi=0; yi<img->dimy; yi++) {
    dy1=dy2=(0.5+(double)yi)*img->sizey - cy;
    dy1-=0.25*img->sizey; dy2+=0.25*img->sizey;
    for(xi=0; xi<img->dimx; xi++) {
      dx1=dx2=(0.5+(double)xi)*img->sizex - cx;
      dx1-=0.25*img->sizex; dx2+=0.25*img->sizex;
      d1=dx1*dx1+dy1*dy1; d2=dx1*dx1+dy2*dy2;
      d3=dx2*dx2+dy1*dy1; d4=dx2*dx2+dy2*dy2;
      n=0; if(d1<=r2) n++; if(d2<=r2) n++; if(d3<=r2) n++; if(d4<=r2) n++;
      for(fi=0; fi<img->dimt; fi++)
        img->m[zi][yi][xi][fi]=0.25*((double)n*cblo[fi]+(double)(4-n)*cbkg[fi]);
    }
  }
 
  return(0);
}

idiSimulateTubeImgPlane()

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

extern

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

simmet	Simulation method: 0=Germano, -1 or 1=No PVE.
img	Pointer to allocated dynamic image; volume must contain pixel sizes and dimensions, and the same time frames as the TAC data.
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r	Radius of vessel (mm).
FWHM	FWHM (mm).
cbkg	Array of background activities; size must be the same as img frame number.
cblo	Array 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.

  {
  int xi, yi, fi;
 
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(zi<0 || zi>=img->dimz) return(2);
  if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
  if(cbkg==NULL || cblo==NULL) return(4);
 
  double d, s=FWHM/2.354820; // 2*sqrt(2*ln(2))
 
  if(simmet==0) { // Germano
    double v, w1, w2, w3, w4, dy, dx;
    for(yi=0; yi<img->dimy; yi++) {
      dy=(0.5+(double)yi)*img->sizey - cy;
      for(xi=0; xi<img->dimx; xi++) {
        dx=(0.5+(double)xi)*img->sizex - cx;
        d=hypot(dx-0.25*img->sizex, dy-0.25*img->sizey);
        w1=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
        d=hypot(dx-0.25*img->sizex, dy+0.25*img->sizey);
        w2=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
        d=hypot(dx+0.25*img->sizex, dy-0.25*img->sizey);
        w3=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
        d=hypot(dx+0.25*img->sizex, dy+0.25*img->sizey);
        w4=0.5*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
        for(fi=0; fi<img->dimt; fi++) {
          v=cblo[fi]-cbkg[fi];
          img->m[zi][yi][xi][fi]= 0.25*v*(w1+w2+w3+w4) + cbkg[fi];
        }
      }
    }
  } else { // Vessel without PVE
    double dx1, dy1, dx2, dy2, d1, d2, d3, d4, r2;
    int n;
    /* Simulate the circle */
    r2=r*r;
    for(yi=0; yi<img->dimy; yi++) {
      dy1=dy2=(0.5+(double)yi)*img->sizey - cy;
      dy1-=0.25*img->sizey; dy2+=0.25*img->sizey;
      for(xi=0; xi<img->dimx; xi++) {
        dx1=dx2=(0.5+(double)xi)*img->sizex - cx;
        dx1-=0.25*img->sizex; dx2+=0.25*img->sizex;
        d1=dx1*dx1+dy1*dy1; d2=dx1*dx1+dy2*dy2;
        d3=dx2*dx2+dy1*dy1; d4=dx2*dx2+dy2*dy2;
        n=0; if(d1<=r2) n++; if(d2<=r2) n++; if(d3<=r2) n++; if(d4<=r2) n++;
        for(fi=0; fi<img->dimt; fi++)
          img->m[zi][yi][xi][fi]=0.25*((double)n*cblo[fi]+(double)(4-n)*cbkg[fi]);
      }
    }
  }
 
  return(0);
}

idiSimulateTubeVol()

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

extern

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

vol	Pointer to allocated image volume; volume must contain pixel sizes and dimensions
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r	Radius of vessel (mm).
FWHM	FWHM (mm).
cbkg	Background activity.
cblo	Vessel activity; can be higher or lower than background.

Definition at line 26 of file vessel.c.

  {
  double s, dx, dy, d, v;
  int xi, yi;
 
  if(vol->status!=IMG_STATUS_OCCUPIED) return(1);
  if(zi<0 || zi>=vol->dimz) return(2);
  if(vol->sizey<=0.0 || vol->sizex<=0.0) return(3);
 
  s=FWHM/2.354820; // 2*sqrt(2*ln(2))
  v=cblo-cbkg;
  for(yi=0; yi<vol->dimy; yi++) {
    dy=(0.5+(double)yi)*vol->sizey - cy;
    for(xi=0; xi<vol->dimx; xi++) {
      dx=(0.5+(double)xi)*vol->sizex - cx;
      vol->v[zi][yi][xi]=0.0;
      d=hypot(dx-0.25*vol->sizex, dy-0.25*vol->sizey);
      vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
      d=hypot(dx-0.25*vol->sizex, dy+0.25*vol->sizey);
      vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
      d=hypot(dx+0.25*vol->sizex, dy-0.25*vol->sizey);
      vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
      d=hypot(dx+0.25*vol->sizex, dy+0.25*vol->sizey);
      vol->v[zi][yi][xi]+=(v/2.0)*( erf((d+r)/(M_SQRT2*s)) - erf((d-r)/(M_SQRT2*s)) );
      vol->v[zi][yi][xi]*=0.25;
      vol->v[zi][yi][xi]+= cbkg;
    }
  }
 
  return(0);
}

imgCircleMask()

int imgCircleMask ( IMG * img, int zi, double cx, double cy, double r, double mv, double * smv, int verbose )

extern

Simulate a mask image of circle with specified radius.

The applied method is only approximate at pixel borders (pixel is divided into 5x5 subpixels).

See also

imgRingMask, imgSimulateSphere, idiSimulateTubeVol, imgMaskInvert

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated static or dynamic image; image must contain pixel sizes and dimensions; mask values are added to pixel values, thus you may need to set pixel values to zero before calling this function.
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r	Radius of circle (mm).
mv	Mask value; this value is added to each pixel value that fits inside the radius; the pixels that are partially inside the radius will get fraction of mask value.
smv	Pointer to value where the sum of added mask pixel values is written; enter NULL if not needed. Use this to validate the results.
verbose	Verbose level; set to <=0 to prevent all prints to stdout.

Definition at line 16 of file circle.c.

  {
  if(verbose>0) printf("%s(img, %d, %g, %g, %g, %g, ...)\n", __func__, zi, cx, cy, r, mv);
 
  int n, xi, yi, fi, i, j;
  double dx[5], dy[5], r2, v;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(zi<0 || zi>=img->dimz) return(2);
  if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
  if(img->sizey!=img->sizex) return(4);
  if(img->dimt<1) return(5);
 
  if(smv!=NULL) *smv=0.0;
  /* Compare r^2 to sum of squared distances instead of square root */
  r2=r*r;
  for(yi=0; yi<img->dimy; yi++) {
    dy[0]=(0.1+(double)yi)*img->sizey - cy;
    for(j=1; j<5; j++) dy[j]=dy[j-1]+0.2*img->sizey;
    for(xi=0; xi<img->dimx; xi++) {
      dx[0]=(0.1+(double)xi)*img->sizex - cx;
      for(i=1; i<5; i++) dx[i]=dx[i-1]+0.2*img->sizex;
      for(i=0, n=0; i<5; i++) for(j=0; j<5; j++) if((dx[i]*dx[i]+dy[j]*dy[j])<r2) n++;
      if(n==0) continue;
      v=(double)n*mv/25.0;
      for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]+=v;
      if(smv!=NULL) *smv+=v;
    }
  }
 
  return(0);
}

Referenced by imgRingMask().

imgGetConcWeightedPeakPos()

int imgGetConcWeightedPeakPos ( IMG * img, float thrs, IMG_PIXEL * pos, int verbose )

extern

Calculates concentration weighted peak position from a sum image.

See also

imgMaskPixelTACs, imgCircleMask, idiSimulateTubeVol

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated image from which the peak position is searched. If IMG struct contains more than one times frame, then only the first one is used.
thrs	Threshold fraction (0<thrs<1) between image min and max.
pos	Pointer to pixel position structure.
verbose	Verbose level; set to <=0 to prevent all prints to stdout and stderr.

Definition at line 14 of file peak.c.

  {
  if(verbose>0) printf("%s(img, %g, ...)\n", __func__, thrs);
 
  int xi, yi, zi, fi=0;
  float d, fmax, fmin, thrslev;
  float px, py, pz, fsum;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimz<1 || img->dimx<2 || img->dimy<2) return(2);
  if(img->dimt<1) return(3);
  if(thrs<=0.0 || thrs>=1.0) return(4);
  if(pos==NULL) return(5);
 
 
  /* Search min and max value inside the image volume */
  zi=yi=xi=0; fmax=fmin=img->m[zi][yi][xi][fi];
  for(zi=0; zi<img->dimz; zi++)
    for(yi=0; yi<img->dimy; yi++)
      for(xi=0; xi<img->dimx; xi++) {
        if(img->m[zi][yi][xi][fi] > fmax) fmax=img->m[zi][yi][xi][fi];
        else if(img->m[zi][yi][xi][fi] < fmin) fmin=img->m[zi][yi][xi][fi];
      }
  if(verbose>1) printf("fmin := %g\nfmax := %g\n", fmin, fmax);
 
  /* Calculate threshold level from min and max */
  d=fmax-fmin; if(d<=0.0) return(10);
  thrslev=d*thrs+fmin; if(verbose>1) printf("thrslev := %g\n", thrslev);
 
  /* Calculate concentration weighted peak position */
  px=py=pz=0.0; fsum=0.0;
  for(zi=0; zi<img->dimz; zi++)
    for(yi=0; yi<img->dimy; yi++)
      for(xi=0; xi<img->dimx; xi++) {
        d=img->m[zi][yi][xi][fi]-thrslev;
        if(d>0.0) {
          px+=(float)(xi+1)*d;
          py+=(float)(yi+1)*d;
          pz+=(float)(zi+1)*d;
          fsum+=d;
        }
      }
  if(fsum<=0.0) return 5;
  px/=fsum; py/=fsum; pz/=fsum;
  if(verbose>1) printf("peak_pos := (%g,%g,%g)\n", px, py, pz);
  pos->x=roundf(px); pos->y=roundf(py); pos->z=roundf(pz); pos->f=1;
 
  return(0);
}

imgMaskPixelTACs()

int imgMaskPixelTACs ( IMG * img, IMG * mask, double thrs, DFT * dft, int verbose )

extern

Extract TACs of every image voxel which has a value > thrs in the mask image.

Mask value is saved as region 'size' in DFT. Mainly for testing IDI methods.

See also

imgGetConcWeightedPeakPos, imgCircleMask, imgSimulateSphere

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated image from which the weighted TAC is calculated.
mask	Pointer to mask image; x, y, and z dimensions must be the same as in the image to which the mask is applied.
thrs	Mask threshold: pixels with mask values above this are included in TACs.
dft	Pointer to initiated DFT struct where pixel TACs will be written; any previous contents are deleted.
verbose	Verbose level; set to <=0 to prevent all prints to stdout .

Definition at line 17 of file idimask.c.

  {
  if(verbose>0) printf("%s()\n", __func__);
 
  int xi, yi, zi, fi, ri, ret;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(mask->status<IMG_STATUS_OCCUPIED) return(2);
  if(mask->dimz!=img->dimz) return(3);
  if(mask->dimy!=img->dimy) return(4);
  if(mask->dimx!=img->dimx) return(5);
  if(img->dimt<1 || mask->dimt<1) return(6);
  if(dft==NULL) return(7);
 
  /* Delete any previous TAC data; do it first so that caller can check
     the nr of extracted TACs inside DFT */
  dftEmpty(dft);
 
  /* Calculate the nr of voxels to extract */
  long long nr=0;
  for(zi=0; zi<mask->dimz; zi++)
    for(yi=0; yi<mask->dimy; yi++)
      for(xi=0; xi<mask->dimx; xi++)
        if(mask->m[zi][yi][xi][0]>thrs) nr++;
  if(verbose>1) printf("mask_pixel_nr := %lld\n", nr);
  /* Return, if no voxels were marked in mask image */
  if(nr==0) return(0);
 
  /* Allocate memory for pixel TACs */
  ret=dftAllocateWithIMG(dft, nr, img);
  if(ret!=0) {
    if(verbose>0) fprintf(stderr, "Error: cannot allocate memory for %lld pixel TACs.\n", nr);
    return(10);
  }
  strcpy(dft->studynr, img->studyNr);
 
  /* Save the pixel TACs in DFT */
  ri=0;
  for(zi=0; zi<mask->dimz; zi++)
    for(yi=0; yi<mask->dimy; yi++)
      for(xi=0; xi<mask->dimx; xi++) if(mask->m[zi][yi][xi][0]>thrs)
  {
    sprintf(dft->voi[ri].voiname, "%d", 1+xi);
    sprintf(dft->voi[ri].hemisphere, "%d", 1+yi);
    sprintf(dft->voi[ri].place, "%d", 1+zi);
    sprintf(dft->voi[ri].name, "%d %d %d", 1+xi, 1+yi, 1+zi);
    for(fi=0; fi<img->dimt; fi++) dft->voi[ri].y[fi]=img->m[zi][yi][xi][fi];
    dft->voi[ri].size=mask->m[zi][yi][xi][0];
    ri++;
  }
  dft->voiNr=ri;
 
  return(0);
}

imgRingMask()

int imgRingMask ( IMG * img, int zi, double cx, double cy, double r1, double r2, double mv, double * smv, int verbose )

extern

Simulate a mask image of ring with specified inner and outer radius.

The applied method is only approximate at pixel borders (pixel is divided into 5x5 subpixels).

See also

imgCircleMask, imgSimulateRing, imgSimulateSphere, imgMaskInvert

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated static or dynamic image; image must contain pixel sizes and dimensions; mask values are added to pixel values, thus you may need to set pixel values to zero before calling this function.
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r1	Inner radius of circle (mm).
r2	Outer radius of circle (mm).
mv	Mask value; this value is added to each pixel value that fits inside the radius; the pixels that are partially inside the radius will get fraction of mask value.
smv	Pointer to value where the sum of added mask pixel values is written; enter NULL if not needed. Use this to validate the results.
verbose	Verbose level; set to <=0 to prevent all prints to stdout.

Definition at line 77 of file circle.c.

  {
  int ret;
  double mv1, mv2, d;
 
  if(verbose>0) printf("%s(img, %d, %g, %g, %g, %g, %g, ...)\n", __func__, zi, cx, cy, r1, r2, mv);
  if(r1>=r2 || r1<0.0) {
    if(verbose>0) fprintf(stderr, "Error: invalid radius.\n");
    return(1);
  }
 
  /* Simulate circle with the larger radius */
  ret=imgCircleMask(img, zi, cx, cy, r2, mv, &mv1, verbose);
  if(ret!=0) {
    if(verbose>0) fprintf(stderr, "Error: cannot simulate circle 1.\n");
    return(ret);
  }
  /* Subtract circle with the smaller radius */
  ret=imgCircleMask(img, zi, cx, cy, r1, -mv, &mv2, verbose);
  if(ret!=0) {
    if(verbose>0) fprintf(stderr, "Error: cannot simulate circle 2.\n");
    return(ret);
  }
  /* Calculate the sum of mask pixel values */
  d=mv1+mv2; if(smv!=NULL) *smv=d;
  if(d==0.0 && verbose>0) {fprintf(stderr, "Warning: empty mask image created.\n"); fflush(stderr);}
  return(0);
}

imgSimulateRing()

int imgSimulateRing ( IMG * img, int fi, int zi, double cx, double cy, double r1, double r2, double vr, double vi, double vo, int verbose )

extern

Simulate an image of ring with specified inner and outer radius and activity in the ring and inside and outside of it. The applied method is only approximate at pixel borders (pixel is divided into 5x5 subpixels). PET resolution effects are not simulated.

See also

imgRingMask, imgSimulateSphere, imgGaussianFIRFilter, idiSimulateTubeVol

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated static or dynamic image; image must contain pixel sizes and dimensions; values are added to any existing pixel values, thus you may need to set pixel values to zero before calling this function.
fi	Frame index [0..dimt-1].
zi	Plane index [0..dimz-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
r1	Inner radius of circle (mm).
r2	Outer radius of circle (mm).
vr	Ring value; this value is added to each pixel value that fits inside the radius's; the pixels that are partially inside the ring will get fraction of the value.
vi	Inside value; this value is added to each pixel value that fits inside the inner radius; the pixels that are partially inside the radius will get fraction of the value.
vo	Outside value; this value is added to each pixel value that fits outside the outer radius; the pixels that are partially outside the radius will get fraction of the value.
verbose	Verbose level; set to <=0 to prevent all prints to stdout.

Definition at line 17 of file heart.c.

  {
  int nr, ni, no, xi, yi, i, j;
  double dx[5], dy[5], v, d;
  double ri2, ro2;
 
  if(verbose>0)
     printf("%s(img, %d, %d, %g, %g, %g, %g, %g, %g, %g, %d)\n",
            __func__, fi, zi, cx, cy, r1, r2, vr, vi, vo, verbose);
  if(r1>=r2 || r1<0.0) {
    if(verbose>0) fprintf(stderr, "Error: invalid radius.\n");
    return(1);
  }
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(zi<0 || zi>=img->dimz) return(2);
  if(img->sizey<=0.0 || img->sizex<=0.0) return(3);
  if(img->sizey!=img->sizex) return(4);
  if(img->dimt<1 || fi>=img->dimt) return(5);
  if(r1<0.0 || r2<r1) return(6);
 
  /* Compare radius ^2 to sum of squared distances instead of square roots */
  ri2=r1*r1; ro2=r2*r2;
  for(yi=0; yi<img->dimy; yi++) {
    dy[0]=(0.1+(double)yi)*img->sizey - cy;
    for(j=1; j<5; j++) dy[j]=dy[j-1]+0.2*img->sizey;
    for(xi=0; xi<img->dimx; xi++) {
      dx[0]=(0.1+(double)xi)*img->sizex - cx;
      for(i=1; i<5; i++) dx[i]=dx[i-1]+0.2*img->sizex;
      nr=ni=no=0;
      for(i=0; i<5; i++) for(j=0; j<5; j++) {
        d=dx[i]*dx[i]+dy[j]*dy[j];
        if(d<ri2) ni++; else if(d<ro2) nr++; else no++;
      }
      v=(double)nr*vr/25.0;
      v+=(double)ni*vi/25.0;
      v+=(double)no*vo/25.0;
      img->m[zi][yi][xi][fi]+=v;
    }
  }
  return(0);
}

imgSimulateSphere()

int imgSimulateSphere ( IMG * img, int fi, double cx, double cy, double cz, double r1, double r2, double vr, double vi, double vo, int verbose )

extern

Simulate a 3D image of circle with specified inner and outer radius and activity in the circle and inside and outside of it.

The applied method is only approximate at pixel borders (pixel is divided into 5x5 subpixels). PET resolution effects are not simulated.

See also

imgRingMask, imgCircleMask, imgGaussianFIRFilter, idiSimulateTubeVol

Returns

Returns 0 if successful.

Parameters

img	Pointer to allocated static or dynamic image; image must contain pixel sizes and dimensions; values are added to any existing pixel values, thus you may need to set pixel values to zero before calling this function.
fi	Frame index [0..dimt-1].
cx	X distance of circle centre (mm) from the upper left corner of the image.
cy	Y distance of circle centre (mm) from the upper left corner of the image.
cz	Z distance of circle centre (mm) from the upper left corner of the image.
r1	Inner radius of circle (mm).
r2	Outer radius of circle (mm).
vr	Sphere wall value; this value is added to each pixel value that fits inside the radius's; the pixels that are partially inside the wall will get fraction of the value.
vi	Inside value; this value is added to each pixel value that fits inside the inner radius; the pixels that are partially inside the radius will get fraction of the value.
vo	Outside value; this value is added to each pixel value that fits outside the outer radius; the pixels that are partially outside the radius will get fraction of the value.
verbose	Verbose level; set to <=0 to prevent all prints to stdout.

Definition at line 154 of file heart.c.

  {
  int nr, ni, no, xi, yi, zi, i, j, k;
  double dx[5], dy[5], dz[5], v, d;
  double ri2, ro2;
 
  if(verbose>0)
     printf("%s(img, %d, %g, %g, %g, %g, %g, %g, %g, %g, %d)\n",
            __func__, fi, cx, cy, cz, r1, r2, vr, vi, vo, verbose);
  if(r1>=r2 || r1<0.0) {
    if(verbose>0) fprintf(stderr, "Error: invalid radius.\n");
    return(1);
  }
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->sizez<=0.0 || img->sizey<=0.0 || img->sizex<=0.0) return(2);
  if(img->sizey!=img->sizex || img->sizez!=img->sizex) return(3);
  if(img->dimt<1 || fi>=img->dimt) return(4);
  if(r1<0.0 || r2<r1) return(5);
 
  /* Compare radius ^2 to sum of squared distances instead of square roots */
  ri2=r1*r1; ro2=r2*r2;
  for(zi=0; zi<img->dimz; zi++) {
    dz[0]=(0.1+(double)zi)*img->sizez - cz;
    for(k=1; k<5; k++) dz[k]=dz[k-1]+0.2*img->sizez;
    for(yi=0; yi<img->dimy; yi++) {
      dy[0]=(0.1+(double)yi)*img->sizey - cy;
      for(j=1; j<5; j++) dy[j]=dy[j-1]+0.2*img->sizey;
      for(xi=0; xi<img->dimx; xi++) {
        dx[0]=(0.1+(double)xi)*img->sizex - cx;
        for(i=1; i<5; i++) dx[i]=dx[i-1]+0.2*img->sizex;
        nr=ni=no=0;
        for(i=0; i<5; i++) for(j=0; j<5; j++) for(k=0; k<5; k++) {
          d=dx[i]*dx[i]+dy[j]*dy[j]+dz[k]*dz[k];
          if(d<ri2) ni++; else if(d<ro2) nr++; else no++;
        }
        v=(double)nr*vr/125.0;
        v+=(double)ni*vi/125.0;
        v+=(double)no*vo/125.0;
        img->m[zi][yi][xi][fi]+=v;
      }
    }
  }
  return(0);
}

rcPeakPET()

double rcPeakPET ( double FWHM, double R )

extern

Calculate the traditional recovery coefficient (RC) of peak value of a circular radioactive object in 2D PET image, assuming that the object length in 3D is relatively long, and assuming that activity of background is zero.

True object activity A0 can be calculated using RC, background activity Abkg, and measured peak activity Aexp using equation A0 = Abkg + (Aexp - Abkg)/RC , or if Abkg=0, A0 = Aexp/RC.

References: Germano et al. JNM 1992; 33: 613-620 and Brix et al. Nuklearmedizin 2002;41:184-190.

See also

idiSimulateTubeVol, heartRecoverySpilloverCorrectionFactors

Returns

Returns the recovery coefficient [0-1].

Parameters

FWHM	Full-width half-maximum value.
R	Radius of the object.

Definition at line 24 of file recovery.c.

  {
  double s, rc;
  s=FWHM/2.354820; // 2*sqrt(2*ln(2))
  rc=1.0-exp(-R*R/(2.0*s*s));
  return(rc);
}

simMyocDiameterCurve()

int simMyocDiameterCurve ( DFT * dft, double t1, double t2, double hbr, double maxdiam, double mindiam )

extern

Calculate the inner LV cavity diameter as a function of time for simulations.

10 samples per heart beat will be calculated. Fractional increase/decrease rate is fixed and coded in here.

See also

imgSimulateRing, imgRingMask, idiSimulateTubeVol

Returns

Returns 0 if successful.

Parameters

dft	Pointer to initiated DFT struct where diameter curve will be written.
t1	Start time (s).
t2	Stop time (s).
hbr	Heart rate (beats/min).
maxdiam	Maximum inner diameter (mm).
mindiam	Minimum inner diameter (mm).

Definition at line 96 of file heart.c.

  {
  int n, i, ret;
  double s;
  double rfract[10]={1.0, 0.6, 0.25, 0.0, 0.1, 0.25, 0.6, 0.8, 0.9, 0.95};
 
  /* check input */
  if(dft==NULL) return(1);
  if(t2<t1) return(1);
  if(maxdiam<mindiam) return(1);
 
  /* delete any previous curve contents */
  dftEmpty(dft);
 
  /* convert heart rate to per sec */
  if(hbr>0.0) hbr/=60.0; else hbr=0.0;
 
  /* calculate the nr of samples */
  if(hbr>0.0 && (t2-t1)>0.0) n=10.0*(t2-t1)/hbr; else n=1;
 
  /* Set up curve struct */
  ret=dftSetmem(dft, n, 1); if(ret!=0) return(3);
  dft->frameNr=n; dft->voiNr=1; dft->timetype=DFT_TIME_STARTEND;
  dft->timeunit=TUNIT_SEC;
 
  /* Fill curve */
  dft->x1[0]=t1; if(hbr>0.0) s=0.1*hbr; else s=t2-t1; 
  for(i=n=0; i<dft->frameNr; i++) {
    if(i>0) dft->x1[i]=dft->x2[i-1];
    dft->x2[i]=dft->x1[i]+s; dft->x[i]=0.5*(dft->x1[i]+dft->x2[i]);
    dft->voi[0].y[i]=mindiam+rfract[n]*(maxdiam-mindiam); n++; if(n==10) n=0;
  }
 
  return 0;
}