imgafind.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Trial program to find the locations of arteries in dynamic PET image.",
  " ",
  "Usage: @P [Options] imgfile template",
  " ",
  "Options:",
  " -peak=<filename>",
  "     Pixel peak values are saved as image.",
  " -time=<filename>",
  "     Pixel peak times (1/peaktime) are saved as image.",
  " -ratio=<filename>",
  "     Pixel peak value/time ratios are saved as image.",
  " -lthr=<percentage value>",
  "     Lower threshold value as percentage of maximum.",
  " -uthr=<percentage value>",
  "     Upper threshold value as percentage of maximum.",
  " -maxpeaktime=<time in sec>",
  "     Pixel TAC values after given time are ignored when finding peak.",
  " -minpeaktime=<time in sec>",
  "     Pixel TAC values before given time are ignored when finding peak.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: eabaort, imgprofi, imgthrs, imgfsegm, imgdysmo, imgmax, img2tif",
  " ",
  "Keywords: image, input, blood, aorta",
  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, fi, fj, zi, yi, xi;
  long long  m, n;
  char      *cptr, imgfile[FILENAME_MAX];
  char       peakfile[FILENAME_MAX], timefile[FILENAME_MAX],
             ratiofile[FILENAME_MAX], templfile[FILENAME_MAX];
  float      lower_threshold=0.0, upper_threshold=1.0;
  float      max_peak_time=0.0;
  float      min_peak_time=0.0;
  float      local_mean, best_mean, f1, f2, best_two_mean,
             peak_time, peak_value_per_time;
  int        best_frame;
  IMG        img, template_img, peak_img, time_img, ratio_img;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=(char)0;
  peakfile[0]=timefile[0]=ratiofile[0]=templfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "PEAK=", 5)==0) {
      strlcpy(peakfile, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "TIME=", 5)==0) {
      strlcpy(timefile, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "RATIO=", 6)==0) {
      strlcpy(ratiofile, cptr+6, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "LTHR=", 5)==0) {
      double v;
      if(atof_with_check(cptr+5, &v)==0) {
        if(v>=0.0 && v<100.0) {lower_threshold=0.01*v; continue;}
      }
    } else if(strncasecmp(cptr, "UTHR=", 5)==0) {
      double v;
      if(atof_with_check(cptr+5, &v)==0) {
        if(v>0.0 && v<=100.0) {upper_threshold=0.01*v; continue;}
      }
    } else if(strncasecmp(cptr, "MAXPEAKTIME=", 12)==0) {
      double v;
      if(atof_with_check(cptr+12, &v)==0) {
        if(v>0.0) {max_peak_time=v; continue;}
      }
    } else if(strncasecmp(cptr, "MINPEAKTIME=", 12)==0) {
      double v;
      if(atof_with_check(cptr+12, &v)==0) {
        if(v>0.0) {min_peak_time=v; 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(imgfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(templfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!templfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(lower_threshold>=upper_threshold) {
    fprintf(stderr, "Error: invalid thresholds.\n");
    return(1);
  }
  if(strcasecmp(imgfile, templfile)==0) {
    fprintf(stderr, "Error: check the output filename.\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("peakfile := %s\n", peakfile);
    printf("timefile := %s\n", timefile);
    printf("ratiofile := %s\n", ratiofile);
    printf("templfile := %s\n", templfile);
    printf("lower_threshold := %g\n", lower_threshold);
    printf("upper_threshold := %g\n", upper_threshold);
    printf("max_peak_time := %g s\n", max_peak_time);
    printf("min_peak_time := %g s\n", min_peak_time);
  }
 
 
  /*
   *  Read dynamic image
   */
  if(verbose>0) fprintf(stdout, "reading dynamic image %s\n", imgfile);
  imgInit(&img);
  ret=imgRead(imgfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret := %d\n", ret);
    return(2);
  }
  /* Check if PET data is raw or image */
  if(img.type!=IMG_TYPE_IMAGE) {
    fprintf(stderr, "Error: %s is not an image.\n", imgfile);
    imgEmpty(&img); return(2);
  }
  if(img.dimt<4) {
    fprintf(stderr, "Error: %s contains only %d time frame(s).\n", imgfile, img.dimt);
    imgEmpty(&img); return(2);
  }
  if(verbose>0) fprintf(stdout, "  image contains %d frames and %d planes.\n", img.dimt, img.dimz);
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
 
  /*
   *  Allocate result images
   */
  if(verbose>1) fprintf(stdout, "allocating memory for parametric images\n");
  imgInit(&peak_img); imgInit(&time_img); imgInit(&ratio_img);
  ret=0;
  if(ret==0 && peakfile[0])
    ret=imgAllocateWithHeader(&peak_img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 && timefile[0])
    ret=imgAllocateWithHeader(&time_img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret==0 /* && ratiofile[0]*/) // needed even if not saved
    ret=imgAllocateWithHeader(&ratio_img, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: out of memory.\n");
    imgEmpty(&img);
    imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
    return(4);
  }
  /* and set the units */
  if(timefile[0]) time_img.unit=IMGUNIT_UNITLESS; //IMGUNIT_PER_MIN;
  /*if(ratiofile[0])*/ ratio_img.unit=IMGUNIT_UNITLESS;
  /* and frame times */
  ratio_img.start[0]=0.0; ratio_img.end[0]=1.0;
 
 
  /* Remove outlier pixels */
  n=imgOutlierFilter(&img, 10.0);
  if(verbose>0) {
    printf("%lld pixel(s) filtered out as outliers from dynamic image.\n", n);
  }
 
 
  /*
   *  Thresholding
   */
  if(verbose>1 || (verbose>0 && (lower_threshold>0.0 || upper_threshold<1.0)))
    fprintf(stdout, "thresholding with limits %g-%g%% of AUC max.\n",
      100.*lower_threshold, 100.*upper_threshold);
  imgInit(&template_img);
  ret=imgThresholdingLowHigh(&img, lower_threshold, upper_threshold, 
                             &template_img, &n, &m);
  if(ret)
    fprintf(stderr, "Warning: error %d in thresholding.\n", ret);
  else if(verbose>0 && (n>0 || m>0)) {
    fprintf(stdout, "  %lld pixels below threshold\n", n);
    fprintf(stdout, "  %lld pixels above threshold\n", m);
  }
  /*ret=imgThresholdByTemplate(&img, &template_img, 0.0);
  if(ret)
    fprintf(stderr, "Warning: error %d in thresholding.\n", ret);*/
 
 
  /*
   *  Compute pixel-by-pixel
   */
  if(verbose>0) fprintf(stdout, "computing pixel-by-pixel\n");
  for(zi=0; zi<img.dimz; zi++) {
    if(img.dimz>1 && verbose==1) {fprintf(stdout, "."); fflush(stdout);}
    else if(verbose>1) printf("processing plane %d\n", 1+zi);
    for(yi=0; yi<img.dimy; yi++) {
      for(xi=0; xi<img.dimx; xi++) {
 
        ratio_img.m[zi][yi][xi][0]=0.0;
        if(template_img.m[zi][yi][xi][0]<=0.0) continue; 
 
        /* Search the four consecutive time frames which have the highest mean */
        best_mean=0.0; best_frame=0;
        for(fj=3; fj<img.dimt; fj++) {
          if(img.mid[fj-2]<min_peak_time) continue;
          local_mean=0.0;
          for(fi=fj-3; fi<=fj; fi++) {
            local_mean+=img.m[zi][yi][xi][fi];
            //printf(" %d", fi);
          }
          if(local_mean>best_mean) {
            best_mean=local_mean;
            best_frame=fj-3;
          }
          //printf(" -> %f\n", local_mean/4.0);
        }
        best_mean/=4.0;
        //printf("best_mean=%f\n", best_mean);
        //printf("best_mean_frames=%d-%d\n", best_frame, best_frame+2);
 
        /* Which two means is higher in this range? */
        f1=img.m[zi][yi][xi][best_frame]+img.m[zi][yi][xi][best_frame+1];
        f2=img.m[zi][yi][xi][best_frame+1]+img.m[zi][yi][xi][best_frame+2];
        if(f1>f2) best_two_mean=f1; else best_two_mean=f2;
        best_two_mean*=0.5;
        //printf("best_two_mean=%f\n", best_two_mean);
        if(peakfile[0]) peak_img.m[zi][yi][xi][0]=best_two_mean;
 
        /* Calculate value-weighted mean of peak mid-frame times */
        peak_time=0.0;
        for(fi=best_frame; fi<=best_frame+3; fi++)
          peak_time+=img.m[zi][yi][xi][fi]*img.mid[fi];
        if(best_mean>0.0) peak_time/=(4.0*best_mean); else peak_time=0.0;
        //printf("peak_time=%f\n", peak_time);
        if(timefile[0]) {
          //time_img.m[zi][yi][xi][0]=best_frame+1;
          if(peak_time>0) time_img.m[zi][yi][xi][0]=1.0/peak_time;
          else time_img.m[zi][yi][xi][0]=0.0;
        }
 
        /* If max peak time was set, and limit is met, then remove pixel
           from template and result images */
        if(max_peak_time>0.0 && peak_time>max_peak_time) {
          template_img.m[zi][yi][xi][0]=0;
          if(peakfile[0]) peak_img.m[zi][yi][xi][0]=0.0;
          if(timefile[0]) time_img.m[zi][yi][xi][0]=0.0;
          continue;
        }
 
        /* Calculate peak value/time ratio */      
        if(peak_time>0.0 && best_two_mean>0.0) {
          peak_value_per_time=(best_two_mean)/(peak_time);
        } else peak_value_per_time=0.0;
        //printf("peak_value_per_time=%f\n", peak_value_per_time);
        /*if(ratiofile[0])*/ ratio_img.m[zi][yi][xi][0]=peak_value_per_time;
      } // next column
    } // next row
  } // next plane
  if(verbose==1) fprintf(stdout, "\ndone.\n");
 
  /* Leave only very peaky pixels into the template */
  ret=imgThresholdingLowHigh(&ratio_img, 0.50, 1.0, 
                             &template_img, &n, &m);
  if(ret)
    fprintf(stderr, "Warning: error %d in thresholding.\n", ret);
  else if(verbose>0 && (n>0 || m>0)) {
    fprintf(stdout, "  %lld pixels below threshold\n", n);
    fprintf(stdout, "  %lld pixels above threshold\n", m);
  }
 
 
  /* Region growing and erode boundaries */
 
  /* Calculate TACs from final regions */
 
  /* If required, save 4D image where arterial regions are set to zero */
 
 
  /*
   *  Save result image(s)
   */
  if(verbose>1) printf("writing result image(s)\n");
  if(templfile[0]) {
    ret=imgWrite(templfile, &template_img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", template_img.statmsg);
      imgEmpty(&img); imgEmpty(&template_img);
      imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
      return(11);
    }
    if(verbose>0) fprintf(stdout, "Template image %s saved.\n", templfile);
  }
  if(peakfile[0]) {
    ret=imgWrite(peakfile, &peak_img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", peak_img.statmsg);
      imgEmpty(&img); imgEmpty(&template_img);
      imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
      return(21);
    }
    if(verbose>0) fprintf(stdout, "Peak value image %s saved.\n", peakfile);
  }
  if(timefile[0]) {
    ret=imgWrite(timefile, &time_img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", time_img.statmsg);
      imgEmpty(&img); imgEmpty(&template_img);
      imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
      return(22);
    }
    if(verbose>0) fprintf(stdout, "Peak time image %s saved.\n", timefile);
  }
  if(ratiofile[0]) {
    ret=imgWrite(ratiofile, &ratio_img);
    if(ret) {
      fprintf(stderr, "Error: %s\n", ratio_img.statmsg);
      imgEmpty(&img); imgEmpty(&template_img);
      imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
      return(23);
    }
    if(verbose>0) fprintf(stdout, "Peak value/time image %s saved.\n", ratiofile);
  }
 
 
  imgEmpty(&img); imgEmpty(&template_img);
  imgEmpty(&peak_img); imgEmpty(&time_img); imgEmpty(&ratio_img);
  if(verbose>0) printf("done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/