imgslope.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 "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate the slope of pixel TACs in dynamic PET image.",
  "Slope is calculated between given start and end times (min).",
  "Static image containing the slope value of each pixel TAC is written.",
  " ",
  "Usage: @P [Options] imgfile starttime endtime slopefile",
  " ",
  "Options:",
  " -median",
  "     Instead of LSQ regression line, the slope is determined as median",
  "     of slopes between any two data points.",
  " -sign",
  "     Change the sign of slopes.",
  " -noneg",
  "     Pixels with negative slope (optionally after sign change) are set",
  "     to zero.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgcalc, imgpeak, imgaumc, imgledif, imgthrs",
  " ",
  "Keywords: image, dynamics, input, mask",
  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    imgfile[FILENAME_MAX], slopefile[FILENAME_MAX];
  double  v, tstart=0.0, tstop=0.0;
  int     mode=0; // 0=LSQ, 1=median
  int     sign=0; // 0=natural, 1=changed
  int     noneg=0; // 0=negative values are kept, 1=negative values set to zero
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=slopefile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strcasecmp(cptr, "MEDIAN")==0) {
      mode=1; continue;
    } else if(strcasecmp(cptr, "SIGN")==0) {
      sign=1; continue;
    } else if(strncasecmp(cptr, "NONEGATIVES", 5)==0) {
      noneg=1; 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 && 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(slopefile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!slopefile[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("imgfile := %s\n", imgfile);
    printf("tstart := %g\n", tstart);
    printf("tstop := %g\n", tstop);
    printf("slopefile := %s\n", slopefile);
    printf("mode := %d\n", mode);
    printf("noneg := %d\n", noneg);
    fflush(stdout);
  }
 
 
  /*
   *  Read dynamic image
   */
  if(verbose>0) printf("reading dynamic image %s\n", imgfile);
  IMG img; imgInit(&img);
  if(imgRead(imgfile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    return(2);
  }
  if(verbose>1) {
    printf("img_dimx := %d\n", img.dimx);
    printf("img_dimy := %d\n", img.dimy);
    printf("img_dimz := %d\n", img.dimz);
    printf("img_dimt := %d\n", img.dimt);
    fflush(stdout);
  }
  if(img.dimt<2) {
    fprintf(stderr, "Error: %s contains only 1 time frame.\n", imgfile);
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
  /* Check that times are available */
  if(imgExistentTimes(&img)==0) {
    fprintf(stderr, "Error: image does not contain frame times.\n");
    imgEmpty(&img); return(2);
  }
 
  /* Check and set time range */
  if(tstop<1.0E-10) {
    tstop=img.end[img.dimt-1]/60.0;
    if(verbose>0) {printf("tstop set to %g min\n", tstop); fflush(stdout);}
  }
  tstart*=60.0; tstop*=60.0;
  /* Get the time range in frame indices */
  int fstart=img.dimt-1;
  int fstop=0;
  for(int i=img.dimt-1; i>=0; i--) if(img.mid[i]>=tstart) fstart=i; else break;
  for(int i=0; i<img.dimt; i++) if(img.mid[i]<=tstop) fstop=i; else break;
  if(verbose>1) {
    printf("fstart := %d\n", fstart);
    printf("fstop := %d\n", fstop);
  }
  int sampleNr=1+fstop-fstart;
  if(sampleNr<2) {
    fprintf(stderr, "Error: invalid time range.\n");
    imgEmpty(&img); return(2);
  }
 
 
  /*
   *  Allocate IMG for slope image
   */
  if(verbose>1) printf("allocating memory for parametric image\n");
  IMG sout; imgInit(&sout);
  if(imgAllocateWithHeader(&sout, img.dimz, img.dimy, img.dimx, 1, &img)) {
    fprintf(stderr, "Error: cannot allocate IMG data.\n");
    imgEmpty(&img); return(4);
  }
  /* Set the frame time */
  sout.start[0]=img.start[fstart]; sout.end[0]=img.end[fstop];
  /* and set the units */
  sout.unit=CUNIT_UNITLESS;
 
 
  /*
   *  Compute pixel-by-pixel
   */
  {
    double x[sampleNr], y[sampleNr], a, b;
    for(int i=0; i<sampleNr; i++) x[i]=img.mid[i+fstart];
    for(int zi=0; zi<img.dimz; zi++) {
      for(int yi=0; yi<img.dimy; yi++) for(int xi=0; xi<img.dimx; xi++) {
        /* Set results to zero */
        sout.m[zi][yi][xi][0]=0.0;
        /* Line fit */
        for(int i=0; i<sampleNr; i++) y[i]=img.m[zi][yi][xi][i+fstart];
        int ret=0;
        if(mode==0) {
          ret=regr_line(x, y, sampleNr, &a, &b);
        } else {
          ret=medianline(x, y, sampleNr, &a, &b);
        }
        if(sign) a*=-1.0;
        if(noneg && !(a>=0.0)) a=0.0;
        if(ret==0) sout.m[zi][yi][xi][0]=60.0*(float)a;
      }
    }
  }
 
 
  /*
   *  Save the slope data
   */
  {
    if(verbose>1) printf("writing slopes\n");
    if(imgWrite(slopefile, &sout)) {
      fprintf(stderr, "Error: %s\n", sout.statmsg);
      imgEmpty(&img); imgEmpty(&sout);
      return(11);
    }
    if(verbose>0) fprintf(stdout, "Image %s saved.\n", slopefile);
  }
 
  imgEmpty(&img); imgEmpty(&sout);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/