yokoi.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 "tpctac.h"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate K1, k2, and K1/k2 using multiple-time graphical analysis (MTGA)",
  "for reversible PET ligands (Yokoi plot) (1,2,3) from regional PET TACs.",
  " ",
  "Usage: @P [options] tacfile input starttime endtime resultfile",
  " ",
  "Options:",
  " -sd=<y|N>",
  "     Standard deviations are saved (y) or not saved (n, default) in results.",
  " -svg=<Filename>",
  "     Plots are written in specified file in Scalable Vector Graphics (SVG) 1.1",
  "     format; specification in https://www.w3.org/TR/SVG/",
  " -plot=<Filename>",
  "     Data for plots is written in TSV format for easy importing in Excel",
  "     or OpenOffice spreadsheet, where the data can be viewed.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "References:",
  "1. Yokoi T, Iida H, Itoh H, Kanno I. A new graphic plot analysis for cerebral",
  "   blood flow and partition coefficient with iodine-123-iodoamphetamine",
  "   and dynamic SPECT validation studies using oxygen-15-water and PET.",
  "   J Nucl Med. 1993; 34(3): 498-505.",
  "2. Ito H, Yokoi T, Ikoma Y, et al. A new graphic plot analysis for",
  "   determination of neuroreceptor binding in positron emission tomography",
  "   studies. NeuroImage 2010; 49(1): 578-586.",
  "3. Ito H, Ikoma Y, Seki C, et al. Visual evaluation of kinetic",
  "   characteristics of PET probe for neuroreceptors using a two-phase graphic",
  "   plot analysis. Ann Nucl Med. 2017; 31(4): 273-282.",
  " ",
  "See also: logan, patlak, fitk2, lhsol",
  " ",
  "Keywords: TAC, MTGA, modelling, K1, Vt",
  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 tacfile[FILENAME_MAX], inpfile[FILENAME_MAX], resfile[FILENAME_MAX];
  char svgfile[FILENAME_MAX], plotfile[FILENAME_MAX];
  int save_stat=0; // 0=do not save, 1=save SD and/or other fit statistics
  double tstart=nan(""), tstop=nan(""); // fit start and end times
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=inpfile[0]=resfile[0]=svgfile[0]=plotfile[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(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); if(strlen(svgfile)>0) continue;
    } else if(strncasecmp(cptr, "PLOT=", 5)==0) {
      strlcpy(plotfile, cptr+5, FILENAME_MAX); if(strlen(plotfile)>0) continue;
    } else if(strncasecmp(cptr, "SD=", 3)==0) {
      save_stat=tpcYesNo(cptr+3); if(save_stat>=0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // later arguments may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-10;
  
  /* 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(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &tstart)) {
      fprintf(stderr, "Error: invalid start time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) {
    if(atofCheck(argv[ai], &tstop)) {
      fprintf(stderr, "Error: invalid stop time '%s'.\n", argv[ai]); return(1);}
    ai++;
  }
  if(ai<argc) strlcpy(resfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!resfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Is something wrong? */
  if(!isnan(tstart) && !isnan(tstop) && tstop<=tstart) {
    fprintf(stderr, "Error: illegal time range %g-%g\n", tstart, tstop);
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("inpfile := %s\n", inpfile);
    printf("resfile := %s\n", resfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(plotfile[0]) printf("plotfile := %s\n", plotfile);
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("save_stat := %d\n", save_stat);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>0) printf("reading tissue and input data\n");
  TAC tac, input; tacInit(&tac); tacInit(&input);
  int fitSampleNr;
  double fitdur;
  if(tstop>0.01) fitdur=tstop; else fitdur=1.0E+10;
 
  if(tacReadModelingData(tacfile, inpfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &tac, &input, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&input); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("tac.sampleNr := %d\n", tac.sampleNr);
    printf("input.sampleNr := %d\n", input.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
 
  /* Get the indices of line fit range */
  int istart=doubleGEIndex(tac.x, fitSampleNr, tstart);
  int istop=doubleGTIndex(tac.x, fitSampleNr, tstop) - 1;
  if(verbose>3) {
    printf("istart := %d\n", istart);
    printf("istop := %d\n", istop);
  }
  int fitNr=1+istop-istart;
  if(fitNr<2) {
    fprintf(stderr, "Error: invalid line fit range.\n");
    tacFree(&tac); tacFree(&input); return(2);
  }
 
 
  /*
   *  Interpolate and integrate PTAC to TTAC times
   */
  if(verbose>1) printf("integrating and interpolating input TAC\n");
  TAC input0; tacInit(&input0); // interpolated PTAC
  TAC input1; tacInit(&input1); // 1st integral 
  if(tacInterpolate(&input, &tac, &input0, &input1, NULL, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&input); tacFree(&input0); tacFree(&input1);
    return(3);
  }
  /* Original input TAC should not be needed any more */
  tacFree(&input);
 
  /* Integrate TTAC */
  if(verbose>1) printf("integrating TTAC\n");
  TAC tac1; tacInit(&tac1); // 1st integral 
  if(tacInterpolate(&tac, &tac, NULL, &tac1, NULL, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1);
    return(3);
  }
 
 
 
  /*
   *  Prepare the room for parameters.
   */
  if(verbose>1) {printf("initializing parameter data\n"); fflush(stdout);}
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 3, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1);
    return(4);
  }
  iftFree(&par.h); // remove stupid header info
  /* Copy titles & file names */
  {
    int i;
    char buf[256];
    time_t t=time(NULL);
    /* set program name */
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    /* set file names */
    iftPut(&par.h, "datafile", tacfile, 0, NULL);
    iftPut(&par.h, "inputfile", inpfile, 0, NULL);
    /* Study number */
    if(tacGetHeaderStudynr(&tac.h, buf, NULL)==TPCERROR_OK)
      iftPut(&par.h, "study_number", buf, 0, NULL);
    /* Model */
    iftPut(&par.h, "model", "Yokoi plot", 0, NULL);
    /* Set current time to results */
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    /* Set the parameter names and units */
    par.parNr=3;
    i=0; strcpy(par.n[i].name,"K1"); par.n[i].unit=UNIT_ML_PER_ML_MIN;
    i++; strcpy(par.n[i].name,"k2"); par.n[i].unit=UNIT_PER_MIN;
    i++; strcpy(par.n[i].name,"K1/k2");
  }
 
  /*
   *  Prepare the room for plots
   */
  MTAC plots; mtacInit(&plots);
  if(mtacAllocate(&plots, tac.tacNr)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate space for plots.\n");
    tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par);
    return(5);
  }
 
 
 
  /*
   *  Analyse each regional TTAC
   */
  if(verbose>1) {printf("\ncalculating one TTAC at a time...\n"); fflush(stdout);}
  int okNr=0;
  for(int ti=0; ti<tac.tacNr; ti++) {
    if(verbose>2 && tac.tacNr>1) {printf("Region %d %s\n", 1+ti, tac.c[ti].name); fflush(stdout);}
    /* Calculate plot x axis */
    double plotx[fitSampleNr];
    for(int i=0; i<fitSampleNr; i++) plotx[i]=tac1.c[ti].y[i]/input1.c[0].y[i];
    /* Calculate plot y axis */
    double ploty[fitSampleNr];
    for(int i=0; i<fitSampleNr; i++) ploty[i]=tac.c[ti].y[i]/input1.c[0].y[i];
    if(verbose>4) {
      printf("\tX\tY\tCt\tiCt\tiCp\tT\n");
      for(int i=0; i<fitSampleNr; i++) if(i>=istart && i<=istop)
        printf("\t%.3g\t%.3g\t%.3g\t%.3g\t%.3g\t%g\n", plotx[i], ploty[i], tac.c[ti].y[i], 
               tac1.c[ti].y[i], input1.c[0].y[i], tac.x[i]);
    }
    /* Fit regression line */
    double slope, ic, slope_sd, ic_sd, xic, xic_sd, y_sd;
    int n;
    if(save_stat==0)
      n=fitLine(plotx+istart, ploty+istart, fitNr, &slope, &ic);
    else
      n=fitLinePearson(plotx+istart, ploty+istart, fitNr,
                       &slope, &slope_sd, &ic, &ic_sd, &xic, &xic_sd, NULL, &y_sd);
    if(n==0 || isnan(slope) || isnan(ic)) {
      fprintf(stderr, "Error: cannot fit %s\n", tac.c[ti].name); fflush(stderr); 
      continue;
    }
    par.r[ti].p[0]=ic;
    par.r[ti].p[1]=-slope;
    par.r[ti].p[2]=-ic/slope;
    if(save_stat!=0) {
      par.r[ti].sd[0]=ic_sd;
      par.r[ti].sd[1]=slope_sd;
      par.r[ti].sd[2]=xic_sd;
    }
    par.r[ti].dataNr=n;
    par.r[ti].start=tac.x[istart];
    par.r[ti].end=tac.x[istop];
    par.r[ti].fitNr=2;
 
    /* Store the data for plots */
    if(tacAllocate(&plots.tac[plots.nr], fitSampleNr, 3)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot copy plot data.\n"); fflush(stderr); 
      continue;
    }
    plots.tac[plots.nr].sampleNr=fitSampleNr;
    plots.tac[plots.nr].tacNr=3;
    for(int i=0; i<fitSampleNr; i++)
      if(isfinite(plotx[i])) plots.tac[plots.nr].x[i]=plotx[i];
      else plots.tac[plots.nr].x[i]=nan("");
    strcpy(plots.tac[plots.nr].c[0].name, tac.c[ti].name); 
    strcpy(plots.tac[plots.nr].c[1].name, tac.c[ti].name); 
    for(int i=0; i<fitSampleNr; i++) {
      plots.tac[plots.nr].c[0].y[i]=plots.tac[plots.nr].c[1].y[i]=nan("");
      if(isfinite(ploty[i])) {
        if(i>=istart && i<=istop) plots.tac[plots.nr].c[0].y[i]=ploty[i];
        else plots.tac[plots.nr].c[1].y[i]=ploty[i];
      }
      /* Regression line point */
      plots.tac[plots.nr].c[2].y[i]=ic+slope*plotx[i];
    }
    plots.nr++;
    okNr++;
  } // next TAC
  if(okNr==0) {
    fprintf(stderr, "Error: cannot calculate Yokoi MTGA.\n");
    tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par); 
    mtacFree(&plots); return(8);
  }
  if(verbose>1) {printf("... done.\n"); fflush(stdout);}
 
  /*
   *  Print parameters on screen
   */
  if(verbose>0 && par.tacNr<80) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, &status);
 
 
  /*
   *  Save CM results
   */
  {
    if(verbose>1) {printf("writing %s\n", resfile); fflush(stdout);}
    int ret=TPCERROR_FAIL;
    par.format=parFormatFromExtension(resfile);
    if(verbose>2) printf("result file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
    FILE *fp; fp=fopen(resfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing parameters.\n");
    } else {
      ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
      fclose(fp);
      if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    }
    if(ret!=TPCERROR_OK) {
      tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par);
      mtacFree(&plots); return(11);
    }
    if(verbose>0) {printf("Results saved in %s.\n", resfile); fflush(stdout);}
  }
 
 
  /*
   *  Write plot data, if requested
   */
  if(plotfile[0]) {
    if(verbose>1) {printf("writing %s\n", plotfile); fflush(stdout);}
    FILE *fp; fp=fopen(plotfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing plot data.\n");
      tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par);
      mtacFree(&plots); return(12);
    }
    fprintf(fp, "Yokoi plot");
    char buf[256];
    if(tacGetHeaderStudynr(&par.h, buf, NULL)==TPCERROR_OK) fprintf(fp, "\t%s", buf);
    fprintf(fp, "\n");
    for(int ti=0; ti<plots.nr; ti++) if(plots.tac[ti].sampleNr>0) {
      fprintf(fp, "\n%s\n", plots.tac[ti].c[0].name);
      fprintf(fp, "X\tFitted\tIgnored\tRegression\n");
      for(int i=0; i<plots.tac[ti].sampleNr; i++) {
        if(isfinite(plots.tac[ti].x[i])) fprintf(fp, "%g", plots.tac[ti].x[i]);
        fprintf(fp, "\t");
        if(isfinite(plots.tac[ti].c[0].y[i])) fprintf(fp, "%g", plots.tac[ti].c[0].y[i]);
        fprintf(fp, "\t");
        if(isfinite(plots.tac[ti].c[1].y[i])) fprintf(fp, "%g", plots.tac[ti].c[1].y[i]);
        fprintf(fp, "\t");
        if(isfinite(plots.tac[ti].c[2].y[i])) fprintf(fp, "%g", plots.tac[ti].c[2].y[i]);
        fprintf(fp, "\n");
      }
    }
    fclose(fp);
    if(verbose>0) {printf("Plot data saved in %s.\n", plotfile); fflush(stdout);}
  }
 
 
  /*
   *  Write SVG plot, if requested
   */
  if(svgfile[0]) {
    if(verbose>1) {printf("writing %s\n", svgfile); fflush(stdout);}
    char mtitle[256]; strcpy(mtitle, "Yokoi plot");
    char buf[256];
    if(tacGetHeaderStudynr(&par.h, buf, NULL)==TPCERROR_OK) {
      strcat(mtitle, " "); strcat(mtitle, buf);
    }
    if(mtgaPlotSVG(&plots, mtitle, svgfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par);
      mtacFree(&plots);
      return(13);
    }
    if(verbose>0) {printf("Plot saved in %s.\n", svgfile); fflush(stdout);}
  }
 
 
  tacFree(&tac); tacFree(&input0); tacFree(&input1); tacFree(&tac1); parFree(&par);
  mtacFree(&plots);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/