MedianFilter.h

#ifndef __MEDIAN_FILTER_H__
#define __MEDIAN_FILTER_H__

#include <vuVector.h>
#include <Utilities.h>
#include <vuML.h>
#include <Dataset.h>

#include <VVector.h>
#include <VVectorFun.h>
using namespace math;

#include <assert.h>
#include <cmath>

#include <string>
#include <algorithm>
using namespace std;

class Neighborhood {

 public: 

  Neighborhood(): m_size(0) {}

  void print() const {
/*     int i; */
/*     for ( i = 0; i < m_size; i++ ) { */
/*       cout << "Neighbor " << i << " is " << m_locations[i] << endl; */
/*     } */
    cout << "Neighborhood size is " << m_size << endl;
  }

  int getSize() const { return m_size; }

  Vector3i& operator[](int i) { 
    assert ( i >= 0 && i < MAX_NUM );
    return m_locations[i]; 
  }

  const Vector3i& operator[](int i) const { 
    assert ( i >= 0 && i < MAX_NUM );
    return m_locations[i]; 
  }


  void setSize(int s) { m_size = s; }

  void addNeighbor(const Vector3i& l) {
    m_locations[m_size] = l;
    m_size++;
  }

 private:
  //cannot have more neighbors than the number of samples.
  Vector3i m_locations[MAX_NUM];
  int m_size;
};

class SphericalFilter {

 public:

  SphericalFilter(): m_radius(0.0) {}

  void filter(const Dataset& in, Dataset& out) {
    int numSamples = in.getNumSamples();

    out.copyAttribs(in);

    int i;
    for ( i = 0; i < numSamples; i++ ) {
      out[i] = filter(in, i);
    }
  }

  int getNeighborhoodSize() const { return m_neighborhood.getSize(); }

  void print() const {
    cout << "Filter radius = " << m_radius << endl;
    int numNeighbors = m_neighborhood.getSize();
    cout << "Neighorhood size = " << numNeighbors << endl;
    int i;
    for ( i = 0; i < numNeighbors; i++ ) {
      cout << "Neighbor " << i << " = " << m_neighborhood[i]
           << " has weight " << m_weights[i] << endl;
    }
  }

  static double gaussianFunc(const double& sigma,
                             const double& x,
                             const double& y,
                             const double& z) {
    return exp( -(x*x + y*y + z*z) / (2.0 * sigma*sigma));
  }

  static double gaussianFunc(const double& sigma, const Vector3i& p) {
    double x = double(p[0]);
    double y = double(p[1]);
    double z = double(p[2]);

    return gaussianFunc(sigma, x, y, z);
  }

  double getRadius() const { return m_radius; }

  virtual void setRadius(const double& r) = 0;

 protected:

  virtual double filter(const Dataset& in, int i) = 0;

  Neighborhood m_neighborhood; 


  double m_weights[MAX_NUM]; //normalized


  double m_radius;
};

class SphericalFilterCC : public SphericalFilter {
 public:

  /*virtual*/ void setRadius(const double& r) {
    m_neighborhood.setSize(0); //clear existing neighborhood

    m_radius = 0;

    //construct integral bounding box of neighborhood candidates
    int radius_int = int(r); 
    Vector3i min(-radius_int);
    Vector3i max(radius_int);

    int x, y, z;
    for ( x = min[0]; x <= max[0]; x++ ) {
      for ( y = min[1]; y <= max[1]; y++ ) {
        for ( z = min[2]; z <= max[2]; z++ ) {
          int p_arr[] = {x,y,z};
          Vector3i p(p_arr);
          double norm_p = norm(p);
          if ( norm_p <= r ) {
            m_neighborhood.addNeighbor(p);

            //track actual radius of discrete neighborhood
            if ( norm_p > m_radius ) {
              m_radius = norm_p;
            }
          }
        }
      }
    }

    printf("Radius is %lf, ", m_radius); 
    m_neighborhood.print();
  } 

};

class MedianSphericalFilterCC : public SphericalFilterCC { 
 private:

  /*virtual*/ double filter(const Dataset& in, int index) {
    int numNeighbors = m_neighborhood.getSize();
    Vector3i p = in.mindex_cc(index);

    assert(numNeighbors > 0);

    int i;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = p + m_neighborhood[i];
      int neighbor_index = in.lindex_cc(neighbor);
      m_weights[i] = in[neighbor_index];
    }

    sort(m_weights, m_weights+numNeighbors);

    if ( ODD(numNeighbors) ) { //median is unique
      int median = (numNeighbors - 1) / 2;
      //i.e. numNeighbors = 5, median = 2, where indices = {0,1,2,3,4}
      return m_weights[median];
    }
 
    //else, take average of the two medians
    int m1 = numNeighbors / 2;
    int m2 = m1 - 1;
    //i.e. numNeighbors = 6, median = 2,3, where indices = {0,1,2,3,4,5}
    return 0.5*(m_weights[m1] + m_weights[m2]);
  }
};

class GaussianSphericalFilterCC : public SphericalFilterCC { 
 public:

