imgbkgrm.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"
#include "libtpcmodel.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Threshold dynamic PET image file in ECAT, NIfTI, or Analyze format.",
  "AIC is used to determine whether the TAC of a single pixel consists of",
  "noise only, or if it resembles the scaled mean of all image pixels.",
  "If no significant signal is found, that pixel is set to zero.",
  "After this operation, PET data can usually be compressed to about",
  "one-third of the original file size.",
  " ",
  "If file name for mask file is given, then original image file is not",
  "modified but a mask file is written, with 0 marking the noise pixels and",
  "1 marking the pixels with signal. Otherwise the original image file is",
  "modified.",
  " ",
  "Usage: @P [Options] imgfile [maskfile]",
  " ",
  "Options:",
  " -keepnegat",
  "     By default, all pixels with AUC<0 are set to zero; with this option",
  "     those are tested as usual.",
  " -bkg=<zero|border>",
  "     Use zero values as background (zero, default), or mean of pixels at",
  "     image x,y-border (border).",
  " -bordertac=<filename>",
  "     Save TAC of average pixel values at image x,y-borders.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P b123dy1.v",
  " ",
  "See also: imgmask, imgthrs, imginteg, imgposv, imgslim",
  " ",
  "Keywords: image, threshold, mask, noise",
  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 imgBorderAverageTAC(
  IMG *img, 
  float *tac,
  float *sd,
  float *med
) {
  if(img==NULL || img->status<IMG_STATUS_OCCUPIED) return(1);
  if(tac==NULL && sd==NULL && med==NULL) return(0);
 
  /* compute the TAC */
  int pxlNr=img->dimz*(2*img->dimx+2*(img->dimy-2));
  float a[pxlNr], fm, fs;
  for(int fi=0; fi<img->dimt; fi++) {
    int i=0;
    for(int zi=0; zi<img->dimz; zi++) {
      int yi, xi;
      yi=0;
      for(xi=0; xi<img->dimx; xi++) a[i++]=img->m[zi][yi][xi][fi];
      yi=img->dimy-1;
      for(xi=0; xi<img->dimx; xi++) a[i++]=img->m[zi][yi][xi][fi];
      xi=0;
      for(yi=1; yi<img->dimy-1; yi++) a[i++]=img->m[zi][yi][xi][fi];
      xi=img->dimx-1;
      for(yi=1; yi<img->dimy-1; yi++) a[i++]=img->m[zi][yi][xi][fi];
    }
    if(tac!=NULL || sd!=NULL) {
      fm=fmean(a, pxlNr, &fs);
      if(tac!=NULL) tac[fi]=fm;
      if(sd!=NULL) sd[fi]=fs;
    }
    if(med!=NULL) med[fi]=fmedian(a, pxlNr);
  }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgNoiseTemplate(
  IMG *img,
  int keepNegat,
  int border,
  IMG *timg,
  int verbose,
  char *status
) {
  if(verbose>0) printf("imgNoiseTemplate(*img, %d, %d, *timg, %d, status)\n",
                       keepNegat, border, verbose);
  if(status!=NULL) strcpy(status, "fault in calling routine");
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(timg==NULL) return(1);
  if(timg->status==IMG_STATUS_OCCUPIED) {
    if(img->dimx!=timg->dimx) return(2);
    if(img->dimy!=timg->dimy) return(2);
    if(img->dimz!=timg->dimz) return(2);
  }
  if(img->dimt<2) {
    if(status!=NULL) strcpy(status, "method works only for dynamic data");
    return(3);
  }
 
  int ret, n;
 
  /*
   *  Compute the mean TAC of all pixels
   */
  n=img->dimt;
  float fdur[n], meany[n], bkgtac[n];
 
  if(verbose>1) fprintf(stdout, "computing mean of pixels\n");
  if(imgExistentTimes(img)) {
    for(int i=0; i<n; i++) fdur[i]=img->end[i]-img->start[i];
  } else {
    for(int i=0; i<n; i++) fdur[i]=1.0;
  }
  ret=imgAverageTAC(img, meany);
  if(ret) {
    if(status!=NULL) strcpy(status, "cannot calculate average TAC");
    if(verbose>1) printf("ret := %d\n", ret);
    return(4);
  }
 
  if(verbose>1) fprintf(stdout, "computing AUC of the mean TAC\n");
  float mean_tac_auc=0.0;
  for(int i=0; i<n; i++) mean_tac_auc+=fdur[i]*meany[i];
  if(verbose>2) printf("mean_tac_auc := %g\n", mean_tac_auc);
  if(fabs(mean_tac_auc)<1.0E-10) {
    if(status!=NULL) strcpy(status, "very small TAC average");
    return(4);
  }
 
  /* If requested, compute the background TAC as the mean of x,y-border pixels */
  double bkg_tac_auc=0.0;
  if(border==0) {
    /* or just zeroes */
    for(int i=0; i<n; i++) bkgtac[i]=0.0;
  } else {
    if(verbose>1) fprintf(stdout, "computing background from border pixels\n");
    if(imgBorderAverageTAC(img, NULL, NULL, bkgtac)!=0) {
      if(status!=NULL) strcpy(status, "cannot calculate background TAC");
      return(5);
    }
    /* Compute AUC for background TAC */
    for(int i=0; i<n; i++) bkg_tac_auc+=fdur[i]*bkgtac[i];
  }
  if(verbose>2) printf("bkg_tac_auc := %g\n", bkg_tac_auc);
 
  /*
   *  Allocate the mask image, if not allocated already
   */
  int maskExisted=0;
  if(timg->status==IMG_STATUS_OCCUPIED) {
    /* Pre-existing mask contents are used later */ 
    maskExisted=1;
  } else {
    ret=imgAllocateWithHeader(timg, img->dimz, img->dimy, img->dimx, 1, img);
    if(ret) {
      if(status!=NULL) strcpy(status, "cannot allocate memory for mask data");
      return(6);
    }
  }
 
  /*
   *  Go through all the pixels
   */
  if(verbose>1) fprintf(stdout, "processing pixels...\n");
  float wsum=0.0; for(int i=0; i<n; i++) wsum+=fdur[i]; if(verbose>4) printf("wsum=%g\n", wsum);
  int zi, yi, xi;
  for(zi=0; zi<img->dimz; zi++) {
    for(yi=0; yi<img->dimy; yi++) {
      for(xi=0; xi<img->dimx; xi++) {
        /* If pixel was already masked out, then keep it that way */
        if(maskExisted && timg->m[zi][yi][xi][0]<0.1) {
          timg->m[zi][yi][xi][0]=0.0; continue;
        }
        /* Compute AUC for this pixel */
        double pxl_tac_auc=0.0;
        for(int i=0; i<n; i++) pxl_tac_auc+=fdur[i]*img->m[zi][yi][xi][i];
        /* If tacAUC<AUCbkg and user did not want to keep negatives, then remove it */
        if(keepNegat==0 && pxl_tac_auc<=bkg_tac_auc) {
          timg->m[zi][yi][xi][0]=0.0; continue;
        }
 
        /* Calculate WSS between this TAC and zero */
        double wss0=0.0;
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i];
          wss0+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and background TAC */
        double wss1=0.0;
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-bkgtac[i];
          wss1+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and scaled background TAC */
        double wss2=0.0; double sf2=1.0;
        if(bkg_tac_auc>1.0E-10) sf2=pxl_tac_auc/bkg_tac_auc;
        if(sf2>2.0) sf2=2.0; // Limit for the scale
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-sf2*bkgtac[i];
          wss2+=fdur[i]*g*g;
        }
 
        /* Calculate WSS between this TAC and scaled mean TAC */
        double wss3=0.0; double sf3=pxl_tac_auc/mean_tac_auc;
        if(sf3<0.5) sf3=0.5; // Limit for the scale
        for(int i=0; i<n; i++) {
          double g=img->m[zi][yi][xi][i]-sf3*meany[i];
          wss3+=fdur[i]*g*g;
        }
 
        /* Compute AICs */
        double aic0=aicSS(wss0, n, 0);
        double aic1=aicSS(wss1, n, 0);
        double aic2=aicSS(wss2, n, 1);
        double aic3=aicSS(wss3, n, 1);
        /* Conserve pixel value if aic3 is better than others */
        if(aic3<aic0 && aic3<aic1 && aic3<aic2) {
          timg->m[zi][yi][xi][0]=1.0;
        } else  {
          timg->m[zi][yi][xi][0]=0.0;
        }
        if(verbose>10) {
          printf(" p[%d][%d][%d]: wss0=%g wss1=%g wss2=%g wss3=%g;  aic0=%g aic1=%g aic2=%g aic3=%g\n",
                 zi, yi, xi, wss0, wss1, wss2, wss3, aic0, aic1, aic2, aic3);
          printf("   -> %g\n", timg->m[zi][yi][xi][0]);
        } else if(verbose>1 && yi==img->dimy-1 && xi==img->dimx-1 && timg->m[zi][yi][xi][0]==1.0) {
          printf(" p[%d][%d][%d]: wss0=%g wss1=%g wss2=%g wss3=%g;  aic0=%g aic1=%g aic2=%g aic3=%g\n",
                 zi, yi, xi, wss0, wss1, wss2, wss3, aic0, aic1, aic2, aic3);
          printf("   sf2=%g sf3=%g\n", sf2, sf3);
          printf("   -> %g\n", timg->m[zi][yi][xi][0]);
        }
      }
    }
  } // next zi
  if(verbose>1) fprintf(stdout, "... done.\n");
 
 
  if(status!=NULL) strcpy(status, "ok");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX], maskfile[FILENAME_MAX], borderfile[FILENAME_MAX];
  int      keepNegat=0; // 0=off; 1=on
  int      bkgMethod=0; // 0=zeroes; 1=x,y-border
  int      ret=0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=maskfile[0]=borderfile[0]=(char)0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "KEEPNEGAT", 5)==0) {
      keepNegat=1; continue;
    } else if(strncasecmp(cptr, "BORDERTAC=", 10)==0) {
      strlcpy(borderfile, cptr+10, FILENAME_MAX); if(borderfile[0]) continue;
    } else if(strncasecmp(cptr, "BKG=", 4)==0) {
      cptr+=4; if(strncasecmp(cptr, "ZERO", 1)==0) {bkgMethod=0; continue;}
      if(strncasecmp(cptr, "BORDER", 1)==0) {bkgMethod=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(petfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(maskfile, 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(!petfile[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("petfile := %s\n", petfile);
    if(maskfile[0]) printf("maskfile := %s\n", maskfile);
    if(borderfile[0]) printf("borderfile := %s\n", borderfile);
    printf("keepNegat := %d\n", keepNegat);
  }
 
 
  /*
   *  Read the contents of the PET file to img data structure
   */
  if(verbose>0) {printf("reading %s\n", petfile); fflush(stdout);}
  IMG img; imgInit(&img);
  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>1) {
    printf("dimt := %d\n", img.dimt);
    printf("dimx := %d\n", img.dimx);
    printf("dimy := %d\n", img.dimy);
    printf("dimz := %d\n", img.dimz);
    printf("type := %d\n", img.type);
  }
  /* Check that we have dynamic data */
  if(img.dimt<2) {
    fprintf(stderr, "Error: method works only for dynamic data.\n");
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Check that frame times are available */
  if(!imgExistentTimes(&img)) {
    if(img.dimt>1) fprintf(stderr, "Warning: unknown frame times.\n");
  }
 
  /* For raw data, divide counts by frame duration */
  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);
  }
 
 
  /*
   *  Calculate border TAC, if requested
   */
  if(borderfile[0]) {
    if(verbose>0) fprintf(stdout, "calculating the TAC of border pixels\n");
    float bv[img.dimt], bs[img.dimt], bm[img.dimt];
    if(imgBorderAverageTAC(&img, bv, bs, bm)!=0) {
      fprintf(stderr, "Error: cannot calculate TAC of border pixels.\n");
      imgEmpty(&img); return(3);
    }
    DFT dft; dftInit(&dft);
    if(dftAllocateWithIMG(&dft, 3, &img)!=0) {
      fprintf(stderr, "Error: cannot make TAC of border pixels.\n");
      imgEmpty(&img); return(3);
    }
    for(int i=0; i<dft.frameNr; i++) dft.voi[0].y[i]=bv[i];
    for(int i=0; i<dft.frameNr; i++) dft.voi[1].y[i]=bs[i];
    for(int i=0; i<dft.frameNr; i++) dft.voi[2].y[i]=bm[i];
    strcpy(dft.voi[0].voiname, "Mean"); strcpy(dft.voi[0].name, dft.voi[0].voiname);
    strcpy(dft.voi[1].voiname, "SD"); strcpy(dft.voi[1].name, dft.voi[1].voiname);
    strcpy(dft.voi[2].voiname, "Median"); strcpy(dft.voi[2].name, dft.voi[2].voiname);
    ret=dftWrite(&dft, borderfile);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", borderfile, dfterrmsg);
      imgEmpty(&img); dftEmpty(&dft); return(11);
    }
    if(verbose>0) printf("TAC of border pixels written in %s\n", borderfile);
    dftEmpty(&dft);
  }
 
 
  /*
   *  Make the mask
   */
  if(verbose>0) fprintf(stdout, "calculating the mask\n");
  IMG mask; imgInit(&mask);
  char buf[128];
  ret=imgNoiseTemplate(&img, keepNegat, bkgMethod, &mask, verbose-1, buf);
  if(ret) {
    fprintf(stderr, "Error: %s.\n", buf);
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&mask); return(4);
  }
  if(imgMaskCount(&mask)<1) {
    fprintf(stderr, "Error: whole image would be considered as noise.\n");
    imgEmpty(&img); imgEmpty(&mask); return(4);
  }
 
 
  /*
   *  If user gave file name for mask, then save it, and quit
   */
  if(maskfile[0]) {
    imgEmpty(&img);
    if(verbose>1) fprintf(stdout, "writing mask image...\n");
    ret=imgWrite(maskfile, &mask);
    if(ret) {
      fprintf(stderr, "Error: %s\n", mask.statmsg);
      imgEmpty(&mask); return(11);
    }
    imgEmpty(&mask);
    if(verbose>0) printf("Written %s\n", maskfile);
    return(0);
  }
 
 
  /*
   *  Otherwise apply the mask to the original PET data and save it.
   */
  if(verbose>1) fprintf(stdout, "applying the mask\n");
  ret=imgThresholdByMask(&img, &mask, 0.0);
  imgEmpty(&mask);
  if(ret) {
    fprintf(stderr, "Error: cannot apply mask to PET data.\n");
    imgEmpty(&img); return(9);
  }
  if(verbose>1) fprintf(stdout, "writing masked image...\n");
  ret=imgWrite(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg);
    imgEmpty(&img); return(12);
  }
  imgEmpty(&img);
  if(verbose>0) printf("Written modified %s\n", petfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/