Skip to content

Commit

Permalink
[#27] QHull DOD refactoring
Browse files Browse the repository at this point in the history
  • Loading branch information
@ilia-glushchenko
ilia-glushchenko committed Nov 5, 2018
1 parent 640b830 commit 6ffaa0b
Show file tree
Hide file tree
Showing 7 changed files with 386 additions and 185 deletions.
1 change: 0 additions & 1 deletion Epona/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ set(EPONA_SOURCES
sources/Analysis.cpp
sources/HalfEdgeDataStructure.cpp
sources/HyperPlane.cpp
sources/JacobiEigenvalue.cpp
)

if (MSVC)
Expand Down
71 changes: 59 additions & 12 deletions Epona/include/Epona/HalfEdgeDataStructure.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#ifndef EPONA_HEDS_HPP
#define EPONA_HEDS_HPP

#include <Epona/HyperPlane.hpp>
#include <iterator>
#include <unordered_map>
#include <list>
Expand All @@ -21,10 +22,36 @@ class HalfEdgeDataStructure
struct HalfEdge;
class Face;

using Faces = std::list<Face>;
using HalfEdges = std::list<HalfEdge>;
using face_iterator = Faces::iterator;
using const_face_iterator = Faces::const_iterator;
using face_iterator = std::list<Face>::iterator;
using const_face_iterator = std::list<Face>::const_iterator;

HalfEdgeDataStructure() = default;

HalfEdgeDataStructure(HalfEdgeDataStructure const&) = delete;

HalfEdgeDataStructure& operator=(HalfEdgeDataStructure const& heds) = delete;

HalfEdgeDataStructure(HalfEdgeDataStructure&& heds)
{
m_halfEdgeList = std::move(heds.m_halfEdgeList);
m_halfEdgePointerIteratorMap = std::move(heds.m_halfEdgePointerIteratorMap);
m_halfEdgeVerticesIteratorMap = std::move(heds.m_halfEdgeVerticesIteratorMap);

m_facesList = std::move(heds.m_facesList);
m_faceIteratorMap = std::move(heds.m_faceIteratorMap);
m_faceVerticesIteratorMap = std::move(heds.m_faceVerticesIteratorMap);
}

HalfEdgeDataStructure& operator=(HalfEdgeDataStructure&& heds)
{
m_halfEdgeList = std::move(heds.m_halfEdgeList);
m_halfEdgePointerIteratorMap = std::move(heds.m_halfEdgePointerIteratorMap);
m_halfEdgeVerticesIteratorMap = std::move(heds.m_halfEdgeVerticesIteratorMap);

m_facesList = std::move(heds.m_facesList);
m_faceIteratorMap = std::move(heds.m_faceIteratorMap);
m_faceVerticesIteratorMap = std::move(heds.m_faceVerticesIteratorMap);
}

/**
* @brief HEDS Face data container
Expand Down Expand Up @@ -321,6 +348,9 @@ class HalfEdgeDataStructure
*/
const_face_iterator GetAdjacentFaceIterator() const;

//! Stores face hyperplane
HyperPlane hyperPlane;

private:
HalfEdge* m_halfEdge;
};
Expand Down Expand Up @@ -360,13 +390,30 @@ class HalfEdgeDataStructure
uint64_t vertexIndex;
};

void MakeFace(uint64_t a, uint64_t b, uint64_t c)
{
MakeFace(a, b, c, {});
}

