Merge pull request #65 from eggrobin/affine-spaces

Affine maps
mockingbirdnest · May 25, 2014 · 18d2ce5 · 18d2ce5
2 parents 5918814 + 9682920
commit 18d2ce5
Show file tree

Hide file tree

Showing 12 changed files with 395 additions and 5 deletions.
diff --git a/geometry/affine_map.hpp b/geometry/affine_map.hpp
@@ -0,0 +1,45 @@
+#pragma once
+
+#include "geometry/point.hpp"
+#include "geometry/grassmann.hpp"
+
+namespace principia {
+namespace geometry {
+
+template<typename FromFrame, typename ToFrame, typename Scalar,
+         template<typename, typename> class LinearMap>
+class AffineMap {
+ public:
+  typedef Vector<Scalar, FromFrame> FromVector;
+  typedef Vector<Scalar, ToFrame> ToVector;
+  // The identity map.
+  AffineMap() = default;
+  AffineMap(Point<FromVector> const& from_origin,
+            Point<ToVector> const& to_origin,
+            LinearMap<FromFrame, ToFrame> const& linear_map);
+
+  AffineMap<ToFrame, FromFrame, Scalar, LinearMap> Inverse() const;
+  Point<ToVector> operator()(Point<FromVector> const& point) const;
+
+ private:
+  // The map is internally represented as x -> linear_map_(x) + translation_.
+  ToVector translation_;
+  LinearMap<FromFrame, ToFrame> linear_map_;
+
+  template<typename FromFrame, typename ThroughFrame, typename ToFrame,
+           typename Scalar, template<typename, typename> class LinearMap>
+  friend AffineMap<FromFrame, ToFrame, Scalar, LinearMap> operator*(
+      AffineMap<ThroughFrame, ToFrame, Scalar, LinearMap> const& left,
+      AffineMap<FromFrame, ToFrame, Scalar, LinearMap> const& right);
+};
+
+template<typename FromFrame, typename ThroughFrame, typename ToFrame,
+         typename Scalar, template<typename, typename> class LinearMap>
+AffineMap<FromFrame, ToFrame, Scalar, LinearMap> operator*(
+    AffineMap<ThroughFrame, ToFrame, Scalar, LinearMap> const& left,
+    AffineMap<FromFrame, ToFrame, Scalar, LinearMap> const& right);
+
+}  // namespace geometry
+}  // namespace principia
+
+#include "geometry/affine_map_body.hpp"
diff --git a/geometry/affine_map_body.hpp b/geometry/affine_map_body.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+#include "geometry/point.hpp"
+#include "geometry/grassmann.hpp"
+
+namespace principia {
+namespace geometry {
+
+// Since the map is represented as x -> linear_map_(x) + translation_ and since
+// it is x -> linear_map(x - from_origin) + to_origin, we get
+// linear_map_ = linear_map, translation_ = to_origin - linear_map(from_origin).
+template<typename FromFrame, typename ToFrame, typename Scalar,
+         template<typename, typename> class LinearMap>
+AffineMap<FromFrame, ToFrame, Scalar, LinearMap>::AffineMap(
+    Point<FromVector> const& from_origin,
+    Point<ToVector> const& to_origin,
+    LinearMap<FromFrame, ToFrame> const& linear_map)
+    : linear_map_(linear_map),
+      translation_(to_origin.coordinates_ -
+                       linear_map(from_origin.coordinates_)) {}
+
+template<typename FromFrame, typename ToFrame, typename Scalar,
+         template<typename, typename> class LinearMap>
+AffineMap<ToFrame, FromFrame, Scalar, LinearMap>
+AffineMap<FromFrame, ToFrame, Scalar, LinearMap>::Inverse() const {
+  AffineMap<ToFrame, FromFrame, Scalar, LinearMap> result;
+  result.linear_map_ = linear_map_.Inverse();
+  result.translation_ = -result.linear_map_(translation_);
+  return result;
+}
+
+template<typename FromFrame, typename ToFrame, typename Scalar,
+         template<typename, typename> class LinearMap>
+Point<typename AffineMap<FromFrame, ToFrame, Scalar, LinearMap>::ToVector>
+AffineMap<FromFrame, ToFrame, Scalar, LinearMap>::operator()(
+    Point<FromVector> const& point) const {
+  return Point<
+      typename AffineMap<FromFrame, ToFrame, Scalar, LinearMap>::ToVector>(
+          linear_map_(point.coordinates_) + translation_);
+}
+
+template<typename FromFrame, typename ThroughFrame, typename ToFrame,
+         typename Scalar, template<typename, typename> class LinearMap>
+AffineMap<FromFrame, ToFrame, Scalar, LinearMap> operator*(
+    AffineMap<ThroughFrame, ToFrame, Scalar, LinearMap> const& left,
+    AffineMap<FromFrame, ToFrame, Scalar, LinearMap> const& right) {
+  AffineMap<ToFrame, FromFrame, Scalar, LinearMap> result;
+  result.linear_map_ = left.linear_map_ * right.linear_map_;
+  result.translation_ = left.linear_map_(right.translation_) +
+                        left.translation_;
+  return result;
+}
+
+}  // namespace geometry
+}  // namespace principia
diff --git a/geometry/affine_map_test.cpp b/geometry/affine_map_test.cpp
@@ -0,0 +1,109 @@
+#include "geometry/affine_map.hpp"
+
+#include <cfloat>
+
+#include <vector>
+
+#include "geometry/grassmann.hpp"
+#include "geometry/orthogonal_map.hpp"
+#include "geometry/permutation.hpp"
+#include "geometry/point.hpp"
+#include "geometry/rotation.hpp"
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+#include "quantities/numbers.hpp"
+#include "quantities/si.hpp"
+#include "testing_utilities/almost_equals.hpp"
+#include "testing_utilities/numerics.hpp"
+
+namespace principia {
+namespace geometry {
+
+using quantities::Length;
+using si::Metre;
+using si::Radian;
+using testing::Contains;
+using testing::Eq;
+using testing::Lt;
+using testing_utilities::AlmostEquals;
+using testing_utilities::RelativeError;
+
+class AffineMapTest : public testing::Test {
+ protected:
+  struct World;
+  typedef OrthogonalMap<World, World> Orth;
+  typedef Permutation<World, World> Perm;
+  typedef Rotation<World, World> Rot;
+  typedef Vector<Length, World> Displacement;
+  typedef Point<Displacement> Position;
+  typedef AffineMap<World, World, Length, Rotation> RigidTransformation;
+
+  void SetUp() override {
+    zero_     = Displacement({0 * Metre, 0 * Metre, 0 * Metre});
+    forward_  = Displacement({1 * Metre, 0 * Metre, 0 * Metre});
+    leftward_ = Displacement({0 * Metre, 1 * Metre, 0 * Metre});
+    upward_   = Displacement({0 * Metre, 0 * Metre, 1 * Metre});
+
+    origin_ = Position(Displacement({0 * Metre, 0 * Metre, 0 * Metre}));
+
+    back_right_bottom_  = Position(Displacement({3 * Metre,
+                                                 4 * Metre,
+                                                 5 * Metre}));
+    front_right_bottom_ = back_right_bottom_ + forward_;
+    back_left_bottom_   = back_right_bottom_ + leftward_;
+    front_left_bottom_  = back_left_bottom_ + forward_;
+    back_right_top_     = back_right_bottom_ + upward_;
+    front_right_top_    = back_right_top_ + forward_;
+    back_left_top_      = back_right_top_ + leftward_;
+    front_left_top_     = back_left_top_ + forward_;
+    vertices_ = {back_left_bottom_, front_left_bottom_, back_right_bottom_,
+                 front_right_bottom_, back_left_top_, front_left_top_,
+                 back_right_top_, front_right_top_};
+    originated_vertices_ = std::vector<Displacement>(vertices_.size());
+    for (std::size_t i = 0; i < vertices_.size(); ++i) {
+      originated_vertices_[i] = vertices_[i] - origin_;
+    }
+  }
+
+  Vector<Length, World> zero_;
+  Vector<Length, World> forward_;
+  Vector<Length, World> leftward_;
+  Vector<Length, World> upward_;
+  Position origin_;
+  Position back_left_bottom_;
+  Position front_left_bottom_;
+  Position back_right_bottom_;
+  Position front_right_bottom_;
+  Position back_left_top_;
+  Position front_left_top_;
+  Position back_right_top_;
+  Position front_right_top_;
+  std::vector<Position> vertices_;
+  std::vector<Displacement> originated_vertices_;
+};
+
+TEST_F(AffineMapTest, Cube) {
+  Rot rotate_left(π / 2 * Radian, upward_);
+  EXPECT_THAT(RelativeError(leftward_, rotate_left(forward_)),
+              Lt(2 * DBL_EPSILON));
+  RigidTransformation map = RigidTransformation(back_right_bottom_,
+                                                front_right_bottom_,
+                                                rotate_left);
+  EXPECT_THAT(map(back_right_bottom_) - origin_,
+              Eq(front_right_bottom_ - origin_));
+  EXPECT_THAT(map(front_left_top_) - origin_,
+              AlmostEquals(back_left_top_ - origin_));
+  // Check that |map| is an isometry of the cube whose vertices are |vertices_|.
+  for (auto const& point : vertices_) {
+    EXPECT_THAT(originated_vertices_,
+                Contains(AlmostEquals(map(point) - origin_)));
+  }
+  // Test that |map.Inverse() * map| acts as the identity on that cube.
+  for (std::size_t i = 0; i < vertices_.size(); ++i) {
+    EXPECT_THAT(originated_vertices_[i],
+                AlmostEquals((map.Inverse() * map)(vertices_[i]) - origin_));
+  }
+}
+
+}  // namespace geometry
+}  // namespace principia
diff --git a/geometry/geometry.vcxproj b/geometry/geometry.vcxproj
@@ -83,6 +83,7 @@
   </ItemDefinitionGroup>
   <ItemGroup>
     <ClInclude Include="affine_map.hpp" />
