vessel.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
int idiSimulateTubeVol(
  VOL *vol,
  int zi,
  double cx,
  double cy, 
  double r,
  double FWHM,
  double cbkg,
  double cblo
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int idiSimulateTubeImg(
  IMG *img,
  int zi,
  double cx,
  double cy, 
  double r,
  double *cbkg,
  double *cblo
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int idiSimulateTubeImgPlane(
  int simmet,
  IMG *img,
  int zi,
  double cx,
  double cy, 
  double r,
  double FWHM,
  double *cbkg,
  double *cblo
) {
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/