1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#include "main.h"
12#include <Eigen/Geometry>
13#include <Eigen/LU>
14#include <Eigen/QR>
15
16template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane)
17{
18  /* this test covers the following files:
19     Hyperplane.h
20  */
21  typedef typename HyperplaneType::Index Index;
22  const Index dim = _plane.dim();
23  enum { Options = HyperplaneType::Options };
24  typedef typename HyperplaneType::Scalar Scalar;
25  typedef typename NumTraits<Scalar>::Real RealScalar;
26  typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime, 1> VectorType;
27  typedef Matrix<Scalar, HyperplaneType::AmbientDimAtCompileTime,
28                         HyperplaneType::AmbientDimAtCompileTime> MatrixType;
29
30  VectorType p0 = VectorType::Random(dim);
31  VectorType p1 = VectorType::Random(dim);
32
33  VectorType n0 = VectorType::Random(dim).normalized();
34  VectorType n1 = VectorType::Random(dim).normalized();
35
36  HyperplaneType pl0(n0, p0);
37  HyperplaneType pl1(n1, p1);
38  HyperplaneType pl2 = pl1;
39
40  Scalar s0 = internal::random<Scalar>();
41  Scalar s1 = internal::random<Scalar>();
42
43  VERIFY_IS_APPROX( n1.dot(n1), Scalar(1) );
44
45  VERIFY_IS_MUCH_SMALLER_THAN( pl0.absDistance(p0), Scalar(1) );
46  VERIFY_IS_APPROX( pl1.signedDistance(p1 + n1 * s0), s0 );
47  VERIFY_IS_MUCH_SMALLER_THAN( pl1.signedDistance(pl1.projection(p0)), Scalar(1) );
48  VERIFY_IS_MUCH_SMALLER_THAN( pl1.absDistance(p1 +  pl1.normal().unitOrthogonal() * s1), Scalar(1) );
49
50  // transform
51  if (!NumTraits<Scalar>::IsComplex)
52  {
53    MatrixType rot = MatrixType::Random(dim,dim).householderQr().householderQ();
54    DiagonalMatrix<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random());
55    Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random());
56
57    pl2 = pl1;
58    VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot).absDistance(rot * p1), Scalar(1) );
59    pl2 = pl1;
60    VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot,Isometry).absDistance(rot * p1), Scalar(1) );
61    pl2 = pl1;
62    VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling).absDistance((rot*scaling) * p1), Scalar(1) );
63    pl2 = pl1;
64    VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*scaling*translation)
65                                 .absDistance((rot*scaling*translation) * p1), Scalar(1) );
66    pl2 = pl1;
67    VERIFY_IS_MUCH_SMALLER_THAN( pl2.transform(rot*translation,Isometry)
68                                 .absDistance((rot*translation) * p1), Scalar(1) );
69  }
70
71  // casting
72  const int Dim = HyperplaneType::AmbientDimAtCompileTime;
73  typedef typename GetDifferentType<Scalar>::type OtherScalar;
74  Hyperplane<OtherScalar,Dim,Options> hp1f = pl1.template cast<OtherScalar>();
75  VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),pl1);
76  Hyperplane<Scalar,Dim,Options> hp1d = pl1.template cast<Scalar>();
77  VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),pl1);
78}
79
80template<typename Scalar> void lines()
81{
82  typedef Hyperplane<Scalar, 2> HLine;
83  typedef ParametrizedLine<Scalar, 2> PLine;
84  typedef Matrix<Scalar,2,1> Vector;
85  typedef Matrix<Scalar,3,1> CoeffsType;
86
87  for(int i = 0; i < 10; i++)
88  {
89    Vector center = Vector::Random();
90    Vector u = Vector::Random();
91    Vector v = Vector::Random();
92    Scalar a = internal::random<Scalar>();
93    while (internal::abs(a-1) < 1e-4) a = internal::random<Scalar>();
94    while (u.norm() < 1e-4) u = Vector::Random();
95    while (v.norm() < 1e-4) v = Vector::Random();
96
97    HLine line_u = HLine::Through(center + u, center + a*u);
98    HLine line_v = HLine::Through(center + v, center + a*v);
99
100    // the line equations should be normalized so that a^2+b^2=1
101    VERIFY_IS_APPROX(line_u.normal().norm(), Scalar(1));
102    VERIFY_IS_APPROX(line_v.normal().norm(), Scalar(1));
103
104    Vector result = line_u.intersection(line_v);
105
106    // the lines should intersect at the point we called "center"
107    VERIFY_IS_APPROX(result, center);
108
109    // check conversions between two types of lines
110    PLine pl(line_u); // gcc 3.3 will commit suicide if we don't name this variable
111    CoeffsType converted_coeffs = HLine(pl).coeffs();
112    converted_coeffs *= (line_u.coeffs()[0])/(converted_coeffs[0]);
113    VERIFY(line_u.coeffs().isApprox(converted_coeffs));
114  }
115}
116
117template<typename Scalar> void hyperplane_alignment()
118{
119  typedef Hyperplane<Scalar,3,AutoAlign> Plane3a;
120  typedef Hyperplane<Scalar,3,DontAlign> Plane3u;
121
122  EIGEN_ALIGN16 Scalar array1[4];
123  EIGEN_ALIGN16 Scalar array2[4];
124  EIGEN_ALIGN16 Scalar array3[4+1];
125  Scalar* array3u = array3+1;
126
127  Plane3a *p1 = ::new(reinterpret_cast<void*>(array1)) Plane3a;
128  Plane3u *p2 = ::new(reinterpret_cast<void*>(array2)) Plane3u;
129  Plane3u *p3 = ::new(reinterpret_cast<void*>(array3u)) Plane3u;
130
131  p1->coeffs().setRandom();
132  *p2 = *p1;
133  *p3 = *p1;
134
135  VERIFY_IS_APPROX(p1->coeffs(), p2->coeffs());
136  VERIFY_IS_APPROX(p1->coeffs(), p3->coeffs());
137
138  #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
139  if(internal::packet_traits<Scalar>::Vectorizable)
140    VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Plane3a));
141  #endif
142}
143
144
145void test_geo_hyperplane()
146{
147  for(int i = 0; i < g_repeat; i++) {
148    CALL_SUBTEST_1( hyperplane(Hyperplane<float,2>()) );
149    CALL_SUBTEST_2( hyperplane(Hyperplane<float,3>()) );
150    CALL_SUBTEST_2( hyperplane(Hyperplane<float,3,DontAlign>()) );
151    CALL_SUBTEST_2( hyperplane_alignment<float>() );
152    CALL_SUBTEST_3( hyperplane(Hyperplane<double,4>()) );
153    CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) );
154    CALL_SUBTEST_1( lines<float>() );
155    CALL_SUBTEST_3( lines<double>() );
156  }
157}
158