Fix Hyperplane::Through(a,b,c) when points are aligned or identical. We use the stratgey as in Quaternion::setFromTwoVectors.

This commit is contained in:
Gael Guennebaud 2014-11-26 15:01:53 +01:00
parent da584912b6
commit 8518ba0bbc
2 changed files with 40 additions and 1 deletions

View File

@ -100,7 +100,17 @@ public:
{
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 3)
Hyperplane result(p0.size());
result.normal() = (p2 - p0).cross(p1 - p0).normalized();
VectorType v0(p2 - p0), v1(p1 - p0);
result.normal() = v0.cross(v1);
RealScalar norm = result.normal().norm();
if(norm <= v0.norm() * v1.norm() * NumTraits<RealScalar>::epsilon())
{
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
result.normal() = svd.matrixV().col(2);
}
else
result.normal() /= norm;
result.offset() = -p0.dot(result.normal());
return result;
}

View File

@ -124,6 +124,33 @@ template<typename Scalar> void lines()
}
}
template<typename Scalar> void planes()
{
using std::abs;
typedef Hyperplane<Scalar, 3> Plane;
typedef Matrix<Scalar,3,1> Vector;
typedef Matrix<Scalar,4,1> CoeffsType;
for(int i = 0; i < 10; i++)
{
Vector v0 = Vector::Random();
Vector v1(v0), v2(v0);
if(internal::random<double>(0,1)>0.25)
v1 += Vector::Random();
if(internal::random<double>(0,1)>0.25)
v2 += v1 * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16));
if(internal::random<double>(0,1)>0.25)
v2 += Vector::Random() * std::pow(internal::random<Scalar>(0,1),internal::random<int>(1,16));
Plane p0 = Plane::Through(v0, v1, v2);
VERIFY_IS_APPROX(p0.normal().norm(), Scalar(1));
VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v0), Scalar(1));
VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v1), Scalar(1));
VERIFY_IS_MUCH_SMALLER_THAN(p0.absDistance(v2), Scalar(1));
}
}
template<typename Scalar> void hyperplane_alignment()
{
typedef Hyperplane<Scalar,3,AutoAlign> Plane3a;
@ -163,5 +190,7 @@ void test_geo_hyperplane()
CALL_SUBTEST_4( hyperplane(Hyperplane<std::complex<double>,5>()) );
CALL_SUBTEST_1( lines<float>() );
CALL_SUBTEST_3( lines<double>() );
CALL_SUBTEST_2( planes<float>() );
CALL_SUBTEST_5( planes<double>() );
}
}