liverpv.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 "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Input for liver models consists of two parts: portal vein and hepatic artery.",
  "This program generates portal vein TAC from arterial TAC using the model",
  "proposed by Munk et al. (2003).",
  "In a pig study, Winterdahl et al. (2011) measured beta values for",
  "[F-18]FDG (0.50), [O-15]H2O (2.17), [O-15]CO (0.10), [C-11]methylglucose",
  "(0.57), and [F-18]FDGal (0.82).",
  " ",
  "Usage: @P [Options] arteryfile beta outputfile",
  " ",
  "Options:",
  " -i=<Interval>",
  "     Sample time interval in convolution; by default the shortest interval",
  "     in the data, but restricted to range 0.0001*beta - 0.02*beta;",
  "     too long interval leads to bias.",
  " -si",
  "     Output is written with sample intervals used in convolution; by default",
  "     data is interpolated back to the sample times of the artery data.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Beta must be given in minutes, and optional interval in seconds.",
  " ",
  "Example:",
  "  @P ue65ap.bld 0.50 ue65pv.bld",
  " ",
  "References:",
  "1. Munk OL, Keiding S, Bass L. Impulse-response function of splanchnic",
  "   circulation with model-independent constraints: theory and experimental",
  "   validation. Am J Physiol Gastrointest Liver Physiol. 2003;285:G671-G680.",
  "2. Winterdahl M, Keiding S, Sorensen M, Mortensen FV, Alstrup AKO, Munk AL.",
  "   Tracer input for kinetic modelling of liver physiology determined without",
  "   sampling portal venous blood in pigs.",
  "   Eur J Nucl Med Mol Imaging 2011;38:263-270.",
  " ",
  "See also: liverinp, taccalc, taccat, b2plasma",
  " ",
  "Keywords: input, blood, TAC, modelling, simulation, liver",
  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   pfile[FILENAME_MAX], afile[FILENAME_MAX];
  double beta=nan("");
  double interval=nan("");
  int    output_sampling=0; // 0=artery, 1=convolution
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  pfile[0]=afile[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, "I=", 2)==0) {
      interval=atofVerified(cptr+2); if(interval>0.0) continue;
    } else if(strcasecmp(cptr, "SI")==0) {
      output_sampling=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);}
 
  /* Arguments */
  if(ai<argc) {strlcpy(afile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    if(atofCheck(argv[ai], &beta)!=0 || beta<=0.0) {
      fprintf(stderr, "Error: invalid beta value.\n");
      return(1);
    }
    ai++;
  }
  if(ai<argc) {strlcpy(pfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!pfile[0]) { 
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]); 
    printf("\n");
    printf("afile := %s\n", afile);
    printf("beta := %g\n", beta);
    printf("pfile := %s\n", pfile);
    printf("output_sampling := %d\n", output_sampling);
    if(!isnan(interval)) printf("interval := %g\n", interval);
    fflush(stdout);
  }
 
 
  /*
   *  Read TAC file
   */
 
  /* Read arterial data */
  if(verbose>0) printf("reading %s\n", afile);
  TAC atac; tacInit(&atac);
  if(tacRead(&atac, afile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(atac.format));
    printf("tacNr := %d\n", atac.tacNr);
    printf("sampleNr := %d\n", atac.sampleNr);
    printf("xunit := %s\n", unitName(atac.tunit));
    printf("yunit := %s\n", unitName(atac.cunit));
  }
  if(atac.tacNr>1) {
    fprintf(stderr, "Warning: only first TAC in arterial data is used.\n");
    atac.tacNr=1;
  }
  /* Check for missing sample times */
  if(tacXNaNs(&atac)>0) {
    fprintf(stderr, "Error: missing frame times.\n");
    tacFree(&atac); return(2);
  }
  /* Check for missing concentrations */
  if(tacYNaNs(&atac, -1)>0) {
    fprintf(stderr, "Error: missing concentrations.\n");
    tacFree(&atac); return(2);
  }
  if(atac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples in plasma data.\n");
    tacFree(&atac); return(2);
  }
 
  /* Sort by sample times */
  if(tacSortByTime(&atac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&atac); return(2);
  }
  if(atac.isframe!=0 && verbose>0) {
    if(tacSetX(&atac, &status)!=TPCERROR_OK) { // make sure that frame middle times are set
      fprintf(stderr, "Error: invalid sample times.\n");
      tacFree(&atac); return(2);
    }
    fprintf(stderr, "Warning: frame durations are ignored.\n");
  }
 
 
 
  /* Beta is in minutes; convert it to seconds if necessary */
  /* Interval is in sec; convert to minutes if necessary */
  if(atac.tunit==UNIT_SEC) {
    beta*=60.0;
    if(verbose>1) printf("beta is converted to %g s.\n", beta);
  } else if(atac.tunit==UNIT_MIN) {
    if(!isnan(interval) && interval>0.0) {
      interval/=60.0;
      if(verbose>1) printf("interval is converted to %g min.\n", interval);
    }
  } else {
    fprintf(stderr, "Warning: time units are assumed to be in minutes.\n");
    if(!isnan(interval) && interval>0.0) {
      interval/=60.0;
      if(verbose>1) printf("interval is converted to %g min.\n", interval);
    }
  }
 
  /* Check user-defined interval time */
  if(!isnan(interval) && interval>0.0) {
    if(interval>0.01*atac.x[atac.sampleNr-1] || interval>0.2*beta) {
      fprintf(stderr, "Error: too long interval time.\n");
      tacFree(&atac); return(2);
    }
    if(interval>0.001*atac.x[atac.sampleNr-1] || interval>0.05*beta) {
      if(verbose>0) fprintf(stderr, "Warning: interval time may be too long.\n");
    }
  }
 
  /* Calculate kernel function integral at the end of data */
  if(verbose>2) {
    double kernel_integral=atac.x[atac.sampleNr-1]/(atac.x[atac.sampleNr-1]+beta);
    if(verbose>1) printf("kernel integral at end of data: %g\n", kernel_integral);
  }
  /* Calculate the time when 99.9% of kernel function integral is reached */
  if(verbose>2) {
    double kernel_time=999.*beta;
    double kernel_integral=kernel_time/(kernel_time+beta);
    if(verbose>1) printf("kernel integral at %g: %g\n", kernel_time, kernel_integral);
  }
  
  /* Interpolate data with even sample intervals for convolution */
  TAC itac; tacInit(&itac);
  int ret=0;
  if(!isnan(interval)) {
    /* user-defined interpolation sample frequency */
    ret=tacInterpolateToEqualLengthFrames(&atac, interval, interval, &itac, &status);
  } else {
    /* automatic interpolation sample frequency */
    ret=tacInterpolateToEqualLengthFrames(&atac, 0.0001*beta, 0.02*beta, &itac, &status);
  }
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&itac); return(3);
  }
  /* Get the sample interval in interpolated data */
  double freq=itac.x2[0]-itac.x1[0]; 
  if(verbose>1) {
    printf("sample_intervals_in_convolution := %g\n", freq);
    printf("interpolated data range: %g - %g\n", itac.x1[0], itac.x2[itac.sampleNr-1]);
  }
 
 
  /*
   *  Calculate the response function for convolution.
   *
   *  Function is f(x)=beta/(beta+x)^2 ; 
   *  its indefinite integral is Constant - beta/(beta+x),
   *  definite integral from 0 to +inf is 1,
   *  definite integral from 0 to a is a/(a+beta),
   *  and definite integral from t to t+dt is beta*dt/((beta+t-dt/2)*(beta+t+dt/2)).
   */
  if(verbose>0) fprintf(stdout, "allocate memory for kernel...\n");
  double *kernel=(double*)malloc(2*itac.sampleNr*sizeof(double));
  if(kernel==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    tacFree(&atac); tacFree(&itac); return(4);
  }
  double *cy=kernel+itac.sampleNr;
 
  if(verbose>0) fprintf(stdout, "computing the kernel...\n");
  double ksum=0.0;
  if(verbose>6) printf("\nData\tKernel:\n");
  for(int i=0; i<itac.sampleNr; i++) {
    kernel[i]=beta*freq/((beta+itac.x1[i])*(beta+itac.x2[i])); // integral, not mean
    ksum+=kernel[i]; // sum of stepwise integrals
    if(verbose>6) printf("%g\t%g\n", itac.c[0].y[i], kernel[i]);
  }
  if(verbose>4) printf("last_kernel_value := %g\n", kernel[itac.sampleNr-1]);
  if(verbose>2) printf("Sum of response function := %g\n", ksum);
  if(!isnormal(ksum)) {
    fprintf(stderr, "Error: invalid kernel contents.\n");
    tacFree(&atac); tacFree(&itac); free(kernel); return(5);
  }
 
 
  /*
   *  Convolution
   */
  if(verbose>0) {fprintf(stdout, "convolution...\n"); fflush(stdout);}
  ret=convolve1D(itac.c[0].y, itac.sampleNr, kernel, itac.sampleNr, cy);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot convolve the data.\n");
    tacFree(&atac); tacFree(&itac); free(kernel); return(6);
  }
  if(verbose>6) {
    printf("\nData x\ty\tKernel\tConvolved\n");
    for(int i=0; i<itac.sampleNr; i++)
      printf("%g\t%g\t%g\t%g\n", itac.x[i], itac.c[0].y[i], kernel[i], cy[i]);
  }
 
  /* Copy convoluted curve over interpolated curve */
  for(int i=0; i<itac.sampleNr; i++) itac.c[0].y[i]=cy[i];
  /* No need for kernel or working space */
  free(kernel);
 
  /*
   *  If user wants to save data with sample intervals used in convolution, then do that and exit
   */
  if(output_sampling==1) {
 
    /* No need for arterial data */
    tacFree(&atac);
 
    /* Change TAC name */
    if( strcasestr(itac.c[0].name, "art")!=NULL || strcasestr(itac.c[0].name, "aort")!=NULL || 
        strcasestr(itac.c[0].name, "LV")!=NULL || strcasestr(itac.c[0].name, "cav")!=NULL )
    {
      strcpy(itac.c[0].name, "PV");
    } else if(strcasestr(itac.c[0].name, "CA")!=NULL) {
      strcpy(itac.c[0].name, "Cpv");
    }
 
    if(verbose>1) printf("writing portal vein data in %s\n", pfile);
    FILE *fp; fp=fopen(pfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", pfile);
      tacFree(&itac); return(11);
    }
    ret=tacWrite(&itac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp); tacFree(&itac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>=0) printf("Portal vein TAC saved in %s\n", pfile);
    return(0);
  }
 
 
  /*
   *  Interpolate convolved data into original sample times
   */
  if(verbose>1) fprintf(stdout, "interpolating to original sample times...\n");
  /* Transform sample interval data into dot-to-dot data */
  TAC dtac; tacInit(&dtac);
  ret=tacFramesToSteps(&itac, &dtac, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&itac); return(10);
  }
  /* Interpolate to original times */
  if(atac.isframe==0)
    ret=liInterpolate(dtac.x, dtac.c[0].y, dtac.sampleNr, atac.x, atac.c[0].y, NULL, NULL, 
                      atac.sampleNr, 0, 0, verbose-10);
  else
    ret=liInterpolateForPET(dtac.x, dtac.c[0].y, dtac.sampleNr, atac.x1, atac.x2, atac.c[0].y, 
                            NULL, NULL, atac.sampleNr, 0, 1, verbose-10);
  tacFree(&dtac);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot interpolate back.\n");
    tacFree(&atac); tacFree(&itac); return(7);
  }
 
  /* No need for interpolated data */
  tacFree(&itac);
 
  /* Change TAC name */
  if( strcasestr(atac.c[0].name, "art")!=NULL || strcasestr(atac.c[0].name, "aort")!=NULL || 
      strcasestr(atac.c[0].name, "LV")!=NULL || strcasestr(atac.c[0].name, "cav")!=NULL )
  {
    strcpy(atac.c[0].name, "PV");
  } else if(strcasestr(atac.c[0].name, "CA")!=NULL) {
    strcpy(atac.c[0].name, "Cpv");
  }
 
  /*
   *  Write the file
   */
  if(verbose>1) printf("writing portal vein data in %s\n", pfile);
  FILE *fp; fp=fopen(pfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", pfile);
    tacFree(&atac); return(11);
  }
  ret=tacWrite(&atac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp); tacFree(&atac);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(12);
  }
  if(verbose>=0) printf("Portal vein TAC saved in %s\n", pfile);
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/