imgthreshold.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "libtpcimgp.h"
/*****************************************************************************/
 
/*****************************************************************************/
int imgThresholding(
  IMG *img,
  float threshold_level,
  long long *thr_nr
) {
  int zi, yi, xi, fi, ret;
  long long nr=0;
  IMG aucimg;
  float maxauc, thr_limit;
 
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(img->dimt>1) { // dynamic image
    /* Calculate AUC image */
    imgInit(&aucimg);
    ret=imgFrameIntegral(img, 0, img->dimt-1, &aucimg, 0);
    if(ret) return(ret);
    /* Search the max AUC */
    ret=imgMax(&aucimg, &maxauc); if(ret) {imgEmpty(&aucimg); return(ret);}
    /* Calculate the limit value */
    thr_limit=threshold_level*maxauc;
    /* Set to zero the pixels with AUC less than this limit */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          if(aucimg.m[zi][yi][xi][0]<thr_limit) {
            for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            nr++;
          }
    imgEmpty(&aucimg);
  } else { // static image
    /* Search the max value */
    ret=imgMax(img, &maxauc); if(ret) return(ret);
    /* Calculate the limit value */
    thr_limit=threshold_level*maxauc;
    /* Set to zero the pixels with AUC less than this limit */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          if(img->m[zi][yi][xi][0]<thr_limit) {
            for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            nr++;
          }
  }
  if(thr_nr!=NULL) *thr_nr=nr;
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgThresholdingLowHigh(
  IMG *img,
  float lower_threshold_level,
  float upper_threshold_level,
  IMG *timg,
  long long *lower_thr_nr,
  long long *upper_thr_nr
) {
  int zi, yi, xi, fi, ret;
  long long lnr=0, unr=0;
  IMG aucimg;
  float maxauc, lower_thr_limit, upper_thr_limit;
  
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
 
  /* Check if mask image is given and if it exists */
  if(timg!=NULL && timg->status!=IMG_STATUS_OCCUPIED) {
    /* Allocate memory for the mask data */
    ret=imgAllocate(timg, img->dimz, img->dimy, img->dimx, 1);
    if(ret) return(ret);
    /* set mask image header information */
    imgCopyhdr(img, timg);
    timg->start[0]=img->start[0]; timg->end[0]=img->end[img->dimt-1];
    timg->mid[0]=(timg->start[0]+timg->end[0])/2.0;
    /* Initiate mask to 1 */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          timg->m[zi][yi][xi][0]=1.0;
  } else if(timg!=NULL) {
    /* Check that mask dimensions are ok */
    if(timg->dimz!=img->dimz || timg->dimy!=img->dimy || timg->dimx!=img->dimx) return 2;
    if(timg->dimt<1) return 3;
  }
 
  /* Threshold */
  if(img->dimt>1) { // dynamic image
    /* Calculate AUC image */
    imgInit(&aucimg);
    ret=imgFrameIntegral(img, 0, img->dimt-1, &aucimg, 0);
    if(ret) return(ret);
    /* Search the max AUC */
    ret=imgMax(&aucimg, &maxauc); if(ret) {imgEmpty(&aucimg); return(ret);}
    /* Calculate the limit values */
    lower_thr_limit=lower_threshold_level*maxauc;
    upper_thr_limit=upper_threshold_level*maxauc;
    /* Set to zero the pixels with AUC less/more than these limits */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          if(aucimg.m[zi][yi][xi][0]<lower_thr_limit) {
            if(timg==NULL) for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            else timg->m[zi][yi][xi][0]=0.0;
            lnr++;
          } else if(aucimg.m[zi][yi][xi][0]>upper_thr_limit) {
            if(timg==NULL) for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            else timg->m[zi][yi][xi][0]=0.0;
            unr++;
          }
    imgEmpty(&aucimg);
  } else { // static image
    /* Search the max value */
    ret=imgMax(img, &maxauc); if(ret) return(ret);
    /* Calculate the limit values */
    lower_thr_limit=lower_threshold_level*maxauc;
    upper_thr_limit=upper_threshold_level*maxauc;
    /* Set to zero the pixels with AUC less/more than these limits */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          if(img->m[zi][yi][xi][0]<lower_thr_limit) {
            if(timg==NULL) for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            else timg->m[zi][yi][xi][0]=0.0;
            lnr++;
          } else if(img->m[zi][yi][xi][0]>upper_thr_limit) {
            if(timg==NULL) for(fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=0.0;
            else timg->m[zi][yi][xi][0]=0.0;
            unr++;
          }
  }
 
  if(lower_thr_nr!=NULL) *lower_thr_nr=lnr;
  if(upper_thr_nr!=NULL) *upper_thr_nr=unr;
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgThresholdMaskCount(
  IMG *img,
  float minValue,
  float maxValue,
  IMG *timg,
  long long *count
) {
  int fi=0, zi, xi, yi, ret;
 
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(count!=NULL) *count=0;
  /* Check if mask image exists */
  if(timg->status!=IMG_STATUS_OCCUPIED) {
    /* Allocate memory for the template data */
    ret=imgAllocate(timg, img->dimz, img->dimy, img->dimx, 1);
    if(ret) return(ret);
    /* set mask image header information */
    imgCopyhdr(img, timg);
    timg->start[0]=img->start[0]; timg->end[0]=img->end[img->dimt-1];
    timg->mid[0]=(timg->start[0]+timg->end[0])/2.0;
    /* Initiate template to 1 */
    for(zi=0; zi<img->dimz; zi++)
      for(yi=0; yi<img->dimy; yi++)
        for(xi=0; xi<img->dimx; xi++)
          timg->m[zi][yi][xi][0]=1.0;
  } else {
    /* Check that mask image dimensions are ok */
    if(timg->dimz!=img->dimz || timg->dimy!=img->dimy || timg->dimx!=img->dimx)
      return 2;
    if(timg->dimt<1) return(3);
  }
  /* Set mask contents */
  fi=0; 
  long long n=0;
  for(zi=0; zi<img->dimz; zi++)
    for(yi=0; yi<img->dimy; yi++)
      for(xi=0; xi<img->dimx; xi++) {
        /*printf("pxl=%g (%g - %g)\n", img->m[zi][yi][xi][fi], minValue, maxValue);*/
        if(img->m[zi][yi][xi][fi]<minValue || img->m[zi][yi][xi][fi]>maxValue) 
          timg->m[zi][yi][xi][0]=0.0;
        else n++;
      }
  if(count!=NULL) *count=n;
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgThresholdMask(
  IMG *img,
  float minValue,
  float maxValue,
  IMG *timg
) {
  return(imgThresholdMaskCount(img, minValue, maxValue, timg, NULL));
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgThresholdByMask(
  IMG *img,
  IMG *templt,
  float thrValue
) {
  if(img->status!=IMG_STATUS_OCCUPIED) return(1);
  if(templt->status!=IMG_STATUS_OCCUPIED) return(2);
  for(int zi=0; zi<img->dimz; zi++)
    for(int yi=0; yi<img->dimy; yi++) for(int xi=0; xi<img->dimx; xi++)
      if(templt->m[zi][yi][xi][0]==0.0) {  
        for(int fi=0; fi<img->dimt; fi++) img->m[zi][yi][xi][fi]=thrValue;
      }
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
void imgCutoff(
  IMG *image,
  float cutoff,
  int mode
) {
  int zi, yi, xi, ti;
 
  if(mode==0) {
    for(zi=0; zi<image->dimz; zi++)
      for(yi=0; yi<image->dimy; yi++)
        for(xi=0; xi<image->dimx; xi++)
          for(ti=0; ti<image->dimt; ti++)
            if(image->m[zi][yi][xi][ti]>cutoff) image->m[zi][yi][xi][ti]=cutoff;
  } else {
    for(zi=0; zi<image->dimz; zi++)
      for(yi=0; yi<image->dimy; yi++)
        for(xi=0; xi<image->dimx; xi++)
          for(ti=0; ti<image->dimt; ti++)
            if(image->m[zi][yi][xi][ti]<cutoff) image->m[zi][yi][xi][ti]=cutoff;
  }
  return;
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgOutlierFilter(
  IMG *img,
  float limit
) {
  int zi, yi, xi, fi;
  long long nr=0;
  float f;
  IMG tmp;
 
  if(img->status!=IMG_STATUS_OCCUPIED || img->dimt<1) return(-1);
  imgInit(&tmp);
  if(imgAllocate(&tmp, img->dimz, img->dimy, img->dimx, 1)!=0) return(-2);
  /* Process one frame at a time */  
  for(fi=0; fi<img->dimt; fi++) {
    /* Copy one frame to safety */
    for(zi=0; zi<tmp.dimz; zi++)
      for(yi=0; yi<tmp.dimy; yi++)
        for(xi=0; xi<tmp.dimx; xi++)
          tmp.m[zi][yi][xi][0]=img->m[zi][yi][xi][fi];
    /* Go through it */
    for(zi=0; zi<img->dimz; zi++) {
      for(yi=1; yi<img->dimy-1; yi++) for(xi=1; xi<img->dimx-1; xi++) {
        f=tmp.m[zi][yi][xi-1][0]+
          tmp.m[zi][yi][xi+1][0]+
          tmp.m[zi][yi-1][xi][0]+
          tmp.m[zi][yi+1][xi][0]+
          tmp.m[zi][yi+1][xi-1][0]+
          tmp.m[zi][yi+1][xi+1][0]+
          tmp.m[zi][yi-1][xi-1][0]+
          tmp.m[zi][yi-1][xi+1][0];
        f/=8.0;
        if(img->m[zi][yi][xi][fi]>(limit*f)) {img->m[zi][yi][xi][fi]=f; nr++;}
      }
    } /* next plane */
  } /* next frame */
  imgEmpty(&tmp);
 
  return(nr);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgRegionGrowingByThreshold(
  IMG *img,
  const int sz,
  const int sy,
  const int sx,
  float lthr,
  float uthr,
  IMG *mask,
  int verbose
) {
  if(verbose>0) {
    printf("imgRegionGrowingByThreshold(img, %d, %d, %d, %g, %g, mask, %d)\n",
       sz, sy, sx, lthr, uthr, verbose);
    fflush(stdout);
  }
  if(img==NULL || img->status!=IMG_STATUS_OCCUPIED) return(2);
  if(mask==NULL || mask->status!=IMG_STATUS_OCCUPIED) return(3);
  if(img->dimz!=mask->dimz || img->dimy!=mask->dimy || img->dimx!=mask->dimx) return(4);
  int dimz=img->dimz, dimy=img->dimy, dimx=img->dimx;
 
  /* Check that seed pixel resides inside image volume */
  if(sz<0 || sy<0 || sx<0 || sz>=dimz || sy>=dimy || sx>=dimx) {
    if(verbose>1) printf("pixels does not reside inside image\n");
    return(1);
  }
 
  /* Check that seed pixel is not already part of mask */
  if(mask->m[sz][sy][sx][0]>0.5) {
    if(verbose>1) printf("seed already belongs to mask\n");
    return(1);
  }
 
  /* Threshold */
  if(img->m[sz][sy][sx][0]<lthr || img->m[sz][sy][sx][0]>uthr) {
    if(verbose>1) {
      printf("  [%d][%d][%d] outside thresholds\n", sz, sy, sx);
      fflush(stdout);
    }
    mask->m[sz][sy][sx][0]=0.0;
    return(1);
  }
  if(verbose>1) {
    printf("  [%d][%d][%d] added to mask\n", sz, sy, sx);
    fflush(stdout);
  }
  mask->m[sz][sy][sx][0]=1.0;
 
  /* Check the neighbouring pixels */
  int nsz, nsy, nsx;
  nsz=sz; nsy=sy; nsx=sx-1;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
  nsz=sz; nsy=sy; nsx=sx+1;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
  nsz=sz; nsy=sy-1; nsx=sx;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
  nsz=sz; nsy=sy+1; nsx=sx;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
  nsz=sz-1; nsy=sy; nsx=sx;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
  nsz=sz+1; nsy=sy; nsx=sx;
  imgRegionGrowingByThreshold(img, nsz, nsy, nsx, lthr, uthr, mask, verbose);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/