taclkup.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 "tpctacmod.h"
/*****************************************************************************/
#define LKUP_HIST_NR 25
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Replace the y values in TAC file with the values from a given",
  "look-up table.",
  "The look-up table must contain two columns: program looks from the first",
  "column a matching value for the TAC y value, and replaces the y",
  "value with the value from the second column of the table.",
  "Look-up table must be sorted in ascending order.",
  " ",
  "Usage: @P [options] TAC lkupfile outputfile",
  " ",
  "Output file is written in TAC file format, unless filename",
  "extension is *.res or *.par.",
  " ",
  "Options:",
  " -c[losest]",
  "     If exact match in look-up table is not found, the closest value",
  "     is selected; by default, value is interpolated from the table.",
  " -his=<Filename>",
  "     Histogram of the results is written in TAC format in specified file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: arlkup, imglkup, dftinteg, taccalc, tacunit, tacsety",
  " ",
  "Keywords: look-up table, autoradiography, ARG, perfusion, blood flow, radiowater",
  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], lkupfile[FILENAME_MAX], resfile[FILENAME_MAX],
          hisfile[FILENAME_MAX];
  char   *cptr;
  int     takeClosest=0;
  int     ret, i, j, n;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=lkupfile[0]=resfile[0]=hisfile[0]=(char)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=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "CLOSEST", 1)==0) {
      takeClosest=1; continue;
    } else if(strncasecmp(cptr, "HIS=", 4)==0) {
      strlcpy(hisfile, cptr+4, FILENAME_MAX); if(strlen(hisfile)>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(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(lkupfile, 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("tacfile := %s\n", tacfile);
    printf("lkupfile := %s\n", lkupfile);
    printf("resfile := %s\n", resfile);
    if(hisfile[0]) printf("hisfile := %s\n", hisfile);
    printf("takeClosest := %d\n", takeClosest);
  }
 
 
 
  /*
   *  Read the look-up table
   */
  if(verbose>1) printf("reading %s\n", lkupfile);
  TAC lkup; tacInit(&lkup);
  ret=tacRead(&lkup, lkupfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&lkup); return(2);
  }
  if(verbose>2) {
    printf("lkup.fileformat := %s\n", tacFormattxt(lkup.format));
    printf("lkup.tacNr := %d\n", lkup.tacNr);
    printf("lkup.sampleNr := %d\n", lkup.sampleNr);
    printf("lkup.xunit := %s\n", unitName(lkup.tunit));
    printf("lkup.yunit := %s\n", unitName(lkup.cunit));
  }
  if(lkup.tacNr>1 || lkup.sampleNr<2) {
    fprintf(stderr, "Error: invalid look-up table %s\n", lkupfile);
    tacFree(&lkup); return(2);
  }
  /* Check for missing values */
  if(tacNaNs(&lkup)>0) {
    fprintf(stderr, "Error: missing look-up table values.\n");
    tacFree(&lkup); return(2);
  }
  /* Set simpler pointers to look-up data */
  int tableSize=lkup.sampleNr;
  double *table1=lkup.x;
  double *table2=lkup.c[0].y;
  if(verbose>1) {
    printf("Look-up y: %g - %g\n", table2[0], table2[tableSize-1]);
    printf("Look-up x: %g - %g\n", table1[0], table1[tableSize-1]);
  }
 
 
  /*
   *  Read the TAC 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: %s\n", errorMsg(status.error));
    tacFree(&lkup); tacFree(&tac); return(3);
  }
  if(verbose>2) {
    printf("tac.fileformat := %s\n", tacFormattxt(tac.format));
    printf("tac.tacNr := %d\n", tac.tacNr);
    printf("tac.sampleNr := %d\n", tac.sampleNr);
    printf("tac.xunit := %s\n", unitName(tac.tunit));
    printf("tac.yunit := %s\n", unitName(tac.cunit));
  }
 
 
  /*
   *  Prepare histogram data
   */
  if(verbose>1) printf("allocating memory for the histogram\n");
  TAC hist; tacInit(&hist);
  ret=tacAllocate(&hist, LKUP_HIST_NR, 1);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory for histogram.\n");
    tacFree(&lkup); tacFree(&tac); return(4);
  }
  hist.tacNr=1; hist.sampleNr=LKUP_HIST_NR;
  hist.format=TAC_FORMAT_PMOD; hist.isframe=1;
  hist.tunit=lkup.cunit; hist.cunit=UNIT_UNITLESS;
  strcpy(hist.c[0].name, "n");
  /* Calculate histogram bins and initiate the histogram values */
  double histogram_bin=(table2[tableSize-1]-table2[0])/(double)LKUP_HIST_NR;
  if(verbose>2) printf("histogram_bin_size := %g\n", histogram_bin);
  hist.x1[0]=table2[0]; hist.x2[0]=hist.x1[0]+histogram_bin;
  for(i=1; i<hist.sampleNr; i++) {
    hist.x1[i]=hist.x2[i-1];
    hist.x2[i]=hist.x1[i]+histogram_bin;
  }
  for(i=0; i<hist.sampleNr; i++) hist.c[0].y[i]=0.0;
 
 
  /*
   *  Replace TAC y values with numbers from the look-up table
   */
  if(verbose>1) fprintf(stdout, "using look-up table\n");
  double v;
  int low, high, mid;
  for(i=0; i<tac.tacNr; i++) for(j=0; j<tac.sampleNr; j++) {
    v=tac.c[i].y[j]; //printf("%g ->\n", v);
    if(isnan(v)) continue;
    if(v<table1[0]) { /* check if below table values */
      v=table2[0];
      //hist.c[0].y[0]+=1.0;
    } else if(v>table1[tableSize-1]) { /* check if over table values */
      v=table2[tableSize-1];
      //hist.c[0].y[hist.sampleNr-1]+=1.0;
    } else { /* binary search */
      low=mid=0; high=tableSize-1; ret=1;
      while(low<=high) {
        mid=(low+high)/2;
        if(v<table1[mid]) high=mid-1;
        else if(v>table1[mid]) low=mid+1;
        else {v=table2[mid]; ret=0; break;}
      }
      /* If not found, get the closest/interpolated value */
      if(ret) {
        n=high; high=low; low=n;
        if(takeClosest) {
          if(v-table1[low]<table1[high]-v) v=table2[low];
          else v=table2[high];
        } else {
          v=table2[low]+
            (v-table1[low])*(table2[high]-table2[low])
            /(table1[high]-table1[low]);
        }
      }
      /* Add to histogram data */
      for(n=0; n<hist.sampleNr; n++) 
        if(v>=hist.x1[n] && v<=hist.x2[n]) {
          hist.c[0].y[n]+=1.0; break;
        }
    }
    tac.c[i].y[j]=v; //printf("       %g\n", v);
  } /* next y */
 
  /* Copy the unit */
  tac.cunit=lkup.cunit;
 
  /* No need for look-up table anymore */
  tacFree(&lkup);
 
 
  /*
   *  Save the modified file
   */
  if(verbose>1) printf("writing %s\n", resfile);
  /* If filename extension is .res or .par, then save in RES/PAR format,
     otherwise in original TAC format */
  int format=0;
  cptr=strrchr(resfile, '.'); if(cptr!=NULL) cptr++;
  if(strcasecmp(cptr, "RES")==0) format=PAR_FORMAT_RES;
  else if(strcasecmp(cptr, "PAR")==0) format=PAR_FORMAT_CSV_UK;
  else if(strncasecmp(cptr, "HTML", 3)==0) format=PAR_FORMAT_HTML;
  /* Open file for writing */
  FILE *fp; fp=fopen(resfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", resfile);
    tacFree(&tac); tacFree(&hist); return(11);
  }
  if(format==0) {
    /* Write as TAC file */
    ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 0, &status);
  } else {
    /* Copy TAC data into a new parameter struct */
    PAR par; parInit(&par);
    ret=tacToPAR(&tac, &par, &status);
    if(ret==TPCERROR_OK)
      ret=parWrite(&par, fp, format, 1, &status);
    parFree(&par);
  }
  fclose(fp); tacFree(&tac);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&hist); return(12);
  }
  if(verbose>0) printf("Converted data saved in %s.\n", resfile);
 
 
  /*
   *  Save or print the histogram, if required
   */
  if(!hisfile[0] && verbose<=1) {
    tacFree(&hist); return(0);
  }
  if(hisfile[0]) {
    if(verbose>1) printf("writing %s\n", hisfile);
    fp=fopen(hisfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", hisfile);
      tacFree(&hist); return(21);
    }
  } else {
    fp=stdout; printf("\nHistogram:\n"); fflush(stdout);
  }
  /* Write histogram as TAC file */
  ret=tacWrite(&hist, fp, TAC_FORMAT_UNKNOWN, 0, &status);
  tacFree(&hist); if(hisfile[0]) fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(22);
  }
  if(verbose>0 && hisfile[0])  printf("Histogram saved in %s.\n", hisfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/