litac.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 "tpcli.h"
#include "tpctac.h"
/*****************************************************************************/
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacIntegrate(
  TAC *inp,
  TAC *out,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || out==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  tacFree(out);
  int ret=TPCERROR_OK;
  ret=tacDuplicate(inp, out);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(inp->isframe==0) {
    for(int i=0; i<inp->tacNr; i++) {
      ret=liIntegrate(inp->x, inp->c[i].y, inp->sampleNr, out->c[i].y, 3, 0);
      if(ret!=TPCERROR_OK) break;
    }
  } else {
    for(int i=0; i<inp->tacNr; i++) {
      ret=liIntegratePET(inp->x1, inp->x2, inp->c[i].y, inp->sampleNr, out->c[i].y, NULL, 0);
      if(ret!=TPCERROR_OK) break;
    }
  }
  if(ret!=TPCERROR_OK) tacFree(out);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  return(ret);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacInterpolate(
  TAC *inp,
  TAC *xinp,
  TAC *tac,
  TAC *itac,
  TAC *iitac,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || xinp==NULL || xinp==tac || xinp==itac || xinp==iitac) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tac==NULL && itac==NULL && iitac==NULL) { // there's nothing to do
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
  /* Check the function input */
  if(inp->sampleNr<1 || inp->tacNr<1 || xinp->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check that input does not have any missing values */
  if(tacNaNs(inp)>0 || tacXNaNs(xinp)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
 
  /* If input data have similar x values as requested already, then
     copy x1 and x2 into input if available */
  if(xinp->isframe && tacXMatch(inp, xinp, verbose-2) && inp->sampleNr<=xinp->sampleNr) {
    for(int i=0; i<inp->sampleNr; i++) {
      inp->x[i]=xinp->x[i]; inp->x1[i]=xinp->x1[i]; inp->x2[i]=xinp->x2[i];
    }
    inp->isframe=xinp->isframe;
  }
 
 
  /* Check that there is no need for excess extrapolation */
  {
    double start1, end1, start2, end2, range2;
    if(tacXRange(inp, &start1, &end1) || tacXRange(xinp, &start2, &end2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    range2=end2-start2;
    if(verbose>2) printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
    if(start1>start2+0.2*range2 || end1<end2-0.3*range2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
 
  /* Delete any previous contents in output structs */
  if(tac!=NULL) tacFree(tac);
  if(itac!=NULL) tacFree(itac);
  if(iitac!=NULL) tacFree(iitac);
 
  /* Allocate memory for output data */
  int ret=TPCERROR_OK;
  if(tac!=NULL) ret=tacAllocate(tac, xinp->sampleNr, inp->tacNr); 
  if(ret==TPCERROR_OK && itac!=NULL) 
    ret=tacAllocate(itac, xinp->sampleNr, inp->tacNr);
  if(ret==TPCERROR_OK && iitac!=NULL) 
    ret=tacAllocate(iitac, xinp->sampleNr, inp->tacNr);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  if(tac!=NULL) {tac->tacNr=inp->tacNr; tac->sampleNr=xinp->sampleNr;}
  if(itac!=NULL) {itac->tacNr=inp->tacNr; itac->sampleNr=xinp->sampleNr;}
  if(iitac!=NULL) {iitac->tacNr=inp->tacNr; iitac->sampleNr=xinp->sampleNr;}
 
 
  /* Copy header information */
  if(tac!=NULL) {
    tacCopyHdr(inp, tac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+i);
  }
  if(itac!=NULL) {
    tacCopyHdr(inp, itac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, itac->c+i);
  }
  if(iitac!=NULL) {
    tacCopyHdr(inp, iitac);
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+i);
  }
  /* Copy X (time) information */
  if(tac!=NULL) {
    tac->isframe=xinp->isframe; tacXCopy(xinp, tac, 0, tac->sampleNr-1);
  }
  if(itac!=NULL) {
    itac->isframe=xinp->isframe; tacXCopy(xinp, itac, 0, itac->sampleNr-1);
  }
  if(iitac!=NULL) {
    iitac->isframe=xinp->isframe; tacXCopy(xinp, iitac, 0, iitac->sampleNr-1);
  }
 
  /* Check if TACs have the same frame times already */
  if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=xinp->sampleNr) {
    if(verbose>2) printf("  same frame times\n");
    int ret=0, ri, fi;
    double *i1, *i2;
    /* copy the values directly and integrate in place, if requested */
    for(ri=0; ri<inp->tacNr && ret==0; ri++) {
      if(tac!=NULL) 
        for(fi=0; fi<tac->sampleNr; fi++)
          tac->c[ri].y[fi]=inp->c[ri].y[fi];
      /* If integrals not requested then that's it for this TAC */
      if(itac==NULL && iitac==NULL) continue;
      /* otherwise set pointers for integrals and then integrate */
      if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;      
      if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
      if(inp->isframe)
        ret=liIntegratePET(inp->x1, inp->x2, inp->c[ri].y,
                           xinp->sampleNr, i1, i2, verbose-2);
      else
        ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                          xinp->x, NULL, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
    }
    if(ret) {
      if(tac!=NULL) tacFree(tac); 
      if(itac!=NULL) tacFree(itac);
      if(iitac!=NULL) tacFree(iitac);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK); // that's it then
  }
 
  /* Interpolate and/or integrate inp data to tac sample times,
     taking into account tac frame lengths if available */
  if(verbose>2) printf("  different frame times\n");
  ret=0;
  double *i0, *i1, *i2;
  for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
    /* set pointers for output y arrays */
    if(tac!=NULL) i0=tac->c[ri].y; else i0=NULL;      
    if(itac!=NULL) i1=itac->c[ri].y; else i1=NULL;      
    if(iitac!=NULL) i2=iitac->c[ri].y; else i2=NULL;
    if(xinp->isframe)
      ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
                xinp->x1, xinp->x2, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
    else
      ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                xinp->x, i0, i1, i2, xinp->sampleNr, 4, 1, verbose-2);
  }
  if(ret) {
    if(tac!=NULL) tacFree(tac); 
    if(itac!=NULL) tacFree(itac);
    if(iitac!=NULL) tacFree(iitac);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacInterpolateInto(
  TAC *inp,
  TAC *tac,
  TAC *itac,
  TAC *iitac,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  
  /* Check the function input */
  if(tac->sampleNr<1 || inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(itac!=NULL && itac->sampleNr!=tac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
  if(iitac!=NULL && iitac->sampleNr!=tac->sampleNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
    return TPCERROR_INVALID_X;
  }
 
 
  /* Check that input does not have any missing values */
  /* Check that there are no missing sample times in output TACs */
  if(tacNaNs(inp)>0 || tacXNaNs(tac)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
  
  /* Check that there is no need for excess extrapolation */
  {
    double start1, end1, start2, end2, range2;
    if(tacXRange(inp, &start1, &end1) || tacXRange(tac, &start2, &end2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    range2=end2-start2;
    if(verbose>2) 
      printf("time ranges: %g - %g vs %g - %g\n", start1, end1, start2, end2);
    if(start1>(start2+0.2*range2) || end1<(end2-0.3*range2)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
  }
 
  /* Allocate memory for interpolated data */
  if(tacAllocateMore(tac, inp->tacNr)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
  /* Copy header information */
  for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, tac->c+tac->tacNr+i);
  /* Same for the integrals */
  if(itac!=NULL) {
    if(tacAllocateMore(itac, inp->tacNr)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      return TPCERROR_OUT_OF_MEMORY;
    }
    for(int i=0; i<inp->tacNr; i++) 
      tacCopyTacchdr(inp->c+i, itac->c+itac->tacNr+i);
  }
  if(iitac!=NULL) {
    if(tacAllocateMore(iitac, inp->tacNr)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      return TPCERROR_OUT_OF_MEMORY;
    }
    for(int i=0; i<inp->tacNr; i++) tacCopyTacchdr(inp->c+i, iitac->c+iitac->tacNr+i);
  }
 
  /* Check if TACs have the same frame times already */
  if(tacXMatch(inp, tac, verbose-2) && inp->sampleNr>=tac->sampleNr) {
    if(verbose>2) printf("same frame times\n");
    int ret=0, ri, fi;
    double *i1, *i2;
    /* copy the values directly and integrate in place, if requested */
    for(ri=0; ri<inp->tacNr && ret==0; ri++) {
      for(fi=0; fi<tac->sampleNr; fi++)
        tac->c[tac->tacNr+ri].y[fi]=inp->c[ri].y[fi];
      /* If integrals not requested then that's it for this TAC */
      if(itac==NULL && iitac==NULL) continue;
      /* otherwise set pointers for integrals and then integrate */
      if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;      
      if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
      if(tac->isframe)
        ret=liIntegratePET(tac->x1, tac->x2, tac->c[tac->tacNr+ri].y,
                           tac->sampleNr, i1, i2, verbose-2);
      else
        ret=liInterpolate(tac->x, tac->c[tac->tacNr+ri].y, tac->sampleNr,
               tac->x, NULL, i1, i2, tac->sampleNr, 4, 1, verbose-2);
    }
    if(ret) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
    tac->tacNr+=inp->tacNr;
    if(itac!=NULL) itac->tacNr+=inp->tacNr;
    if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK); // that's it then
  }
 
  /* Interpolate and, if requested, integrate inp data to tac sample times,
     taking into account tac frame lengths if available */
  int ret=0;
  double *i1, *i2;
  for(int ri=0; ri<inp->tacNr && ret==0; ri++) {
    /* set pointers for integrals */
    if(itac!=NULL) i1=itac->c[itac->tacNr+ri].y; else i1=NULL;
    if(iitac!=NULL) i2=iitac->c[iitac->tacNr+ri].y; else i2=NULL;
    if(tac->isframe)
      ret=liInterpolateForPET(inp->x, inp->c[ri].y, inp->sampleNr,
                tac->x1, tac->x2, tac->c[tac->tacNr+ri].y, 
                i1, i2, tac->sampleNr, 4, 1, verbose-2);
    else
      ret=liInterpolate(inp->x, inp->c[ri].y, inp->sampleNr,
                tac->x, tac->c[tac->tacNr+ri].y, 
                i1, i2, tac->sampleNr, 4, 1, verbose-2);
  }
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  tac->tacNr+=inp->tacNr;
  if(itac!=NULL) itac->tacNr+=inp->tacNr;
  if(iitac!=NULL) iitac->tacNr+=inp->tacNr;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
double tacAUC(
  TAC *tac,
  int ti,
  double t1,
  double t2,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, %d, %g, %g)\n", __func__, ti, t1, t2);
  /* Check the function input */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(nan(""));
  }
  if(ti<0 || ti>=tac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(nan(""));
  }
  if(t1>=t2) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return(nan(""));
  }
 
  /* Get and check the x range */
  double xmin=nan(""), xmax=nan("");
  if(tacXRange(tac, &xmin, &xmax)!=0 || t1>=xmax || t2<=xmin) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return(nan(""));
  }
 
  /* Make temporary x and yi array */
  double tx[2]; tx[0]=t1; tx[1]=t2;
  double tyi[2];
 
  /* If TAC contains just frame mid times, then simple dot-to-dot integration */
  if(tac->isframe==0) {
    if(liInterpolate(tac->x, tac->c[ti].y, tac->sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(nan(""));
    }
    /* That's it then */
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(tyi[1]-tyi[0]);
  }
 
  /* We have frame start and end times */
  /* Convert the data into stepwise dot-to-dot data */
  TAC stac; tacInit(&stac);
  {
    int ret=tacFramesToSteps(tac, &stac, NULL);
    if(ret!=TPCERROR_OK) {
      tacFree(&stac);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(nan(""));
    }
  }
  /* Then simple dot-to-dot integration */
  if(liInterpolate(stac.x, stac.c[ti].y, stac.sampleNr, tx, NULL, tyi, NULL, 2, 3, 1, verbose-10))
  {
    tacFree(&stac);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return(nan(""));
  }
  tacFree(&stac);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(tyi[1]-tyi[0]);
}
/*****************************************************************************/
 
/*****************************************************************************/
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);
}
/*****************************************************************************/
 
/*****************************************************************************/