From 4ba32517575d19a1e285576e7e73347455896fcd Mon Sep 17 00:00:00 2001 From: Ilya Glushchenko Date: Sun, 25 Nov 2018 20:08:11 +0300 Subject: [PATCH] [#27] QHull optimization and 2019 copyrights update --- Epona/include/Epona/Analysis.hpp | 4 +- Epona/include/Epona/Debug.hpp | 2 +- Epona/include/Epona/FloatingPoint.hpp | 2 +- Epona/include/Epona/HalfEdgeDataStructure.hpp | 102 ++-- Epona/include/Epona/HyperPlane.hpp | 2 +- Epona/include/Epona/JacobiEigenvalue.hpp | 8 +- Epona/include/Epona/QuickhullConvexHull.hpp | 493 +++++++++--------- Epona/include/Epona/debug/DebugDummy.hpp | 2 +- .../Epona/debug/DebugImplementation.hpp | 2 +- Epona/sources/Analysis.cpp | 2 +- Epona/sources/HalfEdgeDataStructure.cpp | 20 +- Epona/sources/HyperPlane.cpp | 2 +- demo/Main.cpp | 2 +- test/AnalysisTests.cpp | 2 +- 14 files changed, 317 insertions(+), 328 deletions(-) diff --git a/Epona/include/Epona/Analysis.hpp b/Epona/include/Epona/Analysis.hpp index a0b41d1..9b6a7e4 100644 --- a/Epona/include/Epona/Analysis.hpp +++ b/Epona/include/Epona/Analysis.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ @@ -377,7 +377,7 @@ glm::vec3 CalculateBarycentricCoordinates(glm::vec3 p, glm::vec3 a, glm::vec3 b, /* * @brief Calculate area of the triangle and returns true if area is equal to zero - * + * * @param a triangle's point * @param b triangle's point * @param c triangle's point diff --git a/Epona/include/Epona/Debug.hpp b/Epona/include/Epona/Debug.hpp index 7433606..90704ed 100644 --- a/Epona/include/Epona/Debug.hpp +++ b/Epona/include/Epona/Debug.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/include/Epona/FloatingPoint.hpp b/Epona/include/Epona/FloatingPoint.hpp index 9c97c15..ae8a9aa 100644 --- a/Epona/include/Epona/FloatingPoint.hpp +++ b/Epona/include/Epona/FloatingPoint.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/include/Epona/HalfEdgeDataStructure.hpp b/Epona/include/Epona/HalfEdgeDataStructure.hpp index 2edec70..50eff83 100644 --- a/Epona/include/Epona/HalfEdgeDataStructure.hpp +++ b/Epona/include/Epona/HalfEdgeDataStructure.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ @@ -28,9 +28,9 @@ class HalfEdgeDataStructure HalfEdgeDataStructure() = default; HalfEdgeDataStructure(HalfEdgeDataStructure const&) = delete; - + HalfEdgeDataStructure& operator=(HalfEdgeDataStructure const& heds) = delete; - + HalfEdgeDataStructure(HalfEdgeDataStructure&& heds) { m_halfEdgeList = std::move(heds.m_halfEdgeList); @@ -52,7 +52,7 @@ class HalfEdgeDataStructure m_faceIteratorMap = std::move(heds.m_faceIteratorMap); m_faceVerticesIteratorMap = std::move(heds.m_faceVerticesIteratorMap); - return *this; + return *this; } /** @@ -315,11 +315,11 @@ class HalfEdgeDataStructure { return AdjacentFaceCirculator(*this); } - - operator HalfEdge*() const - { - return m_currentHalfEdge; - } + + operator HalfEdge*() const + { + return m_currentHalfEdge; + } private: HalfEdge* m_currentHalfEdge; @@ -355,11 +355,8 @@ class HalfEdgeDataStructure */ const_face_iterator GetAdjacentFaceIterator() const; - //! Stores face hyperplane - HyperPlane hyperPlane; - - //! Stores face index - uint32_t index; + //! Stores face index + uint32_t index; private: HalfEdge* m_halfEdge; @@ -400,48 +397,48 @@ class HalfEdgeDataStructure uint64_t vertexIndex; }; - /** - * @brief Inserts face into the current data structure - * @param[in] a vertex index - * @param[in] b vertex index - * @param[in] c vertex index - */ - void MakeFace(uint64_t a, uint64_t b, uint64_t c) - { - MakeFace(a, b, c, {}, 0); - } + /** + * @brief Inserts face into the current data structure + * @param[in] a vertex index + * @param[in] b vertex index + * @param[in] c vertex index + */ + void MakeFace(uint64_t a, uint64_t b, uint64_t c) + { + MakeFace(a, b, c, 0); + } /** * @brief Inserts face into the current data structure * @param a vertex index * @param b vertex index * @param c vertex index - * @param hp hyperplane containing input vertices - * @param index outside face index + * @param hp hyperplane containing input vertices + * @param index outside face index */ - void MakeFace(uint64_t a, uint64_t b, uint64_t c, HyperPlane hp, uint32_t index); - - /** - * @brief Inserts face into the current data structure - * @type T array-like type with overaloded operator[] - * @param index face indices array of size 3 - */ - template < typename T > - decltype(auto) GetFace(T index) - { - return GetFace(index[0], index[1], index[2]); - } - - /** - * @brief Inserts face into the current data structure - * @type T array-like type with overaloded operator[] - * @param index face indices array of size 3 - */ - template < typename T > - decltype(auto) GetFace(T index) const - { - return GetFace(index[0], index[1], index[2]); - } + void MakeFace(uint64_t a, uint64_t b, uint64_t c, uint32_t index); + + /** + * @brief Inserts face into the current data structure + * @type T array-like type with overaloded operator[] + * @param index face indices array of size 3 + */ + template < typename T > + decltype(auto) GetFace(T index) + { + return GetFace(index[0], index[1], index[2]); + } + + /** + * @brief Inserts face into the current data structure + * @type T array-like type with overaloded operator[] + * @param index face indices array of size 3 + */ + template < typename T > + decltype(auto) GetFace(T index) const + { + return GetFace(index[0], index[1], index[2]); + } /** * @brief Returns a face iterator @@ -538,12 +535,5 @@ class HalfEdgeDataStructure ); }; -template < typename T > -std::tuple GetRidge(HalfEdgeDataStructure::face_iterator it) -{ - - - return {}; -} } // namespace epona #endif // EPONA_HEDS_HPP diff --git a/Epona/include/Epona/HyperPlane.hpp b/Epona/include/Epona/HyperPlane.hpp index 58b4b63..1d6f966 100644 --- a/Epona/include/Epona/HyperPlane.hpp +++ b/Epona/include/Epona/HyperPlane.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/include/Epona/JacobiEigenvalue.hpp b/Epona/include/Epona/JacobiEigenvalue.hpp index daba690..c5664cb 100644 --- a/Epona/include/Epona/JacobiEigenvalue.hpp +++ b/Epona/include/Epona/JacobiEigenvalue.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ @@ -105,7 +105,7 @@ inline std::tuple CalculateJacobiEigenvectorsEigenvalue( uint8_t j; uint16_t iterations = 0; - do + do { std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix); @@ -114,7 +114,7 @@ inline std::tuple CalculateJacobiEigenvectorsEigenvalue( symmetricMatrix = (glm::transpose(rotationMatrix) * symmetricMatrix) * rotationMatrix; std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix); - } + } while ((++iterations < maxIterations) && (glm::abs(symmetricMatrix[i][j]) > coverageThreshold)); return std::tuple{ eigenvectors, glm::vec3{ symmetricMatrix[0][0], symmetricMatrix[1][1], symmetricMatrix[2][2] } }; @@ -160,7 +160,7 @@ inline glm::mat3 CalculateJacobiEigenvectors( glm::mat3 symmetricMatrix, float coverageThreshold = epona::fp::g_floatingPointThreshold, uint32_t maxIterations = 100 ) { - glm::mat3 eigenvectors(1); + glm::mat3 eigenvectors(1); glm::mat3 rotationMatrix(1); uint8_t i = 0; uint8_t j = 0; diff --git a/Epona/include/Epona/QuickhullConvexHull.hpp b/Epona/include/Epona/QuickhullConvexHull.hpp index 8ecad32..115dae3 100644 --- a/Epona/include/Epona/QuickhullConvexHull.hpp +++ b/Epona/include/Epona/QuickhullConvexHull.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ @@ -28,8 +28,9 @@ namespace epona struct ConvexHull { - std::vector indices; + std::vector indices; std::vector faces; + std::vector planes; HalfEdgeDataStructure heds; }; @@ -40,23 +41,23 @@ namespace { struct FaceOutliers { - std::vector> outliers; - std::vector extremeIndices; - std::vector extremeSignedDistances; + std::vector> outliers; + std::vector extremeIndices; + std::vector extremeSignedDistances; }; struct Ridge { - uint32_t aVertex; - uint32_t bVertex; + uint32_t aVertex; + uint32_t bVertex; - bool operator==(Ridge other) const { - return (aVertex == other.aVertex && bVertex == other.bVertex) - || (aVertex == other.bVertex && bVertex == other.aVertex); - } + bool operator==(Ridge other) const { + return (aVertex == other.aVertex && bVertex == other.bVertex) + || (aVertex == other.bVertex && bVertex == other.aVertex); + } - bool operator==(uint32_t index) const { - return aVertex == index || bVertex == index; - } + bool operator==(uint32_t index) const { + return aVertex == index || bVertex == index; + } }; /** @@ -122,34 +123,36 @@ epona::ConvexHull CalculateConvexHullInitialTetrahedron(std::vector c } } - glm::vec3 mean(0); - for (uint32_t index : indices) - mean += vertices[index]; - mean /= 4; - - assert( !epona::OneLine(vertices[indices[0]], vertices[indices[1]], vertices[indices[2]]) ); - assert( !epona::OneLine(vertices[indices[0]], vertices[indices[1]], vertices[indices[3]]) ); - assert( !epona::OneLine(vertices[indices[0]], vertices[indices[2]], vertices[indices[3]]) ); - assert( !epona::OneLine(vertices[indices[1]], vertices[indices[2]], vertices[indices[3]]) ); - - epona::HalfEdgeDataStructure heds; - heds.MakeFace(indices[0], indices[1], indices[2], - epona::HyperPlane{vertices[indices[0]], vertices[indices[1]], vertices[indices[2]], &mean }, 0 ); - heds.MakeFace(indices[0], indices[1], indices[3], - epona::HyperPlane{vertices[indices[0]], vertices[indices[1]], vertices[indices[3]], &mean }, 1 ); - heds.MakeFace(indices[0], indices[2], indices[3], - epona::HyperPlane{vertices[indices[0]], vertices[indices[2]], vertices[indices[3]], &mean }, 2 ); - heds.MakeFace(indices[1], indices[2], indices[3], - epona::HyperPlane{vertices[indices[1]], vertices[indices[2]], vertices[indices[3]], &mean }, 3 ); + glm::vec3 mean(0); + for (uint32_t index : indices) + mean += vertices[index]; + mean /= 4; + + assert( !epona::OneLine(vertices[indices[0]], vertices[indices[1]], vertices[indices[2]]) ); + assert( !epona::OneLine(vertices[indices[0]], vertices[indices[1]], vertices[indices[3]]) ); + assert( !epona::OneLine(vertices[indices[0]], vertices[indices[2]], vertices[indices[3]]) ); + assert( !epona::OneLine(vertices[indices[1]], vertices[indices[2]], vertices[indices[3]]) ); + + epona::HalfEdgeDataStructure heds; + heds.MakeFace(indices[0], indices[1], indices[2], 0 ); + heds.MakeFace(indices[0], indices[1], indices[3], 1 ); + heds.MakeFace(indices[0], indices[2], indices[3], 2 ); + heds.MakeFace(indices[1], indices[2], indices[3], 3 ); //Create initial convex hull tetrahedron return {{ indices[0], indices[1], indices[2], indices[3] }, - { glm::u32vec3{indices[0], indices[1], indices[2]}, - glm::u32vec3{indices[0], indices[1], indices[3]}, - glm::u32vec3{indices[0], indices[2], indices[3]}, - glm::u32vec3{indices[1], indices[2], indices[3]} }, - std::move(heds) - }; + { glm::u32vec3{indices[0], indices[1], indices[2]}, + glm::u32vec3{indices[0], indices[1], indices[3]}, + glm::u32vec3{indices[0], indices[2], indices[3]}, + glm::u32vec3{indices[1], indices[2], indices[3]} }, + { + epona::HyperPlane{vertices[indices[0]], vertices[indices[1]], vertices[indices[2]], &mean }, + epona::HyperPlane{vertices[indices[0]], vertices[indices[1]], vertices[indices[3]], &mean }, + epona::HyperPlane{vertices[indices[0]], vertices[indices[2]], vertices[indices[3]], &mean }, + epona::HyperPlane{vertices[indices[1]], vertices[indices[2]], vertices[indices[3]], &mean }, + }, + std::move(heds) + }; } } // namespace :: @@ -159,216 +162,220 @@ namespace epona inline bool RecalculateConvexHull(ConvexHull& hull, std::vector const& vertices, uint32_t maxIterations = 100, float distanceThreshold = 1e-3f) { - bool changed = false; - - std::vector indices; - indices.reserve(vertices.size()); - for (uint32_t v = 0; v < vertices.size(); ++v) - if (std::find(hull.indices.begin(), hull.indices.end(), v) == hull.indices.end()) - indices.push_back(v); - std::sort(hull.indices.begin(), hull.indices.end(), std::greater()); - - uint32_t totalOutliersCount = 0; - ::FaceOutliers faceOutliers; - std::vector visibleFaces; - std::vector<::Ridge> boundary; - std::vector temporaryOutliers; - - do { - glm::vec3 mean(0); - for (auto i : hull.indices) - mean += vertices[i]; - mean /= hull.indices.size(); - - //Find outliers and extreme vertices - uint32_t newOutliersCount = 0; - faceOutliers.outliers.resize(hull.faces.size()); - faceOutliers.extremeIndices.resize(hull.faces.size(), UINT32_MAX); - faceOutliers.extremeSignedDistances.resize(hull.faces.size(), -FLT_MAX); - for (uint32_t v : indices) - { - for (uint32_t f = 0; f < hull.faces.size(); ++f) - { - float const signedDistance = hull.heds.GetFace(hull.faces[f])->hyperPlane.SignedDistance(vertices[v]); - if (signedDistance > distanceThreshold) { - faceOutliers.outliers[f].push_back(v); - newOutliersCount++; - if (signedDistance > faceOutliers.extremeSignedDistances[f]) { - faceOutliers.extremeSignedDistances[f] = signedDistance; - faceOutliers.extremeIndices[f] = v; - } - break; - } - } - } - - //Find visible faces from the first extreme outlier - uint32_t extremeOutlierIndex = UINT32_MAX; - visibleFaces.clear(); - visibleFaces.reserve(hull.faces.size() / 2); - for (uint32_t f = 0; f < hull.faces.size(); ++f) - { - if (!faceOutliers.outliers[f].empty()) - { - visibleFaces.push_back(f); - std::vector visitedFaces(1, f); - visitedFaces.reserve(hull.faces.size() / 2); - glm::vec3 const extremeVertex = vertices[faceOutliers.extremeIndices[f]]; - - for (uint32_t vf = 0; vf < visibleFaces.size(); ++vf) - { - auto adjFaceIt = hull.heds.GetFace(hull.faces[visibleFaces[vf]])->GetAdjacentFaceIterator(); - - for (uint8_t i = 0; i < 3; ++i, ++adjFaceIt) { - //NOTE: will work faster then open address map till ~ visitedFaces.size() < 100 - if (std::find(visitedFaces.begin(), visitedFaces.end(), adjFaceIt->index) == visitedFaces.end()) - { - visitedFaces.push_back(adjFaceIt->index); - float const signedDistance = adjFaceIt->hyperPlane.SignedDistance(extremeVertex); - if (fp::IsGreater(signedDistance, 0.0f)) - visibleFaces.push_back(adjFaceIt->index); - } - } - } - extremeOutlierIndex = faceOutliers.extremeIndices[f]; - break; - } - } - std::sort(visibleFaces.begin(), visibleFaces.end()); - - if (visibleFaces.empty()) - { - indices.clear(); - } - else - { - changed = true; - assert(extremeOutlierIndex != UINT32_MAX); - - //Find boundary ridges - boundary.clear(); - boundary.reserve(visibleFaces.size()); - - for (uint32_t vf : visibleFaces) - { - auto adjFaceIt = hull.heds.GetFace(hull.faces[vf])->GetAdjacentFaceIterator(); - for (uint8_t i = 0; i < 3; ++i, ++adjFaceIt) - { - if (std::find(visibleFaces.begin(), visibleFaces.end(), adjFaceIt->index) == visibleFaces.end()) - { - auto halfEdge = static_cast(adjFaceIt); - ::Ridge const ridge{ static_cast(halfEdge->vertexIndex), static_cast(halfEdge->prev->vertexIndex) }; - //NOTE: will work faster then open address map till ~ boundary.size() < 100 - if (std::find(boundary.begin(), boundary.end(), ridge) == boundary.end()) - { - boundary.push_back(ridge); - } - } - } - } - //FixMe: Checks for degenerate case, that should be prevented by threshold check, is there a way to make it better? - if (std::find(boundary.begin(), boundary.end(), extremeOutlierIndex) != boundary.end()) - return true; - - //Collect inliers as outliers and remove them from hull indices - temporaryOutliers.clear(); - temporaryOutliers.reserve(visibleFaces.size() * 3 + newOutliersCount); - for (uint32_t vf : visibleFaces) { - auto face = hull.faces[vf]; - for (uint8_t i = 0; i < 3; ++i) { - //NOTE: will work faster then open address map till ~ (boundary/outliers).size() < 100 - if (std::find(boundary.begin(), boundary.end(), face[i]) == boundary.end() - && std::find(temporaryOutliers.begin(), temporaryOutliers.end(), face[i]) == temporaryOutliers.end()) - { - temporaryOutliers.push_back(face[i]); - //NOTE: will work faster then open address map till ~ hull.indieces.size() < 100 - auto it = std::find(hull.indices.begin(), hull.indices.end(), face[i]); - if (it != hull.indices.end()) - hull.indices.erase(it); - } - } - } - - //Remove visible faces and collect unassigned outliers - for (uint32_t vf : visibleFaces) - { - hull.heds.RemoveFace(hull.heds.GetFace(hull.faces[vf])); - temporaryOutliers.insert(temporaryOutliers.end(), - std::make_move_iterator(faceOutliers.outliers[vf].begin()), - std::make_move_iterator(faceOutliers.outliers[vf].end()) - ); - } - indices.swap(temporaryOutliers); - temporaryOutliers.clear(); - - //Create new faces and assign outliers - for (::Ridge ridge : boundary) - { - uint32_t faceIndex = static_cast(hull.faces.size()); - if (!visibleFaces.empty()) - { - faceIndex = visibleFaces.back(); - hull.faces[faceIndex] = { extremeOutlierIndex, ridge.aVertex, ridge.bVertex }; - visibleFaces.pop_back(); - - faceOutliers.extremeIndices[faceIndex] = UINT32_MAX; - faceOutliers.extremeSignedDistances[faceIndex] = -FLT_MAX; - faceOutliers.outliers[faceIndex].clear(); - } - else - { - hull.faces.emplace_back(extremeOutlierIndex, ridge.aVertex, ridge.bVertex); - } - hull.heds.MakeFace(extremeOutlierIndex, ridge.aVertex, ridge.bVertex, - HyperPlane{ vertices[extremeOutlierIndex], vertices[ridge.aVertex], vertices[ridge.bVertex], &mean }, - faceIndex - ); - //NOTE: will work faster then open address map till ~ hull.indieces.size() < 100 - if (std::find(hull.indices.begin(), hull.indices.end(), extremeOutlierIndex) == hull.indices.end()) - hull.indices.push_back(extremeOutlierIndex); - } - - //Restore structure - if (!visibleFaces.empty()) - { - std::sort(visibleFaces.begin(), visibleFaces.end()); - std::vector faces; - FaceOutliers tmpFaceOutliers; - tmpFaceOutliers.extremeIndices.reserve(faceOutliers.extremeIndices.size()); - tmpFaceOutliers.extremeSignedDistances.reserve(faceOutliers.extremeSignedDistances.size()); - tmpFaceOutliers.outliers.reserve(faceOutliers.outliers.size()); - for (uint32_t i = 0, j = 0; i < hull.faces.size(); ++i) - { - if (j == visibleFaces.size() || i != visibleFaces[j]) - { - faces.push_back(hull.faces[i]); - hull.heds.GetFace(hull.faces[i])->index = static_cast(faces.size() - 1); - tmpFaceOutliers.extremeIndices.push_back(faceOutliers.extremeIndices[i]); - tmpFaceOutliers.extremeSignedDistances.push_back(faceOutliers.extremeSignedDistances[i]); - tmpFaceOutliers.outliers.push_back(std::move(faceOutliers.outliers[i])); - } - else { ++j; } - } - faceOutliers.extremeIndices.swap(tmpFaceOutliers.extremeIndices); - faceOutliers.extremeSignedDistances.swap(tmpFaceOutliers.extremeSignedDistances); - faceOutliers.outliers.swap(tmpFaceOutliers.outliers); - hull.faces.swap(faces); - } - } - - totalOutliersCount = newOutliersCount + static_cast(indices.size()); - for (auto& outliers : faceOutliers.outliers) - totalOutliersCount += static_cast(outliers.size()); - } while (totalOutliersCount && --maxIterations); - - return changed; + bool changed = false; + + std::vector indices; + indices.reserve(vertices.size()); + for (uint32_t v = 0; v < vertices.size(); ++v) + if (std::find(hull.indices.begin(), hull.indices.end(), v) == hull.indices.end()) + indices.push_back(v); + std::sort(hull.indices.begin(), hull.indices.end(), std::greater()); + + uint32_t totalOutliersCount = 0; + ::FaceOutliers faceOutliers; + std::vector visibleFaces; + std::vector<::Ridge> boundary; + std::vector temporaryOutliers; + + do { + glm::vec3 mean(0); + for (auto i : hull.indices) + mean += vertices[i]; + mean /= hull.indices.size(); + + //Find outliers and extreme vertices + uint32_t newOutliersCount = 0; + faceOutliers.outliers.resize(hull.faces.size()); + faceOutliers.extremeIndices.resize(hull.faces.size(), UINT32_MAX); + faceOutliers.extremeSignedDistances.resize(hull.faces.size(), -FLT_MAX); + for (uint32_t v : indices) + { + for (uint32_t f = 0; f < hull.faces.size(); ++f) + { + float const signedDistance = hull.planes[f].SignedDistance(vertices[v]); + if (signedDistance > distanceThreshold) { + faceOutliers.outliers[f].push_back(v); + newOutliersCount++; + if (signedDistance > faceOutliers.extremeSignedDistances[f]) { + faceOutliers.extremeSignedDistances[f] = signedDistance; + faceOutliers.extremeIndices[f] = v; + } + break; + } + } + } + + //Find visible faces from the first extreme outlier + uint32_t extremeOutlierIndex = UINT32_MAX; + visibleFaces.clear(); + visibleFaces.reserve(hull.faces.size() / 2); + for (uint32_t f = 0; f < hull.faces.size(); ++f) + { + if (!faceOutliers.outliers[f].empty()) + { + visibleFaces.push_back(f); + std::vector visitedFaces(1, f); + visitedFaces.reserve(hull.faces.size() / 2); + glm::vec3 const extremeVertex = vertices[faceOutliers.extremeIndices[f]]; + + for (uint32_t vf = 0; vf < visibleFaces.size(); ++vf) + { + auto adjFaceIt = hull.heds.GetFace(hull.faces[visibleFaces[vf]])->GetAdjacentFaceIterator(); + + for (uint8_t i = 0; i < 3; ++i, ++adjFaceIt) { + //NOTE: will work faster then open address map till ~ visitedFaces.size() < 100 + if (std::find(visitedFaces.begin(), visitedFaces.end(), adjFaceIt->index) == visitedFaces.end()) + { + visitedFaces.push_back(adjFaceIt->index); + float const signedDistance = hull.planes[adjFaceIt->index].SignedDistance(extremeVertex); + if (fp::IsGreater(signedDistance, 0.0f)) + visibleFaces.push_back(adjFaceIt->index); + } + } + } + extremeOutlierIndex = faceOutliers.extremeIndices[f]; + break; + } + } + std::sort(visibleFaces.begin(), visibleFaces.end()); + + if (visibleFaces.empty()) + { + indices.clear(); + } + else + { + changed = true; + assert(extremeOutlierIndex != UINT32_MAX); + + //Find boundary ridges + boundary.clear(); + boundary.reserve(visibleFaces.size()); + + for (uint32_t vf : visibleFaces) + { + auto adjFaceIt = hull.heds.GetFace(hull.faces[vf])->GetAdjacentFaceIterator(); + for (uint8_t i = 0; i < 3; ++i, ++adjFaceIt) + { + if (std::find(visibleFaces.begin(), visibleFaces.end(), adjFaceIt->index) == visibleFaces.end()) + { + auto halfEdge = static_cast(adjFaceIt); + ::Ridge const ridge{ static_cast(halfEdge->vertexIndex), static_cast(halfEdge->prev->vertexIndex) }; + //NOTE: will work faster then open address map till ~ boundary.size() < 100 + if (std::find(boundary.begin(), boundary.end(), ridge) == boundary.end()) + { + boundary.push_back(ridge); + } + } + } + } + //FixMe: Checks for degenerate case, that should be prevented by threshold check, is there a way to make it better? + if (std::find(boundary.begin(), boundary.end(), extremeOutlierIndex) != boundary.end()) + return true; + + //Collect inliers as outliers and remove them from hull indices + temporaryOutliers.clear(); + temporaryOutliers.reserve(visibleFaces.size() * 3 + newOutliersCount); + for (uint32_t vf : visibleFaces) { + auto face = hull.faces[vf]; + for (uint8_t i = 0; i < 3; ++i) { + //NOTE: will work faster then open address map till ~ (boundary/outliers).size() < 100 + if (std::find(boundary.begin(), boundary.end(), face[i]) == boundary.end() + && std::find(temporaryOutliers.begin(), temporaryOutliers.end(), face[i]) == temporaryOutliers.end()) + { + temporaryOutliers.push_back(face[i]); + //NOTE: will work faster then open address map till ~ hull.indieces.size() < 100 + auto it = std::find(hull.indices.begin(), hull.indices.end(), face[i]); + if (it != hull.indices.end()) + hull.indices.erase(it); + } + } + } + + //Remove visible faces and collect unassigned outliers + for (uint32_t vf : visibleFaces) + { + hull.heds.RemoveFace(hull.heds.GetFace(hull.faces[vf])); + temporaryOutliers.insert(temporaryOutliers.end(), + std::make_move_iterator(faceOutliers.outliers[vf].begin()), + std::make_move_iterator(faceOutliers.outliers[vf].end()) + ); + } + indices.swap(temporaryOutliers); + temporaryOutliers.clear(); + + //Create new faces and assign outliers + for (::Ridge ridge : boundary) + { + uint32_t faceIndex = static_cast(hull.faces.size()); + if (!visibleFaces.empty()) + { + faceIndex = visibleFaces.back(); + hull.faces[faceIndex] = { extremeOutlierIndex, ridge.aVertex, ridge.bVertex }; + hull.planes[faceIndex] = { vertices[extremeOutlierIndex], vertices[ridge.aVertex], vertices[ridge.bVertex], &mean }; + visibleFaces.pop_back(); + + faceOutliers.extremeIndices[faceIndex] = UINT32_MAX; + faceOutliers.extremeSignedDistances[faceIndex] = -FLT_MAX; + faceOutliers.outliers[faceIndex].clear(); + } + else + { + hull.faces.emplace_back(extremeOutlierIndex, ridge.aVertex, ridge.bVertex); + hull.planes.emplace_back(vertices[extremeOutlierIndex], vertices[ridge.aVertex], vertices[ridge.bVertex], &mean); + } + hull.heds.MakeFace(extremeOutlierIndex, ridge.aVertex, ridge.bVertex, faceIndex); + //NOTE: will work faster then open address map till ~ hull.indieces.size() < 100 + if (std::find(hull.indices.begin(), hull.indices.end(), extremeOutlierIndex) == hull.indices.end()) + hull.indices.push_back(extremeOutlierIndex); + } + + //Restore structure + if (!visibleFaces.empty()) + { + std::sort(visibleFaces.begin(), visibleFaces.end()); + std::vector faces; + std::vector planes; + faces.reserve(hull.faces.size() - visibleFaces.size()); + planes.reserve(hull.faces.size() - visibleFaces.size()); + FaceOutliers tmpFaceOutliers; + tmpFaceOutliers.extremeIndices.reserve(faceOutliers.extremeIndices.size()); + tmpFaceOutliers.extremeSignedDistances.reserve(faceOutliers.extremeSignedDistances.size()); + tmpFaceOutliers.outliers.reserve(faceOutliers.outliers.size()); + for (uint32_t i = 0, j = 0; i < hull.faces.size(); ++i) + { + if (j == visibleFaces.size() || i != visibleFaces[j]) + { + faces.push_back(hull.faces[i]); + planes.push_back(hull.planes[i]); + hull.heds.GetFace(hull.faces[i])->index = static_cast(faces.size() - 1); + tmpFaceOutliers.extremeIndices.push_back(faceOutliers.extremeIndices[i]); + tmpFaceOutliers.extremeSignedDistances.push_back(faceOutliers.extremeSignedDistances[i]); + tmpFaceOutliers.outliers.push_back(std::move(faceOutliers.outliers[i])); + } + else { ++j; } + } + faceOutliers.extremeIndices.swap(tmpFaceOutliers.extremeIndices); + faceOutliers.extremeSignedDistances.swap(tmpFaceOutliers.extremeSignedDistances); + faceOutliers.outliers.swap(tmpFaceOutliers.outliers); + hull.faces.swap(faces); + hull.planes.swap(planes); + } + } + + totalOutliersCount = newOutliersCount + static_cast(indices.size()); + for (auto& outliers : faceOutliers.outliers) + totalOutliersCount += static_cast(outliers.size()); + } while (totalOutliersCount && --maxIterations); + + return changed; } inline ConvexHull CalculateConvexHull(std::vector const& vertices, uint32_t maxIterations = 100, float distanceThreshold = 1e-3f) { ConvexHull hull = ::CalculateConvexHullInitialTetrahedron(vertices); - RecalculateConvexHull(hull, vertices, maxIterations); + RecalculateConvexHull(hull, vertices, maxIterations); return hull; } diff --git a/Epona/include/Epona/debug/DebugDummy.hpp b/Epona/include/Epona/debug/DebugDummy.hpp index 89e293d..bff2c6b 100644 --- a/Epona/include/Epona/debug/DebugDummy.hpp +++ b/Epona/include/Epona/debug/DebugDummy.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/include/Epona/debug/DebugImplementation.hpp b/Epona/include/Epona/debug/DebugImplementation.hpp index 0945325..4a3beb0 100644 --- a/Epona/include/Epona/debug/DebugImplementation.hpp +++ b/Epona/include/Epona/debug/DebugImplementation.hpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/sources/Analysis.cpp b/Epona/sources/Analysis.cpp index 3914b7f..76057f3 100644 --- a/Epona/sources/Analysis.cpp +++ b/Epona/sources/Analysis.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/Epona/sources/HalfEdgeDataStructure.cpp b/Epona/sources/HalfEdgeDataStructure.cpp index 2742dd8..77b4065 100644 --- a/Epona/sources/HalfEdgeDataStructure.cpp +++ b/Epona/sources/HalfEdgeDataStructure.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ @@ -36,24 +36,17 @@ HalfEdgeDataStructure::Face::const_face_iterator HalfEdgeDataStructure::Face::Ge bool HalfEdgeDataStructure::FaceVertices::operator==(FaceVertices const& other) const { - std::array pointers = {{a, b, c}}; - std::sort(pointers.begin(), pointers.end()); - std::array otherPointers = {{other.a, other.b, other.c}}; - std::sort(otherPointers.begin(), otherPointers.end()); - - return pointers[0] == otherPointers[0] - && pointers[1] == otherPointers[1] - && pointers[2] == otherPointers[2]; + return (a == other.a && b == other.b && c == other.c) + || (a == other.c && b == other.a && c == other.b) + || (a == other.b && b == other.c && c == other.a); } size_t HalfEdgeDataStructure::FaceVertices::Hasher::operator()(FaceVertices const& face) const { - return std::hash{}(face.a) - ^ std::hash{}(face.b) - ^ std::hash{}(face.c); + return static_cast(face.a ^ face.b ^ face.c); } -void HalfEdgeDataStructure::MakeFace(uint64_t a, uint64_t b, uint64_t c, HyperPlane hp, uint32_t index) +void HalfEdgeDataStructure::MakeFace(uint64_t a, uint64_t b, uint64_t c, uint32_t index) { FaceVertices const faceVerticesKey{a, b, c}; if (m_faceVerticesIteratorMap.find(faceVerticesKey) == m_faceVerticesIteratorMap.end()) @@ -69,7 +62,6 @@ void HalfEdgeDataStructure::MakeFace(uint64_t a, uint64_t b, uint64_t c, HyperPl auto const backFaceIterator = m_facesList.emplace(m_facesList.end(), *newHalfEdges.back()); m_faceIteratorMap[&*backFaceIterator] = backFaceIterator; m_faceVerticesIteratorMap[faceVerticesKey] = backFaceIterator; - backFaceIterator->hyperPlane = hp; backFaceIterator->index = index; //Initialize half edges diff --git a/Epona/sources/HyperPlane.cpp b/Epona/sources/HyperPlane.cpp index c5bfc62..086a8a9 100644 --- a/Epona/sources/HyperPlane.cpp +++ b/Epona/sources/HyperPlane.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2017 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/demo/Main.cpp b/demo/Main.cpp index df85a2c..a8fbea0 100644 --- a/demo/Main.cpp +++ b/demo/Main.cpp @@ -1,5 +1,5 @@ /* -* Copyright (C) 2018 by Godlike +* Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */ diff --git a/test/AnalysisTests.cpp b/test/AnalysisTests.cpp index 7b38fd0..af1fe55 100644 --- a/test/AnalysisTests.cpp +++ b/test/AnalysisTests.cpp @@ -1,5 +1,5 @@ /* - * Copyright (C) 2017 by Godlike + * Copyright (C) 2019 by Godlike * This code is licensed under the MIT license (MIT) * (http://opensource.org/licenses/MIT) */