  /*virtual*/ void setRadius(const double& r) {
    //need to fill Gaussian filter weights in addition to
    //setRadius function of parent

    SphericalFilterCC::setRadius(r);

    //sigma below cannot be 0, else div by 0 in exponent of gaussian
    if ( EQUALS(m_radius, 0.0) ) {
      m_weights[0] = 1.0; 
      return;
    }

    //now that we know the radius, use it to compute sigma
    double N = 3.0;
    double sigma = m_radius / N; 

    //compute weights
    int i;
    int numNeighbors = m_neighborhood.getSize();
    double sum = 0;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = m_neighborhood[i];
      m_weights[i] = gaussianFunc(sigma, neighbor);
      sum += m_weights[i];
    }
    
    //normalize weights so they add up to 1
    for ( i = 0; i < numNeighbors; i++ ) {
      m_weights[i] /= sum;
    }
  }

 private:

  /*virtual*/ double filter(const Dataset& in, int index) {
    int numNeighbors = m_neighborhood.getSize();
    Vector3i p = in.mindex_cc(index);

    assert(numNeighbors > 0);

    double filtered_val = 0.0;
    int i;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = p + m_neighborhood[i];
      int neighbor_index = in.lindex_cc(neighbor);
      filtered_val += in[neighbor_index] * m_weights[i];
    }

    return filtered_val;
  }
};

class SphericalFilterBCC : public SphericalFilter {
 public:
  /*virtual*/ void setRadius(const double& r) { 
    m_neighborhood.setSize(0); //clear existing neighborhood

    m_radius = 0;

    //construct integral bounding box of neighborhood candidates
    int radius_int = int(r); 
    Vector3i min(-radius_int);
    Vector3i max(radius_int);

    int x, y, z;
    for ( x = min[0]; x <= max[0]; x++ ) {
      for ( y = min[1]; y <= max[1]; y++ ) {
        for ( z = min[2]; z <= max[2]; z++ ) {
          int p_arr[] = {x,y,z};
          Vector3i p(p_arr);
          double norm_p = norm(p);
          if ( norm_p <= r ) {
            //test parity for BCC grid
            if ( ( ODD(x) && ODD(y) && ODD(z) ) ||
                 ( EVEN(x) && EVEN(y) && EVEN(z) ) ) {
              m_neighborhood.addNeighbor(p);

              //track actual radius of discrete neighborhood
              if ( norm_p > m_radius ) {
                m_radius = norm_p;
              }
            }
          }
        }
      }
    }

    printf("Radius is %lf, ", m_radius); 
    m_neighborhood.print();
  }   
};

class MedianSphericalFilterBCC : public SphericalFilterBCC {
 private:

  /*virtual*/ double filter(const Dataset& in, int index) {
    int numNeighbors = m_neighborhood.getSize();
    Vector3i p = in.mindex_bcc(index);

    assert(numNeighbors > 0);

    int i;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = p + m_neighborhood[i];
      int neighbor_index = in.lindex_bcc(neighbor);
      m_weights[i] = in[neighbor_index];
    }

    sort(m_weights, m_weights+numNeighbors);

    if ( ODD(numNeighbors) ) { //median is unique
      int median = (numNeighbors - 1) / 2;
      //i.e. numNeighbors = 5, median = 2, where indices = {0,1,2,3,4}
      return m_weights[median];
    }
 
    //else, take average of the two medians
    int m1 = numNeighbors / 2;
    int m2 = m1 - 1;
    //i.e. numNeighbors = 6, median = 2,3, where indices = {0,1,2,3,4,5}
    return 0.5*(m_weights[m1] + m_weights[m2]);
  }
};

class GaussianSphericalFilterBCC : public SphericalFilterBCC { 
 public:

  /*virtual*/ void setRadius(const double& r) {
    //need to fill Gaussian filter weights in addition to
    //setRadius function of parent

    SphericalFilterBCC::setRadius(r);

    //sigma below cannot be 0, else div by 0 in exponent of gaussian
    if ( EQUALS(m_radius, 0.0) ) {
      m_weights[0] = 1.0; 
      return;
    }

    //now that we know the radius, use it to compute sigma
    double N = 3.0;
    double sigma = m_radius / N; 

    //compute weights
    int i;
    int numNeighbors = m_neighborhood.getSize();
    double sum = 0;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = m_neighborhood[i];
      m_weights[i] = gaussianFunc(sigma, neighbor);
      sum += m_weights[i];
    }
    
    //normalize weights so they add up to 1
    for ( i = 0; i < numNeighbors; i++ ) {
      m_weights[i] /= sum;
    }
  }

 private:

  /*virtual*/ double filter(const Dataset& in, int index) {
    int numNeighbors = m_neighborhood.getSize();
    Vector3i p = in.mindex_bcc(index);

    assert(numNeighbors > 0);

    double filtered_val = 0.0;
    int i;
    for ( i = 0; i < numNeighbors; i++ ) {
      Vector3i neighbor = p + m_neighborhood[i];
      int neighbor_index = in.lindex_bcc(neighbor);
      filtered_val += in[neighbor_index] * m_weights[i];
    }

    return filtered_val;
  }
};

#endif

---