simcmdk.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcfileutil.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpctacmod.h"
#include "tpccm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Simulation of PET tissue time-radioactivity concentration curve (TAC)",
  "based on two-tissue compartment model. Vascular volume is not considered.",
  " ",
  "  dC1(t)/dt = K1*Ca(t) - (k2+k3)*C1(t) + k2*C2(t)",
  "  dC2(t)/dt = k3*C1(t) - k4*C2(t)",
  "  Ct(t) = C1(t) + C2(t) ",
  " ",
  "Usage: @P [options] inputfile simfile",
  " ",
  "Options:",
  " -fa | -fe",
  "     If input file contains frame start and end times, simulated TAC is",
  "     either saved with frame start and end times and frame concentration",
  "     average (option -fa, default), or with start/end times and concentrations",
  "     at those times (option -fe).",
  " -sub",
  "     Save simulated C1 and C2 curves; by default only their sum is saved.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Model parameters are specified as additional columns of the input file.",
  "The first one or two columns must contain the times or frame start and end",
  "times, and after that the input function concentrations.",
  " ",
  "See also: sim_3tcm, p2t_3c, tacadd, simframe, taccbv",
  " ",
  "Keywords: TAC, simulation, compartmental model",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int    ai, help=0, version=0, verbose=1;
  char   inpfile[FILENAME_MAX], simfile[FILENAME_MAX];
  int    framemode=0; // 0=frame start and end times with average concentration;
                      // 1=frame start/end times with concentration of the moment.
  int    savesub=0;
 
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=simfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "FA")==0) {
      framemode=0; continue;
    } else if(strcasecmp(cptr, "FE")==0) {
      framemode=1; continue;
    } else if(strcasecmp(cptr, "SUB")==0) {
      savesub=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(simfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!simfile[0]) {fprintf(stderr, "Error: missing file name.\n"); return(1);}
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("simfile := %s\n", simfile);
    printf("framemode := %d\n", framemode);
    printf("savesub := %d\n", savesub);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read input TAC
   */
  if(verbose>1) fprintf(stdout, "reading input TAC\n");
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, inpfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("frames := %d\n", tac.isframe);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in input TAC.\n");
    tacFree(&tac); return(2);
  }
  if(tac.tacNr<2) {
    fprintf(stderr, "Error: model parameters missing from input TAC.\n");
    tacFree(&tac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: missing values in %s.\n", inpfile);
    tacFree(&tac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Fix small frame gaps and overlaps */
  if(tac.isframe==1) {
    if(tacCorrectFrameOverlap(&tac, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: times frames not consecutive in %s.\n", inpfile);
      tacFree(&tac); return(2);
    }
  }
 
 
  /*
   *  Find the model parameter columns
   */
  int ik1, ik2, ik3, ik4;
  ik1=ik2=ik3=ik4=0;
  if(tacSelectTACs(&tac, "K1", 1, NULL)==1) ik1=tacFirstSelected(&tac);
  if(tacSelectTACs(&tac, "k2", 1, NULL)==1) ik2=tacFirstSelected(&tac);
  if(tacSelectTACs(&tac, "k3", 1, NULL)==1) ik3=tacFirstSelected(&tac);
  if(tacSelectTACs(&tac, "k4", 1, NULL)==1) ik4=tacFirstSelected(&tac);
  if(ik1==0) {
    fprintf(stderr, "Error: missing K1 in %s.\n", inpfile);
    tacFree(&tac); return(3);
  }
  if(ik3==0) ik4=0;
  if(verbose>1) {
    printf("K1 column := %d\n", 1+ik1);
    if(ik2>0) printf("k2 column := %d\n", 1+ik2);
    if(ik3>0) printf("k3 column := %d\n", 1+ik3);
    if(ik4>0) printf("k4 column := %d\n", 1+ik4);
  }
 
  /* Allocate space for simulated data */
  if((status.error=tacAllocateMore(&tac, 3))!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  int iCt=tac.tacNr;
  int iC1=tac.tacNr+1;
  int iC2=tac.tacNr+2;
  tac.tacNr+=3;
  strcpy(tac.c[iCt].name, "Ct");
  strcpy(tac.c[iC1].name, "C1");
  strcpy(tac.c[iC2].name, "C2");
 
 
  /*
   *  Simulate
   */
  if(tac.isframe==0) {
    for(int i=0; i<tac.sampleNr; i++) {
      double K1=0.0, k2=0.0, k3=0.0, k4=0.0;
      if(i==0) {
        tac.c[iCt].y[0]=tac.c[iC1].y[0]=tac.c[iC2].y[0]=0.0;
        continue;
      }
      K1=tac.c[ik1].y[i-1];
      if(ik2>0) k2=tac.c[ik2].y[i-1];
      if(ik3>0) k3=tac.c[ik3].y[i-1];
      if(ik4>0) k4=tac.c[ik4].y[i-1];
      double dt=(tac.x[i]-tac.x[i-1]);
      if(!(dt>0.0)) {
        tac.c[iCt].y[i]=tac.c[iCt].y[i-1];
        tac.c[iC1].y[i]=tac.c[iC1].y[i-1];
        tac.c[iC2].y[i]=tac.c[iC2].y[i-1];
        continue;
      }
      double dt2=0.5*dt;
      tac.c[iC1].y[i]= (
         K1*(1.0+k4*dt2)*dt2*(tac.c[0].y[i-1]+tac.c[0].y[i]) +
         (1.0 - dt2*(k2+k3-k4+k2*k4*dt2))*tac.c[iC1].y[i-1] + 
         k4*dt*tac.c[iC2].y[i-1]
        ) / (
         1.0 + dt2*(k2+k3+k4+k2*k4*dt2)
        );
      tac.c[iC2].y[i]= (
        k3*dt2*(tac.c[iC1].y[i-1]+tac.c[iC1].y[i]) + (1.0 - k4*dt2)*tac.c[iC2].y[i-1]
        ) / (
         1.0 + dt2*k4
        );
      tac.c[iCt].y[i]=tac.c[iC1].y[i]; if(finite(tac.c[iC2].y[i])) tac.c[iCt].y[i]+=tac.c[iC2].y[i];
      if(verbose>3 && i==tac.sampleNr-1) printf("K1=%g k2=%g k3=%g k4=%g\n", K1, k2, k3, k4);
    }
  } else {
    for(int i=0; i<tac.sampleNr; i++) {
      double K1=0.0, k2=0.0, k3=0.0, k4=0.0;
      K1=tac.c[ik1].y[i];
      if(ik2>0) k2=tac.c[ik2].y[i];
      if(ik3>0) k3=tac.c[ik3].y[i];
      if(ik4>0) k4=tac.c[ik4].y[i];
      double dt=tac.x2[i]-tac.x1[i];
      if(!(dt>0.0)) {
        if(i==0) {
          tac.c[iCt].y[0]=tac.c[iC1].y[0]=tac.c[iC2].y[0]=0.0;
        } else {
          tac.c[iCt].y[i]=tac.c[iCt].y[i-1];
          tac.c[iC1].y[i]=tac.c[iC1].y[i-1];
          tac.c[iC2].y[i]=tac.c[iC2].y[i-1];
        }
        continue;
      }
      double dt2=0.5*dt;
      if(i==0) {
        tac.c[iC1].y[i]= (
           K1*(1.0+k4*dt2)*dt*tac.c[0].y[i]
          ) / (
           1.0 + dt2*(k2+k3+k4+k2*k4*dt2)
          );
        tac.c[iC2].y[i]= (
           k3*dt2*(tac.c[iC1].y[i])
          ) / (
           1.0 + dt2*k4
          );
      } else {
        tac.c[iC1].y[i]= (
           K1*(1.0+k4*dt2)*dt*tac.c[0].y[i] +
           (1.0 - dt2*(k2+k3-k4+k2*k4*dt2))*tac.c[iC1].y[i-1] +
           k4*dt*tac.c[iC2].y[i-1]
          ) / (
           1.0 + dt2*(k2+k3+k4+k2*k4*dt2)
          );
        tac.c[iC2].y[i]= (
          k3*dt2*(tac.c[iC1].y[i-1]+tac.c[iC1].y[i]) + (1.0 - k4*dt2)*tac.c[iC2].y[i-1]
          ) / (
           1.0 + dt2*k4
          );
      }
      tac.c[iCt].y[i]=tac.c[iC1].y[i]; if(finite(tac.c[iC2].y[i])) tac.c[iCt].y[i]+=tac.c[iC2].y[i];
      if(verbose>3 && i==tac.sampleNr-1) printf("K1=%g k2=%g k3=%g k4=%g\n", K1, k2, k3, k4);
    }
  }
  if(tac.isframe) {
    if(framemode==0) {
      /* Calculate frame averages, if requested */
      for(int i=tac.sampleNr-1; i>0; i--) {
        tac.c[iCt].y[i]=0.5*(tac.c[iCt].y[i-1]+tac.c[iCt].y[i]);
        tac.c[iC1].y[i]=0.5*(tac.c[iC1].y[i-1]+tac.c[iC1].y[i]);
        tac.c[iC2].y[i]=0.5*(tac.c[iC2].y[i-1]+tac.c[iC2].y[i]);
      }
      tac.c[iCt].y[0]*=0.5;
      tac.c[iC1].y[0]*=0.5;
      tac.c[iC2].y[0]*=0.5;
    } else {
      /* Save concentrations at frame end times, like they were computed,
         and set the times to match that. */
      for(int i=0; i<tac.sampleNr; i++) tac.x[i]=tac.x2[i];
      tac.isframe=0;
      /* Add also the start time of the first frame. Function tacAddZeroSample() does not add zero 
         time, in the unlikely case it already exists, and in that case we are better of without. */
      int prevNr=tac.sampleNr;
      tacAddZeroSample(&tac, &status); 
      if(tac.sampleNr>prevNr) tac.x[0]=tac.x1[1];
    }
  }
 
 
  /*
   *  Save simulated data 
   */
  if(verbose>1) printf("writing %s\n", simfile);
  tacMoveTACC(&tac, iCt, 0);
  tac.tacNr=1;
  if(savesub) {
    tacMoveTACC(&tac, iC1, 1);
    tacMoveTACC(&tac, iC2, 2);
    tac.tacNr=3;
  }
  FILE *fp; fp=fopen(simfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
    tacFree(&tac); return(11);
  }
  tacWrite(&tac, fp, TAC_FORMAT_PMOD, 1, &status);
  fclose(fp); tacFree(&tac);
  if(status.error!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) {printf("%s saved.\n", simfile); fflush(stdout);}
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/