imgprofi.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "List the x and y profiles of PET image at the pixel of maximum intensity",
  "inside the 4D image. X and y profiles are saved separately in TAC format",
  "files; from dynamic image the profiles from every time frame are saved.",
  " ",
  "Usage: @P [Options] imgfile xprofile yprofile [zprofile]",
  " ",
  "Options:",
  " -tif=<TIFF filename>",
  "     Image matrix where profiles are calculated is saved as TIFF image",
  "     showing the profile lines.",
  " -pxl=<x,y,z> or -pxlfile=<filename>",
  "     Instead of maximum, calculate profiles at pixel (x,y,z), where",
  "     x=column (starting from left, 1..width), y=row (starting from top,",
  "     1..height), and z=plane (1..depth).",
  "     Option -pxlfile can be used to give pixel coordinates in file.",
  " -focus",
  "     Profiles are centred on the selected pixel.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example #1:",
  "  @P us3456dy.v profilex.tac profiley.tac",
  " ",
  "See also: imgmaxp, imgpeak, pxl2tac, img2tif, imgthrs, imgdim, mask2pxl",
  " ",
  "Keywords: image, profile, max, input, 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;
  char      imgfile[FILENAME_MAX], tiffile[FILENAME_MAX], pxlfile[FILENAME_MAX];
  char      prxfile[FILENAME_MAX], pryfile[FILENAME_MAX], przfile[FILENAME_MAX];
  int       focus=0; // 0=do not focus profiles to max; 1=do focus
  IMG_PIXEL pos;
  int       ret;
  char     *cptr;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  imgfile[0]=pxlfile[0]=tiffile[0]=prxfile[0]=pryfile[0]=przfile[0]=(char)0;
  pos.x=pos.y=pos.z=pos.f=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, "TIF=", 4)==0) {
      if(strlcpy(tiffile, cptr+4, FILENAME_MAX)>1) continue;
    } else if(strncasecmp(cptr, "PXL=", 4)==0 && strlen(cptr)>4) {
      ret=string_to_xyz(cptr+4, &pos.x, &pos.y, &pos.z);
      if(ret==0 && pos.x>0 && pos.y>0 && pos.z>0) continue;
    } else if(strncasecmp(cptr, "PXLFILE=", 8)==0) {
      if(strlcpy(pxlfile, cptr+8, FILENAME_MAX)>1) continue;
    } else if(strcasecmp(cptr, "FOCUS")==0) {
      focus=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) {strlcpy(prxfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(pryfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(przfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!pryfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(pxlfile[0] && pos.x>0) {
    fprintf(stderr, "Error: do not use -pxl and -pxlfile together.\n");
    return(1);
  }
 
  /* In verbose mode print options */
  if(verbose>1) {
    printf("imgfile := %s\n", imgfile);
    printf("prxfile := %s\n", prxfile);
    printf("pryfile := %s\n", pryfile);
    if(przfile[0]) printf("przfile := %s\n", przfile);
    if(tiffile[0]) printf("tiffile := %s\n", tiffile);
    if(pxlfile[0]) printf("pxlfile := %s\n", pxlfile);
    printf("focus := %d\n", focus);
  }
 
  /*
   *  Read pixel definition file
   */
  if(pxlfile[0]) {
    if(verbose==1) printf("reading %s\n", pxlfile);
    IMG_PIXELS pxl;
    pxlInit(&pxl);
    char buf[256];
    ret=pxlRead(&pxl, pxlfile, buf);
    if(ret) {
      fprintf(stderr, "Error: %s\n", buf);
      fprintf(stderr, "Error: cannot read pixel definition file.\n");
      if(verbose>2) printf("ret := %d\n", ret);
      return(2);
    }
    if(pxl.pxlNr>1) fprintf(stderr, "Warning: using only the first pixel definition.\n");
    pos.x=pxl.p[0].x; pos.y=pxl.p[0].y; pos.z=pxl.p[0].z; 
    pxlFree(&pxl);
  }
  if(verbose>1) {
    if(pos.x>0) printf("pixel := %d,%d,%d\n", pos.x, pos.y, pos.z);
  }
 
 
  /*
   *  Read PET data
   */
  if(verbose>0) printf("reading PET data %s\n", imgfile);
  IMG img; imgInit(&img);
  ret=imgRead(imgfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret := %d\n", ret);
    return(3);
  }
  if(verbose>0) fprintf(stdout, "  image contains %d frames and %d planes.\n", img.dimt, img.dimz);
  if(przfile[0] && img.dimz<3) {
    fprintf(stderr, "Warning: option for z profile ignored for this data.\n");
    przfile[0]=(char)0;
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
 
  /*
   *  If pixel was not specified, then search for max
   */
  double maxv;
  if(pos.x<1) {
    float f;
    if(verbose>1) fprintf(stdout, "  searching for image max\n");
    ret=imgSmoothMax(&img, &f, &pos);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot determine maximum in image.\n");
      imgEmpty(&img);
      return(4);
    }
    maxv=f;
    if(verbose>0) {
      printf("  max_pixel := %d,%d,%d,%d\n", pos.x, pos.y, pos.z, pos.f);
      printf("  max_smoothed_value := %g\n", maxv);
    }
  } else {
    /* maximum value is needed anyway */
    float f;
    if(verbose>1) fprintf(stdout, "  searching for image max value\n");
    ret=imgSmoothMax(&img, &f, NULL);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot determine maximum in image.\n");
      imgEmpty(&img);
      return(4);
    }
    maxv=f;
    if(verbose>0) {printf("  max_smoothed_value := %g\n", maxv);}
  }
 
 
  /*
   *  Make x profile file
   */
  double d;
  DFT dft; dftInit(&dft);
  /* Create DFT */
  if(verbose>1) fprintf(stdout, "creating x profiles\n");
  ret=dftSetmem(&dft, img.dimx, img.dimt); /* frames as columns */
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(&img);
    return(5);
  }
  /* Set DFT contents */
  dft._type=DFT_FORMAT_STANDARD;
  dft.frameNr=img.dimx; dft.voiNr=img.dimt;
  /* set studynumber */
  if(strlen(img.studyNr)>0) strcpy(dft.studynr, img.studyNr);
  else studynr_from_fname(imgfile, dft.studynr);
  strcpy(dft.unit, imgUnit(img.unit)); /* Set calibration unit */
  /* set profile x values */
  dft.timetype=DFT_TIME_MIDDLE; dft.timeunit=TUNIT_MM; 
  if(img.sizex<1.0E-10) img.sizex=1.0;
  d=-(0.5*((double)img.dimx-1.0));
  for(int fi=0; fi<img.dimx; fi++) {
    dft.x[fi]=d*img.sizex; d+=1.0;
  }
  if(focus) {
    double c=(((double)pos.x-1.0) - 0.5*((double)img.dimx-1.0)) * img.sizex;
    if(verbose>2) printf("pos.x=%d cx=%g\n", pos.x, c);
    for(int fi=0; fi<img.dimx; fi++) dft.x[fi]-=c;
  }
  /* set "voi" titles */
  if(dft.voiNr>1) {
    for(int ri=0; ri<dft.voiNr; ri++) {
      char buf[128]; 
      sprintf(buf, "Fr%04d", ri+1); strlcpy(dft.voi[ri].voiname, buf, MAX_REGIONSUBNAME_LEN);
      snprintf(buf, 128, "%g", img.mid[ri]); buf[6]=(char)0;
      snprintf(dft.voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1, "%.6s", buf);
      sprintf(buf, "Pl%04d", pos.z); strlcpy(dft.voi[ri].place, buf, MAX_REGIONSUBNAME_LEN);
      sprintf(dft.voi[ri].name, "%s %s %s",
        dft.voi[ri].voiname, dft.voi[ri].hemisphere, dft.voi[ri].place);
    }
  } else {
    sprintf(dft.voi[0].voiname, "x%04d", pos.x);
    sprintf(dft.voi[0].hemisphere, "y%04d", pos.y);
    sprintf(dft.voi[0].place, "z%04d", pos.z);
    sprintf(dft.voi[0].name, "%s %s %s",
      dft.voi[0].voiname, dft.voi[0].hemisphere, dft.voi[0].place);
  }
  /* copy x profile */
  for(int xi=0; xi<img.dimx; xi++)
    for(int fi=0; fi<img.dimt; fi++)
      dft.voi[fi].y[xi]=img.m[pos.z-1][pos.y-1][xi][fi];
  /* Save DFT */
  ret=dftWrite(&dft, prxfile);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", prxfile, dfterrmsg);
    imgEmpty(&img); dftEmpty(&dft);
    return(11);
  }
  if(verbose>0) fprintf(stdout, "x profile(s) written in %s\n", prxfile);
  dftEmpty(&dft);
 
 
  /*
   *  Make y profile file
   */
  /* Create DFT */
  if(verbose>1) fprintf(stdout, "creating y profiles\n");
  ret=dftSetmem(&dft, img.dimy, img.dimt); /* frames as columns */
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(&img);
    return(5);
  }
  /* Set DFT contents */
  dft._type=DFT_FORMAT_STANDARD;
  dft.frameNr=img.dimy; dft.voiNr=img.dimt;
  /* set studynumber */
  if(strlen(img.studyNr)>0) strcpy(dft.studynr, img.studyNr);
  else studynr_from_fname(imgfile, dft.studynr);
  strcpy(dft.unit, imgUnit(img.unit)); /* Set calibration unit */
  /* set profile x values */
  dft.timetype=DFT_TIME_MIDDLE; dft.timeunit=TUNIT_MM; 
  if(img.sizey<1.0E-10) img.sizey=1.0;
  d=-(0.5*((double)img.dimy-1.0));
  for(int fi=0; fi<img.dimy; fi++) {
    dft.x[fi]=d*img.sizey; d+=1.0;
  }
  if(focus) {
    double c=(((double)pos.y-1.0) - 0.5*((double)img.dimx-1.0)) * img.sizex;
    if(verbose>2) printf("pos.y=%d cy=%g\n", pos.y, c);
    for(int fi=0; fi<img.dimy; fi++) dft.x[fi]-=c;
  }
  /* set "voi" titles */
  if(dft.voiNr>1) {
    for(int ri=0; ri<dft.voiNr; ri++) {
      char buf[128]; sprintf(buf, "Fr%04d", ri+1); strlcpy(dft.voi[ri].voiname, buf, MAX_REGIONSUBNAME_LEN);
      snprintf(buf, 128, "%g", img.mid[ri]); buf[6]=(char)0;
      snprintf(dft.voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1, "%.6s", buf);
      sprintf(buf, "Pl%04d", pos.z); strlcpy(dft.voi[ri].place, buf, MAX_REGIONSUBNAME_LEN);
      sprintf(dft.voi[ri].name, "%s %s %s",
        dft.voi[ri].voiname, dft.voi[ri].hemisphere, dft.voi[ri].place);
    }
  } else {
    sprintf(dft.voi[0].voiname, "x%04d", pos.x);
    sprintf(dft.voi[0].hemisphere, "y%04d", pos.y);
    sprintf(dft.voi[0].place, "z%04d", pos.z);
    sprintf(dft.voi[0].name, "%s %s %s",
      dft.voi[0].voiname, dft.voi[0].hemisphere, dft.voi[0].place);
  }
  /* copy y profile */
  for(int yi=0; yi<img.dimy; yi++)
    for(int fi=0; fi<img.dimt; fi++)
      dft.voi[fi].y[yi]=img.m[pos.z-1][yi][pos.x-1][fi];
  /* Save DFT */
  ret=dftWrite(&dft, pryfile);
  if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", pryfile, dfterrmsg);
    imgEmpty(&img); dftEmpty(&dft);
    return(12);
  }
  if(verbose>0)
    fprintf(stdout, "y profile(s) written in %s\n", pryfile);
  dftEmpty(&dft);
 
 
  /*
   *  Make z profile file if necessary
   */
  if(przfile[0] && img.dimz>2) {
    /* Create DFT */
    if(verbose>1) fprintf(stdout, "creating z profiles\n");
    ret=dftSetmem(&dft, img.dimz, img.dimt); /* frames as columns */
    if(ret) {
      fprintf(stderr, "Error: cannot allocate memory.\n");
      imgEmpty(&img);
      return(5);
    }
    /* Set DFT contents */
    dft._type=DFT_FORMAT_STANDARD;
    dft.frameNr=img.dimz; dft.voiNr=img.dimt;
    /* set studynumber */
    if(strlen(img.studyNr)>0) strcpy(dft.studynr, img.studyNr);
    else studynr_from_fname(imgfile, dft.studynr);
    strcpy(dft.unit, imgUnit(img.unit)); /* Set calibration unit */
    /* set profile x values */
    dft.timetype=DFT_TIME_MIDDLE; dft.timeunit=TUNIT_MM; 
    if(img.sizez<1.0E-10) img.sizez=1.0;
    d=-(0.5*((double)img.dimz-1.0));
    for(int fi=0; fi<img.dimz; fi++) {
      dft.x[fi]=d*img.sizez; d+=1.0;
    }
    if(focus) {
      double c=(((double)pos.z-1.0) - 0.5*((double)img.dimz-1.0)) * img.sizez;
      if(verbose>2) printf("pos.z=%d cy=%g\n", pos.z, c);
      for(int fi=0; fi<img.dimz; fi++) dft.x[fi]-=c;
    }
    /* set "voi" titles */
    if(dft.voiNr>1) {
      for(int ri=0; ri<dft.voiNr; ri++) {
        char buf[128]; sprintf(buf, "Fr%04d", ri+1); strlcpy(dft.voi[ri].voiname, buf, MAX_REGIONSUBNAME_LEN);
        snprintf(buf, 128, "%g", img.mid[ri]); buf[6]=(char)0;
        snprintf(dft.voi[ri].hemisphere, MAX_REGIONSUBNAME_LEN+1, "%.6s", buf);
        sprintf(buf, "Pl%04d", pos.z); strlcpy(dft.voi[ri].place, buf, MAX_REGIONSUBNAME_LEN);
        sprintf(dft.voi[ri].name, "%s %s %s",
          dft.voi[ri].voiname, dft.voi[ri].hemisphere, dft.voi[ri].place);
      }
    } else {
      sprintf(dft.voi[0].voiname, "x%04d", pos.x);
      sprintf(dft.voi[0].hemisphere, "y%04d", pos.y);
      sprintf(dft.voi[0].place, "z%04d", pos.z);
      sprintf(dft.voi[0].name, "%s %s %s",
        dft.voi[0].voiname, dft.voi[0].hemisphere, dft.voi[0].place);
    }
    /* copy z profile */
    for(int zi=0; zi<img.dimz; zi++)
      for(int fi=0; fi<img.dimt; fi++)
        dft.voi[fi].y[zi]=img.m[zi][pos.y-1][pos.x-1][fi];
    /* Save DFT */
    ret=dftWrite(&dft, przfile);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", przfile, dfterrmsg);
      imgEmpty(&img); dftEmpty(&dft);
      return(13);
    }
    if(verbose>0)
      fprintf(stdout, "z profile(s) written in %s\n", przfile);
    dftEmpty(&dft);
  }
 
  /*
   *  Create and save TIFF image if necessary.
   *  Do this as the last step since this messes with the image data.
   */
  if(tiffile[0]) {
    if(verbose>1) fprintf(stdout, "creating TIFF image(s)\n");
    char tmp[FILENAME_MAX+16], extension[FILENAME_MAX];
    int fi, xj, yj;
    float f;
    for(fi=0; fi<img.dimt; fi++) { /* One frame at a time */
      /* Add frame nr to filename, if necessary */
      if(img.dimt>1) {
        if(fi==0) {
          strcpy(tmp, tiffile); cptr=strrchr(tmp, '.');
          if(cptr!=NULL) {strcpy(extension, cptr); *cptr=(char)0;}
          else {strcpy(extension, ".tif");}
        }
        snprintf(tiffile, FILENAME_MAX, "%s_f%d%s", tmp, fi+1, extension);
      }
      /* 'Invert' image values on profile lines */
      for(xj=0; xj<img.dimx; xj++) {
        if(img.m[pos.z-1][pos.y-1][xj][fi]>0.5*maxv)
          img.m[pos.z-1][pos.y-1][xj][fi]=0.0;
        else
          img.m[pos.z-1][pos.y-1][xj][fi]=maxv;
      }
      for(yj=0; yj<img.dimy; yj++) {
        if(img.m[pos.z-1][yj][pos.x-1][fi]>0.5*maxv)
          img.m[pos.z-1][yj][pos.x-1][fi]=0.0;
        else
          img.m[pos.z-1][yj][pos.x-1][fi]=maxv;
      }
      f=-1.0;
      ret=tiffWriteImg(&img, pos.z-1, fi, &f, PET_RAINBOW, tiffile,
                       0, 0, verbose-5, NULL);
      if(ret) fprintf(stderr, "Warning: cannot write %s (%d)\n", tiffile, ret);
      else if(verbose>0) fprintf(stdout, "  %s written.\n", tiffile);
 
    } /* next frame */
  }
 
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/