+    <ClInclude Include="affine_map_body.hpp" />
     <ClInclude Include="point.hpp" />
     <ClInclude Include="point_body.hpp" />
     <ClInclude Include="grassmann_body.hpp" />
@@ -104,6 +105,7 @@
   <ItemGroup>
     <ClCompile Include="grassmann_test.cpp" />
     <ClCompile Include="point_test.cpp" />
+    <ClCompile Include="affine_map_test.cpp" />
     <ClCompile Include="orthogonal_map_test.cpp" />
     <ClCompile Include="permutation_test.cpp" />
     <ClCompile Include="quaternion_test.cpp" />

diff --git a/geometry/geometry.vcxproj.filters b/geometry/geometry.vcxproj.filters
@@ -68,6 +68,9 @@
     <ClInclude Include="point_body.hpp">
       <Filter>Source Files</Filter>
     </ClInclude>
+    <ClInclude Include="affine_map_body.hpp">
+      <Filter>Source Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="sign_test.cpp">
@@ -94,5 +97,8 @@
     <ClCompile Include="r3_element_test.cpp">
       <Filter>Test Files</Filter>
     </ClCompile>
+    <ClCompile Include="affine_map_test.cpp">
+      <Filter>Test Files</Filter>
+    </ClCompile>
   </ItemGroup>
 </Project>
