dft2img.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcroi.h"
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Creating a PET image file in ECAT 6 or 7 format, with contents from",
  "the user-specified TAC file. Extension of output image filename (.img or .v)",
  "determines whether ECAT 6 or 7 format is used.",
  "Image area is divided into equally sized rectangles, containing as pixel",
  "values the regional TAC values.",
  "Only one image plane is written to image (z dimension will be 1).",
  "Frame times are written as specified in the TAC file.",
  " ",
  "Usage: @P [Options] tacfile imagefile [roifile]",
  " ",
  "If roifile is given, the rectangular regions are saved in it.",
  "Not effective with option -tpl.",
  " ",
  "Options:",
  " -scanner=<Advance|931|HR+|HRRT>",
  "     Scanner, by default GE Advance.",
  " -zoom=<zoom>",
  "     Reconstruction zoom factor, minimum and default is 1.41421356.",
  " -tpl=<filename>",
  "     Image area is not divided into rectangles, but instead the template",
  "     image containing pixel values [1,TAC nr] is used; ROIs are not saved",
  "     with this option.",
  " -dim=<image dimension>",
  "     Image x,y-size in pixels; default 256; not effective with template.",
  " -nodecay",
  "     If TAC data is not corrected for decay, use this option; otherwise",
  "     image headers will tell that image is corrected for decay.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: flat2img, imgadd, img2tif, img2dft, tac2nii, pxl2mask, simframe",
  " ",
  "Keywords: image, ECAT, simulation, software testing",
  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   fi, ri, zi, rowNr, colNr, w, h, row, col;
  DFT   dft;
  char  dfile[FILENAME_MAX], ifile[FILENAME_MAX], rfile[FILENAME_MAX],
        tfile[FILENAME_MAX], *cptr;
  IMG   img;
  int   image_dimx=256, image_dimy=256, image_dimz=1;
  int   scanner_type=SCANNER_ADVANCE;
  float zoom=1.41421356;
  int   decay=IMG_DC_CORRECTED;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dfile[0]=ifile[0]=rfile[0]=tfile[0]=(char)0;
  imgInit(&img); dftInit(&dft);
  /* 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, "SCANNER=", 8)==0) {
      cptr+=8;
      if(*cptr=='9') {scanner_type=SCANNER_ECAT931; continue;}
      if(*cptr=='A' || *cptr=='a') {scanner_type=SCANNER_ADVANCE; continue;}
      if(strcasecmp(cptr, "HRRT")==0) {scanner_type=SCANNER_HRRT; continue;}
      if(strncmp(cptr, "HR", 2)==0 || strncmp(cptr, "hr", 2)==0) {
        scanner_type=SCANNER_HRPLUS; continue;}
    } else if(strncasecmp(cptr, "DIM=", 4)==0) {
      cptr+=4; image_dimx=image_dimy=atoi(cptr); if(image_dimx>1) continue;
    } else if(strncasecmp(cptr, "ZOOM=", 5)==0) {
      cptr+=5; zoom=atof_dpi(cptr); if(zoom>1.41421356) continue;
    } else if(strncasecmp(cptr, "NODECAY", 1)==0) {
      decay=IMG_DC_NONCORRECTED; continue;
    } else if(strncasecmp(cptr, "TPL=", 4)==0) {
      cptr+=4; strlcpy(tfile, cptr, FILENAME_MAX); if(strlen(tfile)>0) 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(dfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(ifile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(rfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
 
  /* Is something missing? */
  if(!ifile[0]) {
    fprintf(stderr, "Error: missing command-line argument; try %s --help\n", argv[0]);
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile := %s\n", dfile);
    printf("ifile := %s\n", ifile);
    printf("rfile := %s\n", rfile);
    printf("tfile := %s\n", tfile);
    printf("scanner_type := %d\n", scanner_type);
    if(!tfile[0]) { // these will be reset later if template is given
      printf("image_dimx := %d\n", image_dimx);
      printf("image_dimy := %d\n", image_dimy);
      printf("image_dimz := %d\n", image_dimz);
    }
    printf("zoom := %g\n", zoom);
    printf("decay := %d\n", decay);
    fflush(stdout);
  }
 
 
  /*
   *  Read regional data
   */
  if(verbose>1) printf("Reading TAC file %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
  if(verbose>2) printf("roiNr := %d\nframeNr := %d\n", dft.voiNr, dft.frameNr);
 
 
  /*
   *  If template is given, then read image dimensions etc from there
   */
  if(tfile[0]) {
    /* Read template image header */
    if(verbose>1) fprintf(stdout, "reading header from template %s\n", tfile);
    IMG timg; imgInit(&timg);
    if(imgReadHeader(tfile, &timg, timg._fileFormat)!=0) {
      fprintf(stderr, "Error: %s\n", timg.statmsg); 
      if(verbose>3) imgInfo(&timg);
      dftEmpty(&dft); imgEmpty(&timg);
      return(3);
    }
    if(verbose>19) imgInfo(&timg);
    image_dimx=timg.dimx;
    image_dimy=timg.dimy;
    image_dimz=timg.dimz;
    zoom=timg.zoom;
    imgEmpty(&timg); // not needed for awhile
    if(verbose>2) {
      printf("image_dimx := %d\n", image_dimx);
      printf("image_dimy := %d\n", image_dimy);
      printf("image_dimz := %d\n", image_dimz);
    }
  }
 
 
  /*
   *  Create the image struct
   */
  /* Allocate memory for it */
  if(verbose>1) printf("allocating memory for image data\n");
  if(imgAllocate(&img, image_dimz, image_dimy, image_dimx, dft.frameNr)) {
    fprintf(stderr, "Error in memory allocation.\n");
    dftEmpty(&dft); imgEmpty(&img); return(4);
  }
  /* Set header information */
  if(verbose>1) printf("setting image headers\n");
  img.type=IMG_TYPE_IMAGE;
  imgSetUnit(&img, dft.unit);
  strcpy(img.studyNr, dft.studynr);
  img.zoom=zoom; 
  img.decayCorrection=decay;
  (void)imgSetScanner(&img, scanner_type);
  img._fileFormat=IMG_UNKNOWN;
  img.planeNumber[0]=1;
  for(fi=0; fi<dft.frameNr; fi++) {
    img.start[fi]=dft.x1[fi]; img.end[fi]=dft.x2[fi]; img.mid[fi]=dft.x[fi];
    if(dft.timeunit==TUNIT_MIN) {
      img.start[fi]*=60.; img.end[fi]*=60.; img.mid[fi]*=60.;}
  }
  /* Scan start time */
  {
    struct tm stm;
    if(get_datetime(dft.scanStartTime, &stm, 0)==0) {
      img.scanStart=timegm(&stm);
    } else {
      img.scanStart=time(NULL);
    }
  }
 
 
  /*
   *  Fill the image data, either square areas, or using template
   */
  if(!tfile[0]) { // no template
 
    ROI_list roilist;
    roiInit(&roilist);
 
    /*
     *  Set the regions
     */
    if(verbose>1) printf("Setting ROIs\n");
    /* Allocate memory */
    if((roilist.roi=(ROI*)calloc(dft.voiNr, sizeof(ROI)))==NULL) {
      fprintf(stderr, "Error in memory allocation.\n");
      dftEmpty(&dft); imgEmpty(&img); return(5);
    }
    for(ri=0; ri<dft.voiNr; ri++) {
      roilist.roi[ri].point_nr=5;
      roilist.roi[ri].x=(int*)calloc(roilist.roi[ri].point_nr, sizeof(int));
      roilist.roi[ri].y=(int*)calloc(roilist.roi[ri].point_nr, sizeof(int));
    }
    roilist.nr=dft.voiNr;
    /* Set ROI headers */
    for(ri=0; ri<dft.voiNr; ri++) {
      strcpy(roilist.roi[ri].imgfile, ifile);
      roilist.roi[ri].zoom=1.0; roilist.roi[ri].recon_zoom=zoom;
      roilist.roi[ri].type=ROI_TRACE;
      sprintf(roilist.roi[ri].roiname, "%s %s",
        dft.voi[ri].voiname, dft.voi[ri].hemisphere);
      roilist.roi[ri].matnum=mat_numcod(1, 1, 1, 0, 0);
      roilist.roi[ri].roi=1+ri;
    }
    img.zoom=roilist.roi[0].recon_zoom;
    /* Calculate the number of ROI rectangles inside image matrix */
    colNr=ceil(sqrt((double)dft.voiNr));
    rowNr=ceil((double)dft.voiNr/(double)colNr);
    /* Calculate the width and height of the ROIs in pixels */
    w=image_dimx/colNr; h=image_dimy/rowNr; if(w<2 || h<2) {
      fprintf(stderr, "Error: too many TACs.\n");
      dftEmpty(&dft); imgEmpty(&img); roiEmpty(&roilist); return(5);
    }
    if(verbose>=0) 
      printf("setting %dx%d ROIs of width %d and height %d pixels with %d frames\n",
             colNr, rowNr, w, h, dft.frameNr);
    /* Set region definitions */
    for(ri=col=row=0; ri<dft.voiNr; ri++) {
      if(verbose>1) printf("Setting ROI (%d,%d) %s:\n", col, row, roilist.roi[ri].roiname);
      roilist.roi[ri].pos_x=col*w; roilist.roi[ri].pos_y=row*h;
      roilist.roi[ri].x[0]=roilist.roi[ri].y[0]=0;
      roilist.roi[ri].x[1]=w-1; roilist.roi[ri].y[1]=0;
      roilist.roi[ri].x[2]=w-1; roilist.roi[ri].y[2]=h-1;
      roilist.roi[ri].x[3]=0; roilist.roi[ri].y[3]=h-1;
      roilist.roi[ri].x[4]=0; roilist.roi[ri].y[4]=0;
      col++; if(col==colNr) {col=0; row++;}
      if(verbose>10) roiPrint(roilist.roi+ri);
    }
    /* Save ROI file, if required */
    if(rfile[0]) {
      if(verbose>=0) printf("Writing ROI file %s\n", rfile);
      if(roiSave(rfile, &roilist)) {
        fprintf(stderr, "Error in writing ROI file.\n");
        dftEmpty(&dft); roiEmpty(&roilist); imgEmpty(&img); return(11);
      }
    }
    /* Set contents of image matrices */
    if(verbose>1) printf("setting image data contents\n");
    for(ri=0; ri<dft.voiNr; ri++) {
      if(verbose>2) printf("  roi=%d\n", ri);
      for(row=roilist.roi[ri].pos_y;
          row<=roilist.roi[ri].pos_y+roilist.roi[ri].y[2]; row++)
        for(col=roilist.roi[ri].pos_x;
            col<=roilist.roi[ri].pos_x+roilist.roi[ri].x[2]; col++)
          for(fi=0; fi<dft.frameNr; fi++) {
            if(verbose>30) printf("    m[%d][%d][%d][%d]=%g\n", 0, row, col, fi, dft.voi[ri].y[fi]);
            for(zi=0; zi<img.dimz; zi++)
              img.m[zi][row][col][fi]=dft.voi[ri].y[fi];
          }
    }
    if(verbose>20) imgInfo(&img);
    roiEmpty(&roilist);
 
  } else { // use template image
 
    int xi, yi, zi, notset;
    IMG timg;
    imgInit(&timg);
 
    /* Read template image */
    if(verbose>1) fprintf(stdout, "reading template %s\n", tfile);
    if(imgRead(tfile, &timg)!=0) {
      fprintf(stderr, "Error: %s\n", timg.statmsg); 
      if(verbose>0) imgInfo(&timg);
      dftEmpty(&dft); imgEmpty(&img); imgEmpty(&timg);
      return(6);
    }
    if(verbose>19) imgInfo(&timg);
 
    /* Check template image dimensions */
    if(timg.dimt>1) {
      fprintf(stderr, "Error: template contains more than one frame.\n");
      dftEmpty(&dft); imgEmpty(&img); imgEmpty(&timg);
      return(6);
    }
    if(timg.dimx!=img.dimx || timg.dimy!=img.dimy || timg.dimz!=img.dimz) {
      fprintf(stderr, "Error: wrong template dimensions.\n");
      dftEmpty(&dft); imgEmpty(&img); imgEmpty(&timg);
      return(6);
    }
 
    /* Set zoom and pixel sizes based on template */
    img.scanner=timg.scanner; img.modality=timg.modality;
    img.zoom=timg.zoom;
    img.sizex=timg.sizex; img.sizey=timg.sizey; img.sizez=timg.sizez;
    img.gapx=timg.gapx; img.gapy=timg.gapy; img.gapz=timg.gapz;
    img.resolutionx=timg.resolutionx; img.resolutiony=timg.resolutiony; img.resolutionz=timg.resolutionz;
    img.axialFOV=timg.axialFOV; img.transaxialFOV=timg.transaxialFOV;
    img.sampleDistance=timg.sampleDistance;
 
 
    /* Set contents of image matrix */
    for(ri=0; ri<dft.voiNr; ri++) dft.voi[ri].size=0.0;
    notset=0;
    for(zi=0; zi<img.dimz; zi++) {
      for(yi=0; yi<img.dimy; yi++) {
        for(xi=0; xi<img.dimx; xi++) {
          ri=temp_roundf(timg.m[zi][yi][xi][0]) - 1;
          if(ri<0 || ri>=dft.voiNr) {
            for(fi=0; fi<dft.frameNr; fi++) img.m[zi][yi][xi][fi]=0.0;
            notset++; continue;
          }
          for(fi=0; fi<dft.frameNr; fi++)
            img.m[zi][yi][xi][fi]=dft.voi[ri].y[fi];
          dft.voi[ri].size+=1.0;
        }
      }
    }
    if(verbose>0) {
      printf("%d pixels were not assigned value from TAC.\n", notset);
      for(ri=0; ri<dft.voiNr; ri++) {
        printf("TAC %d was written in %d pixels.\n", 1+ri, (int)temp_roundf(dft.voi[ri].size));
      }
    }
 
    imgEmpty(&timg);
 
  }
  dftEmpty(&dft);
 
 
  /*
   *  Save image
   */
  if(verbose>1) printf("Writing image file %s ...\n", ifile);
  if(imgWrite(ifile, &img)) {
    fprintf(stderr, "Error: %s\n", img.statmsg); 
    imgEmpty(&img); return(13);
  }
  if(verbose>=0) printf("%s written\n", ifile);
 
  /* Free memory */
  imgEmpty(&img);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/