Dataset.h

#ifndef __DATASET_H__
#define __DATASET_H__

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

#include <assert.h>

#include <string>
using namespace std;

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

#define MAX_NUM 2000000

enum GridType {CC, BCC};

class Dataset {
public:

  Dataset() {
    m_gridType = CC;
    m_dim[0] = 0;
    m_dim[1] = 0;
    m_dim[2] = 0;
    m_scaling = vuVector(1.);
  }

  Dataset(GridType g) {
    m_gridType = g;
    m_dim[0] = 0;
    m_dim[1] = 0;
    m_dim[2] = 0;
    m_scaling = vuVector(1.);
  }

  int getNumSamples() const {
    return m_dim[0] * m_dim[1] * m_dim[2];
  }

  double& operator[](int i) {
    assert(i >= 0 && i < getNumSamples());
    return m_data[i];
  }

  const double& operator[](int i) const {
    assert(i >= 0 && i < getNumSamples());
    return m_data[i];
  }

  int getXDim() const { return m_dim[0]; }

  int getYDim() const { return m_dim[1]; }

  int getZDim() const { return m_dim[2]; }

  Vector3i getDims() const { return Vector3i(m_dim); }

  vuVector getScale() const { return m_scaling; }

  GridType getGridType() const { return m_gridType; }

  int lindex_cc(int x, int y, int z) const {
    //index = x + y*xdim + z*xdim*ydim  
    int xdim = getXDim();
    int ydim = getYDim();
    int zdim = getZDim();

    //clamp invalid values
    if ( x >= xdim ) x = xdim - 1; 
    if ( x < 0 ) x = 0;
    if ( y >= ydim ) y = ydim - 1;
    if ( y < 0 ) y = 0;
    if ( z >= zdim ) z = zdim - 1;
    if ( z < 0 ) z = 0;

    return x + y*xdim + z*xdim*ydim;
  }

  int lindex_cc(const Vector3i& p) const { return this->lindex_cc(p[0], p[1], p[2]); }

  Vector3i mindex_cc(int i) const {
    Vector3i mi;
    int x = getXDim();
    int y = getYDim();
    mi[0] = i % (x * y) % x;
    mi[1] = i % (x * y) / x;
    mi[2] = i / (x * y); 
    return mi;
  }

  static Vector3i index_bcc_to_cc(int x, int y, int z) {
    // suppose step size in x,y direction is 2, in z direction is 1
    // for BCC lattice
    //assert points are all of the same parity.
    assert( (ODD(x) && ODD(y) && ODD(z)) ||
            (EVEN(x) && EVEN(y) && EVEN(z)) );

    //First observe z, if it is odd, shift x,y back down
    if ( ODD(z) ) {
      x--;
      y--;
    }

    //Now, both x & y must be even, so, divide each by two so step
    //size becomes 1.
    assert(EVEN(x) && EVEN(y));
    x >>= 1; //Shift right 1 bit == divide by 2
    y >>= 1;

    return makeVec<int>(x, y, z);
  }

  int lindex_bcc(const Vector3i& index_bcc) const {
    Vector3i index_cc = Dataset::index_bcc_to_cc(index_bcc[0], index_bcc[1], index_bcc[2]);
    return lindex_cc(index_cc);
  }

  static Vector3i index_cc_to_bcc(int x, int y, int z) {
    x <<= 1; //Shift left 1 bit == mult by 2
    y <<= 1;

    //observe z, if it is odd, shift x,y up
    if ( ODD(z) ) {
      x++;
      y++;
    }

    assert( (ODD(x) && ODD(y) && ODD(z)) ||
            (EVEN(x) && EVEN(y) && EVEN(z)) );

    return makeVec<int>(x, y, z);
  }

  Vector3i mindex_bcc(int index) const {
    Vector3i index_cc = mindex_cc(index);
    return Dataset::index_cc_to_bcc(index_cc[0], index_cc[1], index_cc[2]);
  }

  void print() const;

  void read(const string& fullname);

  void write(const string& outDir, const string& filename, const string& ext) const;