diff --git a/geometry/grassmann.hpp b/geometry/grassmann.hpp
@@ -25,6 +25,7 @@ class Multivector<Scalar, Frame, 1> {
   ~Multivector() = default;
 
   R3Element<Scalar> coordinates() const;
+  Scalar Norm() const;
 
  private:
   R3Element<Scalar> coordinates_;
@@ -37,6 +38,7 @@ struct Multivector<Scalar, Frame, 2> {
   ~Multivector() = default;
 
   R3Element<Scalar> coordinates() const;
+  Scalar Norm() const;
 
  private:
   R3Element<Scalar> coordinates_;
@@ -49,6 +51,7 @@ struct Multivector<Scalar, Frame, 3> {
   ~Multivector() = default;
 
   Scalar coordinates() const;
+  Scalar Norm() const;
 
  private:
   Scalar coordinates_;

diff --git a/geometry/grassmann_body.hpp b/geometry/grassmann_body.hpp
@@ -30,6 +30,21 @@ inline Scalar Multivector<Scalar, Frame, 3>::coordinates() const {
   return coordinates_;
 }
 
+template<typename Scalar, typename Frame>
+inline Scalar Multivector<Scalar, Frame, 1>::Norm() const {
+  return coordinates_.Norm();
+}
+
+template<typename Scalar, typename Frame>
+inline Scalar Multivector<Scalar, Frame, 2>::Norm() const {
+  return coordinates_.Norm();
+}
+
+template<typename Scalar, typename Frame>
+inline Scalar Multivector<Scalar, Frame, 3>::Norm() const {
+  return quantities::Abs(coordinates_);
+}
+
 template<typename LScalar, typename RScalar, typename Frame>
 inline quantities::Product<LScalar, RScalar> InnerProduct(
     Vector<LScalar, Frame> const& left,

diff --git a/geometry/grassmann_test.cpp b/geometry/grassmann_test.cpp
@@ -223,5 +223,16 @@ TEST_F(GrassmannTest, Actions) {
 
 #pragma warning(default: 4566)
 
+TEST_F(GrassmannTest, Norm) {
+  Vector<Length, World> const v({-3 * 4 * Metre, 4 * 4 * Metre, 5 * 3 * Metre});
+  EXPECT_THAT(v.Norm(), Eq(5 * 5 * Metre));
+  Bivector<Length, World> const w({+20 * 21 * Metre,
+                                   -21 * 21 * Metre,
+                                   +29 * 20 * Metre});
+  EXPECT_THAT(w.Norm(), Eq(29 * 29 * Metre));
+  Trivector<Length, World> const u(-4 * Furlong);
+  EXPECT_THAT(u.Norm(), Eq(4 * Furlong));
+}
+
 }  // namespace geometry
 }  // namespace principia
diff --git a/geometry/point.hpp b/geometry/point.hpp
@@ -24,6 +24,13 @@ class Point {
   bool operator==(Point const& right) const;
   bool operator!=(Point const& right) const;
 
+ private:
+  Vector coordinates_;
+
+  template<typename FromFrame, typename ToFrame, typename Scalar,
+           template<typename, typename> class LinearMap>
+  friend class AffineMap;
+
   template<typename Vector>
   friend Point<Vector> operator+(Vector const& translation,
                                  Point<Vector> const& point);
@@ -33,9 +40,6 @@ class Point {
                                   Weight const& left_weight,
                                   Point<Vector> const& right,
                                   Weight const& right_weight);
-
- private:
-  Vector coordinates_;
 };
 
 template<typename Vector>

diff --git a/testing_utilities/numerics.hpp b/testing_utilities/numerics.hpp
@@ -33,6 +33,11 @@ template<typename Scalar>
 Scalar AbsoluteError(geometry::R3Element<Scalar> const& expected,
                      geometry::R3Element<Scalar> const& actual);
 
+// Uses Multivector.Norm().
+template<typename Scalar, typename Frame, unsigned int Rank>
+Scalar AbsoluteError(
+    geometry::Multivector<Scalar, Frame, Rank> const& expected,
+    geometry::Multivector<Scalar, Frame, Rank> const& actual);
 
 template<typename T, typename Norm>
 quantities::Dimensionless RelativeError(T const& expected, T const& actual,
@@ -55,6 +60,12 @@ quantities::Dimensionless RelativeError(
     geometry::R3Element<Scalar> const& expected,
     geometry::R3Element<Scalar> const& actual);
 
+// Uses Multivector.Norm().
+template<typename Scalar, typename Frame, unsigned int Rank>
+quantities::Dimensionless RelativeError(
+    geometry::Multivector<Scalar, Frame, Rank> const& expected,
+    geometry::Multivector<Scalar, Frame, Rank> const& actual);
+
 std::int64_t ULPDistance(double const x, double const y);
 
 }  // namespace testing_utilities

diff --git a/testing_utilities/numerics_body.hpp b/testing_utilities/numerics_body.hpp
@@ -43,6 +43,20 @@ Scalar AbsoluteError(geometry::R3Element<Scalar> const& expected,
               [](geometry::R3Element<Scalar> const& v) { return v.Norm(); });
 }
 
+template<typename Scalar, typename Frame, unsigned int Rank>
+Scalar AbsoluteError(
+    geometry::Multivector<Scalar, Frame, Rank> const& expected,
+    geometry::Multivector<Scalar, Frame, Rank> const& actual) {
+  return AbsoluteError<
+      geometry::Multivector<Scalar, Frame, Rank>,
+      std::function<Scalar(geometry::Multivector<Scalar, Frame, Rank>)>,
+      Scalar>(expected,
+              actual,
+              [](geometry::Multivector<Scalar, Frame, Rank> v) {
+                return v.Norm();
+              });
+}
+
 template<typename T, typename Norm>
 quantities::Dimensionless RelativeError(T const& expected, T const& actual,
                                         Norm const norm) {
@@ -70,6 +84,16 @@ quantities::Dimensionless RelativeError(
                        [](geometry::R3Element<Scalar> v) { return v.Norm(); });
 }
 
+template<typename Scalar, typename Frame, unsigned int Rank>
+quantities::Dimensionless RelativeError(
+    geometry::Multivector<Scalar, Frame, Rank> const& expected,
+    geometry::Multivector<Scalar, Frame, Rank> const& actual) {
+  return RelativeError(
+      expected,
+      actual,
+      [](geometry::Multivector<Scalar, Frame, Rank> v) { return v.Norm(); });
+}
+
 union Qword {
   double double_value;
   std::int64_t long_value;