From: Mikko Rasa <tdb@tdb.fi>
Date: Wed, 15 May 2013 07:27:06 +0000 (+0300)
Subject: Beginnings of a geometry library
X-Git-Url: http://git.tdb.fi/?a=commitdiff_plain;h=6ff13022b53830d35283905d562c2ef3af198cc1;p=libs%2Fmath.git

Beginnings of a geometry library
---

diff --git a/Build b/Build
index 517bb2d..62e7375 100644
--- a/Build
+++ b/Build
@@ -3,6 +3,7 @@ package "mspmath"
 	library "mspmath"
 	{
 		source "source/linal";
+		source "source/geometry";
 		install true;
 		install_map
 		{
diff --git a/source/geometry/affinetransformation.h b/source/geometry/affinetransformation.h
new file mode 100644
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
diff --git a/source/geometry/angle.h b/source/geometry/angle.h
new file mode 100644
index 0000000..edd7521
--- /dev/null
+++ b/source/geometry/angle.h
@@ -0,0 +1,221 @@
+#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
diff --git a/source/geometry/box.h b/source/geometry/box.h
new file mode 100644
index 0000000..bfcfa2f
--- /dev/null
+++ b/source/geometry/box.h
@@ -0,0 +1,26 @@
+#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
diff --git a/source/geometry/circle.h b/source/geometry/circle.h
new file mode 100644
index 0000000..ebc1c41
--- /dev/null
+++ b/source/geometry/circle.h
@@ -0,0 +1,20 @@
+#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
diff --git a/source/geometry/dummy.cpp b/source/geometry/dummy.cpp
new file mode 100644
index 0000000..885bc5a
--- /dev/null
+++ b/source/geometry/dummy.cpp
@@ -0,0 +1,9 @@
+#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"
diff --git a/source/geometry/hyperbox.h b/source/geometry/hyperbox.h
new file mode 100644
index 0000000..2ade28a
--- /dev/null
+++ b/source/geometry/hyperbox.h
@@ -0,0 +1,79 @@
+#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
diff --git a/source/geometry/hypersphere.h b/source/geometry/hypersphere.h
new file mode 100644
index 0000000..8b7d0a4
--- /dev/null
+++ b/source/geometry/hypersphere.h
@@ -0,0 +1,71 @@
+#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
diff --git a/source/geometry/ray.h b/source/geometry/ray.h
new file mode 100644
index 0000000..4825cfe
--- /dev/null
+++ b/source/geometry/ray.h
@@ -0,0 +1,39 @@
+#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
diff --git a/source/geometry/rectangle.h b/source/geometry/rectangle.h
new file mode 100644
index 0000000..296ba6c
--- /dev/null
+++ b/source/geometry/rectangle.h
@@ -0,0 +1,25 @@
+#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
diff --git a/source/geometry/shape.h b/source/geometry/shape.h
new file mode 100644
index 0000000..528a77c
--- /dev/null
+++ b/source/geometry/shape.h
@@ -0,0 +1,30 @@
+#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
diff --git a/source/geometry/transformedshape.h b/source/geometry/transformedshape.h
new file mode 100644
index 0000000..9a98f31
--- /dev/null
+++ b/source/geometry/transformedshape.h
@@ -0,0 +1,77 @@
+#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
diff --git a/source/linal/matrix4.h b/source/linal/matrix4.h
deleted file mode 100644
index 8b89a65..0000000
--- a/source/linal/matrix4.h
+++ /dev/null
@@ -1,21 +0,0 @@
-#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