imglkup.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
//#include <float.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
//#include "libtpcmodext.h"
/*****************************************************************************/
#define LKUP_HIST_NR 25
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Replace the voxel values in PET image or scan file with the values from",
  "a look-up table.",
  "PET image file must be in ECAT, NIfTI, or Analyze format.",
  "The look-up table must contain two columns: program looks from the first",
  "column a matching value for the voxel value, and replaces the voxel value",
  "with the value from the second column of the table.",
  "Look-up table must be sorted in ascending order.",
  " ",
  "Usage: @P [Options] petfile lkupfile outputfile",
  " ",
  "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.",
  " -dyn[amic]",
  "     By default, dynamic PET data is considered as an error, because",
  "     this program is mostly used in autoradiography method with static or",
  "     or integrated PET data; this option enables usage of dynamic data.",
  " -his=<Filename>",
  "     Histogram of the results is written in TAC format in specified file.",
/*
  " -u=<unit id>",
  "     Set image unit to a specified value: 0=unknown, 1=cnts/sec,",
  "     2=counts, 3=kBq/mL, 4=sec*kBq/mL, 5=1/sec, 6=1/min, 7=mL/mL,",
  "     8=mL/dL, 9=mL/(mL*min), 10=mL/(dL*min), 11=unitless, 12=nCi/mL,",
  "     13=MBq/mL, 14=Bq/cc, 15=uCi/cc, 16=umol/(min*100g), 17=mg/(min*100g)",
*/
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: arlkup, imginteg, taclkup, imgcalc, imgunit, img2dft",
  " ",
  "Keywords: image, look-up table, autoradiography, ARG, perfusion, 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;
  int      ret=0;
  char     petfile[FILENAME_MAX], lkupfile[FILENAME_MAX], outfile[FILENAME_MAX];
  char    *cptr, hisfile[FILENAME_MAX];
  int      takeClosest=0;
  int      acceptDynamic=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=outfile[0]=lkupfile[0]=hisfile[0]=(char)0;
  /* Get 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==NULL) continue;
    if(strncasecmp(cptr, "CLOSEST", 1)==0) {
      takeClosest=1; continue;
    } else if(strncasecmp(cptr, "DYNAMIC", 3)==0) {
      acceptDynamic=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;
 
  /* 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(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(lkupfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, 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(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(strcasecmp(petfile, outfile)==0) {
    fprintf(stderr, "Error: same name for input and output file.\n"); 
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("lkupfile := %s\n", lkupfile);
    printf("outfile := %s\n", outfile);
    if(hisfile[0]) printf("hisfile := %s\n", hisfile);
    printf("takeClosest := %d\n", takeClosest);
    printf("acceptDynamic := %d\n", acceptDynamic);
  }
 
 
 
 
  /*
   *  Read the look-up table
   */
  if(verbose>0) printf("reading %s\n", lkupfile);
  DFT table; dftInit(&table);
  if(dftRead(lkupfile, &table)) {
    fprintf(stderr, "Error: %s\n", dfterrmsg);
    return(2);
  }
  if(verbose>20) dftPrint(&table);
  int tableSize=table.frameNr; 
  double *table1=table.x; double *table2=table.voi[0].y;
  if(table.voiNr>1 || tableSize<2) {
    fprintf(stderr, "Error: invalid look-up table %s\n", lkupfile);
    dftEmpty(&table); return(2);
  }
  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 PET file
   */
  if(verbose>0) fprintf(stdout, "reading %s\n", petfile);
  IMG img; imgInit(&img);
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret=%d\n", ret);
    dftEmpty(&table); imgEmpty(&img); return(3);
  }
  if(verbose>10) imgInfo(&img);
  /* Check the frame number */
  if(acceptDynamic==0 && img.dimt>1) {
    fprintf(stderr, "Error: file contains more than one frame.\n");
    dftEmpty(&table); imgEmpty(&img); return(3);
  }
 
 
  /*
   *  Prepare histogram data
   */
  if(verbose>1) printf("allocating memory for the histogram\n");
  DFT hist; dftInit(&hist);
  ret=dftSetmem(&hist, LKUP_HIST_NR, 1);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot allocate memory for histogram.\n");
    dftEmpty(&table); imgEmpty(&img); return(4);
  }
  hist.voiNr=1; hist.frameNr=LKUP_HIST_NR;
  hist._type=DFT_FORMAT_PMOD; hist.timetype=DFT_TIME_STARTEND;
  strcpy(hist.voi[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(int i=1; i<hist.frameNr; i++) {
    hist.x1[i]=hist.x2[i-1];
    hist.x2[i]=hist.x1[i]+histogram_bin;
  }
  for(int i=0; i<hist.frameNr; i++) hist.voi[0].y[i]=0.0;
 
 
 
 
  /*
   *  Replace voxel values with numbers from the look-up table
   */
  if(verbose>0) {fprintf(stdout, "using look-up table\n"); fflush(stdout);}
  double v;
  int low, high, mid, n;
  int zi, yi, xi, ti;
  for(zi=0; zi<img.dimz; zi++) {
    for(yi=0; yi<img.dimy; yi++) {
      for(xi=0; xi<img.dimx; xi++) {
        for(ti=0; ti<img.dimt; ti++) {
          v=(double)img.m[zi][yi][xi][ti]; //printf("%g ->\n", v);
          if(isnan(v)) continue;
          if(v<table1[0]) { /* check if below table values */
            v=table2[0];
            //hist.voi[0].y[0]+=1.0;
          } else if(v>table1[tableSize-1]) { /* check if over table values */
            v=table2[tableSize-1];
            //hist.voi[0].y[hist.frameNr-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.frameNr; n++) 
              if(v>=hist.x1[n] && v<=hist.x2[n]) {
              hist.voi[0].y[n]+=1.0; break;
            }
          }
          img.m[zi][yi][xi][ti]=(float)v; //printf("       %g\n", v);
        } // next frame
      } // next column
    } // next row
  } // next plane
  if(verbose>0 && img.dimz>2) {fprintf(stdout, "\n"); fflush(stdout);}
 
  /* Copy the unit */
  img.unit=petCunitId(table.unit);
 
  /* No need for look-up table anymore */
  dftEmpty(&table);
 
 
 
  /*
   *  Write the modified image
   */
  if(verbose>1) fprintf(stdout, "writing %s\n", outfile);
  if(imgWrite(outfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); 
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img);
    return(11);
  }
  if(verbose>0) printf("Written %s\n", outfile);
  imgEmpty(&img);
 
 
  /*
   *  Save or print the histogram, if required
   */
  if(!hisfile[0] && verbose<=1) {
    dftEmpty(&hist); return(0);
  }
  if(!hisfile[0]) {
    printf("\nHistogram:\n"); fflush(stdout);
    dftPrint(&hist);
    dftEmpty(&hist); return(0);
  }
  if(verbose>1) printf("writing %s\n", hisfile);
  if(dftWrite(&hist, hisfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", hisfile, dfterrmsg);
    dftEmpty(&hist); return(21);
  }
  if(verbose>0) fprintf(stdout, "Histogram written in %s\n", hisfile);
  dftEmpty(&hist);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/