heart.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
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 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);
}
/*****************************************************************************/
 
/*****************************************************************************/
int simMyocDiameterCurve(
  DFT *dft,
  double t1,
  double t2,
  double hbr,
  double maxdiam,
  double mindiam
) {
  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;
}
/*****************************************************************************/
 
/*****************************************************************************/
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 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);
}
/*****************************************************************************/
 
/*****************************************************************************/