lhtest.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"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Method testing.",
  " ",
  "Usage: @P [options] PTAC TTAC results",
  " ",
  "Options:",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted regional TTACs are written in specified file.",
  " -lp=<Filename>",
  "     Parameters of linear model are saved in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  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    ptacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], resfile[FILENAME_MAX],
          fitfile[FILENAME_MAX], svgfile[FILENAME_MAX], lpfile[FILENAME_MAX];
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ptacfile[0]=ttacfile[0]=resfile[0]=fitfile[0]=svgfile[0]=lpfile[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, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); if(strlen(fitfile)>0) continue;
    } else if(strncasecmp(cptr, "LP=", 3)==0) {
      strlcpy(lpfile, cptr+3, FILENAME_MAX); if(strlen(lpfile)>0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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(ptacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ttacfile, argv[ai++], FILENAME_MAX);
  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);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ptacfile := %s\n", ptacfile);
    printf("ttacfile := %s\n", ttacfile);
    printf("resfile := %s\n", resfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("lpfile := %s\n", lpfile);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>1) printf("reading tissue and input data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC ptac, ttac0; tacInit(&ptac); tacInit(&ttac0);
  double fitdur=1.0E+10;
  if(tacReadModelingData(ttacfile, ptacfile, NULL, NULL, &fitdur, 0,
                          NULL, &ttac0, &ptac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0); tacFree(&ptac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac0.format));
    printf("tacNr := %d\n", ttac0.tacNr);
    printf("ttac.sampleNr := %d\n", ttac0.sampleNr);
    printf("ptac.sampleNr := %d\n", ptac.sampleNr);
    printf("xunit := %s\n", unitName(ttac0.tunit));
    printf("yunit := %s\n", unitName(ttac0.cunit));
  }
  if(ttac0.isframe==0) {
    fprintf(stderr, "Error: TTAC file does not contain frame times.\n");
    tacFree(&ttac0); tacFree(&ptac); return(2);
  }
 
 
  /*
   *  Interpolate PTAC to TTAC frame times
   */
  if(verbose>1) printf("interpolating PTAC\n");
  TAC input0; tacInit(&input0); // PTAC mean during PET frames
  if(tacInterpolate(&ptac, &ttac0, &input0, NULL, NULL, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac0); tacFree(&ptac); 
    tacFree(&input0);
    return(3);
  }
 
  /*
   *  Integrate PTAC to TTAC frame end times
   */
  if(verbose>1) printf("integrating PTAC\n");
  TAC input1; tacInit(&input1); // 1st integral 
  tacDuplicate(&input0, &input1);
  if(liIntegrateFE(input0.x1, input0.x2, input0.c[0].y, input0.sampleNr, input1.c[0].y, NULL, 0)!=0) {
    fprintf(stderr, "Error: cannot integrate PTAC.\n");
    tacFree(&ttac0); tacFree(&ptac); 
    tacFree(&input0); tacFree(&input1);
    return(3);
  }
  /* Original PTAC should not be needed any more */
  tacFree(&ptac);
 
  /*
   *  Integrate TTAC to TTAC frame end times
   */
  if(verbose>1) printf("integrating TTAC\n");
  TAC ttac1; tacInit(&ttac1); // 1st integral 
  tacDuplicate(&ttac0, &ttac1);
  for(int i=0; i<ttac0.tacNr; i++) {
    if(liIntegrateFE(ttac0.x1, ttac0.x2, ttac0.c[i].y, ttac0.sampleNr, ttac1.c[i].y, NULL, 0)!=0) {
      fprintf(stderr, "Error: cannot integrate PTAC.\n");
      tacFree(&ttac0); tacFree(&ttac1); 
      tacFree(&input0); tacFree(&input1);
      return(3);
    }
  }
 
 
  /*
   *  Allocate memory required by NNLS
   */
  if(verbose>1) printf("allocating memory for NNLS\n");
  int llsq_m=ttac0.sampleNr;
  int llsq_n=3;
  double *llsq_mat=(double*)malloc((llsq_n*llsq_m)*sizeof(double));
  if(llsq_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    tacFree(&ttac0); tacFree(&ttac1);
    tacFree(&input0); tacFree(&input1);
    return(5);
  }
  double **llsq_a=(double**)malloc(llsq_n*sizeof(double*));
  if(llsq_a==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    tacFree(&ttac0); tacFree(&ttac1);
    tacFree(&input0); tacFree(&input1);
    free(llsq_mat);
    return(5);
  }
  for(int ni=0; ni<llsq_n; ni++) llsq_a[ni]=llsq_mat+ni*llsq_m;
  double r2, llsq_b[llsq_m], llsq_x[llsq_n], llsq_wp[llsq_n], llsq_zz[llsq_m];
  int indexp[llsq_n];
 
 
  /*
   *  Fit each regional TTAC
   */
  for(int ti=0; ti<ttac0.tacNr; ti++) {
 
    if(verbose>1 && ttac0.tacNr>1) {
      printf("Region %d %s\n", 1+ti, ttac0.c[ti].name); fflush(stdout);}
 
    /* Setup data matrix A and vector B */
    for(int mi=0; mi<llsq_m; mi++)
      llsq_b[mi]=ttac0.c[ti].y[mi];
    for(int mi=0; mi<llsq_m; mi++) {
      llsq_mat[mi]=input0.c[0].y[mi];           // Vb
      if(mi>0)                                  // K1+Vb*k2
        llsq_mat[mi+llsq_m]=0.5*(input1.c[0].y[mi]+input1.c[0].y[mi-1]);
      else
        llsq_mat[mi+llsq_m]=0.5*(input1.c[0].y[mi]+0.0);
      if(mi>0)                                  // k2
        llsq_mat[mi+2*llsq_m]=-0.5*(ttac1.c[0].y[mi]+ttac1.c[0].y[mi-1]);
      else
        llsq_mat[mi+2*llsq_m]=-0.5*(ttac1.c[0].y[mi]+0.0);
      }
      if(verbose>5) {
        printf("Matrix A and vector B:\n");
        for(int mi=0; mi<llsq_m; mi++) {
          printf("%.2e", llsq_a[0][mi]);
          for(int ni=1; ni<llsq_n; ni++) printf(", %.2e", llsq_a[ni][mi]);
          printf("; %.3e\n", llsq_b[mi]);
        }
      }
 
      /* Compute NNLS */
      if(verbose>3) printf("starting NNLS...\n");
      int ret=nnls(llsq_a, llsq_m, llsq_n, llsq_b, llsq_x, &r2, llsq_wp, llsq_zz, indexp);
      if(verbose>3) printf("  ... done.\n");
      if(ret>1) {
        fprintf(stderr, "Warning: no NNLS solution for %s\n", ttac0.c[ti].name);
        for(int ni=0; ni<llsq_n; ni++) llsq_x[ni]=0.0;
        r2=0.0;
      } else if(ret==1) {
        fprintf(stderr, "Warning: NNLS iteration max exceeded for %s\n", ttac0.c[ti].name);
      }
      if(verbose>4) {
        printf("solution_vector: %g", llsq_wp[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_wp[ni]);
        printf("\n");
      }
      if(verbose>3) {
        printf("parameter_vector: %g", llsq_x[0]);
        for(int ni=1; ni<llsq_n; ni++) printf(", %g", llsq_x[ni]);
        printf("\n");
      }
 
    } // next TAC
 
    /* Free allocated memory */
    free(llsq_a); free(llsq_mat);
 
 
 
  tacFree(&ttac0); tacFree(&ttac1);
  tacFree(&input0); tacFree(&input1);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/