/**
* @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);
void MakeFace(uint64_t a, uint64_t b, uint64_t c, HyperPlane hp);

template < typename T >
decltype(auto) GetFace(T index)
{
return GetFace(index[0], index[1], index[2]);
}

template < typename T >
decltype(auto) GetFace(T index) const
{
return GetFace(index[0], index[1], index[2]);
}

/**
* @brief Returns a face iterator
Expand Down Expand Up @@ -429,20 +476,20 @@ class HalfEdgeDataStructure
bool operator==(HalfEdgeVertices const& other) const;
};

HalfEdges m_halfEdgeList;
std::list<HalfEdge> m_halfEdgeList;
std::unordered_map<
HalfEdge const*, HalfEdges::iterator
HalfEdge const*, std::list<HalfEdge>::iterator
> m_halfEdgePointerIteratorMap;
std::unordered_map<
HalfEdgeVertices, HalfEdges::iterator, HalfEdgeVertices::Hasher
HalfEdgeVertices, std::list<HalfEdge>::iterator, HalfEdgeVertices::Hasher
> m_halfEdgeVerticesIteratorMap;

Faces m_facesList;
std::list<Face> m_facesList;
std::unordered_map<
Face const*, Faces::iterator
Face const*, std::list<Face>::iterator
> m_faceIteratorMap;
std::unordered_map<
FaceVertices, Faces::iterator, FaceVertices::Hasher
FaceVertices, std::list<Face>::iterator, FaceVertices::Hasher
> m_faceVerticesIteratorMap;

/**
Expand All @@ -455,7 +502,7 @@ class HalfEdgeDataStructure
* @param[in] vertexIndexTo end edge vertex index
*/
void IntializeHalfEdge(
HalfEdges::iterator he,
std::list<HalfEdge>::iterator he,
HalfEdge* next, HalfEdge* prev,
Face* face,
uint64_t vertexIndexFrom,
Expand Down
231 changes: 166 additions & 65 deletions Epona/include/Epona/JacobiEigenvalue.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,77 +7,178 @@
#define EPONA_JACOBI_EIGENVALUE_HPP

#include <Epona/FloatingPoint.hpp>
#define GLM_ENABLE_EXPERIMENTAL
#include <glm/gtx/norm.hpp>
#include <glm/glm.hpp>
#include <tuple>

namespace
{
/**
* @brief Finds maximal absolute value off diagonal matrix element and sets its indices
*
* @param[in] mat matrix to search
* @param[out] i max element row index
* @param[out] j max element column index
*/
inline std::tuple<uint8_t, uint8_t> FindMaxNormOffDiagonal(glm::mat3 mat)
{
uint8_t i = 0;
uint8_t j = 1;

for (uint8_t p = 0; p < 3; ++p)
{
for (uint8_t q = 0; q < 3; ++q)
{
if (p != q)
{
if (glm::abs(mat[i][j]) < glm::abs(mat[p][q]))
{
i = p;
j = q;
}
}
}
}

return { i, j };
}

/**
* @brief Calculates rotation angle for a given matrix and its element
*
* @param[in] mat matrix to rotate
* @param[in] i element row index
* @param[in] j element column index
*
* @return angle in radians
*/
inline float CalculateRotationAngle(glm::mat3 mat, uint8_t i, uint8_t j, float coverageThreshold)
{
if (glm::abs(mat[i][i] - mat[j][j]) < coverageThreshold)
{
return (glm::pi<float>() / 4.0f) * (mat[i][j] > 0 ? 1.0f : -1.0f);
}

return 0.5f * glm::atan(2.0f * mat[i][j] / (mat[i][i] - mat[j][j]));
}

/**
* @brief Makes Givens rotation matrix from the angle and indices
*
* @param[in] theta rotation angle in radians
* @param[in] i row index
* @param[in] j column index
*
* @return Givens rotation matrix
*/
inline glm::mat3 MakeGivensRotationMatrix(float theta, uint8_t i, uint8_t j)
{
glm::mat3 g;

g[i][i] = glm::cos(theta);
g[i][j] = glm::sin(theta);
g[j][i] = -glm::sin(theta);
g[j][j] = glm::cos(theta);

return g;
}
} // namespace ::

namespace epona
{

/**
* @brief Calculates eigenvalues and eigenvectors
*
* @param[in] symmetricMatrix target symmetric matrix
* @param[in] coverageThreshold coverage precision threshold
* @param[in] maxIterations maximum amount of iteration
*/
inline std::tuple<glm::mat3, glm::vec3> CalculateJacobiEigenvectorsEigenvalue(
glm::mat3 symmetricMatrix, float coverageThreshold = epona::fp::g_floatingPointThreshold, uint32_t maxIterations = 100
)
{
glm::mat3 eigenvectors;
glm::mat3 rotationMatrix;
uint8_t i;
uint8_t j;
uint16_t iterations = 0;

do
{
std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);

rotationMatrix = ::MakeGivensRotationMatrix(::CalculateRotationAngle(symmetricMatrix, i, j, coverageThreshold), i, j);
eigenvectors = eigenvectors * rotationMatrix;
symmetricMatrix = (glm::transpose(rotationMatrix) * symmetricMatrix) * rotationMatrix;

std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);
}
while ((++iterations < maxIterations) && (glm::abs(symmetricMatrix[i][j]) > coverageThreshold));

return { eigenvectors, { symmetricMatrix[0][0], symmetricMatrix[1][1], symmetricMatrix[2][2] } };
}

