library "mspmath"
{
source "source/linal";
+ source "source/geometry";
install true;
install_map
{
--- /dev/null
+#ifndef MSP_GEOMETRY_AFFINETRANSFORMATION_H_
+#define MSP_GEOMETRY_AFFINETRANSFORMATION_H_
+
+#include <msp/linal/squarematrix.h>
+#include "angle.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class AffineTransformation;
+
+
+/**
+Helper class to provide specialized operations for AffineTransformation.
+*/
+template<typename T, unsigned D>
+class AffineTransformationOps
+{
+protected:
+ AffineTransformationOps() { }
+};
+
+template<typename T>
+class AffineTransformationOps<T, 2>
+{
+protected:
+ AffineTransformationOps() { }
+
+public:
+ static AffineTransformation<T, 2> rotation(const Angle<T> &);
+};
+
+template<typename T>
+class AffineTransformationOps<T, 3>
+{
+protected:
+ AffineTransformationOps() { }
+
+public:
+ static AffineTransformation<T, 3> rotation(const Angle<T> &, const LinAl::Vector<T, 3> &);
+};
+
+
+/**
+An affine transformation in D dimensions. Affine transformations preserve
+straightness of lines and ratios of distances. Angles and distances themselves
+may change. Internally this is represented by a square matrix of size D+1.
+*/
+template<typename T, unsigned D>
+class AffineTransformation: public AffineTransformationOps<T, D>
+{
+ friend class AffineTransformationOps<T, D>;
+
+private:
+ LinAl::SquareMatrix<T, D+1> matrix;
+
+public:
+ AffineTransformation();
+
+ static AffineTransformation<T, D> translation(const LinAl::Vector<T, D> &);
+ static AffineTransformation<T, D> scaling(const LinAl::Vector<T, D> &);
+ static AffineTransformation<T, D> shear(const LinAl::Vector<T, D> &, const LinAl::Vector<T, D> &);
+
+ const LinAl::SquareMatrix<T, D+1> &get_matrix() const { return matrix; }
+ operator const LinAl::SquareMatrix<T, D+1> &() const { return matrix; }
+
+ LinAl::Vector<T, D> transform(const LinAl::Vector<T, D> &) const;
+ LinAl::Vector<T, D> transform_linear(const LinAl::Vector<T, D> &) const;
+};
+
+template<typename T, unsigned D>
+inline AffineTransformation<T, D>::AffineTransformation()
+{
+ this->matrix = LinAl::SquareMatrix<T, D+1>::identity();
+}
+
+
+template<typename T, unsigned D>
+AffineTransformation<T, D> AffineTransformation<T, D>::translation(const LinAl::Vector<T, D> &v)
+{
+ AffineTransformation<T, D> r;
+ for(unsigned i=0; i<D; ++i)
+ r.matrix(D, i) = v[i];
+ return r;
+}
+
+template<typename T, unsigned D>
+AffineTransformation<T, D> AffineTransformation<T, D>::scaling(const LinAl::Vector<T, D> &factors)
+{
+ AffineTransformation<T, D> r;
+ for(unsigned i=0; i<D; ++i)
+ r.matrix(i, i) = factors[i];
+ return r;
+}
+
+template<typename T, unsigned D>
+AffineTransformation<T, D> AffineTransformation<T, D>::shear(const LinAl::Vector<T, D> &normal, const LinAl::Vector<T, D> &shift)
+{
+ AffineTransformation<T, D> r;
+ for(unsigned i=0; i<D; ++i)
+ for(unsigned j=0; j<D; ++j)
+ r.matrix(i, j) += normal[j]*shift[i];
+ return r;
+}
+
+template<typename T>
+AffineTransformation<T, 2> AffineTransformationOps<T, 2>::rotation(const Angle<T> &angle)
+{
+ AffineTransformation<T, 2> r;
+ T c = cos(angle);
+ T s = sin(angle);
+ r.matrix(0, 0) = c;
+ r.matrix(0, 1) = -s;
+ r.matrix(1, 0) = s;
+ r.matrix(1, 1) = c;
+ return r;
+}
+
+template<typename T>
+AffineTransformation<T, 3> AffineTransformationOps<T, 3>::rotation(const Angle<T> &angle, const LinAl::Vector<T, 3> &axis)
+{
+ AffineTransformation<T, 3> r;
+ LinAl::Vector<T, 3> axn = normalize(axis);
+ T c = cos(angle);
+ T s = sin(angle);
+ // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
+ r.matrix(0, 0) = c+axn.x*axn.x*(1-c);
+ r.matrix(0, 1) = axn.x*axn.y*(1-c)-axn.z*s;
+ r.matrix(0, 2) = axn.x*axn.z*(1-c)+axn.y*s;
+ r.matrix(1, 0) = axn.y*axn.x*(1-c)+axn.z*s;
+ r.matrix(1, 1) = c+axn.y*axn.y*(1-c);
+ r.matrix(1, 2) = axn.y*axn.z*(1-c)-axn.x*s;
+ r.matrix(2, 0) = axn.z*axn.x*(1-c)-axn.y*s;
+ r.matrix(2, 1) = axn.z*axn.y*(1-c)+axn.x*s;
+ r.matrix(2, 2) = c+axn.z*axn.z*(1-c);
+ return r;
+}
+
+
+template<typename T, unsigned N>
+inline LinAl::Vector<T, N+1> augment_vector(const LinAl::Vector<T, N> &v, T s)
+{
+ LinAl::Vector<T, N+1> r;
+ for(unsigned i=0; i<N; ++i)
+ r[i] = v[i];
+ r[N] = s;
+ return r;
+}
+
+template<typename T, unsigned N>
+inline LinAl::Vector<T, N-1> reduce_vector(const LinAl::Vector<T, N> &v)
+{
+ LinAl::Vector<T, N-1> r;
+ for(unsigned i=0; i<N-1; ++i)
+ r[i] = v[i];
+ return r;
+}
+
+template<typename T, unsigned N>
+inline LinAl::Vector<T, N-1> divide_vector(const LinAl::Vector<T, N> &v)
+{
+ LinAl::Vector<T, N-1> r;
+ for(unsigned i=0; i<N-1; ++i)
+ r[i] = v[i]/v[N-1];
+ return r;
+}
+
+
+template<typename T, unsigned D>
+inline LinAl::Vector<T, D> AffineTransformation<T, D>::transform(const LinAl::Vector<T, D> &v) const
+{
+ return reduce_vector(matrix*augment_vector(v, T(1)));
+}
+
+template<typename T, unsigned D>
+inline LinAl::Vector<T, D> AffineTransformation<T, D>::transform_linear(const LinAl::Vector<T, D> &v) const
+{
+ return reduce_vector(matrix*augment_vector(v, T(0)));
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_ANGLE_H_
+#define MSP_GEOMETRY_ANGLE_H_
+
+#include <cmath>
+
+namespace Msp {
+namespace Geometry {
+
+/**
+A planar angle. Creating an Angle from a raw value or extracting it requires
+specifying whether degrees or radians are used. This eliminates a common
+source of errors.
+*/
+template<typename T>
+class Angle
+{
+private:
+ T value;
+
+ explicit Angle(T);
+public:
+ template<typename U>
+ Angle(const Angle<U> &);
+
+ static Angle from_degrees(T);
+ static Angle from_radians(T);
+
+ static Angle right();
+ static Angle quarter_circle();
+ static Angle straight();
+ static Angle half_circle();
+ static Angle full_circle();
+
+ T degrees() const;
+ T radians() const;
+
+ Angle &operator+=(const Angle &);
+ Angle &operator-=(const Angle &);
+ Angle &operator*=(T);
+ Angle &operator/=(T);
+};
+
+template<typename T>
+inline Angle<T>::Angle(T v):
+ value(v)
+{ }
+
+template<typename T>
+template<typename U>
+inline Angle<T>::Angle(const Angle<U> &other):
+ value(other.value)
+{ }
+
+template<typename T>
+inline Angle<T> Angle<T>::from_degrees(T d)
+{
+ return Angle<T>(d*M_PI/180);
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::from_radians(T r)
+{
+ return Angle<T>(r);
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::right()
+{
+ return from_radians(M_PI/2);
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::quarter_circle()
+{
+ return right();
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::straight()
+{
+ return from_radians(M_PI);
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::half_circle()
+{
+ return straight();
+}
+
+template<typename T>
+inline Angle<T> Angle<T>::full_circle()
+{
+ return from_radians(M_PI*2);
+}
+
+template<typename T>
+inline T Angle<T>::degrees() const
+{
+ return value*180/M_PI;
+}
+
+template<typename T>
+inline T Angle<T>::radians() const
+{
+ return value;
+}
+
+template<typename T>
+inline Angle<T> &Angle<T>::operator+=(const Angle &a)
+{
+ value += a.value;
+ return *this;
+}
+
+template<typename T>
+inline Angle<T> operator+(const Angle<T> &a1, const Angle<T> &a2)
+{
+ Angle<T> r(a1);
+ return r += a2;
+}
+
+template<typename T>
+inline Angle<T> &Angle<T>::operator-=(const Angle &a)
+{
+ value -= a.value;
+ return *this;
+}
+
+template<typename T>
+inline Angle<T> operator-(const Angle<T> &a1, const Angle<T> &a2)
+{
+ Angle<T> r(a1);
+ return r -= a2;
+}
+
+template<typename T>
+inline Angle<T> &Angle<T>::operator*=(T s)
+{
+ value *= s;
+ return *this;
+}
+
+template<typename T>
+inline Angle<T> operator*(const Angle<T> &a, T s)
+{
+ Angle<T> r(a);
+ return r *= s;
+}
+
+template<typename T>
+inline Angle<T> operator*(T s, const Angle<T> &a)
+{
+ return a*s;
+}
+
+template<typename T>
+inline Angle<T> &Angle<T>::operator/=(T s)
+{
+ value /= s;
+ return *this;
+}
+
+template<typename T>
+inline Angle<T> operator/(const Angle<T> &a, T s)
+{
+ Angle<T> r(a);
+ return r /= s;
+}
+
+template<typename T>
+inline Angle<T> operator/(T s, const Angle<T> &a)
+{
+ return a/s;
+}
+
+template<typename T>
+inline T sin(const Angle<T> &angle)
+{
+ return std::sin(angle.radians());
+}
+
+template<typename T>
+inline Angle<T> asin(T s)
+{
+ return Angle<T>::from_radians(std::asin(s));
+}
+
+template<typename T>
+inline T cos(const Angle<T> &angle)
+{
+ return std::cos(angle.radians());
+}
+
+template<typename T>
+inline Angle<T> acos(T s)
+{
+ return Angle<T>::from_radians(std::acos(s));
+}
+
+template<typename T>
+inline T tan(const Angle<T> &angle)
+{
+ return std::tan(angle.radians());
+}
+
+template<typename T>
+inline Angle<T> atan(T s)
+{
+ return Angle<T>::from_radians(std::atan(s));
+}
+
+template<typename T>
+inline Angle<T> atan2(T y, T x)
+{
+ return Angle<T>::from_radians(std::atan2(y, x));
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_BOX_H_
+#define MSP_GEOMETRY_BOX_H_
+
+#include <msp/linal/vector3.h>
+#include "hyperbox.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T>
+class Box: public HyperBox<T, 3>
+{
+public:
+ Box() { }
+ explicit Box(const LinAl::Vector<T, 3> &d): HyperBox<T, 3>(d) { }
+ Box(T w, T h, T d): HyperBox<T, 3>(LinAl::Vector3<T>(w, h, d)) { }
+
+ T get_width() const { return this->get_dimension(0); }
+ T get_height() const { return this->get_dimension(1); }
+ T get_depth() const { return this->get_dimension(2); }
+};
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_CIRCLE_H_
+#define MSP_GEOMETRY_CIRCLE_H_
+
+#include "hypersphere.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T>
+class Circle: public HyperSphere<T, 2>
+{
+public:
+ Circle() { }
+ explicit Circle(T r): HyperSphere<T, 2>(r) { }
+};
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#include "affinetransformation.h"
+#include "angle.h"
+#include "box.h"
+#include "circle.h"
+#include "hyperbox.h"
+#include "hypersphere.h"
+#include "rectangle.h"
+#include "shape.h"
+#include "transformedshape.h"
--- /dev/null
+#ifndef MSP_GEOMETRY_HYPERBOX_H_
+#define MSP_GEOMETRY_HYPERBOX_H_
+
+#include <msp/linal/vector.h>
+#include "ray.h"
+#include "shape.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class HyperBox: public Shape<T, D>
+{
+private:
+ LinAl::Vector<T, D> dimensions;
+
+public:
+ HyperBox();
+ explicit HyperBox(const LinAl::Vector<T, D> &);
+
+ virtual HyperBox *clone() const;
+
+ const LinAl::Vector<T, D> &get_dimensions() const { return dimensions; }
+ T get_dimension(unsigned) const;
+
+ virtual HyperBox<T, D> get_axis_aligned_bounding_box() const { return *this; }
+ virtual bool check_intersection(const Ray<T, D> &) const;
+};
+
+template<typename T, unsigned D>
+inline HyperBox<T, D>::HyperBox()
+{
+ for(unsigned i=0; i<D; ++i)
+ dimensions[i] = 1;
+}
+
+template<typename T, unsigned D>
+inline HyperBox<T, D>::HyperBox(const LinAl::Vector<T, D> &d):
+ dimensions(d)
+{ }
+
+template<typename T, unsigned D>
+inline HyperBox<T, D> *HyperBox<T, D>::clone() const
+{
+ return new HyperBox<T, D>(dimensions);
+}
+
+template<typename T, unsigned D>
+inline T HyperBox<T, D>::get_dimension(unsigned i) const
+{
+ return dimensions[i];
+}
+
+template<typename T, unsigned D>
+inline bool HyperBox<T, D>::check_intersection(const Ray<T, D> &ray) const
+{
+ LinAl::Vector<T, D> half_dim = dimensions/T(2);
+ for(unsigned i=0; i<D; ++i)
+ for(int j=-1; j<=1; j+=2)
+ {
+ T x = (T(j)*half_dim[i]-ray.get_start()[i])/ray.get_direction()[i];
+ if(x>0)
+ {
+ LinAl::Vector<T, D> p = ray.get_start()+ray.get_direction()*x;
+ bool inside = true;
+ for(unsigned k=0; (inside && k<D); ++k)
+ inside = (k==i || (p[k]>=-half_dim[k] && p[k]<half_dim[k]));
+ if(inside)
+ return true;
+ }
+ }
+
+ return false;
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_HYPERSPHERE_H_
+#define MSP_GEOMETRY_HYPERSPHERE_H_
+
+#include <msp/linal/vector.h>
+#include "hyperbox.h"
+#include "ray.h"
+#include "shape.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class HyperSphere: public Shape<T, D>
+{
+private:
+ T radius;
+
+public:
+ HyperSphere();
+ explicit HyperSphere(T);
+
+ virtual HyperSphere *clone() const;
+
+ T get_radius() const { return radius; }
+
+ virtual HyperBox<T, D> get_axis_aligned_bounding_box() const;
+ virtual bool check_intersection(const Ray<T, D> &) const;
+};
+
+template<typename T, unsigned D>
+inline HyperSphere<T, D>::HyperSphere():
+ radius(1)
+{ }
+
+template<typename T, unsigned D>
+inline HyperSphere<T, D>::HyperSphere(T r):
+ radius(r)
+{ }
+
+template<typename T, unsigned D>
+inline HyperSphere<T, D> *HyperSphere<T, D>::clone() const
+{
+ return new HyperSphere<T, D>(radius);
+}
+
+template<typename T, unsigned D>
+inline HyperBox<T, D> HyperSphere<T, D>::get_axis_aligned_bounding_box() const
+{
+ LinAl::Vector<T, D> dimensions;
+ for(unsigned i=0; i<D; ++i)
+ dimensions[i] = radius;
+ return HyperBox<T, D>(dimensions);
+}
+
+template<typename T, unsigned D>
+inline bool HyperSphere<T, D>::check_intersection(const Ray<T, D> &ray) const
+{
+ T x = inner_product(ray.get_direction(), ray.get_start());
+ if(x>0)
+ return inner_product(ray.get_start(), ray.get_start())<=radius*radius;
+ else
+ {
+ LinAl::Vector<T, D> nearest = ray.get_start()-ray.get_direction()*x;
+ return inner_product(nearest, nearest)<=radius*radius;
+ }
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_RAY_H_
+#define MSP_GEOMETRY_RAY_H_
+
+#include <msp/linal/vector.h>
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class Ray
+{
+private:
+ LinAl::Vector<T, D> start;
+ LinAl::Vector<T, D> direction;
+
+public:
+ Ray();
+ Ray(const LinAl::Vector<T, D> &, const LinAl::Vector<T, D> &);
+
+ const LinAl::Vector<T, D> &get_start() const { return start; }
+ const LinAl::Vector<T, D> &get_direction() const { return direction; }
+};
+
+template<typename T, unsigned D>
+Ray<T, D>::Ray()
+{
+ direction[0] = 1;
+}
+
+template<typename T, unsigned D>
+Ray<T, D>::Ray(const LinAl::Vector<T, D> &s, const LinAl::Vector<T, D> &d):
+ start(s),
+ direction(normalize(d))
+{ }
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_RECTANGLE_H_
+#define MSP_GEOMETRY_RECTANGLE_H_
+
+#include <msp/linal/vector2.h>
+#include "hyperbox.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T>
+class Rectangle: public HyperBox<T, 2>
+{
+public:
+ Rectangle() { }
+ explicit Rectangle(const LinAl::Vector<T, 2> &d): HyperBox<T, 2>(d) { }
+ Rectangle(T w, T h): HyperBox<T, 2>(LinAl::Vector2<T>(w, h)) { }
+
+ T get_width() const { return this->get_dimension(0); }
+ T get_height() const { return this->get_dimension(1); }
+};
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_SHAPE_H_
+#define MSP_GEOMETRY_SHAPE_H_
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class HyperBox;
+
+template<typename T, unsigned D>
+class Ray;
+
+template<typename T, unsigned D>
+class Shape
+{
+protected:
+ Shape() { }
+public:
+ virtual ~Shape() { }
+
+ virtual Shape *clone() const = 0;
+
+ virtual HyperBox<T, D> get_axis_aligned_bounding_box() const = 0;
+ virtual bool check_intersection(const Ray<T, D> &) const = 0;
+};
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
--- /dev/null
+#ifndef MSP_GEOMETRY_TRANSFORMEDSHAPE_H_
+#define MSP_GEOMETRY_TRANSFORMEDSHAPE_H_
+
+#include "affinetransformation.h"
+#include "ray.h"
+#include "shape.h"
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class TransformedShape
+{
+private:
+ Shape<T, D> *shape;
+ AffineTransformation<T, D> transformation;
+
+public:
+ TransformedShape(const Shape<T, D> &, const AffineTransformation<T, D> &);
+ TransformedShape(const TransformedShape &);
+ TransformedShape &operator=(const TransformedShape &);
+ ~TransformedShape();
+
+ virtual TransformedShape *clone() const;
+
+ const Shape<T, D> &get_shape() const { return *shape; }
+ const AffineTransformation<T, D> &get_transformation() const { return transformation; }
+
+ virtual bool check_intersection(const Ray<T, D> &) const;
+};
+
+template<typename T, unsigned D>
+inline TransformedShape<T, D>::TransformedShape(const Shape<T, D> &s, const AffineTransformation<T, D> &t):
+ shape(s.clone()),
+ transformation(t)
+{ }
+
+template<typename T, unsigned D>
+inline TransformedShape<T, D>::TransformedShape(const TransformedShape &other):
+ shape(other.shape->clone()),
+ transformation(other.transformation)
+{ }
+
+template<typename T, unsigned D>
+inline TransformedShape<T, D> &TransformedShape<T, D>::operator=(const TransformedShape<T, D> &other)
+{
+ delete shape;
+ shape = other.shape->clone();
+ transformation = other.transformation();
+}
+
+template<typename T, unsigned D>
+inline TransformedShape<T, D>::~TransformedShape()
+{
+ delete shape;
+}
+
+template<typename T, unsigned D>
+inline TransformedShape<T, D> *TransformedShape<T, D>::clone() const
+{
+ return new TransformedShape<T, D>(*this);
+}
+
+template<typename T, unsigned D>
+inline bool TransformedShape<T, D>::check_intersection(const Ray<T, D> &ray) const
+{
+ // TODO cache the inverse transformation for performance
+ LinAl::SquareMatrix<T, D+1> inverse_trans = LinAl::invert(transformation.get_matrix());
+ Ray<T, D> trans_ray(reduce_vector(inverse_trans*augment_vector(ray.get_start(), T(1))),
+ reduce_vector(inverse_trans*augment_vector(ray.get_direction(), T(0))));
+ return shape->check_intersection(trans_ray);
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif
+++ /dev/null
-#ifndef MSP_LINAL_MATRIX4_H_
-#define MSP_LINAL_MATRIX4_H_
-
-#include "squarematrix.h"
-
-namespace Msp {
-namespace LinAl {
-
-/**
-A 4x4 square matrix, capable of expressing affine transformations in a
-three-dimensional vector space.
-*/
-template<typename T>
-class Matrix4: public SquareMatrix<T, 4>
-{
-};
-
-} // namespace LinAl
-} // namespace Msp
-
-#endif
projects / libs / math.git / commitdiff