Beginnings of a geometry library

[libs/math.git] / source / geometry / affinetransformation.h

diff --git a/source/geometry/affinetransformation.h b/source/geometry/affinetransformation.h
new file mode 100644 (file)
index 0000000..e0ca370
--- /dev/null
+++ b/source/geometry/affinetransformation.h
@@ -0,0 +1,185 @@
+#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