diff --git a/src/Magnum/MeshTools/BoundingVolume.cpp b/src/Magnum/MeshTools/BoundingVolume.cpp new file mode 100644 index 000000000..9e1bac909 --- /dev/null +++ b/src/Magnum/MeshTools/BoundingVolume.cpp @@ -0,0 +1,79 @@ +/* + This file is part of Magnum. + + Copyright © 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, + 2020, 2021, 2022 Vladimír Vondruš + Copyright © 2022 Pablo Escobar + + Permission is hereby granted, free of charge, to any person obtaining a + copy of this software and associated documentation files (the "Software"), + to deal in the Software without restriction, including without limitation + the rights to use, copy, modify, merge, publish, distribute, sublicense, + and/or sell copies of the Software, and to permit persons to whom the + Software is furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included + in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + DEALINGS IN THE SOFTWARE. +*/ + +#include "BoundingVolume.h" + +#include + +#include "Magnum/Math/FunctionsBatch.h" + +namespace Magnum { namespace MeshTools { + +Range3D boundingBoxAxisAligned(const Containers::StridedArrayView1D& points) { + return Math::minmax(points); +} + +Containers::Pair boundingSphereBouncingBubble(const Containers::StridedArrayView1D& points) { + /* See comment about radius below, this is done for consistency */ + if(points.isEmpty()) return {{}, Math::TypeTraits::epsilon()}; + + /** @todo Skip NaNs here? To match behaviour of boundingBoxAxisAligned() + which uses minmax(). */ + Vector3 center = points[0]; + /* Radius ends up in the denominator in the first loop so we can't + initialize it to 0. Unfortunately this also means the returned radius is + always above or equal to epsilon. */ + Float radius = Math::TypeTraits::epsilon(); + Float radiusSquared = radius*radius; + + for(int i = 0; i < 2; ++i) { + for(const Vector3& p : points) { + const Float ds = (p - center).dot(); + if(ds > radiusSquared) { + const Float alphaInv = radius/Math::sqrt(ds); + const Float alphaSqInv = alphaInv*alphaInv; + radius = (1.0f/alphaInv + alphaInv)*0.5f*radius; + center = ((1.0f + alphaSqInv)*center + (1.0f - alphaSqInv)*p)*0.5f; + radiusSquared = radius*radius; + } + } + } + + for(const Vector3& p : points) { + const Vector3 diff = p - center; + const Float ds = diff.dot(); + if(ds > radiusSquared) { + const Float d = Math::sqrt(ds); + radius = (radius + d)*0.5f; + center = center + ((d - radius)/d*diff); + radiusSquared = radius*radius; + } + } + + return {center, radius}; +} + +}} diff --git a/src/Magnum/MeshTools/BoundingVolume.h b/src/Magnum/MeshTools/BoundingVolume.h new file mode 100644 index 000000000..c3e1ffe9d --- /dev/null +++ b/src/Magnum/MeshTools/BoundingVolume.h @@ -0,0 +1,72 @@ +#ifndef Magnum_MeshTools_BoundingVolume_h +#define Magnum_MeshTools_BoundingVolume_h +/* + This file is part of Magnum. + + Copyright © 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, + 2020, 2021, 2022 Vladimír Vondruš + Copyright © 2022 Pablo Escobar + + Permission is hereby granted, free of charge, to any person obtaining a + copy of this software and associated documentation files (the "Software"), + to deal in the Software without restriction, including without limitation + the rights to use, copy, modify, merge, publish, distribute, sublicense, + and/or sell copies of the Software, and to permit persons to whom the + Software is furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included + in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + DEALINGS IN THE SOFTWARE. +*/ + +/** @file + * @brief Function @ref Magnum::MeshTools::boundingBoxAxisAligned(), @ref Magnum::MeshTools::boundingSphereBouncingBubble() + * @m_since_latest + */ + +#include + +#include "Magnum/Magnum.h" +#include "Magnum/Math/Range.h" +#include "Magnum/Math/Vector3.h" +#include "Magnum/MeshTools/visibility.h" + +namespace Magnum { namespace MeshTools { + +/** +@brief Calculate an axis-aligned bounding box +@param positions Vertex positions +@return Bounding box +@m_since_latest + +Same as @ref Math::minmax(const Corrade::Containers::StridedArrayView1D&). +*/ +MAGNUM_MESHTOOLS_EXPORT Range3D boundingBoxAxisAligned(const Containers::StridedArrayView1D& positions); + +/** +@brief Calculate an approximate bounding sphere using the Bouncing Bubble + algorithm +@param positions Vertex positions +@return Sphere center and radius +@m_since_latest + +The resulting bounding sphere is not usually minimal. According to the author, +a 1% to 2% error can be expected. Due to the nature of the algorithm, the +radius is never below @ref Math::TypeTraits::epsilon(), even for empty or +entirely overlapping lists of points. NaNs are ignored, unless the +first position is NaN in which case it is propagated. Algorithm used: +* *Bo Tian --- Bouncing Bubble: A fast algorithm for Minimal Enclosing Ball +problem, 2012, https://www.grin.com/document/204869* +*/ +MAGNUM_MESHTOOLS_EXPORT Containers::Pair boundingSphereBouncingBubble(const Containers::StridedArrayView1D& positions); + +}} + +#endif diff --git a/src/Magnum/MeshTools/CMakeLists.txt b/src/Magnum/MeshTools/CMakeLists.txt index 326df3dd5..04a68d66c 100644 --- a/src/Magnum/MeshTools/CMakeLists.txt +++ b/src/Magnum/MeshTools/CMakeLists.txt @@ -33,6 +33,7 @@ set(MagnumMeshTools_SRCS # Files compiled with different flags for main library and unit test library set(MagnumMeshTools_GracefulAssert_SRCS + BoundingVolume.cpp Combine.cpp CompressIndices.cpp Concatenate.cpp @@ -47,6 +48,7 @@ set(MagnumMeshTools_GracefulAssert_SRCS Transform.cpp) set(MagnumMeshTools_HEADERS + BoundingVolume.h Combine.h CompressIndices.h Concatenate.h diff --git a/src/Magnum/MeshTools/Test/BoundingVolumeTest.cpp b/src/Magnum/MeshTools/Test/BoundingVolumeTest.cpp new file mode 100644 index 000000000..57b80c86c --- /dev/null +++ b/src/Magnum/MeshTools/Test/BoundingVolumeTest.cpp @@ -0,0 +1,264 @@ +/* + This file is part of Magnum. + + Copyright © 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, + 2020, 2021, 2022 Vladimír Vondruš + Copyright © 2022 Pablo Escobar + + Permission is hereby granted, free of charge, to any person obtaining a + copy of this software and associated documentation files (the "Software"), + to deal in the Software without restriction, including without limitation + the rights to use, copy, modify, merge, publish, distribute, sublicense, + and/or sell copies of the Software, and to permit persons to whom the + Software is furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included + in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + DEALINGS IN THE SOFTWARE. +*/ + +#include +#include +#include +#include + +#include "Magnum/Math/FunctionsBatch.h" +#include "Magnum/Math/TypeTraits.h" +#include "Magnum/Math/Angle.h" +#include "Magnum/Math/Vector3.h" +#include "Magnum/MeshTools/BoundingVolume.h" +#include "Magnum/MeshTools/Reference.h" +#include "Magnum/MeshTools/Transform.h" +#include "Magnum/Primitives/Capsule.h" +#include "Magnum/Primitives/Cube.h" +#include "Magnum/Primitives/Icosphere.h" +#include "Magnum/Trade/MeshData.h" + +namespace Magnum { namespace MeshTools { namespace Test { namespace { + +struct BoundingVolumeTest: TestSuite::Tester { + explicit BoundingVolumeTest(); + + void boxAxisAligned(); + void boxAxisAlignedNaN(); + + void sphereBouncingBubble(); + void sphereBouncingBubbleNaN(); + + void benchmarkBoxAxisAligned(); + void benchmarkSphereBouncingBubble(); +}; + +BoundingVolumeTest::BoundingVolumeTest() { + addTests({&BoundingVolumeTest::boxAxisAligned, + &BoundingVolumeTest::boxAxisAlignedNaN, + &BoundingVolumeTest::sphereBouncingBubble, + &BoundingVolumeTest::sphereBouncingBubbleNaN}); + + addBenchmarks({&BoundingVolumeTest::benchmarkBoxAxisAligned, + &BoundingVolumeTest::benchmarkSphereBouncingBubble}, 150); +} + +void BoundingVolumeTest::boxAxisAligned() { + /* boundingboxAxisAligned() is just a wrapper around minmax() so only test + that the input and output are forwarded correctly */ + constexpr Float cylinderLength = 7.0f; + const Trade::MeshData cylinderMesh = Primitives::capsule3DSolid(3, 1, 12, cylinderLength*0.5f); + const Range3D box = MeshTools::boundingBoxAxisAligned( + cylinderMesh.attribute(Trade::MeshAttribute::Position)); + + CORRADE_COMPARE(box.center(), (Vector3{})); + CORRADE_COMPARE(box.size(), (Vector3{2.0f, cylinderLength + 2.0f, 2.0f})); +} + +void BoundingVolumeTest::boxAxisAlignedNaN() { + /* NaNs are skipped (unless it's all NaNs), matching minmax() behaviour */ + { + const Range3D box = MeshTools::boundingBoxAxisAligned(Containers::stridedArrayView({ + Vector3{Constants::nan()}, + Vector3{1.0f, 1.0f, 1.0f}, + Vector3{Constants::nan()}, + Vector3{2.0f, 2.0f, 2.0f}, + Vector3{Constants::nan()} + })); + CORRADE_COMPARE(box.min(), (Vector3{1.0f, 1.0f, 1.0f})); + CORRADE_COMPARE(box.max(), (Vector3{2.0f, 2.0f, 2.0f})); + } { + const Range3D box = MeshTools::boundingBoxAxisAligned(Containers::stridedArrayView({ + Vector3{Constants::nan()}, + Vector3{Constants::nan()} + })); + CORRADE_VERIFY(Math::isNan(box.min())); + CORRADE_VERIFY(Math::isNan(box.max())); + } +} + +void BoundingVolumeTest::sphereBouncingBubble() { + /* Empty positions -- produces radius epsilon for consistency with all + all identical positions */ + { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::StridedArrayView1D{}); + CORRADE_COMPARE(sphere.first(), (Vector3{})); + CORRADE_COMPARE(sphere.second(), Math::TypeTraits::epsilon()); + + /* Identical positions -- radius is always >= epsilon due to specifics of + the algorithm */ + } { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{1.0f, 2.0f, 3.0f} + })); + CORRADE_COMPARE(sphere.first(), (Vector3{1.0f, 2.0f, 3.0f})); + CORRADE_COMPARE(sphere.second(), Math::TypeTraits::epsilon()); + } { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{3.0f, 1.0f, 2.0f}, + Vector3{3.0f, 1.0f, 2.0f}, + Vector3{3.0f, 1.0f, 2.0f} + })); + CORRADE_COMPARE(sphere.first(), (Vector3{3.0f, 1.0f, 2.0f})); + CORRADE_COMPARE(sphere.second(), Math::TypeTraits::epsilon()); + + /* Simple cases */ + } { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{1.0f, 1.0f, 1.0f}, + Vector3{2.0f, 2.0f, 2.0f} + })); + CORRADE_COMPARE(sphere.first(), (Vector3{1.5f})); + CORRADE_COMPARE(sphere.second(), Vector3{0.5f}.length()); + } { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{2.0f, 0.0f, 0.0f}, + Vector3{-2.0f, 0.0f, 0.0f}, + Vector3{0.0f, 2.0f, 0.0f}, + Vector3{0.0f, -2.0f, 0.0f} + })); + CORRADE_COMPARE(sphere.first(), (Vector3{0.0f, 0.0f, 0.0f})); + CORRADE_COMPARE(sphere.second(), 2.0f); + + /* Icosphere -- original/translated and scaled */ + } { + const Trade::MeshData sphereMesh = Primitives::icosphereSolid(1); + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(sphereMesh.attribute(Trade::MeshAttribute::Position)); + /* No error */ + CORRADE_COMPARE(sphere.first(), (Vector3{0.0f, 0.0f, 0.0f})); + CORRADE_COMPARE(sphere.second(), 1.0f); + } { + Trade::MeshData sphereMesh = Primitives::icosphereSolid(1); + constexpr Vector3 translation{1.0f, 2.0f, 3.0f}; + constexpr Vector3 scale{0.5f, 1.2f, 2.8f}; + MeshTools::transform3DInPlace(sphereMesh, Matrix4::translation(translation)*Matrix4::scaling(scale)); + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(sphereMesh.attribute(Trade::MeshAttribute::Position)); + /* Noticeable error */ + constexpr Float Delta = 0.04f; + CORRADE_COMPARE_WITH(sphere.first(), translation, TestSuite::Compare::around(Vector3{Delta})); + /* Radius should never be smaller than the ground truth */ + CORRADE_COMPARE_WITH(sphere.second(), scale.max() + Delta, TestSuite::Compare::around(Delta)); + + /* Cube -- translated and scaled */ + } { + Trade::MeshData cubeMesh = MeshTools::owned(Primitives::cubeSolid()); + constexpr Vector3 translation{1.0f, 2.0f, 3.0f}; + constexpr Float scale = 13.2f; + MeshTools::transform3DInPlace(cubeMesh, Matrix4::translation(translation)*Matrix4::scaling(Vector3{scale})); + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(cubeMesh.attribute(Trade::MeshAttribute::Position)); + /* Noticeable error */ + constexpr Float Delta = 0.04f; + CORRADE_COMPARE_WITH(sphere.first(), translation, TestSuite::Compare::around(Vector3{Delta})); + CORRADE_COMPARE_WITH(sphere.second(), Constants::sqrt3()*scale + Delta, + TestSuite::Compare::around(Delta)); + } + + /* Radius is rotationally invariant */ + using namespace Math::Literals; + + for(Deg degrees = 0.0_degf; degrees < 360.0_degf; degrees += 60.0_degf) { + CORRADE_ITERATION(degrees); + Trade::MeshData cubeMesh = MeshTools::owned(Primitives::cubeSolid()); + constexpr Vector3 translation{1.0f, 2.0f, 3.0f}; + MeshTools::transform3DInPlace(cubeMesh, Matrix4::rotationY(degrees)*Matrix4::translation(translation)); + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(cubeMesh.attribute(Trade::MeshAttribute::Position)); + CORRADE_COMPARE(sphere.second(), Constants::sqrt3()); + } +} + +void BoundingVolumeTest::sphereBouncingBubbleNaN() { + /* NaN is ignored except for the first position element */ + { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{1.0f, 1.0f, 1.0f}, + Vector3{Constants::nan()}, + Vector3{2.0f, 2.0f, 2.0f}, + Vector3{Constants::nan()} + })); + CORRADE_COMPARE(sphere.first(), (Vector3{1.5f})); + CORRADE_COMPARE(sphere.second(), Vector3{0.5f}.length()); + + } { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(Containers::stridedArrayView({ + Vector3{Constants::nan()}, + Vector3{1.0f, 1.0f, 1.0f}, + Vector3{2.0f, 2.0f, 2.0f} + })); + { + CORRADE_EXPECT_FAIL("NaN in the first position is not ignored."); + CORRADE_COMPARE(sphere.first(), (Vector3{1.5f})); + CORRADE_COMPARE(sphere.second(), Vector3{0.5f}.length()); + } + CORRADE_VERIFY(Math::isNan(sphere.first())); + CORRADE_COMPARE(sphere.second(), Math::TypeTraits::epsilon()); + } +} + +void BoundingVolumeTest::benchmarkBoxAxisAligned() { + Containers::Array points{NoInit, 500}; + for(size_t i = 0; i < points.size(); ++i) { + points[i] = Vector3{Float(i)*0.01f}; + } + + Float r = 0.0f; + CORRADE_BENCHMARK(50) { + const Range3D box = MeshTools::boundingBoxAxisAligned(points); + r += box.size().x(); + } + + CORRADE_COMPARE_AS(r, 1.0f, TestSuite::Compare::Greater); +} + +void BoundingVolumeTest::benchmarkSphereBouncingBubble() { + Containers::Array points{NoInit, 500}; + for(size_t i = 0; i < points.size(); ++i) { + points[i] = Vector3{Float(i)*0.01f}; + } + + Float r = 0.0f; + CORRADE_BENCHMARK(50) { + const Containers::Pair sphere = + MeshTools::boundingSphereBouncingBubble(points); + r += sphere.second(); + } + + CORRADE_COMPARE_AS(r, 1.0f, TestSuite::Compare::Greater); +} + +}}}} + +CORRADE_TEST_MAIN(Magnum::MeshTools::Test::BoundingVolumeTest) diff --git a/src/Magnum/MeshTools/Test/CMakeLists.txt b/src/Magnum/MeshTools/Test/CMakeLists.txt index a6de05189..a221a5d7e 100644 --- a/src/Magnum/MeshTools/Test/CMakeLists.txt +++ b/src/Magnum/MeshTools/Test/CMakeLists.txt @@ -27,6 +27,7 @@ # property that would have to be set on each target separately. set(CMAKE_FOLDER "Magnum/MeshTools/Test") +corrade_add_test(MeshToolsBoundingVolumeTest BoundingVolumeTest.cpp LIBRARIES MagnumMeshToolsTestLib MagnumPrimitives) corrade_add_test(MeshToolsCombineTest CombineTest.cpp LIBRARIES MagnumMeshToolsTestLib) corrade_add_test(MeshToolsCompressIndicesTest CompressIndicesTest.cpp LIBRARIES MagnumMeshToolsTestLib) corrade_add_test(MeshToolsConcatenateTest ConcatenateTest.cpp LIBRARIES MagnumMeshToolsTestLib)