From 6425e875a1158e1e2a0afcf703105e8ddbfee7bd Mon Sep 17 00:00:00 2001 From: Martin Pecka Date: Mon, 28 Sep 2020 18:06:23 +0000 Subject: [PATCH] Added AlignedBox::transform(AffineTransform). --- Eigen/src/Geometry/AlignedBox.h | 94 ++++++++ test/geo_alignedbox.cpp | 400 +++++++++++++++++++++++++++++++- 2 files changed, 483 insertions(+), 11 deletions(-) diff --git a/Eigen/src/Geometry/AlignedBox.h b/Eigen/src/Geometry/AlignedBox.h index c902d8f0a..f8fb423cb 100644 --- a/Eigen/src/Geometry/AlignedBox.h +++ b/Eigen/src/Geometry/AlignedBox.h @@ -7,6 +7,42 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +// Function void Eigen::AlignedBox::transform(const Transform& transform) +// is provided under the following license agreement: +// +// Software License Agreement (BSD License) +// +// Copyright (c) 2011-2014, Willow Garage, Inc. +// Copyright (c) 2014-2015, Open Source Robotics Foundation +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above +// copyright notice, this list of conditions and the following +// disclaimer in the documentation and/or other materials provided +// with the distribution. +// * Neither the name of Open Source Robotics Foundation nor the names of its +// contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +// COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + #ifndef EIGEN_ALIGNEDBOX_H #define EIGEN_ALIGNEDBOX_H @@ -246,6 +282,15 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) return *this; } + /** \returns a copy of \c *this translated by the vector \a t. */ + template + EIGEN_DEVICE_FUNC inline AlignedBox translated(const MatrixBase& a_t) const + { + AlignedBox result(m_min, m_max); + result.translate(a_t); + return result; + } + /** \returns the squared distance between the point \a p and the box \c *this, * and zero if \a p is inside the box. * \sa exteriorDistance(const MatrixBase&), squaredExteriorDistance(const AlignedBox&) @@ -274,6 +319,55 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) EIGEN_DEVICE_FUNC inline NonInteger exteriorDistance(const AlignedBox& b) const { EIGEN_USING_STD_MATH(sqrt) return sqrt(NonInteger(squaredExteriorDistance(b))); } + /** + * Specialization of transform for pure translation. + */ + template + EIGEN_DEVICE_FUNC inline void transform( + const typename Transform::TranslationType& translation) + { + this->translate(translation); + } + + /** + * Transforms this box by \a transform and recomputes it to + * still be an axis-aligned box. + * + * \note This method is provided under BSD license (see the top of this file). + */ + template + EIGEN_DEVICE_FUNC inline void transform(const Transform& transform) + { + // Projective transform is not (yet) supported + EIGEN_STATIC_ASSERT(Mode != Projective, THIS_METHOD_IS_ONLY_FOR_SPECIFIC_TRANSFORMATIONS); + + // Method adapted from FCL src/shape/geometric_shapes_utility.cpp#computeBV(...) + // https://github.com/flexible-collision-library/fcl/blob/fcl-0.4/src/shape/geometric_shapes_utility.cpp#L292 + // + // Here's a nice explanation why it works: https://zeuxcg.org/2010/10/17/aabb-from-obb-with-component-wise-abs/ + + // two times rotated extent + const VectorType rotated_extent_2 = transform.linear().cwiseAbs() * sizes(); + // two times new center + const VectorType rotated_center_2 = transform.linear() * (this->m_max + this->m_min) + + Scalar(2) * transform.translation(); + + this->m_max = (rotated_center_2 + rotated_extent_2) / Scalar(2); + this->m_min = (rotated_center_2 - rotated_extent_2) / Scalar(2); + } + + /** + * \returns a copy of \c *this transformed by \a transform and recomputed to + * still be an axis-aligned box. + */ + template + EIGEN_DEVICE_FUNC AlignedBox transformed(const Transform& transform) const + { + AlignedBox result(m_min, m_max); + result.transform(transform); + return result; + } + /** \returns \c *this with scalar type casted to \a NewScalarType * * Note that if \a NewScalarType is equal to the current scalar type of \c *this diff --git a/test/geo_alignedbox.cpp b/test/geo_alignedbox.cpp index c6c051ce4..0a866a961 100644 --- a/test/geo_alignedbox.cpp +++ b/test/geo_alignedbox.cpp @@ -9,14 +9,12 @@ #include "main.h" #include -#include -#include #include using namespace std; // NOTE the following workaround was needed on some 32 bits builds to kill extra precision of x87 registers. -// It seems that it os not needed anymore, but let's keep it here, just in case... +// It seems that it is not needed anymore, but let's keep it here, just in case... template EIGEN_DONT_INLINE void kill_extra_precision(T& /* x */) { @@ -34,7 +32,8 @@ template void alignedbox(const BoxType& _box) AlignedBox.h */ typedef typename BoxType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; + typedef NumTraits ScalarTraits; + typedef typename ScalarTraits::Real RealScalar; typedef Matrix VectorType; const Index dim = _box.dim(); @@ -88,7 +87,386 @@ template void alignedbox(const BoxType& _box) } +template void alignedboxTranslatable(const BoxType& _box) +{ + typedef typename BoxType::Scalar Scalar; + typedef Matrix VectorType; + typedef Transform IsometryTransform; + typedef Transform AffineTransform; + alignedbox(_box); + + const VectorType Ones = VectorType::Ones(); + const VectorType UnitX = VectorType::UnitX(); + const Index dim = _box.dim(); + + // box((-1, -1, -1), (1, 1, 1)) + BoxType a(-Ones, Ones); + + VERIFY_IS_APPROX(a.sizes(), Ones * Scalar(2)); + + BoxType b = a; + VectorType translate = Ones; + translate[0] = Scalar(2); + b.translate(translate); + // translate by (2, 1, 1) -> box((1, 0, 0), (3, 2, 2)) + + VERIFY_IS_APPROX(b.sizes(), Ones * Scalar(2)); + VERIFY_IS_APPROX((b.min)(), UnitX); + VERIFY_IS_APPROX((b.max)(), Ones * Scalar(2) + UnitX); + + // Test transform + + IsometryTransform tf = IsometryTransform::Identity(); + tf.translation() = -translate; + + BoxType c = b.transformed(tf); + // translate by (-2, -1, -1) -> box((-1, -1, -1), (1, 1, 1)) + VERIFY_IS_APPROX(c.sizes(), a.sizes()); + VERIFY_IS_APPROX((c.min)(), (a.min)()); + VERIFY_IS_APPROX((c.max)(), (a.max)()); + + c.transform(tf); + // translate by (-2, -1, -1) -> box((-3, -2, -2), (-1, 0, 0)) + VERIFY_IS_APPROX(c.sizes(), a.sizes()); + VERIFY_IS_APPROX((c.min)(), Ones * Scalar(-2) - UnitX); + VERIFY_IS_APPROX((c.max)(), -UnitX); + + // Scaling + + AffineTransform atf = AffineTransform::Identity(); + atf.scale(Scalar(3)); + c.transform(atf); + // scale by 3 -> box((-9, -6, -6), (-3, 0, 0)) + VERIFY_IS_APPROX(c.sizes(), Scalar(3) * a.sizes()); + VERIFY_IS_APPROX((c.min)(), Ones * Scalar(-6) - UnitX * Scalar(3)); + VERIFY_IS_APPROX((c.max)(), UnitX * Scalar(-3)); + + atf = AffineTransform::Identity(); + atf.scale(Scalar(-3)); + c.transform(atf); + // scale by -3 -> box((27, 18, 18), (9, 0, 0)) + VERIFY_IS_APPROX(c.sizes(), Scalar(9) * a.sizes()); + VERIFY_IS_APPROX((c.min)(), UnitX * Scalar(9)); + VERIFY_IS_APPROX((c.max)(), Ones * Scalar(18) + UnitX * Scalar(9)); + + // test for roundoff errors + IsometryTransform identity = IsometryTransform::Identity(); + BoxType transformedC; + transformedC.extend(c.transformed(identity)); + VERIFY(transformedC.contains(c)); + + for (size_t i = 0; i < 10; ++i) + { + VectorType minCorner; + VectorType maxCorner; + for (Index d = 0; d < dim; ++d) + { + minCorner[d] = internal::random(-10,10); + maxCorner[d] = minCorner[d] + internal::random(0, 10); + } + + c = BoxType(minCorner, maxCorner); + + translate = VectorType::Random(); + c.translate(translate); + + VERIFY_IS_APPROX((c.min)(), minCorner + translate); + VERIFY_IS_APPROX((c.max)(), maxCorner + translate); + } +} + +template +Rotation rotate2D(Scalar _angle) { + return Rotation2D(_angle); +} + +template +Rotation rotate2DIntegral(typename NumTraits::NonInteger _angle) { + typedef typename NumTraits::NonInteger NonInteger; + return Rotation2D(_angle).toRotationMatrix(). + template cast(); +} + +template +Rotation rotate3DZAxis(Scalar _angle) { + return AngleAxis(_angle, Matrix(0, 0, 1)); +} + +template +Rotation rotate3DZAxisIntegral(typename NumTraits::NonInteger _angle) { + typedef typename NumTraits::NonInteger NonInteger; + return AngleAxis(_angle, Matrix(0, 0, 1)). + toRotationMatrix().template cast(); +} + +template +Rotation rotate4DZWAxis(Scalar _angle) { + Rotation result = Matrix::Identity(); + result.block(0, 0, 3, 3) = rotate3DZAxis(_angle).toRotationMatrix(); + return result; +} + +template +MatrixType randomRotationMatrix() +{ + // algorithm from + // https://www.isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/III-7/103/2016/isprs-annals-III-7-103-2016.pdf + const MatrixType rand = MatrixType::Random(); + const MatrixType q = rand.householderQr().householderQ(); + const JacobiSVD svd = q.jacobiSvd(ComputeFullU | ComputeFullV); + const typename MatrixType::Scalar det = (svd.matrixU() * svd.matrixV().transpose()).determinant(); + MatrixType diag = rand.Identity(); + diag(MatrixType::RowsAtCompileTime - 1, MatrixType::ColsAtCompileTime - 1) = det; + const MatrixType rotation = svd.matrixU() * diag * svd.matrixV().transpose(); + return rotation; +} + +template +std::vector > boxGetCorners(const Matrix& _min, const Matrix& _max, int dim = Dim) +{ + std::vector > result; + if (dim == 1) + { + result.push_back(_min); + result.push_back(_max); + } + else + { + std::vector > shorter = boxGetCorners(_min, _max, dim - 1); + for (size_t i = 0; i < shorter.size(); ++i) + { + Matrix vec = shorter[i]; + + Matrix vec1 = _min; + vec1.block(Dim - dim, 0, dim - 1, 1) = vec.block(Dim - dim, 0, dim - 1, 1); + result.push_back(vec1); + + Matrix vec2 = _max; + vec2.block(Dim - dim, 0, dim - 1, 1) = vec.block(Dim - dim, 0, dim - 1, 1); + result.push_back(vec2); + } + } + return result; +} + +template void alignedboxRotatable( + const BoxType& _box, + Rotation (*_rotate)(typename NumTraits::NonInteger /*_angle*/)) +{ + alignedboxTranslatable(_box); + + typedef typename BoxType::Scalar Scalar; + typedef typename NumTraits::NonInteger NonInteger; + typedef Matrix VectorType; + typedef Transform IsometryTransform; + typedef Transform AffineTransform; + + const VectorType Zero = VectorType::Zero(); + const VectorType Ones = VectorType::Ones(); + const VectorType UnitX = VectorType::UnitX(); + const VectorType UnitY = VectorType::UnitY(); + // this is vector (0, 0, -1, -1, -1, ...), i.e. with zeros at first and second dimensions + const VectorType UnitZ = Ones - UnitX - UnitY; + + // in this kind of comments the 3D case values will be illustrated + // box((-1, -1, -1), (1, 1, 1)) + BoxType a(-Ones, Ones); + + // to allow templating this test for both 2D and 3D cases, we always set all + // but the first coordinate to the same value; so basically 3D case works as + // if you were looking at the scene from top + + VectorType minPoint = -2 * Ones; + minPoint[0] = -3; + VectorType maxPoint = Zero; + maxPoint[0] = -1; + BoxType c(minPoint, maxPoint); + // box((-3, -2, -2), (-1, 0, 0)) + + IsometryTransform tf2 = IsometryTransform::Identity(); + // for some weird reason the following statement has to be put separate from + // the following rotate call, otherwise precision problems arise... + Rotation rot = _rotate(NonInteger(EIGEN_PI)); + tf2.rotate(rot); + + c.transform(tf2); + // rotate by 180 deg around origin -> box((1, 0, -2), (3, 2, 0)) + + VERIFY_IS_APPROX(c.sizes(), a.sizes()); + VERIFY_IS_APPROX((c.min)(), UnitX - UnitZ * Scalar(2)); + VERIFY_IS_APPROX((c.max)(), UnitX * Scalar(3) + UnitY * Scalar(2)); + + rot = _rotate(NonInteger(EIGEN_PI / 2)); + tf2.setIdentity(); + tf2.rotate(rot); + + c.transform(tf2); + // rotate by 90 deg around origin -> box((-2, 1, -2), (0, 3, 0)) + + VERIFY_IS_APPROX(c.sizes(), a.sizes()); + VERIFY_IS_APPROX((c.min)(), Ones * Scalar(-2) + UnitY * Scalar(3)); + VERIFY_IS_APPROX((c.max)(), UnitY * Scalar(3)); + + // box((-1, -1, -1), (1, 1, 1)) + AffineTransform atf = AffineTransform::Identity(); + atf.linearExt()(0, 1) = Scalar(1); + c = BoxType(-Ones, Ones); + c.transform(atf); + // 45 deg shear in x direction -> box((-2, -1, -1), (2, 1, 1)) + + VERIFY_IS_APPROX(c.sizes(), Ones * Scalar(2) + UnitX * Scalar(2)); + VERIFY_IS_APPROX((c.min)(), -Ones - UnitX); + VERIFY_IS_APPROX((c.max)(), Ones + UnitX); +} + +template void alignedboxNonIntegralRotatable( + const BoxType& _box, + Rotation (*_rotate)(typename NumTraits::NonInteger /*_angle*/)) +{ + alignedboxRotatable(_box, _rotate); + + typedef typename BoxType::Scalar Scalar; + typedef typename NumTraits::NonInteger NonInteger; + typedef Matrix VectorType; + typedef Transform IsometryTransform; + typedef Transform AffineTransform; + + const Index dim = _box.dim(); + const VectorType Zero = VectorType::Zero(); + const VectorType Ones = VectorType::Ones(); + const VectorType UnitX = VectorType::UnitX(); + const VectorType UnitY = VectorType::UnitY(); + // this is vector (0, 0, -1, -1, -1, ...), i.e. with zeros at first and second dimensions + const VectorType UnitZ = Ones - UnitX - UnitY; + + VectorType minPoint = -2 * Ones; + minPoint[1] = 1; + VectorType maxPoint = Zero; + maxPoint[1] = 3; + BoxType c(minPoint, maxPoint); + // ((-2, 1, -2), (0, 3, 0)) + + VectorType cornerBL = (c.min)(); + VectorType cornerTR = (c.max)(); + VectorType cornerBR = (c.min)(); cornerBR[0] = cornerTR[0]; + VectorType cornerTL = (c.max)(); cornerTL[0] = cornerBL[0]; + + NonInteger angle = NonInteger(EIGEN_PI/3); + Rotation rot = _rotate(angle); + IsometryTransform tf2; + tf2.setIdentity(); + tf2.rotate(rot); + + c.transform(tf2); + // rotate by 60 deg -> box((-3.59, -1.23, -2), (-0.86, 1.5, 0)) + + cornerBL = tf2 * cornerBL; + cornerBR = tf2 * cornerBR; + cornerTL = tf2 * cornerTL; + cornerTR = tf2 * cornerTR; + + VectorType minCorner = Ones * Scalar(-2); + VectorType maxCorner = Zero; + minCorner[0] = (min)((min)(cornerBL[0], cornerBR[0]), (min)(cornerTL[0], cornerTR[0])); + maxCorner[0] = (max)((max)(cornerBL[0], cornerBR[0]), (max)(cornerTL[0], cornerTR[0])); + minCorner[1] = (min)((min)(cornerBL[1], cornerBR[1]), (min)(cornerTL[1], cornerTR[1])); + maxCorner[1] = (max)((max)(cornerBL[1], cornerBR[1]), (max)(cornerTL[1], cornerTR[1])); + + for (Index d = 2; d < dim; ++d) + VERIFY_IS_APPROX(c.sizes()[d], Scalar(2)); + + VERIFY_IS_APPROX((c.min)(), minCorner); + VERIFY_IS_APPROX((c.max)(), maxCorner); + + VectorType minCornerValue = Ones * Scalar(-2); + VectorType maxCornerValue = Zero; + minCornerValue[0] = Scalar(Scalar(-sqrt(2*2 + 3*3)) * Scalar(cos(Scalar(atan(2.0/3.0)) - angle/2))); + minCornerValue[1] = Scalar(Scalar(-sqrt(1*1 + 2*2)) * Scalar(sin(Scalar(atan(2.0/1.0)) - angle/2))); + maxCornerValue[0] = Scalar(-sin(angle)); + maxCornerValue[1] = Scalar(3 * cos(angle)); + VERIFY_IS_APPROX((c.min)(), minCornerValue); + VERIFY_IS_APPROX((c.max)(), maxCornerValue); + + // randomized test - translate and rotate the box and compare to a box made of transformed vertices + for (size_t i = 0; i < 10; ++i) + { + for (Index d = 0; d < dim; ++d) + { + minCorner[d] = internal::random(-10,10); + maxCorner[d] = minCorner[d] + internal::random(0, 10); + } + + c = BoxType(minCorner, maxCorner); + + std::vector corners = boxGetCorners(minCorner, maxCorner); + const size_t numCorners = corners.size(); + + typename AffineTransform::LinearMatrixType rotation = + randomRotationMatrix(); + + tf2.setIdentity(); + tf2.rotate(rotation); + tf2.translate(VectorType::Random()); + + c.transform(tf2); + for (size_t corner = 0; corner < numCorners; ++corner) + corners[corner] = tf2 * corners[corner]; + + for (Index d = 0; d < dim; ++d) + { + minCorner[d] = corners[0][d]; + maxCorner[d] = corners[0][d]; + + for (size_t corner = 0; corner < numCorners; ++corner) + { + minCorner[d] = (min)(minCorner[d], corners[corner][d]); + maxCorner[d] = (max)(maxCorner[d], corners[corner][d]); + } + } + + VERIFY_IS_APPROX((c.min)(), minCorner); + VERIFY_IS_APPROX((c.max)(), maxCorner); + } + + // randomized test - transform the box with a random affine matrix and compare to a box made of transformed vertices + for (size_t i = 0; i < 10; ++i) + { + for (Index d = 0; d < dim; ++d) + { + minCorner[d] = internal::random(-10,10); + maxCorner[d] = minCorner[d] + internal::random(0, 10); + } + + c = BoxType(minCorner, maxCorner); + + std::vector corners = boxGetCorners(minCorner, maxCorner); + const size_t numCorners = corners.size(); + + AffineTransform atf = AffineTransform::Identity(); + atf.linearExt() = AffineTransform::LinearPart::Random(); + atf.translate(VectorType::Random()); + + c.transform(atf); + for (size_t corner = 0; corner < numCorners; ++corner) + corners[corner] = atf * corners[corner]; + + for (Index d = 0; d < dim; ++d) + { + minCorner[d] = corners[0][d]; + maxCorner[d] = corners[0][d]; + + for (size_t corner = 0; corner < numCorners; ++corner) + { + minCorner[d] = (min)(minCorner[d], corners[corner][d]); + maxCorner[d] = (max)(maxCorner[d], corners[corner][d]); + } + } + + VERIFY_IS_APPROX((c.min)(), minCorner); + VERIFY_IS_APPROX((c.max)(), maxCorner); + } +} template void alignedboxCastTests(const BoxType& _box) @@ -173,21 +551,21 @@ EIGEN_DECLARE_TEST(geo_alignedbox) { for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_1( alignedbox(AlignedBox2f()) ); + CALL_SUBTEST_1( (alignedboxNonIntegralRotatable(AlignedBox2f(), &rotate2D)) ); CALL_SUBTEST_2( alignedboxCastTests(AlignedBox2f()) ); - CALL_SUBTEST_3( alignedbox(AlignedBox3f()) ); + CALL_SUBTEST_3( (alignedboxNonIntegralRotatable(AlignedBox3f(), &rotate3DZAxis)) ); CALL_SUBTEST_4( alignedboxCastTests(AlignedBox3f()) ); - CALL_SUBTEST_5( alignedbox(AlignedBox4d()) ); + CALL_SUBTEST_5( (alignedboxNonIntegralRotatable(AlignedBox4d(), &rotate4DZWAxis)) ); CALL_SUBTEST_6( alignedboxCastTests(AlignedBox4d()) ); - CALL_SUBTEST_7( alignedbox(AlignedBox1d()) ); + CALL_SUBTEST_7( alignedboxTranslatable(AlignedBox1d()) ); CALL_SUBTEST_8( alignedboxCastTests(AlignedBox1d()) ); - CALL_SUBTEST_9( alignedbox(AlignedBox1i()) ); - CALL_SUBTEST_10( alignedbox(AlignedBox2i()) ); - CALL_SUBTEST_11( alignedbox(AlignedBox3i()) ); + CALL_SUBTEST_9( alignedboxTranslatable(AlignedBox1i()) ); + CALL_SUBTEST_10( (alignedboxRotatable(AlignedBox2i(), &rotate2DIntegral)) ); + CALL_SUBTEST_11( (alignedboxRotatable(AlignedBox3i(), &rotate3DZAxisIntegral)) ); CALL_SUBTEST_14( alignedbox(AlignedBox(4)) ); }