peak.c

Go to the documentation of this file.

/*****************************************************************************/
#include "libtpcidi.h"
/*****************************************************************************/
 
/*****************************************************************************/
int imgGetConcWeightedPeakPos(
  IMG *img,
  float thrs,
  IMG_PIXEL *pos,
  int verbose
) {
  if(verbose>0) printf("%s(img, %g, ...)\n", __func__, thrs);
 
  int xi, yi, zi, fi=0;
  float d, fmax, fmin, thrslev;
  float px, py, pz, fsum;
 
  if(img->status<IMG_STATUS_OCCUPIED) return(1);
  if(img->dimz<1 || img->dimx<2 || img->dimy<2) return(2);
  if(img->dimt<1) return(3);
  if(thrs<=0.0 || thrs>=1.0) return(4);
  if(pos==NULL) return(5);
 
 
  /* Search min and max value inside the image volume */
  zi=yi=xi=0; fmax=fmin=img->m[zi][yi][xi][fi];
  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][fi] > fmax) fmax=img->m[zi][yi][xi][fi];
        else if(img->m[zi][yi][xi][fi] < fmin) fmin=img->m[zi][yi][xi][fi];
      }
  if(verbose>1) printf("fmin := %g\nfmax := %g\n", fmin, fmax);
 
  /* Calculate threshold level from min and max */
  d=fmax-fmin; if(d<=0.0) return(10);
  thrslev=d*thrs+fmin; if(verbose>1) printf("thrslev := %g\n", thrslev);
 
  /* Calculate concentration weighted peak position */
  px=py=pz=0.0; fsum=0.0;
  for(zi=0; zi<img->dimz; zi++)
    for(yi=0; yi<img->dimy; yi++)
      for(xi=0; xi<img->dimx; xi++) {
        d=img->m[zi][yi][xi][fi]-thrslev;
        if(d>0.0) {
          px+=(float)(xi+1)*d;
          py+=(float)(yi+1)*d;
          pz+=(float)(zi+1)*d;
          fsum+=d;
        }
      }
  if(fsum<=0.0) return 5;
  px/=fsum; py/=fsum; pz/=fsum;
  if(verbose>1) printf("peak_pos := (%g,%g,%g)\n", px, py, pz);
  pos->x=roundf(px); pos->y=roundf(py); pos->z=roundf(pz); pos->f=1;
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/