p2t_3c.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Deprecated, use sim_3tcm instead.",
  " ",
  "Simulation of PET tissue time-radioactivity concentration curve (TAC)",
  "from arterial plasma (Ca) TAC based on three-tissue compartmental model,",
  "where the compartments are in series (default):"
  " ",
  "  ____    K1   ____    k3   ____    k5   ____   ",
  " | Ca | ----> | C1 | ----> | C2 | ----> | C3 |  ",
  " |____| <---- |____| <---- |____| <---- |____|  ",
  "          k2           k4           k6          ",
  " ",
  "  dC1(t)/dt = K1*Ca(T) - (k2+k3)*C1(T) + k4*C2(T) ",
  "  dC2(t)/dt = k3*C1(T) - (k4+k5)*C2(T) + k6*C3(T) ",
  "  dC3(t)/dt = k5*C2(T) - k6*C3(T)                 ",
  "  Ct(T) = C1(T) + C2(T) + C3(T)                 ",
  " ",
  ", or, optionally, three-tissue compartmental model, where the 2nd and 3rd",
  "tissue compartments are parallel, often used to represent specific and",
  "non-specific binding:",
  " ",
  "  ____    K1   ____    k3   ____   ",
  " | Ca | ----> | C1 | ----> | C2 |  ",
  " |____| <---- |____| <---- |____|  ",
  "          k2           k4          ",
  "                | ^                ",
  "             k5 | | k6             ",
  "                v |                ",
  "               ____                ",
  "              | C3 |               ",
  "              |____|               ",
  " ",
  "  dC1(t)/dt = K1*Ca(T) - (k2+k3+k5)*C1(T) + k4*C2(T) + k6*C3(T) ",
  "  dC2(t)/dt = k3*C1(T) - k4*C2(T)  ",
  "  dC3(t)/dt = k5*C1(T) - k6*C3(T)  ",
  "  Ct(T) = C1(T) + C2(T) + C3(T)    ",
  " ",
  "Usage: @P [options] plasmafile K1 k2 k3 k4 k5 k6 simfile",
  " ",
  "Options:",
  " -paral[lel]",
  "     Model with parallel compartments C2 and C3 is applied.",
  " -ser[ies]",
  "     Model with compartments C1, C2, and C3 in series is applied (default).",
  " -sub | -nosub",
  "     TACs of sub-compartments (C1, C2 and C3) are written (-sub, default)",
  "     or not written (-nosub) into the simfile.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P -nosub plasma.dat 0.3 0.2 0.1 0.2 0 0 sim.dat",
  " ",
  "See also: sim_3tcm, p2t_v3c, p2t_di, tacadd, tacren, simframe, dft2img",
  " ",
  "Keywords: TAC, simulation, modelling, 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;
  int        ret;
  int        fi, ri;
  int        save_only_total=0;
  int        parallel=0;
  DFT        plasma, sim;
  double     k1, k2, k3, k4, k5, k6, *t, *ca, *ct, *ct1, *ct2, *ct3;
  char       pfile[FILENAME_MAX], sfile[FILENAME_MAX], *cptr;
  char       tmp[FILENAME_MAX+128];
 
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&plasma); dftInit(&sim);
  pfile[0]=sfile[0]=(char)0;
  k1=k2=k3=k4=k5=k6=-1.0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(!*cptr) continue;
    if(strncasecmp(cptr, "NOSUB", 5)==0) {
      save_only_total=1; continue;
    } else if(strncasecmp(cptr, "SUB", 3)==0) {
      save_only_total=0; continue;
    } else if(strncasecmp(cptr, "PARALLEL", 5)==0) {
      parallel=1; continue;
    } else if(strncasecmp(cptr, "SERIES", 3)==0) {
      parallel=0; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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 */
  for(; ai<argc; ai++) {
    if(!pfile[0]) {
      strlcpy(pfile, argv[ai], FILENAME_MAX); continue;
    } else if(k1<0.0) {
      k1=atof_dpi(argv[ai]); if(k1>=0.0) continue;
    } else if(k2<0.0) {
      k2=atof_dpi(argv[ai]); if(k2>=0.0) continue;
    } else if(k3<0.0) {
      k3=atof_dpi(argv[ai]); if(k3>=0.0) continue;
    } else if(k4<0.0) {
      k4=atof_dpi(argv[ai]); if(k4>=0.0) continue;
    } else if(k5<0.0) {
      k5=atof_dpi(argv[ai]); if(k5>=0.0) continue;
    } else if(k6<0.0) {
      k6=atof_dpi(argv[ai]); if(k6>=0.0) continue;
    } else if(!sfile[0]) {
      strlcpy(sfile, argv[ai], FILENAME_MAX); continue;
    } else {
      fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
      return(1);
    }
  }  
 
  /* Is something missing? */
  if(!sfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(k1<=0.0) {fprintf(stderr, "Error: k1 must be > 0.\n"); return(1);}
  if(k1>1.0E+08 || k2>1.0E+08 || k3>1.0E+08 || k4>1.0E+08 ||
     k5>1.0E+08 || k6>1.0E+08) {
    fprintf(stderr, "Error: too high rate constant.\n"); return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("pfile := %s\n", pfile);
    printf("sfile := %s\n", sfile);
    printf("parameters := %g %g %g %g %g %g\n", k1, k2, k3, k4, k5, k6);
    printf("save_only_total := %d\n", save_only_total);
    printf("parallel := %d\n", parallel);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", pfile);
  if(dftRead(pfile, &plasma)) {
    fprintf(stderr, "Error in reading '%s': %s\n", pfile, dfterrmsg);
    dftEmpty(&plasma);
    return(2);
  }
  if(plasma.frameNr<3) {
    fprintf(stderr, "Error: too few samples in plasma data.\n");
    dftEmpty(&plasma); return(2);
  }
  if(plasma.voiNr>1) {
    fprintf(stderr, "Error: plasma data contains more than one curve.\n");
    dftEmpty(&plasma); return(2);
  }
 
  /*
   *  Allocate memory for simulated tissue data.
   */
  if(verbose>1) printf("allocating memory\n");
  /* allocate */
  ret=dftSetmem(&sim, plasma.frameNr, 4);
  if(ret) {
    fprintf(stderr, 
            "Error (%d): cannot allocate memory for simulated TACs.\n", ret);
    dftEmpty(&plasma); return(3);
  }
  sim.frameNr=plasma.frameNr; sim.voiNr=4;
  /* Copy times & header info */
  dftCopymainhdr(&plasma, &sim);
  for(fi=0; fi<sim.frameNr; fi++) {
    sim.x[fi]=plasma.x[fi];
    sim.x1[fi]=plasma.x1[fi]; sim.x2[fi]=plasma.x2[fi];
  }
  /* column headers */
  for(ri=0; ri<sim.voiNr; ri++) {
    sim.voi[ri].sw=0; sim.voi[ri].size=100.0;
    switch(ri) {
      case 0: strcpy(sim.voi[ri].voiname, "Ct"); break;
      case 1: strcpy(sim.voi[ri].voiname, "C1"); break;
      case 2: strcpy(sim.voi[ri].voiname, "C2"); break;
      case 3: strcpy(sim.voi[ri].voiname, "C3"); break;
    }
    strcpy(sim.voi[ri].hemisphere, ""); strcpy(sim.voi[ri].place, "");
    strcpy(sim.voi[ri].name, sim.voi[ri].voiname);
  }
 
 
  /*
   *  Simulate TACs
   */
  if(verbose>1) printf("simulating\n");
  t=sim.x; ca=plasma.voi[0].y; ct=sim.voi[0].y;
  ct1=sim.voi[1].y; ct2=sim.voi[2].y; ct3=sim.voi[3].y;
  if(parallel==0)
    ret=simC3s(t, ca, sim.frameNr, k1, k2, k3, k4, k5, k6, ct, ct1, ct2, ct3);
  else
    ret=simC3p(t, ca, sim.frameNr, k1, k2, k3, k4, k5, k6, ct, ct1, ct2, ct3);
  dftEmpty(&plasma);
  if(ret!=0) {
    fprintf(stderr, "Error (%d) in simulation.\n", ret);
    dftEmpty(&sim); return(8);
  }
 
 
  /*
   *  Save simulated TACs
   */
  if(verbose>1) printf("saving PET curves\n");
  DFT_NR_OF_DECIMALS=6;
  if(save_only_total) sim.voiNr=1;
  /* Set comments */
  dftSetComments(&sim);
  sprintf(tmp, "# plasmafile := %s\n", pfile); strcat(sim.comments, tmp);
  strcpy(tmp, "# model := ");
  if(parallel==0) strcat(tmp, "C3S\n"); else strcat(tmp, "C3P\n");
  strcat(sim.comments, tmp);
  sprintf(tmp, "# K1 := %g\n", k1); strcat(sim.comments, tmp);
  sprintf(tmp, "# k2 := %g\n", k2); strcat(sim.comments, tmp);
  sprintf(tmp, "# k3 := %g\n", k3); strcat(sim.comments, tmp);
  sprintf(tmp, "# k4 := %g\n", k4); strcat(sim.comments, tmp);
  sprintf(tmp, "# k5 := %g\n", k5); strcat(sim.comments, tmp);
  sprintf(tmp, "# k6 := %g\n", k6); strcat(sim.comments, tmp);
  /* Set file format to PMOD, if extension is .tac */
  if(!strcasecmp(filenameGetExtension(sfile), ".tac"))
    sim._type=DFT_FORMAT_PMOD;
  /* Some format has to be set anyway, and simple format would lose information */
  if(sim._type==DFT_FORMAT_PLAIN || sim._type==DFT_FORMAT_UNKNOWN)
    sim._type=DFT_FORMAT_STANDARD;
  /* Write file */
  if(dftWrite(&sim, sfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", sfile, dfterrmsg);
    dftEmpty(&sim); return(11);
  }
  if(verbose>0) fprintf(stdout, "simulated TAC(s) written in %s\n", sfile);
  dftEmpty(&sim);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/