imgssegm.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[] = {
  "Testing image segmentation for static PET image.",
  "This program can be used for method testing purposes only!",
  " ",
  "Segmented image will contain the segment numbers as integers 0, 1, 2, 3, ...",
  " ",
  "Usage: @P [Options] image segmimage",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imginteg, imgaumc, imgthrs, imgfsegm",
  " ",
  "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;
/*****************************************************************************/
 
/*****************************************************************************/
void psegListNr(
  unsigned int *list,
  unsigned int n
) {
  if(n<1) {printf("(empty)\n"); return;}
  if(n==1) {printf("%u\t1\n", list[0]); return;}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) return;
  unsigned int *sn=sg+n;
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        sn[j]++; break;
      }
    if(j==s) { // no
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  for(unsigned int j=0; j<s; j++)
    printf("%u\t%u\n", sg[j], sn[j]);
  free(sg);
}
/*****************************************************************************/
void psegListContinuous(
  unsigned int *list,
  unsigned int n
) {
  if(n<1) {printf("(empty)\n"); return;}
  if(n==1) {list[0]=0; return;}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) return;
  unsigned int *sn=sg+n;
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        sn[j]++; break;
      }
    if(j==s) { // no
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  //printf("Nr of segments: %u\n", s);
 
  //printf("Unsorted list:\n");
  //for(unsigned int j=0; j<s; j++) printf("%u\t%u\n", sg[j], sn[j]);
  // Sort by pixel number
  for(unsigned int j=0; j<s; j++)
    for(unsigned int k=j+1; k<s; k++) if(sn[k]>sn[j]) {
      unsigned int t=sg[j]; sg[j]=sg[k]; sg[k]=t;
      t=sn[j]; sn[j]=sn[k]; sn[k]=t;
    }
  //printf("Sorted list:\n");
  //for(unsigned int j=0; j<s; j++) printf("%u\t%u\n", sg[j], sn[j]);
 
  /* Reset tag numbers */
  for(unsigned int i=0; i<n; i++)
    for(unsigned int j=0; j<s; j++) if(sg[j]==list[i]) {list[i]=j; break;}
 
  free(sg);
}
/*****************************************************************************/
void psegListImg(
  unsigned int *list,
  unsigned int n,
  IMG *img
) {
  if(img==NULL || n<1) {printf("(empty)\n"); return;}
  if((unsigned int )img->dimz*img->dimy*img->dimx!=n) {printf("(error)\n"); return;}
  //if(n==1) {printf("%u\t1\n", list[0]); return;}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) {printf("(error)\n"); return;}
  unsigned int *sn=sg+n;
 
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        //pmean[j]+=img->pixel[i];
        sn[j]++; break;
      }
    if(j==s) { // no
      //pmean[s]+=img->pixel[i];
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  printf("Nr of segments: %u\n", s);
 
  /* Calculate mean and SD for each segment */
  double *pmean=malloc(2*s*sizeof(double)); if(pmean==NULL) {printf("(error)\n"); free(sg); return;}
  double *psd=pmean+s;
  for(unsigned int j=0; j<s; j++) {
    pmean[j]=psd[j]=0.0;
    unsigned int pi=0;
    double *p=malloc(sn[j]*sizeof(double));
    if(p==NULL) {printf("(error)\n"); free(sg); free(pmean); return;}
    for(unsigned int i=1; i<n && pi<sn[j]; i++) if(list[i]==sg[j]) p[pi++]=img->pixel[i];
    pmean[j]=dmean(p, pi, psd+j);
    free(p);
  }
 
  for(unsigned int j=0; j<s; j++)
    printf("%u\t%u\t%g\t%g\n", sg[j], sn[j], pmean[j], psd[j]);
 
  free(sg); free(pmean);
}
/*****************************************************************************/
void psegListCombineEqual(
  unsigned int *list,
  unsigned int n,
  IMG *img
) {
  if(img==NULL || n<1) {printf("(empty)\n"); return;}
  if((unsigned int )img->dimz*img->dimy*img->dimx!=n) {printf("(error)\n"); return;}
  //if(n==1) {printf("%u\t1\n", list[0]); return;}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) {printf("(error)\n"); return;}
  unsigned int *sn=sg+n;
 
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        //pmean[j]+=img->pixel[i];
        sn[j]++; break;
      }
    if(j==s) { // no
      //pmean[s]+=img->pixel[i];
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  //printf("Nr of segments: %u\n", s);
 
  /* Calculate mean and SD for each segment */
  double *pmean=malloc(2*s*sizeof(double)); if(pmean==NULL) {printf("(error)\n"); free(sg); return;}
  double *psd=pmean+s;
  for(unsigned int j=0; j<s; j++) {
    pmean[j]=psd[j]=0.0;
    unsigned int pi=0;
    double *p=malloc(sn[j]*sizeof(double));
    if(p==NULL) {printf("(error)\n"); free(sg); free(pmean); return;}
    for(unsigned int i=1; i<n && pi<sn[j]; i++) if(list[i]==sg[j]) p[pi++]=img->pixel[i];
    pmean[j]=dmean(p, pi, psd+j);
    free(p);
  }
 
  for(unsigned int j=0; j<s; j++) {
    if(sn[j]==0) continue;
    if(fabs(psd[j])>1.0E-08) continue;
    for(unsigned int k=j+1; k<s; k++) {
      if(sn[k]==0) continue;
      if(fabs(psd[k])>1.0E-08) continue;
      if(fabs(pmean[j]-pmean[k])>1.0E-08) continue;
      printf("combining\n");
      printf("  %u\t%u\t%g\t%g\n", sg[j], sn[j], pmean[j], psd[j]);
      printf("  %u\t%u\t%g\t%g\n", sg[k], sn[k], pmean[k], psd[k]);
      sn[k]=0;
      for(unsigned int i=0; i<n; i++) if(list[i]==sg[k]) list[i]=sg[j];
    }
  }
 
  free(sg); free(pmean);
}
/*****************************************************************************/
double psegListCombineMostSimilar(
  unsigned int *list,
  unsigned int n,
  IMG *img
) {
  if(img==NULL || n<1) {printf("(empty)\n"); return(nan(""));}
  if((unsigned int )img->dimz*img->dimy*img->dimx!=n) {printf("(error)\n"); return(nan(""));}
  //if(n==1) {printf("%u\t1\n", list[0]); return;}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) {printf("(error)\n"); return(nan(""));}
  unsigned int *sn=sg+n;
 
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        //pmean[j]+=img->pixel[i];
        sn[j]++; break;
      }
    if(j==s) { // no
      //pmean[s]+=img->pixel[i];
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  //printf("Nr of segments: %u\n", s);
 
  /* Calculate mean and SD for each segment */
  double *pmean=malloc(2*s*sizeof(double)); if(pmean==NULL) {printf("(error)\n"); free(sg); return(nan(""));}
  double *psd=pmean+s;
  for(unsigned int j=0; j<s; j++) {
    pmean[j]=psd[j]=0.0;
    unsigned int pi=0;
    double *p=malloc(sn[j]*sizeof(double));
    if(p==NULL) {printf("(error)\n"); free(sg); free(pmean); return(nan(""));}
    for(unsigned int i=1; i<n && pi<sn[j]; i++) if(list[i]==sg[j]) p[pi++]=img->pixel[i];
    pmean[j]=dmean(p, pi, psd+j);
    free(p);
  }
 
  double min_t=1.0E+10;
  unsigned int min_j=s+1, min_k=s+1;
  for(unsigned int j=0; j<s; j++) {
    if(sn[j]==0) continue;
    for(unsigned int k=j+1; k<s; k++) {
      if(sn[k]==0) continue;
      if(fabs(psd[k])<1.0E-08 && fabs(psd[j])<1.0E-08) continue;
      // t = | mean1 - mean2 | / sqrt( (sd1^2)/n1 + (sd2^2)/n2 )
      double t=fabs(pmean[j]-pmean[k])/sqrt(psd[j]*psd[j]/(double)sn[j] + psd[k]*psd[k]/(double)sn[k]);
      if(t<min_t) {
        min_t=t;
        min_j=j; min_k=k;
      }
    }
  }
  if(min_j>s || min_k>s) {
    printf("nothing to combine.\n");
    free(sg); free(pmean);
    return(nan(""));
  }
  if(min_t>5.0 /*0.2*/) {
    printf("smallest t=%g is too large, nothing to combine.\n", min_t);
    free(sg); free(pmean);
    return(nan(""));
  }
 
  printf("combining, t=%g\n", min_t);
  printf("  %u\t%u\t%g\t%g\n", sg[min_j], sn[min_j], pmean[min_j], psd[min_j]);
  printf("  %u\t%u\t%g\t%g\n", sg[min_k], sn[min_k], pmean[min_k], psd[min_k]);
  sn[min_k]=0;
  for(unsigned int i=0; i<n; i++) if(list[i]==sg[min_k]) list[i]=sg[min_j];
 
  free(sg); free(pmean);
  return(min_t);
}
/*****************************************************************************/
int psegListPxlMove(
  unsigned int *list,
  unsigned int n,
  IMG *img
) {
  if(img==NULL || n<1) {printf("(empty)\n"); return(2);}
  if((unsigned int )img->dimz*img->dimy*img->dimx!=n) {printf("(error)\n"); return(2);}
 
  unsigned int *sg=malloc(2*n*sizeof(unsigned int)); if(sg==NULL) {printf("(error)\n"); return(2);}
  unsigned int *sn=sg+n;
 
  sg[0]=list[0]; sn[0]=1;
  unsigned int s=1;
  for(unsigned int i=1; i<n; i++) {
    unsigned int j;
    /* check if this segment number is listed already */
    for(j=0; j<s; j++)
      if(list[i]==sg[j]) { // yes
        //pmean[j]+=img->pixel[i];
        sn[j]++; break;
      }
    if(j==s) { // no
      //pmean[s]+=img->pixel[i];
      sg[s]=list[i]; sn[s]=1; s++;
    }
  }
  //printf("Nr of segments: %u\n", s);
 
  /* Calculate mean and SD for each segment */
  double *pmean=malloc(2*s*sizeof(double)); if(pmean==NULL) {printf("(error)\n"); free(sg); return(3);}
  double *psd=pmean+s;
  for(unsigned int j=0; j<s; j++) {
    pmean[j]=psd[j]=0.0;
    unsigned int pi=0;
    double *p=malloc(sn[j]*sizeof(double));
    if(p==NULL) {printf("(error)\n"); free(sg); free(pmean); return(4);}
    for(unsigned int i=1; i<n && pi<sn[j]; i++) if(list[i]==sg[j]) p[pi++]=img->pixel[i];
    pmean[j]=dmean(p, pi, psd+j);
    free(p);
  }
 
  /* Find the worst pixel in each segment and check if would fit better into another segment */
  unsigned int moved=0;
  for(unsigned int j=0; j<s; j++) {
    if(fabs(psd[j])<1.0E-08) continue; // all pixels in this segment fit well already
    double max_dif=0.0, pxlv=nan("");
    unsigned int max_i=n+1;
    for(unsigned int i=0; i<n; i++) if(list[i]==sg[j]) {
      double d=fabs(pmean[j]-img->pixel[i]);
      if(d>max_dif) {max_dif=d; max_i=i; pxlv=img->pixel[i];}
    }
    if(max_i>n || max_dif<1.0E-06) continue;
    printf("  group %u largest dif is %g\n", sg[j], max_dif);
    /* Is the mean of another group closer? */
    unsigned int better_k=s+1;
    double min_d=max_dif;
    for(unsigned int k=0; k<s; k++) if(k!=j) {
      double d=fabs(pmean[k]-pxlv);
      if(d<min_d) {min_d=d; better_k=k;} 
    }
    if(better_k>s) continue;
    printf("    group %u would be better match with dif %g\n", sg[better_k], min_d);
    printf("    pixel value is %g; means of own and better groups are %g and %g\n", pxlv, pmean[j], pmean[better_k]);
    list[max_i]=sg[better_k];
    moved++;
  }
 
  free(sg); free(pmean);
  if(moved==0) return(1);
  return(0);
}
/*****************************************************************************/
int imgSegmentStaticTest(
  IMG *img, 
  IMG *out,
  int verbose
) {
  if(verbose>0) {printf("%s(IMG, IMG, ...)\n", __func__); fflush(stdout);}
  if(img==NULL || img->status<IMG_STATUS_OCCUPIED || out==NULL) return(1);
 
  /* Number of pixels */
  unsigned int pxlNr=img->dimz*img->dimy*img->dimx;
  if(verbose>1) printf("pxlNr := %u\n", pxlNr);
  if(pxlNr<54) {
    if(verbose>0) fprintf(stderr, "Error: too few pixels.\n");
    return(2);
  }
 
  /* Number of blocks in each dimension */
  unsigned int bdx=1+(img->dimx-1)/3;
  unsigned int bdy=1+(img->dimy-1)/3;
  unsigned int bdz=1+(img->dimz-1)/3;
  if(verbose>2) printf("block dimensions x y z := %u %u %u\n", bdx, bdy, bdz);
 
 
  /*
   *  Allocate memory for segmented image
   */
  imgEmpty(out);
  if(imgAllocateWithHeader(out, img->dimz, img->dimy, img->dimx, 1, img)) {
    if(verbose>0) fprintf(stderr, "Error: cannot allocate memory for segmented image.\n");
    return(3);
  }
 
 
  /*
   *  Allocate a list of segment numbers for each pixel
   */
  unsigned int *pseg=malloc(sizeof(unsigned int)*pxlNr);
  if(pseg==NULL) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    imgEmpty(out); return(4);
  }
  /* Initiate segmentation by defining each (preferably) 3x3x3 volume as one segment */
  unsigned int s=0, si=0;
  for(int zi=0; zi<img->dimz; zi++)
    for(int yi=0; yi<img->dimy; yi++)
      for(int xi=0; xi<img->dimx; xi++) {
        s=(zi/3)*bdx*bdy + (yi/3)*bdx + (xi/3);
        pseg[si]=s;
        si++;
      }
  //if(verbose>4) psegListNr(pseg, si);
 
  if(verbose>2) psegListImg(pseg, si, img);
 
  psegListCombineEqual(pseg, si, img);
  if(verbose>2) psegListImg(pseg, si, img);
 
  double t=0.0;
  do {
    if(verbose>2) printf("combining most similar\n");
    t=psegListCombineMostSimilar(pseg, si, img);
  } while(!isnan(t));
  if(verbose>2) psegListImg(pseg, si, img);
 
  do {
    if(verbose>2) printf("move pixels between segments\n");
  } while(psegListPxlMove(pseg, si, img)==0);
  if(verbose>2) psegListImg(pseg, si, img);
 
  psegListCombineEqual(pseg, si, img);
  if(verbose>2) psegListImg(pseg, si, img);
 
  do {
    if(verbose>2) printf("combining most similar\n");
    t=psegListCombineMostSimilar(pseg, si, img);
  } while(!isnan(t));
  if(verbose>2) psegListImg(pseg, si, img);
 
  /* Make segment number contiguous */
  psegListContinuous(pseg, si);
 
  if(verbose>2) psegListImg(pseg, si, img);
 
  /* Write segments into output image */
  si=0;
  for(int zi=0; zi<out->dimz; zi++)
    for(int yi=0; yi<out->dimy; yi++)
      for(int xi=0; xi<out->dimx; xi++) {
        out->m[zi][yi][xi][0]=(float)pseg[si];
        si++;
      }
  free(pseg);
 
 
 
 
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     petfile[FILENAME_MAX], segfile[FILENAME_MAX];
  int      ret;
 
  
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  petfile[0]=segfile[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;
    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) {strlcpy(segfile, 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(!segfile[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, segfile)==0) {
    fprintf(stderr, "Error: same file name for input and output data.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("petfile := %s\n", petfile);
    printf("segfile := %s\n", segfile);
  }
  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 only one frame */
  if(img.dimt>1) {
    fprintf(stderr, "Error: image must be static (one frame).\n");
    imgEmpty(&img); return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
 
 
  /*
   *  Segmentation
   */
  if(verbose>0) fprintf(stdout, "segmenting image %s\n", petfile);
  IMG seg; imgInit(&seg);
  ret=imgSegmentStaticTest(&img, &seg, verbose-2);
  if(ret) {
    fprintf(stderr, "Error: segmentation not successful.\n");
    imgEmpty(&img); imgEmpty(&seg); return(3);
  }
 
 
  /*
   *  Write the segmented image
   */
  if(verbose>1) printf("writing segmented image %s\n", segfile);
  if(imgWrite(segfile, &seg)) {
    fprintf(stderr, "Error: %s\n", seg.statmsg); 
    imgEmpty(&img); imgEmpty(&seg);
    return(11);
  }
 
  imgEmpty(&img); imgEmpty(&seg);
  if(verbose>0) printf("done.\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/