-#ifndef MSP_GEOMETRY_AFFINETRANSFORMATION_H_
-#define MSP_GEOMETRY_AFFINETRANSFORMATION_H_
-
-#include <msp/linal/matrix.h>
-#include "angle.h"
-#include "boundingbox.h"
-#include "ray.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::Matrix<T, D+1, 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> &);
-
- AffineTransformation &operator*=(const AffineTransformation &);
- AffineTransformation &invert();
-
- const LinAl::Matrix<T, D+1, D+1> &get_matrix() const { return matrix; }
- operator const LinAl::Matrix<T, D+1, 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;
- Ray<T, D> transform(const Ray<T, D> &) const;
- BoundingBox<T, D> transform(const BoundingBox<T, D> &) const;
-};
-
-template<typename T, unsigned D>
-inline AffineTransformation<T, D>::AffineTransformation()
-{
- this->matrix = LinAl::Matrix<T, D+1, 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(i, D) = 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 D>
-inline AffineTransformation<T, D> &AffineTransformation<T, D>::operator*=(const AffineTransformation<T, D> &other)
-{
- matrix *= other.get_matrix();
- return *this;
-}
-
-template<typename T, unsigned D>
-inline AffineTransformation<T, D> operator*(const AffineTransformation<T, D> &at1, const AffineTransformation<T, D> &at2)
-{
- AffineTransformation<T, D> r = at1;
- return r *= at2;
-}
-
-template<typename T, unsigned D>
-inline AffineTransformation<T, D> &AffineTransformation<T, D>::invert()
-{
- matrix.invert();
- return *this;
-}
-
-template<typename T, unsigned D>
-inline AffineTransformation<T, D> invert(const AffineTransformation<T, D> &at)
-{
- AffineTransformation<T, D> r = at;
- return r.invert();
-}
-
-template<typename T, unsigned D>
-inline LinAl::Vector<T, D> AffineTransformation<T, D>::transform(const LinAl::Vector<T, D> &v) const
-{
- return (matrix*compose(v, T(1))).template slice<D>(0);
-}
-
-template<typename T, unsigned D>
-inline LinAl::Vector<T, D> AffineTransformation<T, D>::transform_linear(const LinAl::Vector<T, D> &v) const
-{
- return (matrix*compose(v, T(0))).template slice<D>(0);
-}
-
-template<typename T, unsigned D>
-inline Ray<T, D> AffineTransformation<T, D>::transform(const Ray<T, D> &ray) const
-{
- LinAl::Vector<T, D> dir = transform_linear(ray.get_direction());
- return Ray<T, D>(transform(ray.get_start()), dir, ray.get_limit()*dir.norm());
-}
-
-template<typename T, unsigned D>
-inline BoundingBox<T, D> AffineTransformation<T, D>::transform(const BoundingBox<T, D> &bbox) const
-{
- LinAl::Vector<T, D> min_pt;
- LinAl::Vector<T, D> max_pt;
- for(unsigned i=0; i<(1<<D); ++i)
- {
- LinAl::Vector<T, D> point;
- for(unsigned j=0; j<D; ++j)
- point[j] = ((i>>j)&1 ? bbox.get_maximum_coordinate(j) : bbox.get_minimum_coordinate(j));
-
- point = transform(point);
-
- if(i==0)
- {
- min_pt = point;
- max_pt = point;
- }
- else
- {
- for(unsigned j=0; j<D; ++j)
- {
- min_pt[j] = std::min(min_pt[j], point[j]);
- max_pt[j] = std::max(max_pt[j], point[j]);
- }
- }
- }
-
- return BoundingBox<T, D>(min_pt, max_pt);
-}
-
-} // namespace Geometry
-} // namespace Msp
-
-#endif