perfrat.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 "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculates perfusion ratio from regional radiowater PET data.",
  "In animal disease models, such as UUO model, one of kidneys is affected",
  "by the disease model and the other kidney serves as reference organ.",
  "Perfusion ratio can be calculated from AUCs of tissue curves without",
  "input function (1, 2). AUCs are limited to the time of the first peak of",
  "the control tissue curve.",
  " ",
  "Usage: @P [options] tacfile reference [outputfile]",
  " ",
  "Options:",
  " -etime=<time>",
  "     Extend or shorten AUC calculation past the peak by specified time",
  "     (sec). Positive time extends, negative time shortens the time range.",
  " -fr",
  "     Force given single reference region for every other TAC.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "As an example, if TAC file rabbit34.tac contains tissue regions with names",
  "UUO cortex, UUO medulla, CTRL cortex, CTRL medulla,",
  "then to calculate AUC ratio between UUO and CTRL cortex, and",
  "between UUO and CTRL medulla, enter command",
  "  @P rabbit34.tac CTRL rabbit34.par",
  "Program calculates AUCs from 0 to the peak time of the individual",
  "reference TAC, and saves AUC ratios in parameter file.",
  " ",
  "References:",
  "1. Xia et al. Hypertension 2008;51(2):466-473.",
  "2. Gulaldi et al. Biomed Res Int. 2013;835859.",
  " ",
  "See also: tacpeak, interpol, taccalc, fit_h2o, dftratio",
  " ",
  "Keywords: TAC, modelling, perfusion, ratio",
  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;
  char tacfile[FILENAME_MAX], outfile[FILENAME_MAX], refname[256];
  double extra_time=nan("");
  int force_ref=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=outfile[0]=refname[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, "DRY")==0) {
      continue;
    } else if(strcasecmp(cptr, "FR")==0) {
      force_ref=1; continue;
    } else if(strncasecmp(cptr, "ETIME=", 6)==0) {
      extra_time=atofVerified(cptr+6); if(!isnan(extra_time)) 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-5;
  
  /* 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(refname, argv[ai++], 256);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!refname[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(strcasecmp(tacfile, outfile)==0) {
    fprintf(stderr, "Error: input file would be overwritten.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("refname := %s\n", refname);
    if(!isnan(extra_time)) printf("extra_time := %g\n", extra_time);
    if(outfile[0]) printf("outfile := %s\n", outfile);
    printf("force_ref := %d\n", force_ref);
    fflush(stdout);
  }
 
 
  /*
   *  Read the file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  ret=tacRead(&tac, tacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, 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("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
  }
  if(verbose>3) iftWrite(&tac.h, stdout, NULL);
  if(tac.tacNr<2 || tac.sampleNr<3) {
    fprintf(stderr, "Error: no data to calculate AUC ratio.\n");
    tacFree(&tac); return(2);
  }
  if(!tacIsX(&tac)) {
    fprintf(stderr, "Error: sample times not available.\n");
    tacFree(&tac); return(2);
  }
  if(tacNaNs(&tac)) {
    fprintf(stderr, "Error: missing values in the TAC file.\n");
    tacFree(&tac); return(2);
  }
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  double x1, x2;
  if(tacSampleXRange(&tac, &x1, &x2)) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) printf("time range: %g - %g\n", x1, x2);
 
  /* If sample times are in minutes, then change extra_time to minutes, too */
  if(!isnan(extra_time)) {
    if(tac.tunit==UNIT_MIN) {
      extra_time/=60.0;
      if(verbose>0) printf("extra_time converted to minutes.\n");
    } else if(tac.tunit!=UNIT_SEC) {
      fprintf(stderr, "Warning: sample time assumed to be in seconds.\n");
    }
  }
 
  /* How many matches for reference name we get? */
  int refNr=tacSelectTACs(&tac, refname, 1, &status);
  if(refNr<=0) {
    fprintf(stderr, "Error: specified reference TAC not found.\n");
    tacFree(&tac); return(3);
  }
  if(refNr==tac.tacNr) {
    fprintf(stderr, "Error: all regions match the reference.\n");
    tacFree(&tac); return(3);
  }
  int refindex=-1;
  if(force_ref==0) {
    if(verbose>2 || (verbose>0 && refNr!=1)) {
      printf("%d tac(s) match name '%s'\n", refNr, refname); fflush(stdout);
    }
  } else {
    if(refNr>1 && verbose>0) {
      printf("%d tac(s) match name '%s'\n", refNr, refname); fflush(stdout);
    }
    refindex=tacSelectBestReference(&tac);
    if(refindex<0) {
      fprintf(stderr, "Error: cannot select best match for '%s'.\n", refname);
      tacFree(&tac); return(3);
    }
  }
 
 
  /*
   *  Pre-calculate AUC from start to each frame mid time for each TAC
   */
  if(verbose>1) {printf("integrating TACs\n"); fflush(stdout);}
  TAC auc; tacInit(&auc);
  if(tacInterpolate(&tac, &tac, NULL, &auc, NULL, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(4);
  }
 
 
  /*
   *  Create a new parameter structure for results.
   */
  if(verbose>2) {printf("allocating place for results\n"); fflush(stdout);}
  PAR par; parInit(&par);
  if(parAllocate(&par, 1, tac.tacNr)) {
    fprintf(stderr, "Error: cannot allocate memory for results.\n");
    parFree(&par); tacFree(&tac); tacFree(&auc); return(5);
  }
  par.tacNr=0; par.parNr=1;
  /* Set results information */
  par.format=parFormatFromExtension(outfile); 
  if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_TSV_UK;
  strcpy(par.n[0].name, "PR"); par.n[0].unit=UNIT_UNITLESS;
  {
    char buf[256];
    tpcProgramName(argv[0], 1, 1, buf, 256);
    iftPut(&par.h, "program", buf, 0, NULL);
    time_t t=time(NULL);
    iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
    char studynr[MAX_STUDYNR_LEN+1];
    if(tacGetHeaderStudynr(&tac.h, studynr, NULL)==TPCERROR_OK)
      studynrFromFilename(tacfile, studynr, 1);
    parSetStudyNr(&par, studynr);
  }
  iftPut(&par.h, "datafile", tacfile, 0, NULL);
  iftPut(&par.h, "reference_name", refname, 0, NULL);
 
 
  /*
   *  Calculate perfusion ratios
   */
  if(verbose>1) {printf("computing...\n"); fflush(stdout);}
  if(force_ref!=0) {
    /* Using the one selected reference for all other TACs */
    if(verbose>0) {printf("reference '%s'\n", tac.c[refindex].name); fflush(stdout);}
    /* Get the ref peak time */
    int maxi;
    if(tacYRange(&tac, refindex, NULL, NULL, NULL, &maxi, NULL, NULL)!=0) {
      fprintf(stderr, "Error: no peak found for %s\n", tac.c[refindex].name);
      tacFree(&tac); tacFree(&auc); parFree(&par); return(8);
    }
    if(verbose>3) {
      printf("  peak at time %g, sample %d, AUC %g\n", tac.x[maxi], 1+maxi, auc.c[refindex].y[maxi]);
      fflush(stdout);
    }
    for(int rj=0; rj<tac.tacNr; rj++) if(rj!=refindex) {
      /* Set result region name */
      strcpy(par.r[par.tacNr].name, tac.c[rj].name);
      /* Set result time range */
      par.r[par.tacNr].start=tac.x[0]; par.r[par.tacNr].end=tac.x[maxi];
      if(par.r[par.tacNr].start>0.0) par.r[par.tacNr].start=0.0;
      /* Calculate AUC ratio */
      if(isnan(extra_time) || fabs(extra_time)<1.0E-06) {
        /* Take AUCs directly from pre-calculated AUC-curves */
        if(verbose>4) {printf("      AUC %g\n", auc.c[rj].y[maxi]); fflush(stdout);}
        par.r[par.tacNr].p[0]=auc.c[rj].y[maxi]/auc.c[refindex].y[maxi];
      } else {
        /* Calculate AUCs from start to peak+extra_time */
        double et=tac.x[maxi]+extra_time;
        if(et>x2) et=x2; else if(et<x1) et=x1;
        par.r[par.tacNr].end=et; // refine result time range
        double newx[1], newy[1], auc1, auc2;
        newx[0]=et;
        ret=liInterpolate(auc.x, auc.c[refindex].y, auc.sampleNr, newx, newy, NULL, NULL, 1, 3, 1, 0);
        auc1=newy[0];
        if(!ret)
          ret=liInterpolate(auc.x, auc.c[rj].y, auc.sampleNr, newx, newy, NULL, NULL, 1, 3, 1, 0);
        auc2=newy[0];
        if(ret || isnan(auc1) || isnan(auc2)) {
          fprintf(stderr, "Error: cannot interpolate.\n");
          tacFree(&tac); tacFree(&auc); parFree(&par); return(8);
        }
        if(verbose>4) {printf("      AUCs %g/%g\n", auc2, auc1); fflush(stdout);}
        par.r[par.tacNr].p[0]=auc2/auc1;
        if(!isfinite(par.r[par.tacNr].p[0])) par.r[par.tacNr].p[0]=0.0;
      }
      /* Convert result time range into minutes, if necessary */
      if(tac.tunit==UNIT_SEC) {par.r[par.tacNr].start/=60.0; par.r[par.tacNr].end/=60.0;}
      /* prepare for next one */
      par.tacNr++; if(par.tacNr>=par._tacNr) break;
    }
 
  } else for(int ri=0; ri<tac.tacNr; ri++) if(tac.c[ri].sw==1) {
    if(verbose>2) {printf("reference '%s'\n", tac.c[ri].name); fflush(stdout);}
    /* Get the peak time */
    int maxi;
    if(tacYRange(&tac, ri, NULL, NULL, NULL, &maxi, NULL, NULL)!=0) {
      if(verbose>0) fprintf(stderr, "Error: no peak found for %s\n", tac.c[ri].name);
      continue;
    }
    if(verbose>3) {
      printf("  peak at time %g, sample %d, AUC %g\n", tac.x[maxi], 1+maxi, auc.c[ri].y[maxi]);
      fflush(stdout);
    }
    /* Extract the parts of region name before and after the given id */
    char part1[FILENAME_MAX], part2[FILENAME_MAX];
    strcpy(part1, tac.c[ri].name); char *cptr=strdelstr(part1, refname); 
    strcpy(part2, cptr); *cptr=(char)0;
    /* Try to find regions with matching name */
    if(verbose>3) {printf("  searching matches for '%s'...'%s'\n", part1, part2); fflush(stdout);}
    int matchNr=0;
    for(int rj=0; rj<tac.tacNr; rj++) if(tac.c[rj].sw==0) { // ... excluding references
      if(verbose>4) {printf("    checking '%s'\n", tac.c[rj].name); fflush(stdout);}
      int m=strlen(tac.c[rj].name); if(m<1) continue;
      int n=strlen(part1); if(m<n) continue;
      if(n>0 && strncasecmp(tac.c[rj].name, part1, n)) continue;
      n=strlen(part2); if(m<n) continue;
      if(n>0 && strcasecmp(tac.c[rj].name+m-n, part2)) continue;
      if(verbose>4) {printf("      match\n"); fflush(stdout);}
      /* Set result region name */
      char *p1=strdup(part1); strCleanSpaces(p1);
      char *p2=strdup(part2); strCleanSpaces(p2);
      if(strlen(p1)>0) strcpy(par.r[par.tacNr].name, p1);
      if(strlen(p1)>0 && strlen(p2)>0) strcat(par.r[par.tacNr].name, " ");
      if(strlen(p2)>0) strcat(par.r[par.tacNr].name, p2);
      free(p1); free(p2);
      if(par.r[par.tacNr].name[0]=='-' || par.r[par.tacNr].name[0]=='_')
        strTrimLeft(par.r[par.tacNr].name, 1);
      if(strlen(par.r[par.tacNr].name)<1) sprintf(par.r[par.tacNr].name, "tac%d", 1+par.tacNr);
      /* Set result time range */
      par.r[par.tacNr].start=tac.x[0]; par.r[par.tacNr].end=tac.x[maxi];
      if(par.r[par.tacNr].start>0.0) par.r[par.tacNr].start=0.0;
      /* Calculate AUC ratio */
      if(isnan(extra_time) || fabs(extra_time)<1.0E-06) {
        /* Take AUCs directly from pre-calculated AUC-curves */
        if(verbose>4) {printf("      AUC %g\n", auc.c[rj].y[maxi]); fflush(stdout);}
        par.r[par.tacNr].p[0]=auc.c[rj].y[maxi]/auc.c[ri].y[maxi];
      } else {
        /* Calculate AUCs from start to peak+extra_time */
        double et=tac.x[maxi]+extra_time;
        if(et>x2) et=x2; else if(et<x1) et=x1;
        par.r[par.tacNr].end=et; // refine result time range
        double newx[1], newy[1], auc1, auc2;
        newx[0]=et;
        ret=liInterpolate(auc.x, auc.c[ri].y, auc.sampleNr, newx, newy, NULL, NULL, 1, 3, 1, 0);
        auc1=newy[0];
        if(!ret)
          ret=liInterpolate(auc.x, auc.c[rj].y, auc.sampleNr, newx, newy, NULL, NULL, 1, 3, 1, 0);
        auc2=newy[0];
        if(ret || isnan(auc1) || isnan(auc2)) {
          fprintf(stderr, "Error: cannot interpolate.\n");
          tacFree(&tac); tacFree(&auc); parFree(&par); return(8);
        }
        if(verbose>4) {printf("      AUCs %g/%g\n", auc2, auc1); fflush(stdout);}
        par.r[par.tacNr].p[0]=auc2/auc1;
        if(!isfinite(par.r[par.tacNr].p[0])) par.r[par.tacNr].p[0]=0.0;
      }
      matchNr++;
      /* Convert result time range into minutes, if necessary */
      if(tac.tunit==UNIT_SEC) {par.r[par.tacNr].start/=60.0; par.r[par.tacNr].end/=60.0;}
      /* prepare for next one */
      par.tacNr++; if(par.tacNr>=par._tacNr) break;
    }
    if(verbose>0) {
      if(matchNr==0)
        fprintf(stderr, "Warning: no matches for '%s'\n", tac.c[ri].name);
      else if(matchNr>1) 
        fprintf(stderr, "Warning: %d matches for '%s'\n", matchNr, tac.c[ri].name);
      fflush(stdout);
    }
  }
 
 
  /* Original data and AUCs no more needed */
  tacFree(&tac); tacFree(&auc);
  /* Check that something could be calculated */
  if(par.tacNr<1) {
    fprintf(stderr, "Error: no data to calculate perfusion ratio.\n");
    parFree(&par); return(10);
  }
 
  /* Print and save the results */
  if(verbose>0 || !outfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, NULL);
  if(outfile[0]) {
    /* Save file */
    if(verbose>1) printf("  saving %s\n", outfile);
    FILE *fp=fopen(outfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      parFree(&par); return(12);
    }
    if(verbose>0) printf("perfusion ratios saved in %s\n", outfile);
  }
 
  parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/