  void setDim(int xdim, int ydim, int zdim) {
    m_dim[0] = xdim;
    m_dim[1] = ydim;
    m_dim[2] = zdim;
    assert(getNumSamples() <= MAX_NUM);
  }

  void setGridType(GridType g) { m_gridType = g; }

  void setScaling(const vuVector& s) { m_scaling = s; }

  void copyAttribs(const Dataset& d);

  static void whiteNoiseError(const Dataset& orig,
                              const Dataset& noisy,
                              double& mean,
                              double& stdDev) {
    int numSamples = orig.getNumSamples();
    assert(numSamples == noisy.getNumSamples());
    assert(numSamples > 0 );
    double* diff = new double[numSamples];

    int i;
    for ( i = 0; i < numSamples; i++ ) {
      diff[i] = orig.m_data[i] - noisy.m_data[i];
    }
    
    mean = sampleMean(diff, numSamples);
    stdDev = sampleStdDev(diff, numSamples);

    delete [] diff;
    diff = NULL;
  }

  static void saltPepperError(const Dataset& orig,
                              const Dataset& noisy,
                              double& percentage) {
    //possible to undercount
    //if orig val = 0, and salted val = 0, not counted
    int numSamples = orig.getNumSamples();
    assert(numSamples == noisy.getNumSamples());
    assert(numSamples > 0 );

    int i;
    int numSalted = 0;
    for ( i = 0; i < numSamples; i++ ) {
      double nv = noisy.m_data[i];
      if ( orig.m_data[i] != nv ) {
        assert(EQUALS(nv, 0) || EQUALS(nv, MAX_INTENSITY));
        numSalted++;
      }
    }
    
    percentage = double(numSalted)/double(numSamples);
  }


  static void denoiseErrorWithBoundaries(const Dataset& orig,
                                         const Dataset& denoised,
                                         double& mean,
                                         double& stdDev) {
    whiteNoiseError(orig, denoised, mean, stdDev);
  }

  void getBoundingBox(Vector3i& min, Vector3i& max) const {
    min = Vector3i(0);
    if ( m_gridType == CC ) {
      max[0] = getXDim() - 1;
      max[1] = getYDim() - 1;
      max[2] = getZDim() - 1;
    }
    else {
      max[0] = getXDim()*2 - 1;
      max[1] = getYDim()*2 - 1;
      max[2] = getZDim() - 1;
    }
  }

  bool isSampleInBoundary( int index, const double& filterRadius ) const {
    Vector3i min, max;
    getBoundingBox(min, max);

    assert( filterRadius >= 0);

    int padding = myCeil(filterRadius);
    min += padding;
    max -= padding;

    Vector3i mindex;
    if ( m_gridType == CC ) {
      mindex = mindex_cc(index);
    }
    else {
      mindex = mindex_bcc(index);
    }

    if ( mindex <= max && mindex >= min ) {
      return false;
    }
    return true;
  }


  static void denoiseErrorWithoutBoundaries(const Dataset& orig,
                                            const Dataset& denoised,
                                            double& mean,
                                            double& stdDev,
                                            const double& filterRadius) {
    int numSamples = orig.getNumSamples();
    assert(numSamples == denoised.getNumSamples());
    assert(numSamples > 0 );
    double* diff = new double[numSamples];

    int i;
    int numSamplesSelected = 0;
    for ( i = 0; i < numSamples; i++ ) {
      if ( !orig.isSampleInBoundary(i, filterRadius) ) {
        diff[numSamplesSelected] = orig.m_data[i] - denoised.m_data[i];
        numSamplesSelected++;
      }
    }

    mean = sampleMean(diff, numSamplesSelected);
    stdDev = sampleStdDev(diff, numSamplesSelected);

    delete [] diff;
    diff = NULL;    
  }

private:

  void writeBytes(const string& outDir, const string& filename, const string& ext) const;
  void writeDoubles(const string& outDir, const string& filename) const;
  void readDoubles(const string& fullname);

  GridType m_gridType;
  double m_data[MAX_NUM];
  int m_dim[3]; //num samples along x, y, z axis
  vuVector m_scaling; //scaling factor
};

#endif

---