/**
* @brief Jacobi eigenvalue calculation algorithm
* @brief Calculates eigenvalues
*
* @param[in] symmetricMatrix target symmetric matrix
* @param[in] coverageThreshold coverage precision threshold
* @param[in] maxIterations maximum amount of iteration
*/
class JacobiEigenvalue
inline glm::vec3 CalculateJacobiEigenvalue(
glm::mat3 symmetricMatrix, float coverageThreshold = epona::fp::g_floatingPointThreshold, uint32_t maxIterations = 100
)
{
public:
/**
* @brief Constructs and calculates eigenvalues and eigenvectors of a given symmetric matrix
* @param[in] symmetricMatrix target symmetric matrix
* @param[in] coverageThreshold coverage precision threshold
* @param[in] maxIterations maximum amount of iteration
*/
explicit JacobiEigenvalue(
glm::mat3 const& symmetricMatrix,
float coverageThreshold = epona::fp::g_floatingPointThreshold,
uint32_t maxIterations = 100
);

/**
* @brief Returns eigenvectors
* @return eigenvectos matrix
*/
glm::mat3 const& GetEigenvectors() const;

/**
* @brief Returns eigenvalues
* @return eigenvalue vector
*/
glm::vec3 const& GetEigenvalues() const;

private:
glm::mat3 const& m_symmetricMatrix;
float const m_coverageThreshold;
uint32_t const m_maxIterations;
glm::mat3 m_eigenvectors;
glm::vec3 m_eigenvalues;

/**
* @brief Finds maximal absolute value off diagonal matrix element and sets its indices
* @param[in] mat matrix to search
* @param[out] i max element row index
* @param[out] j max element column index
*/
static void FindMaxNormOffDiagonal(glm::mat3 const& mat, uint8_t& i, uint8_t& j);

/**
* @brief Calculates rotation angle for a given matrix and its element
* @param[in] mat matrix to rotate
* @param[in] i element row index
* @param[in] j element column index
* @return angle in radians
*/
float CalculateRotationAngle(glm::mat3 const& mat, uint8_t i, uint8_t j) const;

/**
* @brief Makes Givens rotation matrix from the angle and indices
* @param[in] theta rotation angle in radians
* @param[in] i row index
* @param[in] j column index
* @return Givens rotation matrix
*/
static glm::mat3 MakeGivensRotationMatrix(float theta, uint8_t i, uint8_t j);

/**
* @brief Calculates eigenvalues and eigenvectors
*/
void Calculate();
};
glm::mat3 rotationMatrix;
uint8_t i;
uint8_t j;
uint16_t iterations = 0;

do
{
std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);

rotationMatrix = ::MakeGivensRotationMatrix(::CalculateRotationAngle(symmetricMatrix, i, j, coverageThreshold), i, j);
symmetricMatrix = (glm::transpose(rotationMatrix) * symmetricMatrix) * rotationMatrix;

std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);
} while ((++iterations < maxIterations) && (glm::abs(symmetricMatrix[i][j]) > coverageThreshold));

return { symmetricMatrix[0][0], symmetricMatrix[1][1], symmetricMatrix[2][2] };
}

/**
* @brief Calculates eigenvectors
*
* @param[in] symmetricMatrix target symmetric matrix
* @param[in] coverageThreshold coverage precision threshold
* @param[in] maxIterations maximum amount of iteration
*/
inline glm::mat3 CalculateJacobiEigenvectors(
glm::mat3 symmetricMatrix, float coverageThreshold = epona::fp::g_floatingPointThreshold, uint32_t maxIterations = 100
)
{
glm::mat3 eigenvectors;
glm::mat3 rotationMatrix;
uint8_t i;
uint8_t j;
uint16_t iterations = 0;

do
{
std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);

rotationMatrix = ::MakeGivensRotationMatrix(::CalculateRotationAngle(symmetricMatrix, i, j, coverageThreshold), i, j);
eigenvectors = eigenvectors * rotationMatrix;
symmetricMatrix = (glm::transpose(rotationMatrix) * symmetricMatrix) * rotationMatrix;

std::tie(i, j) = ::FindMaxNormOffDiagonal(symmetricMatrix);
} while ((++iterations < maxIterations) && (glm::abs(symmetricMatrix[i][j]) > coverageThreshold));

return eigenvectors;
}

} // namespace epona
#endif // EPONA_JACOBI_EIGENVALUE_HPP
Loading

0 comments on commit 6ffaa0b

Please sign in to comment.