lisim.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcstatist.h"
#include "tpctac.h"
#include "tpcli.h"
/*****************************************************************************/
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacInput2sim(
  TAC *itac,
  TAC *ttac,
  TAC *stac,
  TPCSTATUS *status
) {
  if(status==NULL) return TPCERROR_FAIL;
  int verbose=status->verbose;
  if(verbose>0) printf("%s(itac, ttac, ...)\n", __func__);
 
  /* Check the input */
  if(itac==NULL || ttac==NULL || stac==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(itac->sampleNr<1 || ttac->sampleNr<1 || itac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Delete any previous data */
  tacFree(stac);
 
  /* Make a working copy of input data */
  TAC inp; tacInit(&inp);
  status->error=tacDuplicate(itac, &inp); if(status->error!=TPCERROR_OK) return(status->error);
  /* Make sure that input data contains also frame mid times */
  if(inp.isframe) tacSetX(&inp, NULL);
  /* Convert time units to TTAC units */
  tacXUnitConvert(&inp, ttac->tunit, status);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
 
  /* Sort the input data by sample times */
  tacSortByTime(&inp, status); if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  /* Add zero sample if necessary */
  tacAddZeroSample(&inp, status);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
 
  /* Make a combined list of sample times based on input and tissue TACs */
  unsigned int maxn=inp.sampleNr+ttac->sampleNr;
  if(inp.isframe) maxn+=2*inp.sampleNr;
  if(ttac->isframe) maxn+=ttac->sampleNr;
  double x[maxn];
 
  int n=0;
  for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x[i])) {
    int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x[i])<1.0E-20) break;
    if(j==n) x[n++]=inp.x[i];
  }
  if(inp.isframe) {
    for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x1[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x1[i])<1.0E-20) break;
      if(j==n) x[n++]=inp.x1[i];
    }
    for(int i=0; i<inp.sampleNr; i++) if(isfinite(inp.x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-inp.x2[i])<1.0E-20) break;
      if(j==n) x[n++]=inp.x2[i];
    }
  }
  if(ttac->isframe) {
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x2[i];
    }
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x2[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x2[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x2[i];
    }
  } else {
    for(int i=0; i<ttac->sampleNr; i++) if(isfinite(ttac->x[i])) {
      int j; for(j=0; j<n; j++) if(fabs(x[j]-ttac->x[i])<1.0E-20) break;
      if(j==n) x[n++]=ttac->x[i];
    }
  }
  if(verbose>2) printf("n=%d\n", n);
  statSortDouble(x, n, 0);
 
  /* Put new sample times into output data */
  status->error=tacDuplicate(itac, stac);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  status->error=tacAllocateMoreSamples(stac, n-stac->sampleNr);
  if(status->error!=TPCERROR_OK) {tacFree(&inp); return(status->error);}
  stac->sampleNr=n; stac->tunit=inp.tunit; stac->isframe=0;
  for(int i=0; i<n; i++) stac->x[i]=x[i];
 
  /* Interpolate input data into the new sample times */
  for(int i=0; i<inp.tacNr; i++) {
    status->error=liInterpolate(inp.x, inp.c[i].y, inp.sampleNr,
                            stac->x, stac->c[i].y, NULL, NULL, stac->sampleNr, 4, 1, verbose-10);
    if(status->error!=TPCERROR_OK) break;
  }
  if(status->error!=TPCERROR_OK) {tacFree(&inp); tacFree(stac); return(status->error);}
 
  tacFree(&inp);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacVb(
  TAC *ttac,
  const int i,
  TAC *btac,
  double Vb,
  const int simVb,
  const int petVolume,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(ttac, %d, btac, %g, %d, %d)\n", __func__, i, Vb, simVb, petVolume);
  if(ttac==NULL || btac==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(ttac->tacNr<1 || btac->tacNr<1 || ttac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(ttac->sampleNr>btac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
    return TPCERROR_INCOMPATIBLE_DATA;
  }
  if(!(Vb>=0.0) || !(Vb<=1.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(Vb==1.0 && petVolume==0 && simVb==0) { // combination would lead to divide by zero
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(i>=ttac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
 
  if(simVb==0) {
    if(verbose>1) printf("subtracting blood\n");
    for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
      for(int fi=0; fi<ttac->sampleNr; fi++) {
        ttac->c[ri].y[fi]-=Vb*btac->c[0].y[fi];
        if(petVolume==0) ttac->c[ri].y[fi]/=(1.0-Vb);
      }
    }
  } else {
    if(verbose>1) printf("adding blood\n");
    for(int ri=0; ri<ttac->tacNr; ri++) if(i<0 || i==ri) {
      for(int fi=0; fi<ttac->sampleNr; fi++) {
        if(petVolume==0) ttac->c[ri].y[fi]*=(1.0-Vb);
        ttac->c[ri].y[fi]+=Vb*btac->c[0].y[fi];
      }
    }
  }
 
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacInterpolateToEqualLengthFrames(
  TAC *inp,
  double minfdur,
  double maxfdur,
  TAC *tac,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(tac1, %g, %g, ...)\n", __func__, minfdur, maxfdur); fflush(stdout);}
  if(tac==NULL || inp==NULL || inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(minfdur>0.0 && maxfdur>0.0 && minfdur>maxfdur) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
  if(!tacIsX(inp)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
 
  /* Get the minimum sample interval (that is higher than 0) */
  double freq=nan("");
  int ret=tacGetSampleInterval(inp, 1.0E-08, &freq, NULL);
  //printf("dataset based freq: %g\n", freq);
  if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return(ret);}
  if(minfdur>0.0 && freq<minfdur) freq=minfdur;
  if(maxfdur>0.0 && freq>maxfdur) freq=maxfdur;
  //printf("refined freq: %g\n", freq);
 
  /* Get the x range */
  double xmax=nan("");
  if(tacXRange(inp, NULL, &xmax)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  //printf("xmax: %g\n", xmax);
  if(freq>0.5*xmax) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
 
 
  /* Calculate the number of required interpolated samples */
  int nreq=ceil(xmax/freq); //printf("required_n := %d\n", nreq);
  /* Check that it is not totally stupid */
  if(nreq<1 || nreq>2E+06) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_BIG);
    return TPCERROR_TOO_BIG;
  }
 
  /* Setup the output TAC */
  if(tacAllocate(tac, nreq, inp->tacNr)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  tac->tacNr=inp->tacNr; tac->sampleNr=nreq;
  /* Copy header information */
  tacCopyHdr(inp, tac);
  for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
 
  /* Set X (sample time) information */
  tac->isframe=1;
  for(int i=0; i<nreq; i++) tac->x1[i]=(double)i*freq;
  for(int i=0; i<nreq; i++) tac->x2[i]=(double)(i+1)*freq;
  for(int i=0; i<nreq; i++) tac->x[i]=0.5*(tac->x1[i]+tac->x2[i]);
 
  /* Calculate mean value during each interval */
  if(verbose>2) {printf("  interpolating...\n"); fflush(stdout);}
  for(int i=0; i<inp->tacNr; i++) {
    if(liInterpolateForPET(inp->x, inp->c[i].y, inp->sampleNr, tac->x1, tac->x2, 
                           tac->c[i].y, NULL, NULL, tac->sampleNr, 2, 1, 0)!=0) {
      tacFree(tac);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
  if(verbose>2) {printf("  ... done.\n"); fflush(stdout);}
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/