imgfsegm.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
#include "libtpcmodext.h"
/*****************************************************************************/
unsigned _stklen = 16777216;  /* need a 16MB stack (not necessary on UNIX) */
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Flexible image segmentation (clustering) for dynamic PET image (1).",
  "This program can be used for method testing purposes only!",
  "Default stack size on Windows may be too small, causing program to crash.",
  " ",
  "User must provide threshold percentage to exclude the background (bkgthrs),",
  "maximal CV percentage allowed for pixels in the same cluster (maxcv),",
  "and threshold percentage for correlation coefficient between cluster TACs",
  "(ccthrs).",
  " ",
  "Static cluster image will contain the cluster numbers as integers, where",
  "0 represents pixels below background threshold, and 1, 2, 3, ... represent",
  "clusters in DECREASING order of TAC integrals.",
  " ",
  "Usage: @P [Options] image bkgthrs maxcv ccthrs clusterimage",
  " ",
  "Options:",
  " -ctac=<filename for cluster TACs>",
  "     Cluster TACs are saved in TAC format.",
  " -cf=<filename for correction factor image>",
  "     Save image file that will contain correction factors for individual",
  "     pixels, calculated as (AUC of cluster avg TAC)/(AUC of pixel TAC).",
  " -sm=<filename for smoothed dynamic image>",
  "     TACs in the original dynamic image are replaced by the cluster TACs",
  "     which have been divided by the correction factors.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P -ctac=b123cluster.dat b123dy1.v 30 5 50 b123cluster.v",
  " ",
  "References:",
  "1. Bentourkia M. A flexible image segmentation prior to parametric",
  "   estimation. Comput Med Imaging Graph. 2001;25:501-506.",
  " ",
  "See also: imgthrs, imgdysmo, imgfiltg, imgmask, img2tif, ecat2ana, pxl2tac",
  " ",
  "Keywords: image, smoothing, mask, threshold",
  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     petfile[FILENAME_MAX], clufile[FILENAME_MAX],
           tacfile[FILENAME_MAX], corfile[FILENAME_MAX],
           smofile[FILENAME_MAX];
  double   bkgthrs, maxcv, ccthrs;
  char    *cptr;
  int      ret;
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=clufile[0]=tacfile[0]=corfile[0]=smofile[0]=(char)0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "CTAC=", 5)==0) {
      strlcpy(tacfile, cptr+5, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "CF=", 3)==0) {
      strlcpy(corfile, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "SM=", 3)==0) {
      strlcpy(smofile, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "SMO=", 4)==0) {
      strlcpy(smofile, cptr+4, FILENAME_MAX); 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);}
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &bkgthrs);
    if(ret!=0 || bkgthrs>=100.0) {
      fprintf(stderr, "Error: invalid background threshold.\n");
      return(1);
    }
    if(bkgthrs<0.0) bkgthrs=0.0; else bkgthrs/=100.0;
    ai++;
  }
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &maxcv);
    if(ret!=0 || maxcv>=100.0) {
      fprintf(stderr, "Error: invalid CV limit.\n");
      return(1);
    }
    if(maxcv<0.0) maxcv=0.0; else maxcv/=100.0;
    ai++;
  }
  if(ai<argc) {
    ret=atof_with_check(argv[ai], &ccthrs);
    if(ret!=0 || ccthrs>=100.0) {
      fprintf(stderr, "Error: invalid correlation coefficient threshold.\n");
      return(1);
    }
    if(ccthrs<0.0) ccthrs=0.0; else ccthrs/=100.0;
    ai++;
  }
  if(ai<argc) {strlcpy(clufile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Did we get all the information that we need? */
  if(!clufile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Check that output files do not have the same name as input */
  if(strcasecmp(petfile, clufile)==0 || strcasecmp(petfile, tacfile)==0 ||
     strcasecmp(petfile, corfile)==0 || strcasecmp(petfile, smofile)==0) {
    fprintf(stderr, "Error: same filename for input and output data.\n");
    return(1);    
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("clufile := %s\n", clufile);
    if(tacfile[0]) printf("tacfile:=%s\n", tacfile);
    if(corfile[0]) printf("corfile:=%s\n", corfile);
    if(smofile[0]) printf("smofile:=%s\n", smofile);
    printf("bkgthrs := %g\n", 100.*bkgthrs);
    printf("maxcv := %g\n", 100.*maxcv);
    printf("ccthrs := %g\n", 100.*ccthrs);
  }
  if(verbose>9) IMG_TEST=verbose-9; else IMG_TEST=0;
 
 
  /*
   *  Read image
   */
  if(verbose>0) fprintf(stdout, "reading image %s\n", petfile);
  IMG img; imgInit(&img);
  ret=imgRead(petfile, &img);
  if(ret) {
    fprintf(stderr, "Error: %s\n", img.statmsg); if(verbose>1) imgInfo(&img);
    return(2);
  }
  if(verbose>1) {
    printf("image dimensions: %d %d %d\n", img.dimz, img.dimy, img.dimx);
    printf("image frame nr: %d\n", img.dimt);
  }
  /* Check that at least two frames */
  if(img.dimt<2) {
    fprintf(stderr, "Error: image is not dynamic.\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)) {
    fprintf(stderr, "Error: image does not contain frame times.\n");
    imgEmpty(&img); return(2);
  }
  /* For raw (sinogram) 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 the integral (area-under-curve) image */
  if(verbose>0) printf("calculating AUC image...\n");
  IMG sumimg; imgInit(&sumimg);
  ret=imgFrameIntegral(&img, 0, img.dimt-1, &sumimg, verbose-5);
  if(ret) {
    fprintf(stderr, "Error: cannot calculate AUC image.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); return(4);
  }
  /* Set unit */
  sumimg.unit=CUNIT_UNITLESS;
 
  /* Search for maximum AUC pixel value */
  if(verbose>0) printf("searching for max AUC pixel value...\n");
  float maxauc;
  ret=imgMax(&sumimg, &maxauc);
  if(ret) {
    fprintf(stderr, "Error: cannot find maximum AUC.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&sumimg);
    return(5);
  }
  if(verbose>1) printf("AUCmax := %g\n", maxauc);
 
  /* Make a mask image by thresholding with given background limit */
  float minauc=maxauc*bkgthrs;
  if(verbose>0) { 
    printf("thresholding image (min=%g, max=%g)...\n", minauc, maxauc); 
    fflush(stdout);
  }
  IMG templimg; imgInit(&templimg);
  ret=imgsegmThresholdMask(&sumimg, (float)minauc, (float)maxauc, &templimg);
  if(ret) {
    fprintf(stderr, "Error: cannot do thresholding.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    return(6);
  }
  if(verbose>10) {
    char tmp[FILENAME_MAX+20], buf[128];
    sprintf(tmp, "%s.templimg.tif", petfile);
    float g=-1.0;
    ret=tiffWriteImg(&templimg, -1, -1, &g, PET_GRAYSCALE_INV, tmp, 0, 0, 
    verbose-10, buf);
    if(ret) {fprintf(stderr, "Error: %s\n", buf); fflush(stderr);}
  }
 
  /* Calculate the CVmax */
  double meanauc, CVmax, CVlim;
  meanauc=0.5*(minauc+maxauc);
  CVmax= ((maxauc-meanauc)*(maxauc-meanauc) + 
          (minauc-meanauc)*(minauc-meanauc)) / meanauc;
  CVlim=maxcv*CVmax;
  if(verbose>1) {
    printf("CVmax :=%g\nCVlim := %g\n", CVmax, CVlim);
    fflush(stdout);
  }
 
  /* Convert mask image to cluster image */
  if(verbose>0) printf("conversion of mask to initial cluster image...\n");
  /* 0th cluster is formed by pixels below background threshold, 
     others get temp value -1 */
  ret=imgsegmMaskToCluster(&templimg);
  if(ret) {
    fprintf(stderr, "Error: cannot make cluster template.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    return(7);
  }
 
  /* Search all clusters */
  if(verbose>0) {printf("searching for clusters...\n"); fflush(stdout);}
  int pi, ri, ci, clusterID=0, clusterNr=0;;
  do {
    /* Find max pixel value, excluding pixels belonging to existing clusters */
    if(verbose>4) {
      printf("finding max pixel value for %dth cluster.\n", clusterID+2);
      fflush(stdout);
    }
    ret=imgsegmFindMaxOutsideClusters(&sumimg, &templimg, &maxauc,
                                      &pi, &ri, &ci);  
    if(ret>0) { /* something was very wrong */
      fprintf(stderr, "Error: cannot search max AUC.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      return(8);
    } else if(ret==-1) { /* we are done */
      if(verbose>4) fprintf(stdout, "All pixels included in clusters.\n");
      break;
    } else { /* ok, we found a maximum */
      if(verbose>4) {
        printf("Max %g found at [%d][%d][%d] (%g)\n", maxauc, pi, ri, ci, 
               templimg.m[pi][ri][ci][0]);
        fflush(stdout);
      }
      /* Check whether all of its neighbours already belong to some cluster */
      ret=imgsegmCheckNeighbours(&templimg, pi, ri, ci);
      /* If so, then find the neighbour with best correlation, and add */
      /* this to the same cluster */
      if(ret) {
        ret=imgsegmFindBestNeighbour(&img, &templimg, pi, ri, ci);
        if(verbose>4) {
          printf("All neighbours are part of clusters; this was added to %d.\n",
                 (int)templimg.m[pi][ri][ci][0]);
          fflush(stdout);
        }
        continue;
      }
      /* Otherwise, this is the basis of a new cluster */
      if(verbose>4) {
        printf("At least one neighbour is free.\n"); fflush(stdout);
      }
      clusterID++;
    }
 
    /* Expand the cluster volume around this pixel */
    if(verbose>4) {
      printf("expanding the cluster around pixel [%d][%d][%d]\n", pi, ri, ci);
      fflush(stdout);
    }
    ret=imgsegmClusterExpand(&templimg, &sumimg, &img, clusterID,
          pi, ri, ci, pi, ri, ci, CVlim, ccthrs, verbose-5);
    if(ret>1) {
      fprintf(stderr, "Error: problem in clustering.\n");
      if(verbose>1) printf("ret=%d\n", ret);
      imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      return(9);
    }
    if(verbose>20) {
      char tmp[FILENAME_MAX+20], buf[128];
      sprintf(tmp, "%s.step%d.tif", petfile, clusterID);
      printf("writing %s of temporary cluster situation\n", tmp);
      fflush(stdout);
      float g=(float)clusterID;
      ret=tiffWriteImg(&templimg, -1, -1, &g, PET_GRAYSCALE_INV, tmp, 0, 0, 
                       verbose-10, buf);
      if(ret) {fprintf(stderr, "Error: %s\n", buf); fflush(stderr);}
    }
    if(verbose==1) {fprintf(stdout, "."); fflush(stdout);}
  } while(ret==0);
  if(verbose==1) {fprintf(stdout, "\n"); fflush(stdout);}
  clusterNr=clusterID;
  if(verbose>0) {printf("\n%d clusters found.\n", clusterNr); fflush(stdout);}
 
 
  /*
   *  Write the cluster image
   */
  if(verbose>1) printf("writing cluster image %s\n", clufile);
  if(imgWrite(clufile, &templimg)) {
    fprintf(stderr, "Error: %s\n", templimg.statmsg); 
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    return(11);
  }
 
  /* If user did not ask for any optional data to be saved, then that's it */
  if(!tacfile[0] && !corfile[0] && !smofile[0]) {
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    if(verbose>0) printf("done.\n");
    return(0);
  }
 
 
  /*
   *  Calculate average TACs from each cluster
   */
  if(verbose>0) printf("calculating cluster average TACs...\n");
  DFT tac; dftInit(&tac);
  ret=clusterTACs(&img, &templimg, clusterNr, &tac, verbose-2);
  if(ret) {
    fprintf(stderr, "Error: cannot calculate cluster average TACs.\n");
    if(verbose>1) printf("ret=%d\n", ret);
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    return(12);
  }
  if(tacfile[0]) {
    /* Save TACs */
    if(verbose>1) printf("writing cluster average TACs in %s\n", tacfile);
    ret=dftWrite(&tac, tacfile);
    if(ret) {
      fprintf(stderr, "Error: %s\n", dfterrmsg);
      imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      dftEmpty(&tac); return(13);
    }
  }
 
  /* If user did not ask for any further data to be saved, then that's it */
  if(!corfile[0] && !smofile[0]) {
    imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
    if(verbose>0) printf("done.\n");
    dftEmpty(&tac); return(0);
  }
 
 
  /*
   *  Calculate the correction image
   */
  if(verbose>0) printf("calculating correction image...\n");
  /* Integrate the cluster TACs */
  if(verbose>2) printf("integrating cluster TACs...\n");
  for(int i=0; i<tac.voiNr; i++) {
    ret=petintegrate(tac.x1, tac.x2, tac.voi[i].y, tac.frameNr,
                     tac.voi[i].y2, NULL);
    if(ret) {   
      fprintf(stderr, "Error: cannot integrate cluster TACs.\n");
      imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      dftEmpty(&tac); return(14);
    }
  }
  /* Convert the sum image to correction image */
  if(verbose>2) printf("calculating image...\n");
  for(pi=0; pi<sumimg.dimz; pi++)
    for(ri=0; ri<sumimg.dimy; ri++)
      for(ci=0; ci<sumimg.dimx; ci++) {
        clusterID=(int)templimg.m[pi][ri][ci][0];
        if(clusterID>0 && sumimg.m[pi][ri][ci][0]>1.0E-006) {
          sumimg.m[pi][ri][ci][0]=
            tac.voi[clusterID-1].y2[tac.frameNr-1]/sumimg.m[pi][ri][ci][0];
        } else sumimg.m[pi][ri][ci][0]=1.0;
      }
  if(corfile[0]) {
    if(verbose>1) printf("writing correction image %s\n", corfile);
    if(imgWrite(corfile, &sumimg)) {
      fprintf(stderr, "Error: %s\n", sumimg.statmsg); 
      imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      dftEmpty(&tac); return(15);
    }
  }
 
 
  /*
   *  Calculate and save the smoothed dynamic image, if required
   */
  if(smofile[0]) {
    if(verbose>0) printf("calculating smoothed dynamic image...\n");
    /* Replace pixel TACs with cluster mean TACs times cf */
    int fi;
    for(pi=0; pi<sumimg.dimz; pi++)
      for(ri=0; ri<sumimg.dimy; ri++)
        for(ci=0; ci<sumimg.dimx; ci++) {
          clusterID=(int)templimg.m[pi][ri][ci][0];
            if(clusterID>0) {
              for(fi=0; fi<img.dimt; fi++)
                img.m[pi][ri][ci][fi]=
                  tac.voi[clusterID-1].y[fi]/sumimg.m[pi][ri][ci][0];
            } else {
              for(fi=0; fi<img.dimt; fi++)
                img.m[pi][ri][ci][fi]=0.0;
            }
        }
    /* For raw data (sinogram), multiply counts by frame duration */
    ret=imgRawCountsPerTime(&img, 0);
    if(ret) {
      fprintf(stderr, "Error: cannot correct sinogram counts.\n");
      dftEmpty(&tac); imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      return(17);
    }
    /* Save the smoothed image */
    if(verbose>1) printf("writing smoothed image %s\n", smofile);
    if(imgWrite(smofile, &img)) {
      fprintf(stderr, "Error: %s\n", img.statmsg); 
      dftEmpty(&tac); imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
      return(18);
    }
  }
 
 
  dftEmpty(&tac);
  imgEmpty(&img); imgEmpty(&sumimg); imgEmpty(&templimg);
  if(verbose>0) printf("done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/