imgmaxp.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 "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
int imgRangeWeightedMax(IMG *img, IMG_RANGE *r, IMG_PIXEL *maxp);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Find the position of maximum pixel value in PET image in ECAT 6.3 or 7.x,",
  "NIfTI-1, or Analyze 7.5 format.",
  "Position of the found pixel is written in specified file or stdout as its",
  "coordinates (x,y,z,f; 1..dimension).",
  " ",
  "Usage: @P [Options] imgfile [pxlfile]",
  " ",
  "Options:",
  " -min",
  "     Find minimum instead of maximum.",
  " -vrdfile=<filename>",
  "     Volume range definition (vrd) file is an ASCII text file, which",
  "     contains pixel coordinates (x y z f; 1..dimension) of the two opposite",
  "     corners of the image space to be searched for; for example:",
  "       corner1 := 63 57 26 1",
  "       corner2 := 84 71 44 17",
  " -wmax",
  "     Find pixel value weighted position inside the image.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgmax, imgpeak, imginteg, pxl2tac, pxl2mask, imgprofi, imgmask",
  " ",
  "Keywords: image, pixel, max, min",
  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;
  int     find_what=0; // 0=max, 1=min, 2=wmax
  char    imgfile[FILENAME_MAX], pxlfile[FILENAME_MAX], vrdfile[FILENAME_MAX];
  char   *cptr=NULL;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=pxlfile[0]=vrdfile[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, "VRDFILE=", 8)==0) {
      strlcpy(vrdfile, cptr+8, FILENAME_MAX); if(vrdfile[0]) continue;
    } else if(strncasecmp(cptr, "MINIMUM", 3)==0) {
      find_what=1; continue;
    } else if(strncasecmp(cptr, "MAXIMUM", 3)==0) {
      find_what=0; continue;
    } else if(strncasecmp(cptr, "WMAX", 1)==0) {
      find_what=2; 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(pxlfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Did we get all the information that we need? */
  if(!imgfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    if(vrdfile[0]) printf("vrdfile := %s\n", vrdfile);
    if(pxlfile[0]) printf("pxlfile := %s\n", pxlfile);
    printf("find_what := %d\n", find_what);
    fflush(stdout);
  }
 
  if(verbose>2) IMG_TEST=ECAT63_TEST=ECAT7_TEST=verbose-2;
  else IMG_TEST=ECAT63_TEST=ECAT7_TEST=0;
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  IMG img; imgInit(&img);
  if(verbose>0) printf("reading %s\n", imgfile);
  ret=imgRead(imgfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); 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);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Correct sinogram data for decay */
  if(img.type==IMG_TYPE_RAW) {
    if(verbose>0) fprintf(stdout, "correcting sinogram data for decay\n");
    ret=imgDecayCorrection(&img, 1);
    if(ret) {
      fprintf(stderr, "Error %d: cannot decay correct dynamic sinogram.\n", ret);
      if(ret==1) fprintf(stderr, "      Sinogram contains no isotope.\n");
      if(ret==2) fprintf(stderr, "      Sinogram is already decay corrected.\n");
      imgEmpty(&img);
      return(2);
    }
  }
  /* Correct sinogram data for frame lengths */
  if(img.type==IMG_TYPE_RAW) {
    if(verbose>0) fprintf(stdout, "correcting sinogram data for frame lengths\n");
    ret=imgRawCountsPerTime(&img, 1);
    if(ret) {
      fprintf(stderr, "Error: cannot correct sinogram counts.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); return(2);
    }
  }
 
 
 
  /*
   *  Read image volume range definition, if requested
   */
  IMG_RANGE img_range={0,0,0,0,0,0,0,0};
  if(vrdfile[0]) {
    if(verbose==1) printf("reading image volume range definition file\n");
    if(verbose>1) printf("reading %s\n", vrdfile);
    char tmp[128];
    ret=irdRead(vrdfile, &img_range, tmp);
    if(ret) {
      fprintf(stderr, "Error: could not read image range.\n");
      if(verbose>2) printf("ret := %d\n", ret);
      imgEmpty(&img); return(3);
    }
    if(verbose>1) {
      printf("img_range.x := %d - %d\n", img_range.x1, img_range.x2);
      printf("img_range.y := %d - %d\n", img_range.y1, img_range.y2);
      printf("img_range.z := %d - %d\n", img_range.z1, img_range.z2);
      printf("img_range.f := %d - %d\n", img_range.f1, img_range.f2);
    }
 
    /* Check that image range is inside image dimensions;
       fix zero time frame range if necessary. */
    ret=irdCheck(&img_range, &img);
    if(ret) {
      fprintf(stderr, "Error: invalid image range.\n");
      if(verbose>2) printf("ret := %d\n", ret);
      imgEmpty(&img); return(3);
    }
  }
 
 
  /*
   *  Find image min and max
   */
  IMG_PIXEL maxp, minp;
  float maxv, minv;
  if(find_what==0) {
    if(verbose>0) printf("searching for max in %s\n", imgfile);
    if(vrdfile[0]) ret=imgRangeMinMax(&img, &img_range, &maxp, &maxv, &minp, &minv);
    else ret=imgRangeMinMax(&img, NULL, &maxp, &maxv, &minp, &minv);
  } else if(find_what==1) {
    if(verbose>0) printf("searching for min in %s\n", imgfile);
    if(vrdfile[0]) ret=imgRangeMinMax(&img, &img_range, &maxp, &maxv, &minp, &minv);
    else ret=imgRangeMinMax(&img, NULL, &maxp, &maxv, &minp, &minv);
  } else {
    if(verbose>0) printf("searching for weighted max in %s\n", imgfile);
    if(vrdfile[0]) ret=imgRangeWeightedMax(&img, &img_range, &maxp);
    else ret=imgRangeWeightedMax(&img, NULL, &maxp);
  }
  if(ret) {
    fprintf(stderr, "Error: cannot find min/max in %s\n", imgfile);
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(5);
  }
  if(find_what!=2 && verbose>1) {
    printf("image_min := %g\n", minv);
    printf("image_max := %g\n", maxv);
  }
  imgEmpty(&img);
 
 
  /*
   *  Save or print the pixel coordinates
   */
  IMG_PIXEL *ip;
  if(find_what==0 || find_what==2) ip=&maxp; else ip=&minp;
  if(!pxlfile[0]) {
    fprintf(stdout, "%d,%d,%d,%d\n", ip->x, ip->y, ip->z, ip->f);
    return(0);
  }
  FILE *fp; fp=fopen(pxlfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s).\n", pxlfile);
    return(11);
  }
  fprintf(fp, "%d,%d,%d,%d\n", ip->x, ip->y, ip->z, ip->f);
  fclose(fp);
  if(ret) {
    fprintf(stderr, "Error: cannot write in file %s.\n", pxlfile);
    return(12);
  }
  if(verbose>0) printf("Pixel coordinates written in %s\n", pxlfile);
 
  return(0);
}
/*****************************************************************************/
/*****************************************************************************/
int imgRangeWeightedMax(
  IMG *img,
  IMG_RANGE *r,
  IMG_PIXEL *maxp
) {
  int zi, yi, xi, fi;
  double wsum, wz, wy, wx, wf;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt<1 || img->dimz<1 || img->dimy<1 || img->dimx<1) return(1);
  if(maxp==NULL) return(1);
  wsum=wz=wy=wx=wf=0.0;
 
  if(r!=NULL) {
    if(r->z1<1 || r->y1<1 || r->x1<1 || r->f1<1) return(2);
    if(r->z2<r->z1 || r->y2<r->y1 || r->x2<r->x1 || r->f2<r->f1) return(3);
    if(r->z2>img->dimz || r->y2>img->dimy || r->x2>img->dimx || r->f2>img->dimt) return(4);
 
    maxp->z=r->z1; maxp->y=r->y1; maxp->x=r->x1; maxp->f=r->f1;
    for(zi=r->z1-1; zi<r->z2; zi++) {
      for(yi=r->y1-1; yi<r->y2; yi++) {
        for(xi=r->x1-1; xi<r->x2; xi++) {
          for(fi=r->f1-1; fi<r->f2; fi++) {
            wsum+=img->m[zi][yi][xi][fi];
            wx+=(double)xi*img->m[zi][yi][xi][fi];
            wy+=(double)yi*img->m[zi][yi][xi][fi];
            wz+=(double)zi*img->m[zi][yi][xi][fi];
            wf+=(double)fi*img->m[zi][yi][xi][fi];
          }
        }
      }
    }
  } else {
    maxp->x=maxp->y=maxp->z=maxp->f=1;
    for(zi=0; zi<img->dimz; zi++) {
      for(yi=0; yi<img->dimy; yi++) {
        for(xi=0; xi<img->dimx; xi++) {
          for(fi=0; fi<img->dimt; fi++) {
            wsum+=img->m[zi][yi][xi][fi];
            wx+=(double)xi*img->m[zi][yi][xi][fi];
            wy+=(double)yi*img->m[zi][yi][xi][fi];
            wz+=(double)zi*img->m[zi][yi][xi][fi];
            wf+=(double)fi*img->m[zi][yi][xi][fi];
          }
        }
      }
    }
  }
  if(wsum<1.0E-08) return(5);
  wx/=wsum; wy/=wsum; wz/=wsum; wf/=wsum;
  maxp->x=1+(int)roundf(wx); maxp->y=1+(int)roundf(wy); 
  maxp->z=1+(int)roundf(wz); maxp->f=1+(int)roundf(wf);
  return(0);
}