averager.cc

/********************************************************************\

  Name:         averager.cpp
  Created by:   Stefan Ritt

  Contents:     Robust averager

  $Id: averager.cpp 21210 2013-12-12 11:36:59Z ritt $

\********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>

#include "averager.h"

/*----------------------------------------------------------------*/

Averager::Averager(int nx, int ny, int nz, int dim)
{
   fNx = nx;
   fNy = ny;
   fNz = nz;
   fDim = dim;
   
   int size = sizeof(float)*nx*ny*nz * dim;
   fArray = (float *)malloc(size);
   assert(fArray);
   memset(fArray, 0, size);
   size = sizeof(float)*nx*ny*nz;
   fN = (unsigned short *)malloc(size);
   assert(fN);
   memset(fN, 0, size);
}

/*----------------------------------------------------------------*/

Averager::~Averager()
{
   if (fN)
      free(fN);
   if (fArray)
      free(fArray);
   fN = NULL;
   fArray = NULL;
}

/*----------------------------------------------------------------*/

void Averager::Add(int x, int y, int z, float value)
{
   assert(x < fNx);
   assert(y < fNy);
   assert(z < fNz);
   
   int nIndex = (x*fNy + y)*fNz + z;
   if (fN[nIndex] == fDim - 1) // check if array full
      return;
   
   int aIndex = ((x*fNy + y)*fNz + z) * fDim + fN[nIndex];
   fN[nIndex]++;
   fArray[aIndex] = value;
}

/*----------------------------------------------------------------*/

void Averager::Reset()
{
   int size = sizeof(float)*fNx*fNy*fNz * fDim;
   memset(fArray, 0, size);
   size = sizeof(float)*fNx*fNy*fNz;
   memset(fN, 0, size);
}

/*----------------------------------------------------------------*/

int compar(const void *a, const void *b);

int compar(const void *a, const void *b)
{
   if (*((float *)a) == *((float *)b))
      return 0;
   
   return (*((float *)a) < *((float *)b)) ? -1 : 1;
}
           
double Averager::Average(int x, int y, int z)
{
   assert(x < fNx);
   assert(y < fNy);
   assert(z < fNz);

   double a = 0;
   
   int nIndex = (x*fNy + y)*fNz + z;
   int aIndex = ((x*fNy + y)*fNz + z) * fDim;

   for (int i=0 ; i<fN[nIndex] ; i++)
      a += fArray[aIndex + i];
   
   if (fN[nIndex] > 0)
      a /= fN[nIndex];
   
   return a;
}

/*----------------------------------------------------------------*/

double Averager::Median(int x, int y, int z)
{
   assert(x < fNx);
   assert(y < fNy);
   assert(z < fNz);
   
   double m = 0;
   
   int nIndex = (x*fNy + y)*fNz + z;
   int aIndex = ((x*fNy + y)*fNz + z) * fDim;
   
   qsort(&fArray[aIndex], fN[nIndex], sizeof(float), compar);
   m = fArray[aIndex + fN[nIndex]/2];
   
   return m;
}

/*----------------------------------------------------------------*/

double Averager::RobustAverage(double range, int x, int y, int z)
{
   assert(x < fNx);
   assert(y < fNy);
   assert(z < fNz);
   
   double ra = 0;
   int n = 0;
   double m = Median(x, y, z);
   
   int nIndex = (x*fNy + y)*fNz + z;
   int aIndex = ((x*fNy + y)*fNz + z) * fDim;
   
   for (int i=0 ; i<fN[nIndex] ; i++) {
      if (fArray[aIndex + i] > m - range && fArray[aIndex + i] < m + range) {
         ra += fArray[aIndex + i];
         n++;
      }
   }
   
   if (n > 0)
      ra /= n;

   //if (y == 0 && z == 7 && fN[nIndex] > 10)
   //   printf("%d %lf %lf %lf\n", fN[nIndex], a, m, ra);
   
   return ra;
}

/*----------------------------------------------------------------*/

int Averager::SaveNormalizedDistribution(const char *filename, int x, float range)
{
   assert(x < fNx);
   FILE *f = fopen(filename, "wt");
   
   if (!f)
      return 0;
   
   fprintf(f, "X, Y, Z, Min, Max, Ave, Sigma\n");
   
   for (int y=0 ; y<fNy ; y++)
      for (int z=0 ; z<fNz ; z++) {
         
         int nIndex = (x*fNy + y)*fNz + z;
         int aIndex = ((x*fNy + y)*fNz + z) * fDim;
         
         if (fN[nIndex] > 1) {
            fprintf(f, "%d,%d, %d, ", x, y, z);
            
            double s = 0;
            double s2 = 0;
            double min = 0;
            double max = 0;
            int n = fN[nIndex];
            double m = Median(x, y, z);
            
            for (int i=0 ; i<n ; i++) {
               double v = fArray[aIndex + i] - m;
               s += v;
               s2 += v*v;
               if (v < min)
                  min = v;
               if (v > max)
                  max = v;
            }
            double sigma = sqrt((n * s2 - s * s) / (n * (n-1)));
            double average = s / n;
            
            fprintf(f, "%3.1lf, %3.1lf, %3.1lf, %3.3lf, ", min, max, average, sigma);
            
            if (min < -range || max > range) {
               for (int i=0 ; i<n ; i++)
                  fprintf(f, "%3.1lf,", fArray[aIndex + i] - m);
            }
            
            fprintf(f, "\n");
         }
      }
   
   fclose(f);
   return 1;
}