pxl2tac.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 "libtpccurveio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Extract the TACs of specified pixel(s) in PET image.",
  "Pixels can be specified as mask images, volume range definition files, or",
  "as file containing list of pixel coordinates.",
  " ",
  "Usage: @P [Options] imgfile tacfile pixel(s)",
  " ",
  "Options:",
  " -positives",
  "     Only those pixel TACs which have values over zero are extracted.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Volume range definition (vrd) file is an ASCII text file, which contains",
  "pixel coordinates (x y z; 1..dimension) of the two opposite corners of",
  "the extracted image volume, for example:",
  "     corner1 := 63 57 26",
  "     corner2 := 84 71 44",
  "One or more pixel coordinates (x y z; 1..dimension) can be listed in file,",
  "for example:",
  "     24,52,13",
  "     25,52,14",
  " ",
  "See also: imgthrs, mask2pxl, imgbox, imgmask, imghead, img2dft, imgmaxp",
  " ",
  "Keywords: image, pixel, TAC, 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        ret, fileNr=0, firstfile=0;
  char       petfile[FILENAME_MAX], tacfile[FILENAME_MAX];
  char      *cptr, tmp[128], pxlfile[FILENAME_MAX];
  int        positives=0; // 0=extract all, 1=extract only pixels with value >0
  IMG        img;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=tacfile[0]=pxlfile[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, "POSITIVES", 3)==0) {
      positives=1; continue;
    }
    /* We should not be here */
    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(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  for(; ai<argc; ai++) { // pixel def file(s)
    if(firstfile==0) firstfile=ai;
    fileNr++;
  }
 
  /* Is something missing or wrong? */
  if(!tacfile[0] || fileNr<1) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("tacfile := %s\n", tacfile);
    printf("fileNr := %d\n", fileNr);
    printf("positives := %d\n", positives);
    fflush(stdout);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  imgInit(&img);
  if(verbose>0) printf("reading %s\n", petfile);
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(2);
  }
  if(verbose>0) {
    printf("pet_dimx := %d\n", img.dimx);
    printf("pet_dimy := %d\n", img.dimy);
    printf("pet_dimz := %d\n", img.dimz);
    printf("pet_dimt := %d\n", img.dimt);
  }
  /* 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 pixel definitions
   */
  if(verbose==1) printf("reading pixel positions\n");
  IMG_PIXELS pxl;
  pxlInit(&pxl);
  for(ai=firstfile; ai<argc; ai++) {
    strlcpy(pxlfile, argv[ai], FILENAME_MAX);
    if(verbose>1) printf("reading %s\n", pxlfile);
    /* First, try to read file as pixel list */
    if((ret=pxlRead(&pxl, pxlfile, tmp))==0) {
      continue;
    } else if(verbose>1) {
      printf("could not read as pixel list: %s\n", tmp);
      if(verbose>2) printf("ret := %d\n", ret);
    }
    /* Then try to read as Read Volume Range Definition File */
    VOL_RANGE vol_range={0,0,0,0,0,0};
    IMG_PIXEL p={0,0,0,0};
    if((ret=vrdRead(pxlfile, &vol_range, tmp))==0) {
      vrdReorder(&vol_range);
      if(verbose>1) {
        printf("vol_range.x := %d - %d\n", vol_range.x1, vol_range.x2);
        printf("vol_range.y := %d - %d\n", vol_range.y1, vol_range.y2);
        printf("vol_range.z := %d - %d\n", vol_range.z1, vol_range.z2);
      }
      /* Add range pixels to the list */
      long long n=vrdVxlNr(&vol_range);
      if(n<1) {
        fprintf(stderr, "Warning: no pixels defined in %s\n", pxlfile);
      } else {
        if(pxlAllocateMore(&pxl, n)!=0) {
          fprintf(stderr, "Error: out of memory.\n");
          imgEmpty(&img); pxlFree(&pxl); return(3);
        }
        for(p.z=vol_range.z1; p.z<=vol_range.z2; p.z++)
          for(p.x=vol_range.x1; p.x<=vol_range.x2; p.x++)
            for(p.y=vol_range.y1; p.y<=vol_range.y2; p.y++) 
              pxlAdd(&pxl, &p);
      }
      continue;
    } else if(verbose>1) {
      printf("could not read as vrd file: %s\n", tmp);
      if(verbose>2) printf("ret := %d\n", ret);
    }
    /* Then try to read as mask image */
    if(verbose>1) printf("reading mask image %s\n", pxlfile);
    IMG mask; imgInit(&mask);
    ret=imgRead(pxlfile, &mask);
    if(ret) {
      fprintf(stderr, "Error: %s\n", mask.statmsg);
      if(verbose>1) printf("ret := %d\n", ret);
      imgEmpty(&img); return(4);
    }
    /* Check that dimensions are compatible */
    if(img.dimx!=mask.dimx || img.dimy!=mask.dimy || img.dimz!=mask.dimz) {
      fprintf(stderr, "Error: different image dimensions.\n");
      if(verbose>0) {
        printf("image := %d x %d x %d\n", img.dimx, img.dimy, img.dimz);
        printf("mask := %d x %d x %d\n", mask.dimx, mask.dimy, mask.dimz);
      }
      imgEmpty(&img); imgEmpty(&mask); return(4);
    }
    if(pxlAddFromMask(&pxl, &mask)<1) {
      fprintf(stderr, "Warning: no pixels defined in mask %s\n", pxlfile);
    }
    imgEmpty(&mask);
    continue;
  } // next pixel list, vrd, or mask file
  if(pxl.pxlNr<1) {
    fprintf(stderr, "Error: no pixels to extract.\n");
    imgEmpty(&img); pxlFree(&pxl); return(4);
  }
  if(verbose>6) {
    printf("list of pixels to extract:\n");
    pxlWrite(&pxl, stdout, NULL);
  }
  /* Remove any duplicates */
  ret=pxlRmDuplicates(&pxl);
  if(ret>0 && verbose>0) {
    printf("%d pixel duplicates removed.\n", ret);
  }
  if(verbose>1) {
    printf("list of pixels to extract:\n");
    pxlWrite(&pxl, stdout, NULL);
  }
 
  /* Check that pixels are inside image dimensions */
  ret=0;
  for(long long int i=0; i<pxl.pxlNr; i++) {
    if(pxl.p[i].z<1 || pxl.p[i].z>img.dimz) {
      fprintf(stderr, "Error: pixel outside image z dimension.\n");
      ret++;
    }
    if(pxl.p[i].x<1 || pxl.p[i].x>img.dimx) {
      fprintf(stderr, "Error: pixel outside image x dimension.\n");
      ret++;
    }
    if(pxl.p[i].y<1 || pxl.p[i].y>img.dimy) {
      fprintf(stderr, "Error: pixel outside image x dimension.\n");
      ret++;
    }
  }
  if(ret!=0) {
    imgEmpty(&img); pxlFree(&pxl); return(4);
  }
 
 
  /*
   *  Prepare TAC data
   */
  if(verbose>0) printf("reading %llu pixel TACs.\n", pxl.pxlNr);
  DFT dft; dftInit(&dft);
  /* Allocate memory for TAC */
  ret=dftAllocateWithIMG(&dft, pxl.pxlNr, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for TACs.\n");
    imgEmpty(&img); pxlFree(&pxl);
    return(5);
  }
  if(verbose>13) {
    fflush(stdout);
    printf("  _voidataNr := %d\n", dft._voidataNr);
    printf("  _dataSize := %d\n", dft._dataSize);
    printf("  voiNr := %d\n", dft.voiNr);
    printf("  frameNr := %d\n", dft.frameNr);
    fflush(stdout);
  }
  /* Set DFT information */
  /* File format based on filename */
  {
    char *p; p=strrchr(tacfile, '.');
    if(strcasecmp(p, ".dft")==0) dft._type=DFT_FORMAT_STANDARD;
    else if(strcasecmp(p, ".tac")==0) dft._type=DFT_FORMAT_PMOD;
    else if(strcasecmp(p, ".bld")==0) dft._type=DFT_FORMAT_PMOD;
    else if(strcasecmp(p, ".csv")==0) dft._type=DFT_FORMAT_CSV_UK;
    else dft._type=DFT_FORMAT_PMOD;
  }
  /* Set studynumber */
  if(strlen(dft.studynr)==0) studynr_from_fname2(petfile, dft.studynr, 1);
  if(verbose>1) printf("studynr := %s\n", dft.studynr);
  /* Set scan start time */
  ctime_r_int(&img.scanStart, dft.scanStartTime);
  /* Radiopharmaceutical and isotope */
  strlcpy(dft.radiopharmaceutical, img.radiopharmaceutical, 32);
  dft.decayCorrected=img.decayCorrection;
  if(verbose>10) printf("  isotope := %s\n", imgIsotope(&img));
  strcpy(dft.isotope, imgIsotope(&img));
 
  /* Fill data */
  if(verbose>2) printf("filling TAC data\n");
  int n=0; // TAC
  for(long long i=0; i<pxl.pxlNr; i++) {
    if(verbose>11) {printf("  tac %llu\n", 1+i); fflush(stdout);}
    for(int j=0; j<img.dimt; j++) {
      if(verbose>30) {
        printf("j=%d n=%d z=%d y=%d x=%d\n", j, n, 
               pxl.p[i].z-1, pxl.p[i].y-1, pxl.p[i].x-1);
      }
      dft.voi[n].y[j]=img.m[pxl.p[i].z-1][pxl.p[i].y-1][pxl.p[i].x-1][j];
    }
    /* If required, check that TAC is >0 */
    if(positives!=0) {
      int j;
      for(j=0; j<img.dimt; j++) if(dft.voi[n].y[j]>1.0E-08) break;
      if(j==img.dimt) continue; // all were <= 0
    }
    /* Set "voi" titles */
    if(verbose>20) {printf("  setting ROI name\n"); fflush(stdout);}
    sprintf(dft.voi[n].voiname, "x%04d", pxl.p[i].x);
    sprintf(dft.voi[n].hemisphere, "y%04d", pxl.p[i].y);
    sprintf(dft.voi[n].place, "z%04d", pxl.p[i].z);
    sprintf(dft.voi[n].name, "%s_%s_%s", 
            dft.voi[n].voiname, dft.voi[n].hemisphere, dft.voi[n].place);
    if(verbose>10) printf("   %s\n", dft.voi[n].voiname);
    dft.voi[n].size=img.sizex*img.sizey*img.sizez;
    n++;
  } /* next pixel */
  dft.voiNr=n; dft.frameNr=img.dimt;
  imgEmpty(&img); pxlFree(&pxl);
  if(n<1) {
    fprintf(stderr, "Error: no pixels extracted.\n");
    dftEmpty(&dft); return(9);
  }
  /* Convert times to minutes */
  dftSec2min(&dft);
 
 
  /*
   *  Write the TACs
   */
  if(verbose>1) fprintf(stdout, "writing TACs in %s\n", tacfile);
  if(verbose>3) printf("dft._type := %d\n", dft._type);
  if(verbose>30) dftPrint(&dft);
  if(dftWrite(&dft, tacfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", tacfile, dfterrmsg);
    dftEmpty(&dft); return(11);
  }
  if(verbose>0) fprintf(stdout, "Pixel TACs written in %s\n", tacfile);
  dftEmpty(&dft);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/