ProjectionPlane.h

Go to the documentation of this file.

#ifndef ProjectionPlane_HEADER
#define ProjectionPlane_HEADER

#include <math.h>
#include <boost/array.hpp>
#include "Point.h"
#include "Matrix.h"

namespace Geometry {

template <unsigned ambient, unsigned topological, bool b = (ambient==topological)>
class ProjectionPlane;

/***************************************************************************
 * Easy case: if the two dimensions are the same, there is no projection
 * to do.
 ***************************************************************************/
template <unsigned d> class ProjectionPlane<d, d, true> {
private:
typedef Geometry::Point<d> Point;
public:
template <class iterator> ProjectionPlane(iterator, iterator) { }
template <class iterator> ProjectionPlane(iterator, iterator, int) { }
const Point& project(const Point& p) const { return p; }
Point project1(Point p) const { return p; }
const Point& unproject(const Point& p) const { return p; }
std::string toString() const { return "identity-plane"; }
};



/***************************************************************************
 * The interesting case.
 ***************************************************************************/
template <unsigned ambient, unsigned topological>
class ProjectionPlane<ambient, topological, false> {
public:
typedef Geometry::Point<ambient> Point;
typedef Geometry::Point<topological> PPoint; // projected point

private:
boost::array<Point, topological> basis_; 
Point centroid_;

/***************************************************************************
 * Compute the best projection plane of a set of points.
 ***************************************************************************/
private:
template <class point_iterator>
void init(const point_iterator& begin, const point_iterator& end) {
  assert(ambient > topological);

  unsigned num_pts = 0;

  // Compute the centroid and the moment matrix.
  centroid_.assign(0.0);
  BOUNDED_FOREACH(Point p, begin, end) {
    centroid_ += p;
    num_pts++;
  }

  assert(num_pts >= topological + 1);
  centroid_ /= (double)num_pts; 

  Matrix<ambient, ambient> moment(0.0);
  BOUNDED_FOREACH(Point p, begin, end) {
    p -= centroid_;
    moment += Matrix<ambient, ambient>(p, p);
  }

  // Copy the largest eigenvectors into the basis.
  moment.template eigenvectors<topological>(basis_);

#ifndef NDEBUG
  // Check the basis is actually (almost) orthonormal.
  bool broken = false;
  for(unsigned i = 0; i < topological; ++i) {
    broken |= (fabs(basis_[i].norm2() - 1.0) > 1e-6); // normal?
    for(unsigned j = i + 1; j < topological; ++j) {
      broken |= (fabs(basis_[i].dot(basis_[j])) > 1e-6); // ortho
    }
  }
  if (broken) {
    fprintf(stderr, "Broken projection plane: %s\n", toString().c_str());
    assert(!broken);
  }
#endif
}

public:
template <class point_iterator>
ProjectionPlane(const point_iterator& begin, const point_iterator& end) {
  init(begin, end);
}

/***************************************************************************
 * Project down to the coordinate plane.
 ***************************************************************************/
PPoint project(const Point& p) const {
  PPoint pp;  // uninitialized
  Point p2 = p - centroid_ ;

  for(unsigned i = 0; i < topological; ++i) {
    pp[i] = p2.dot(basis_[i]);
  }
  return pp;
}

// necessary for STL functional stuff
PPoint project1(Point p) const {
  return project(p);
}

/***************************************************************************
 * Project up from the coordinate plane to the real plane.
 ***************************************************************************/
Point unproject(const PPoint& pp) const {
  /*
    pp = centroid + \sum_{i\in topo} ( pp[i] * basis_[i] )
   */

  Point unprojected = centroid_;

  for(unsigned i = 0; i < topological; ++i) {
    unprojected += (basis_[i]) * (pp[i]);
  }

  return unprojected;
}

std::string toString() const {
  char buffer[4096];
  snprintf(buffer, sizeof(buffer), "dim %u, centroid %s, basis << ",
      unsigned(topological), centroid_.toString().c_str());
  std::string s = buffer;
  BOOST_FOREACH(const Point& p, basis_) {
    s += p.toString().c_str();
    s += ' ';
  }
  s += ">>";
  return s;
}

/***************************************************************************
 * FIN
 ***************************************************************************/
};

}

#endif

---