imgratio.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 "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpccurveio.h"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computes an AUC ratio image from dynamic PET image and reference region TTAC.",
  "AUCs are calculated between start and end time given in minutes.",
  "Image data can be in ECAT, Analyze 7.5, or NIfTI-1 format.",
  " ",
  "Usage: @P [Options] imgfile ttacfile starttime endtime outputfile",
  " ",
  "Options:",
  " -Bound",
  "     Instead of VOXEL/REF, (VOXEL-REF)/REF ratio is calculated.",
  " -Thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x ref AUC) are set to zero;",
  "     default is 10%.",
  " -max=<Max value>",
  "     Upper limit for output pixel values.",
  " -noneg",
  "     Negative pixel values are set to zero.", 
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P ua2918dy1.v ua2918cer.tac 60 90 ua2918ratio.v",
  " ",
  "See also: imginteg, imgcalc, dftratio, imgcbv, imgunit, tacunit",
  " ",
  "Keywords: image, modelling",
  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 ecat7CopyHeadersNoQuant(
  const char *file1,
  const char *file2,
  int verbose
) {
  int ii, mi, ret;
  char tmp[1024];
  ECAT_HEADERS ehdr;
 
  if(verbose>1) printf("ecatCopyHeadersNoQuant(%s, %s)\n", file1, file2);  
  if(file1==NULL || file2==NULL) return(1);
  
  
  /* Read headers */
  ehdrInitiate(&ehdr);
  ret=ecat7ReadHeaders(file1, &ehdr, verbose-1); if(ret!=0) return(10);
  /* Remove main header information that should NOT be copied */
  strcpy(tmp, "user_process_code");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "ecat_calibration_factor");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "calibration_units");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "calibration_units_label");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
  strcpy(tmp, "num_frames");
  ii=iftGet(&ehdr.mh, tmp, 0); iftDeleteItem(&ehdr.mh, ii, 0);
 
  /* Remove subheader information that should NOT be copied */
  for(mi=0; mi<ehdr.nr; mi++) {
    strcpy(tmp, "scale_factor");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "image_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_min");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "scan_max");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "frame_duration");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "frame_start_time");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
    strcpy(tmp, "decay_corr_fctr");
    ii=iftGet(&ehdr.m[mi].sh, tmp, 0); iftDeleteItem(&ehdr.m[mi].sh, ii, 0);
  }
  
  /* Write headers */
  ret=ecat7WriteHeaders(file2, &ehdr, verbose-1);
  ehdrEmpty(&ehdr);
  if(ret!=0) return(10+ret);
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  int     ret;
  int     b_per_r=0; // 0 or 1
  char    tacfile[FILENAME_MAX], petfile[FILENAME_MAX], ratfile[FILENAME_MAX];
  char    tmp[512];
  float   calcThreshold=0.10;
  float   upperLimit=-1.0;
  int     nonnegativity_constraint=0;
  DFT     tac;
  IMG     img;
  double  v, tstart, tstop, aucref;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=petfile[0]=ratfile[0]=(char)0;
  imgInit(&img); dftInit(&tac);
  tstart=tstop=nan("");
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "BOUND", 1)==0) {
      b_per_r=1; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v;
      if(atof_with_check(cptr+4, &v)==0) {
        calcThreshold=0.01*v; if(calcThreshold<=2.0) continue;
      }
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      double v; if(atof_with_check(cptr+4, &v)==0 && v>0.0) {upperLimit=(float)v; continue;}
    } else if(strncasecmp(cptr, "NONEG", 3)==0) {
      nonnegativity_constraint=1; continue;
    }
    fprintf(stderr, "Error: unknown 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); ai++;}
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc && atof_with_check(argv[ai], &v)==0) {tstart=v; ai++;}
  if(ai<argc && atof_with_check(argv[ai], &v)==0) {tstop=v; ai++;}
  if(ai<argc) {strlcpy(ratfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing? */
  if(!ratfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(tstop<tstart || tstart<0.0) {
    fprintf(stderr, "Error: invalid time range arguments.\n");
    return(1);
  }
  if(strcasecmp(petfile, ratfile)==0) {
    fprintf(stderr, "Error: check the output filename.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    printf("petfile := %s\n", petfile);
    printf("ratfile := %s\n", ratfile);
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("b_per_r := %d\n", b_per_r);
    printf("calcThreshold := %g%%\n", 100.*calcThreshold);
    printf("max := %f\n", upperLimit);
    printf("nonnegativity_constraint := %d\n", nonnegativity_constraint);
  }
 
 
  /*
   *  Read PET image and reference TTAC
   */
  {
    if(verbose>1) printf("reading data files\n");
    double endtime=-1.0;
    FILE *lfp; if(verbose>1) lfp=stdout; else lfp=NULL;
    ret=imgReadModelingData(
      petfile, NULL, tacfile, NULL, NULL, &endtime, &ai, &img,
      NULL, &tac, 0, lfp, verbose-2, tmp);
    if(ret!=0) {
      fprintf(stderr, "Error in reading data: %s.\n", tmp);
      if(verbose>1) printf("  ret := %d\n", ret);
      return(2);
    }
    if(imgNaNs(&img, 1)>0)
      if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
    if(verbose>2) printf("endtime := %g min\n", endtime);
    /* Convert TTAC sample times to min */
    dftTimeunitConversion(&tac, TUNIT_MIN);
    if(verbose>3) dftPrint(&tac);
  }
  /* Check that PET data is image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", petfile);
    imgEmpty(&img); dftEmpty(&tac); return(3);
  }
 
 
 
  /*
   *  Calculate the reference region AUC in time range
   */
  {
    if(verbose>1) printf("calculating TTAC AUC\n");
    double t[2], auc[2];
    t[0]=tstart; t[1]=tstop;
    ret=interpolate(tac.x, tac.voi[0].y, tac.frameNr, t, NULL, auc, NULL, 2);
    if(ret) {
      fprintf(stderr, "Error in integration of reference region TAC.\n");
      if(verbose>1) printf("ret := %d\n", ret);
      imgEmpty(&img); dftEmpty(&tac); return(4);
    }
    aucref=auc[1]-auc[0];
    if(verbose>1) {
      printf("ref_AUC1 := %g\n", auc[0]);
      printf("ref_AUC2 := %g\n", auc[1]);
    }
    if(verbose>0) printf("ref_AUC := %g\n", aucref);
    if(aucref<=1.0E-10) {
      fprintf(stderr, "Error: invalid reference region AUC; check the data.\n");
      imgEmpty(&img); dftEmpty(&tac); return(5);
    }
  }
 
 
  /*
   *  Allocate IMG for ratio image and fill the header info fields
   */
  if(verbose>1) printf("allocating memory for parametric image\n");
  IMG rout;
  imgInit(&rout);
  ret=imgAllocateWithHeader(&rout, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate IMG data.\n");
    if(verbose>1) fprintf(stderr, "ret := %d\n", ret);
    imgEmpty(&img); dftEmpty(&tac); return(6);
  }
  /* Set the frame time */
  rout.start[0]=60.0*tstart; rout.end[0]=60.0*tstop;
  /* and set the units */
  rout.unit=CUNIT_UNITLESS;
 
 
  /*
   *  Compute pixel-by-pixel
   */
  {
    /* Determine the threshold (AUC unit = sec x concentration) */
    float threshold=calcThreshold*tac.voi[0].y2[tac.frameNr-1];
    if(verbose>2) {
      printf("ref_AUC = %g\n", tac.voi[0].y2[tac.frameNr-1]);
      printf("threshold x AUC = %g\n", threshold);
    }
    if(verbose>0) printf("computing pixel-by-pixel\n");
    long long okNr=0; 
    int zi, yi, xi, fi;
    float peti[img.dimt];
    float t[2], auc[2], aucroi;
    t[0]=60.0*tstart; t[1]=60.0*tstop;
    for(zi=0; zi<img.dimz; zi++) {
      if(verbose>5) printf("computing plane %d\n", img.planeNumber[zi]);
      else if(verbose>0 && img.dimz>2) {fprintf(stdout, "."); fflush(stdout);}
      for(yi=0; yi<img.dimy; yi++) for(xi=0; xi<img.dimx; xi++) {
        /* Set results to zero */
        rout.m[zi][yi][xi][0]=0.0;
        /* Check the threshold */
        ret=fpetintegral(img.start, img.end, img.m[zi][yi][xi], img.dimt, peti, NULL);
        /* if AUC at the end is less than threshold value, then 
           do nothing with this pixel */
        if(peti[img.dimt-1]<threshold) continue;
        /* Subtract reference TTAC, if required */
        if(b_per_r!=0) 
          for(fi=0; fi<img.dimt; fi++) 
            img.m[zi][yi][xi][fi]-=tac.voi[0].y[fi];
        /* Calculate PET TAC integrals at time range start and end */
        ret=finterpolate(img.mid, img.m[zi][yi][xi], img.dimt, t, NULL, auc, NULL, 2);
        auc[0]/=60.0; auc[1]/=60.0; 
        if(ret) aucroi=0.0; else aucroi=auc[1]-auc[0];
        /* Calculate the ratio */
        rout.m[zi][yi][xi][0]=aucroi/aucref;
        if(verbose>4 && zi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
          printf("\nVoxel [%d,%d,%d] AUC: %g - %g = %g\n", zi, yi, xi, auc[1], auc[0], aucroi);
          printf("ratio := %g\n", aucroi/aucref);
        }
        /* Prepare for next pixel */
        okNr++;
      } /* next pixel on this plane */
    } /* next image plane */
    if(verbose>2) printf("okNr := %lld\n", okNr);
  }
 
  /* Remove pixel values that are higher than the user-defined limit */
  if(upperLimit>0.0) {
    if(verbose>0) fprintf(stdout, "filtering out ratio values exceeding %g\n", upperLimit);
    imgCutoff(&rout, upperLimit, 0);
  }
  /* Set negative BP values to zero, if required */
  if(nonnegativity_constraint!=0) {
    if(verbose>1) printf("setting negative ratio values to zero\n");
    imgCutoff(&rout, 0.0, 1);
  }
 
 
  /*
   *  Save the ratio image
   */
  if(verbose>1) printf("writing image %s\n", ratfile);
  ret=imgWrite(ratfile, &rout);
  if(ret) {
    fprintf(stderr, "Error: %s\n", rout.statmsg);
    imgEmpty(&img); dftEmpty(&tac); imgEmpty(&rout); return(11);
  }
  if(verbose==1) printf("Ratio image saved: %s\n", ratfile);
 
  if(rout._fileFormat==IMG_E7 || rout._fileFormat==IMG_E7_2D) {
    if(verbose>1) printf("copying header information\n");
    ret=ecat7CopyHeadersNoQuant(petfile, ratfile, verbose-1);
    if(verbose>0 && ret!=0) printf("copying headers not successful (%d).\n", ret);
  }
 
  imgEmpty(&img); dftEmpty(&tac); imgEmpty(&rout);
  if(verbose>0) printf("done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/