Point.h

Go to the documentation of this file.

#ifndef Point_HEADER
#define Point_HEADER

#include <boost/array.hpp>
#include <math.h> // sqrt
#include <limits> // numeric_limits

namespace Geometry {

template <size_t d>
class Point : public boost::array<double, d> {
  private:
    bool operator==(Point other);

  public:
    typedef boost::array<double, d> super;
    Point() { }
    Point(double v) { assign(v); }
    Point(const super& t) : super(t) { }
    Point(const Point& p) : super(p) { }
    Point(const double *v, size_t
#ifndef NDEBUG
          n
#endif
          ) {
      assert(n >= d);
      for(unsigned i = 0; i < d; i++) {
        (*this)[i] = v[i];
      }
    }

    template <class iterator>
    Point(iterator begin, const iterator&
#ifndef NDEBUG
        end
#endif
        ) {
      assert(end - begin == d);
      for(unsigned i = 0; i < d; i++, begin++) {
        (*this)[i] = *begin;
      }
    }

    /* This constructor is appealing, but then it
       destroys typechecking between points of
       different dimension.  Leave it dead. */
    /* template <class container> Point(const container& c); */

    /* gcc seems to be stupid about inlining and constant-propagating through
       array::assign (which calls fill_n).  So we define it here, where hopefully
       gcc will notice the loop has known length. */
    void assign(double v) {
      for(unsigned i = 0; i < d; ++i) {
        this->elems[i] = v;
      }
    }

    static unsigned dim() { return d; }
    static const unsigned dimension = d;

    //unit vector in the j dimension
    static Point unit(unsigned j) {
      Point rval(0.0);
      rval[j] = 1.0;
      return rval;
    }

    // the all ones vector
    static Point allones() {
      Point rval(1.0);
      return rval;
    }

    // a random unit vector drawn u.a.r. from the cube
    // and normalized.  TODO: This is not u.a.r. from the
    // sphere.
    static Point unitrand() {
      Point p; // uninitialized
      for(unsigned i = 0; i < d; ++i) {
        p[i] = random();
      }
      double len = sqrt(p.norm2());
      p *= 1.0 / len;
      return p;
    }
    static Point randPoint(double max) {
      Point p = unitrand();
      double r = (double)random() / (double)std::numeric_limits<long>::max();
      assert(r < 1.0); // may flip vector; whatever.
      p *= r * max;
      return p;
    }

    void operator+= (const Point& p) {
      for (unsigned i = 0; i < d; ++i) {
        (*this)[i] += p[i];
      }
    }

    Point operator+ (const Point& p) const {
      Point copy = *this; copy += p; return copy;
    }

    void operator-= (const Point& p) {
      for (unsigned i = 0; i < d; ++i) {
        (*this)[i] -= p[i];
      }
    }

    Point operator- (const Point& p) const {
      Point copy = *this; copy -= p; return copy;
    }

    void operator*= (double x) {
      for (unsigned i = 0; i < d; ++i) {
        (*this)[i] *= x;
      }
    }

    Point operator* (double x) const {
      Point copy = *this; copy *= x; return copy;
    }

    void operator/= (double x) {
      for (unsigned i = 0; i < d; ++i) {
        (*this)[i] /= x;
      }
    }

    Point operator/ (double x) const {
      Point copy = *this; copy /= x; return copy;
    }

    double dot(const Point& p) const {
      double sum = 0.0;
      for(unsigned i = 0; i < d; i++) {
        sum += (*this)[i] * p[i];
      }
      return sum;
    }

    double dist2(const Point& q) const {
      double sum = 0.0;
      for(unsigned i = 0; i < d; i++) {
        double diff = (*this)[i] - q[i];
        sum += diff * diff;
      }
      return sum;
    }

    double norm2() const {
      double sum = 0.0;
      for(unsigned i = 0; i < d; i++) {
        sum += (*this)[i] * (*this)[i];
      }
      return sum;
    }

    double norm1() const {
      double sum = 0.0;
      for(unsigned i = 0; i < d; i++) {
        sum += (*this)[i];
      }
      return sum;
    }

    bool bitwiseEqual(const Point& other) const {
      return std::equal(this->begin(), this->end(), other.begin());
    }

    std::string toString() const {
      std::string s = "(";
      for (unsigned i = 0; i < d; i++) {
        char buffer[4000];
        snprintf(buffer, sizeof(buffer), " %g", (*this)[i]);
        s += buffer;
      }
      s += " )";
      return s;
    }
};

inline Point<3> crossProduct(const Point<3>& a, const Point<3>& b) {
  Point<3> p;
  p[0] = a[1]*b[2] - a[2]*b[1];
  p[1] = a[2]*b[0] - a[0]*b[2];
  p[2] = a[0]*b[1] - a[1]*b[0];
  return p;
}


}

#endif

---