From 996bcb18ee4d23ea14b04fc834ae2114c7f276d4 Mon Sep 17 00:00:00 2001 From: "erwin.coumans" Date: Thu, 7 Jul 2011 00:28:15 +0000 Subject: [PATCH] Added HACD library by Khaled Mamou, for improved convex decomposition See also http://sourceforge.net/projects/hacd/ and http://bulletphysics.org/Bullet/phpBB3/viewtopic.php?f=4&t=6619 Thanks to Sujeong Kim for the integration with the Bullet/Demos/ConvexDecompositionDemo --- Demos/ConvexDecompositionDemo/CMakeLists.txt | 6 +- .../ConvexDecompositionDemo.cpp | 118 +- Extras/CMakeLists.txt | 2 +- Extras/HACD/CMakeLists.txt | 51 + Extras/HACD/hacdCircularList.h | 80 ++ Extras/HACD/hacdCircularList.inl | 163 +++ Extras/HACD/hacdGraph.cpp | 240 ++++ Extras/HACD/hacdGraph.h | 117 ++ Extras/HACD/hacdHACD.cpp | 766 +++++++++++++ Extras/HACD/hacdHACD.h | 276 +++++ Extras/HACD/hacdICHull.cpp | 1010 +++++++++++++++++ Extras/HACD/hacdICHull.h | 120 ++ Extras/HACD/hacdManifoldMesh.cpp | 577 ++++++++++ Extras/HACD/hacdManifoldMesh.h | 250 ++++ Extras/HACD/hacdVector.h | 67 ++ Extras/HACD/hacdVector.inl | 178 +++ Extras/HACD/hacdVersion.h | 20 + 17 files changed, 4032 insertions(+), 9 deletions(-) create mode 100644 Extras/HACD/CMakeLists.txt create mode 100644 Extras/HACD/hacdCircularList.h create mode 100644 Extras/HACD/hacdCircularList.inl create mode 100644 Extras/HACD/hacdGraph.cpp create mode 100644 Extras/HACD/hacdGraph.h create mode 100644 Extras/HACD/hacdHACD.cpp create mode 100644 Extras/HACD/hacdHACD.h create mode 100644 Extras/HACD/hacdICHull.cpp create mode 100644 Extras/HACD/hacdICHull.h create mode 100644 Extras/HACD/hacdManifoldMesh.cpp create mode 100644 Extras/HACD/hacdManifoldMesh.h create mode 100644 Extras/HACD/hacdVector.h create mode 100644 Extras/HACD/hacdVector.inl create mode 100644 Extras/HACD/hacdVersion.h diff --git a/Demos/ConvexDecompositionDemo/CMakeLists.txt b/Demos/ConvexDecompositionDemo/CMakeLists.txt index 03ac8189c..f532f16a0 100644 --- a/Demos/ConvexDecompositionDemo/CMakeLists.txt +++ b/Demos/ConvexDecompositionDemo/CMakeLists.txt @@ -14,6 +14,8 @@ SET(GLUT_ROOT ${BULLET_PHYSICS_SOURCE_DIR}/Glut) INCLUDE_DIRECTORIES( ${BULLET_PHYSICS_SOURCE_DIR}/src ${BULLET_PHYSICS_SOURCE_DIR}/Demos/OpenGL + +${BULLET_PHYSICS_SOURCE_DIR}/Extras/HACD ${BULLET_PHYSICS_SOURCE_DIR}/Extras/ConvexDecomposition ${BULLET_PHYSICS_SOURCE_DIR}/Extras/Serialize/BulletFileLoader ${BULLET_PHYSICS_SOURCE_DIR}/Extras/Serialize/BulletWorldImporter @@ -22,7 +24,7 @@ ${BULLET_PHYSICS_SOURCE_DIR}/Extras/Serialize/BulletWorldImporter IF (USE_GLUT) LINK_LIBRARIES( - OpenGLSupport BulletWorldImporter BulletDynamics BulletCollision LinearMath BulletFileLoader ConvexDecomposition ${GLUT_glut_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} + OpenGLSupport BulletWorldImporter BulletDynamics BulletCollision LinearMath BulletFileLoader HACD ConvexDecomposition ${GLUT_glut_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ) ADD_EXECUTABLE(AppConvexDecompositionDemo @@ -50,7 +52,7 @@ IF (USE_GLUT) ELSE (USE_GLUT) LINK_LIBRARIES( - OpenGLSupport BulletWorldImporter BulletDynamics BulletCollision LinearMath BulletFileLoader ConvexDecomposition ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} + OpenGLSupport BulletWorldImporter BulletDynamics BulletCollision LinearMath BulletFileLoader HACD ConvexDecomposition ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ) ADD_EXECUTABLE(AppConvexDecompositionDemo diff --git a/Demos/ConvexDecompositionDemo/ConvexDecompositionDemo.cpp b/Demos/ConvexDecompositionDemo/ConvexDecompositionDemo.cpp index 5c0b88d53..f62b9459a 100644 --- a/Demos/ConvexDecompositionDemo/ConvexDecompositionDemo.cpp +++ b/Demos/ConvexDecompositionDemo/ConvexDecompositionDemo.cpp @@ -13,6 +13,12 @@ subject to the following restrictions: 3. This notice may not be removed or altered from any source distribution. */ +#include "hacdCircularList.h" +#include "hacdVector.h" +#include "hacdICHull.h" +#include "hacdGraph.h" +#include "hacdHACD.h" + #include "cd_wavefront.h" #include "ConvexBuilder.h" @@ -416,6 +422,9 @@ void ConvexDecompositionDemo::initPhysics(const char* filename) if (tcount) { + //----------------------------------- + // Bullet Convex Decomposition + //----------------------------------- char outputFileName[512]; strcpy(outputFileName,filename); @@ -433,7 +442,7 @@ void ConvexDecompositionDemo::initPhysics(const char* filename) printf("WavefrontObj num triangles read %i\n",tcount); ConvexDecomposition::DecompDesc desc; - desc.mVcount = wo.mVertexCount; + desc.mVcount = wo.mVertexCount; desc.mVertices = wo.mVertices; desc.mTcount = wo.mTriCount; desc.mIndices = (unsigned int *)wo.mIndices; @@ -445,13 +454,64 @@ void ConvexDecompositionDemo::initPhysics(const char* filename) MyConvexDecomposition convexDecomposition(outputFile,this); desc.mCallback = &convexDecomposition; - - + + + //----------------------------------------------- + // HACD + //----------------------------------------------- + + std::vector< HACD::Vec3 > points; + std::vector< HACD::Vec3 > triangles; + + for(int i=0; i vertex(wo.mVertices[index], wo.mVertices[index+1],wo.mVertices[index+2]); + points.push_back(vertex); + } + + for(int i=0;i triangle(wo.mIndices[index], wo.mIndices[index+1], wo.mIndices[index+2]); + triangles.push_back(triangle); + } + + + HACD::HACD myHACD; + myHACD.SetPoints(&points[0]); + myHACD.SetNPoints(points.size()); + myHACD.SetTriangles(&triangles[0]); + myHACD.SetNTriangles(triangles.size()); + myHACD.SetCompacityWeight(0.1); + myHACD.SetVolumeWeight(0.0); + + // HACD parameters + // Recommended parameters: 2 100 0 0 0 0 + size_t nClusters = 2; + double concavity = 100; + bool invert = false; + bool addExtraDistPoints = false; + bool addNeighboursDistPoints = false; + bool addFacesPoints = false; + + myHACD.SetNClusters(nClusters); // minimum number of clusters + myHACD.SetNVerticesPerCH(100); // max of 100 vertices per convex-hull + myHACD.SetConcavity(concavity); // maximum concavity + myHACD.SetAddExtraDistPoints(addExtraDistPoints); + myHACD.SetAddNeighboursDistPoints(addNeighboursDistPoints); + myHACD.SetAddFacesPoints(addFacesPoints); + + myHACD.Compute(); + nClusters = myHACD.GetNClusters(); + + myHACD.Save("output.wrl", false); + //convexDecomposition.performConvexDecomposition(desc); - ConvexBuilder cb(desc.mCallback); - cb.process(desc); +// ConvexBuilder cb(desc.mCallback); +// cb.process(desc); //now create some bodies if (1) @@ -461,7 +521,53 @@ void ConvexDecompositionDemo::initPhysics(const char* filename) btTransform trans; trans.setIdentity(); + + for (int c=0;c * pointsCH = new HACD::Vec3[nPoints]; + HACD::Vec3 * trianglesCH = new HACD::Vec3[nTriangles]; + myHACD.GetCH(c, pointsCH, trianglesCH); + + // points + for(size_t v = 0; v < nPoints; v++) + { + vertices[3*v] = pointsCH[v].X(); + vertices[3*v+1] = pointsCH[v].Y(); + vertices[3*v+2] = pointsCH[v].Z(); + } + // triangles + for(size_t f = 0; f < nTriangles; f++) + { + triangles[3*f] = trianglesCH[f].X(); + triangles[3*f+1] = trianglesCH[f].Y(); + triangles[3*f+2] = trianglesCH[f].Z(); + } + + delete [] pointsCH; + delete [] trianglesCH; + + ConvexResult r(nPoints, vertices, nTriangles, triangles); + convexDecomposition.ConvexDecompResult(r); + } + for (int i=0;iaddChildShape(trans,convexShape); + + btRigidBody* body; + body = localCreateRigidBody( 1.0, trans,convexShape); + } +/* for (int i=0;i +#include +namespace HACD +{ + //! CircularListElement class. + template < typename T > class CircularListElement + { + public: + T & GetData() { return m_data; } + const T & GetData() const { return m_data; } + CircularListElement * & GetNext() { return m_next; } + CircularListElement * & GetPrev() { return m_prev; } + const CircularListElement * & GetNext() const { return m_next; } + const CircularListElement * & GetPrev() const { return m_prev; } + //! Constructor + CircularListElement(const T & data) {m_data = data;} + CircularListElement(void){} + //! Destructor + ~CircularListElement(void){} + private: + T m_data; + CircularListElement * m_next; + CircularListElement * m_prev; + + CircularListElement(const CircularListElement & rhs); + }; + + + //! CircularList class. + template < typename T > class CircularList + { + public: + CircularListElement * & GetHead() { return m_head;} + const CircularListElement * GetHead() const { return m_head;} + bool IsEmpty() const { return (m_size == 0);} + size_t GetSize() const { return m_size; } + const T & GetData() const { return m_head->GetData(); } + T & GetData() { return m_head->GetData();} + bool Delete() ; + bool Delete(CircularListElement * element); + CircularListElement * Add(const T * data = 0); + CircularListElement * Add(const T & data); + bool Next(); + bool Prev(); + void Clear() { while(Delete());}; + const CircularList& operator=(const CircularList& rhs); + //! Constructor + CircularList() + { + m_head = 0; + m_size = 0; + } + CircularList(const CircularList& rhs); + //! Destructor + virtual ~CircularList(void) {Clear();}; + private: + CircularListElement * m_head; //!< a pointer to the head of the circular list + size_t m_size; //!< number of element in the circular list + + }; +} +#include +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdCircularList.inl b/Extras/HACD/hacdCircularList.inl new file mode 100644 index 000000000..108f288cc --- /dev/null +++ b/Extras/HACD/hacdCircularList.inl @@ -0,0 +1,163 @@ +#pragma once +#ifndef HACD_CIRCULAR_LIST_INL +#define HACD_CIRCULAR_LIST_INL +#include +#include +namespace HACD +{ + template < typename T > + inline bool CircularList::Delete(CircularListElement * element) + { + if (!element) + { + return false; + } + if (m_size > 1) + { + CircularListElement * next = element->GetNext(); + CircularListElement * prev = element->GetPrev(); + delete element; + m_size--; + if (element == m_head) + { + m_head = next; + } + next->GetPrev() = prev; + prev->GetNext() = next; + return true; + } + else if (m_size == 1) + { + delete m_head; + m_size--; + m_head = 0; + return true; + } + else + { + return false; + } + } + + template < typename T > + inline bool CircularList::Delete() + { + if (m_size > 1) + { + CircularListElement * next = m_head->GetNext(); + CircularListElement * prev = m_head->GetPrev(); + delete m_head; + m_size--; + m_head = next; + next->GetPrev() = prev; + prev->GetNext() = next; + return true; + } + else if (m_size == 1) + { + delete m_head; + m_size--; + m_head = 0; + return true; + } + else + { + return false; + } + } + template < typename T > + inline CircularListElement * CircularList::Add(const T * data) + { + if (m_size == 0) + { + if (data) + { + m_head = new CircularListElement(*data); + } + else + { + m_head = new CircularListElement(); + } + m_head->GetNext() = m_head->GetPrev() = m_head; + } + else + { + CircularListElement * next = m_head->GetNext(); + CircularListElement * element = m_head; + if (data) + { + m_head = new CircularListElement(*data); + } + else + { + m_head = new CircularListElement; + } + m_head->GetNext() = next; + m_head->GetPrev() = element; + element->GetNext() = m_head; + next->GetPrev() = m_head; + } + m_size++; + return m_head; + } + template < typename T > + inline CircularListElement * CircularList::Add(const T & data) + { + const T * pData = &data; + return Add(pData); + } + template < typename T > + inline bool CircularList::Next() + { + if (m_size == 0) + { + return false; + } + m_head = m_head->GetNext(); + return true; + } + template < typename T > + inline bool CircularList::Prev() + { + if (m_size == 0) + { + return false; + } + m_head = m_head->GetPrev(); + return true; + } + template < typename T > + inline CircularList::CircularList(const CircularList& rhs) + { + if (rhs.m_size > 0) + { + CircularListElement * current = rhs.m_head; + do + { + current = current->GetNext(); + Add(current->GetData()); + } + while ( current != rhs.m_head ); + } + } + template < typename T > + inline const CircularList& CircularList::operator=(const CircularList& rhs) + { + if (&rhs != this) + { + Clear(); + if (rhs.m_size > 0) + { + CircularListElement * current = rhs.m_head; + do + { + current = current->GetNext(); + Add(current->GetData()); + } + while ( current != rhs.m_head ); + } + } + return (*this); + } +} +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdGraph.cpp b/Extras/HACD/hacdGraph.cpp new file mode 100644 index 000000000..b69270e18 --- /dev/null +++ b/Extras/HACD/hacdGraph.cpp @@ -0,0 +1,240 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#include +namespace HACD +{ + + GraphEdge::GraphEdge() + { + m_convexHull = 0; + m_v1 = -1; + m_v2 = -1; + m_name = -1; + m_error = 0; + m_surf = 0; + m_perimeter = 0; + m_concavity = 0; + m_volume = 0; + m_deleted = false; + } + + GraphVertex::GraphVertex() + { + m_convexHull = 0; + m_name = -1; + m_error = 0; + m_surf = 0; + m_perimeter = 0; + m_concavity = 0; + m_volume = 0; + m_deleted = false; + } + + bool GraphVertex::DeleteEdge(long name) + { + std::set::iterator it = m_edges.find(name); + if (it != m_edges.end() ) + { + m_edges.erase(it); + return true; + } + return false; + } + + Graph::Graph() + { + m_nV = 0; + m_nE = 0; + } + + Graph::~Graph() + { + } + + void Graph::Allocate(size_t nV, size_t nE) + { + m_nV = nV; + m_edges.reserve(nE); + m_vertices.resize(nV); + for(size_t i = 0; i < nV; i++) + { + m_vertices[i].m_name = static_cast(i); + } + } + + long Graph::AddVertex() + { + size_t name = m_vertices.size(); + m_vertices.resize(name+1); + m_vertices[name].m_name = name; + m_nV++; + return static_cast(name); + } + + long Graph::AddEdge(long v1, long v2) + { + size_t name = m_edges.size(); + m_edges.push_back(GraphEdge()); + m_edges[name].m_name = name; + m_edges[name].m_v1 = v1; + m_edges[name].m_v2 = v2; + m_vertices[v1].AddEdge(name); + m_vertices[v2].AddEdge(name); + m_nE++; + return static_cast(name); + } + + bool Graph::DeleteEdge(long name) + { + if (name < static_cast(m_edges.size())) + { + long v1 = m_edges[name].m_v1; + long v2 = m_edges[name].m_v2; + m_edges[name].m_deleted = true; + m_vertices[v1].DeleteEdge(name); + m_vertices[v2].DeleteEdge(name); + delete m_edges[name].m_convexHull; + m_edges[name].m_distPoints.clear(); + m_edges[name].m_boudaryEdges.clear(); + m_edges[name].m_convexHull = 0; + m_nE--; + return true; + } + return false; + } + bool Graph::DeleteVertex(long name) + { + if (name < static_cast(m_vertices.size())) + { + m_vertices[name].m_deleted = true; + m_vertices[name].m_edges.clear(); + m_vertices[name].m_ancestors = std::vector(); + delete m_vertices[name].m_convexHull; + m_vertices[name].m_distPoints.clear(); + m_vertices[name].m_boudaryEdges.clear(); + m_vertices[name].m_convexHull = 0; + m_nV--; + return true; + } + return false; + } + bool Graph::EdgeCollapse(long v1, long v2) + { + long edgeToDelete = GetEdgeID(v1, v2); + if (edgeToDelete >= 0) + { + // delete the edge (v1, v2) + DeleteEdge(edgeToDelete); + // add v2 to v1 ancestors + m_vertices[v1].m_ancestors.push_back(v2); + // add v2's ancestors to v1's ancestors + m_vertices[v1].m_ancestors.insert(m_vertices[v1].m_ancestors.begin(), + m_vertices[v2].m_ancestors.begin(), + m_vertices[v2].m_ancestors.end()); + // update adjacency information + std::set & v1Edges = m_vertices[v1].m_edges; + std::set::const_iterator ed(m_vertices[v2].m_edges.begin()); + std::set::const_iterator itEnd(m_vertices[v2].m_edges.end()); + long b = -1; + for(; ed != itEnd; ++ed) + { + if (m_edges[*ed].m_v1 == v2) + { + b = m_edges[*ed].m_v2; + } + else + { + b = m_edges[*ed].m_v1; + } + if (GetEdgeID(v1, b) >= 0) + { + m_edges[*ed].m_deleted = true; + m_vertices[b].DeleteEdge(*ed); + m_nE--; + } + else + { + m_edges[*ed].m_v1 = v1; + m_edges[*ed].m_v2 = b; + v1Edges.insert(*ed); + } + } + // delete the vertex v2 + DeleteVertex(v2); + return true; + } + return false; + } + + long Graph::GetEdgeID(long v1, long v2) const + { + if (v1 < static_cast(m_vertices.size()) && !m_vertices[v1].m_deleted) + { + std::set::const_iterator ed(m_vertices[v1].m_edges.begin()); + std::set::const_iterator itEnd(m_vertices[v1].m_edges.end()); + for(; ed != itEnd; ++ed) + { + if ( (m_edges[*ed].m_v1 == v2) || + (m_edges[*ed].m_v2 == v2) ) + { + return m_edges[*ed].m_name; + } + } + } + return -1; + } + + void Graph::Print() const + { + std::cout << "-----------------------------" << std::endl; + std::cout << "vertices (" << m_nV << ")" << std::endl; + for (size_t v = 0; v < m_vertices.size(); ++v) + { + const GraphVertex & currentVertex = m_vertices[v]; + if (!m_vertices[v].m_deleted) + { + + std::cout << currentVertex.m_name << "\t"; + std::set::const_iterator ed(currentVertex.m_edges.begin()); + std::set::const_iterator itEnd(currentVertex.m_edges.end()); + for(; ed != itEnd; ++ed) + { + std::cout << "(" << m_edges[*ed].m_v1 << "," << m_edges[*ed].m_v2 << ") "; + } + std::cout << std::endl; + } + } + + std::cout << "vertices (" << m_nE << ")" << std::endl; + for (size_t e = 0; e < m_edges.size(); ++e) + { + const GraphEdge & currentEdge = m_edges[e]; + if (!m_edges[e].m_deleted) + { + std::cout << currentEdge.m_name << "\t(" + << m_edges[e].m_v1 << "," + << m_edges[e].m_v2 << ") "<< std::endl; + } + } + } + void Graph::Clear() + { + m_vertices.clear(); + m_edges.clear(); + m_nV = 0; + m_nE = 0; + } + +} \ No newline at end of file diff --git a/Extras/HACD/hacdGraph.h b/Extras/HACD/hacdGraph.h new file mode 100644 index 000000000..56421624d --- /dev/null +++ b/Extras/HACD/hacdGraph.h @@ -0,0 +1,117 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_GRAPH_H +#define HACD_GRAPH_H +#include +#include +#include +#include +#include +#include + +namespace HACD +{ + class GraphVertex; + class GraphEdge; + class Graph; + class HACD; + + class GraphVertex + { + public: + bool AddEdge(long name) + { + m_edges.insert(name); + return true; + } + bool DeleteEdge(long name); + GraphVertex(); + ~GraphVertex(){ delete m_convexHull;}; + private: + long m_name; + std::set m_edges; + bool m_deleted; + std::vector m_ancestors; + std::map m_distPoints; + + Real m_error; + double m_surf; + double m_volume; + double m_perimeter; + double m_concavity; + ICHull * m_convexHull; + std::set m_boudaryEdges; + + + friend class GraphEdge; + friend class Graph; + friend class HACD; + }; + + class GraphEdge + { + public: + GraphEdge(); + ~GraphEdge(){delete m_convexHull;}; + private: + long m_name; + long m_v1; + long m_v2; + std::map m_distPoints; + Real m_error; + double m_surf; + double m_volume; + double m_perimeter; + double m_concavity; + ICHull * m_convexHull; + std::set m_boudaryEdges; + bool m_deleted; + + + + friend class GraphVertex; + friend class Graph; + friend class HACD; + }; + + class Graph + { + public: + size_t GetNEdges() const { return m_nE;} + size_t GetNVertices() const { return m_nV;} + bool EdgeCollapse(long v1, long v2); + long AddVertex(); + long AddEdge(long v1, long v2); + bool DeleteEdge(long name); + bool DeleteVertex(long name); + long GetEdgeID(long v1, long v2) const; + void Clear(); + void Print() const; + + Graph(); + virtual ~Graph(); + void Allocate(size_t nV, size_t nE); + + private: + size_t m_nV; + size_t m_nE; + std::vector m_edges; + std::vector m_vertices; + + friend class HACD; + }; +} +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdHACD.cpp b/Extras/HACD/hacdHACD.cpp new file mode 100644 index 000000000..bfe0f9902 --- /dev/null +++ b/Extras/HACD/hacdHACD.cpp @@ -0,0 +1,766 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#define _CRT_SECURE_NO_WARNINGS +#include +#include +#include +#include +#include +#include +namespace HACD +{ + double HACD::Concavity(ICHull & ch, std::map & distPoints) + { + double concavity = 0.0; + double distance = 0.0; + std::map::iterator itDP(distPoints.begin()); + std::map::iterator itDPEnd(distPoints.end()); + for(; itDP != itDPEnd; ++itDP) + { + if (!(itDP->second).m_computed) + { + if (itDP->first >= 0) + { + distance = ch.ComputeDistance(itDP->first, m_points[itDP->first], m_normals[itDP->first], (itDP->second).m_computed, true); + } + else + { + distance = ch.ComputeDistance(itDP->first, m_facePoints[-itDP->first-1], m_faceNormals[-itDP->first-1], (itDP->second).m_computed, true); + } + } + else + { + distance = (itDP->second).m_dist; + } + if (concavity < distance) + { + concavity = distance; + } + } + return concavity; + } + + void HACD::CreateGraph() + { + // vertex to triangle adjacency information + std::vector< std::set > vertexToTriangles; + vertexToTriangles.resize(m_nPoints); + for(size_t t = 0; t < m_nTriangles; ++t) + { + vertexToTriangles[m_triangles[t].X()].insert(t); + vertexToTriangles[m_triangles[t].Y()].insert(t); + vertexToTriangles[m_triangles[t].Z()].insert(t); + } + + m_graph.Clear(); + m_graph.Allocate(m_nTriangles, 3 * m_nTriangles); + unsigned long long tr1[3]; + unsigned long long tr2[3]; + long i1, j1, k1, i2, j2, k2; + long t1, t2; + for (size_t v = 0; v < m_nPoints; v++) + { + std::set::const_iterator it1(vertexToTriangles[v].begin()), itEnd(vertexToTriangles[v].end()); + for(; it1 != itEnd; ++it1) + { + t1 = *it1; + i1 = m_triangles[t1].X(); + j1 = m_triangles[t1].Y(); + k1 = m_triangles[t1].Z(); + tr1[0] = GetEdgeIndex(i1, j1); + tr1[1] = GetEdgeIndex(j1, k1); + tr1[2] = GetEdgeIndex(k1, i1); + std::set::const_iterator it2(it1); + for(++it2; it2 != itEnd; ++it2) + { + t2 = *it2; + i2 = m_triangles[t2].X(); + j2 = m_triangles[t2].Y(); + k2 = m_triangles[t2].Z(); + tr2[0] = GetEdgeIndex(i2, j2); + tr2[1] = GetEdgeIndex(j2, k2); + tr2[2] = GetEdgeIndex(k2, i2); + int shared = 0; + for(int i = 0; i < 3; ++i) + { + for(int j = 0; j < 3; ++j) + { + if (tr1[i] == tr2[j]) + { + shared++; + } + } + } + if (shared == 1) // two triangles are connected if they share exactly one edge + { + m_graph.AddEdge(t1, t2); + } + } + } + } + } + void HACD::InitializeDualGraph() + { + long i, j, k; + Vec3 u, v, w, normal; + delete [] m_normals; + m_normals = new Vec3[m_nPoints]; + if (m_addFacesPoints) + { + delete [] m_facePoints; + delete [] m_faceNormals; + m_facePoints = new Vec3[m_nTriangles]; + m_faceNormals = new Vec3[m_nTriangles]; + } + memset(m_normals, 0, sizeof(Vec3) * m_nPoints); + for(unsigned long f = 0; f < m_nTriangles; f++) + { + i = m_triangles[f].X(); + j = m_triangles[f].Y(); + k = m_triangles[f].Z(); + + m_graph.m_vertices[f].m_distPoints[i].m_distOnly = false; + m_graph.m_vertices[f].m_distPoints[j].m_distOnly = false; + m_graph.m_vertices[f].m_distPoints[k].m_distOnly = false; + + ICHull * ch = new ICHull; + m_graph.m_vertices[f].m_convexHull = ch; + ch->AddPoint(m_points[i], i); + ch->AddPoint(m_points[j], j); + ch->AddPoint(m_points[k], k); + ch->SetDistPoints(&m_graph.m_vertices[f].m_distPoints); + + u = m_points[j] - m_points[i]; + v = m_points[k] - m_points[i]; + w = m_points[k] - m_points[j]; + normal = u ^ v; + + m_normals[i] += normal; + m_normals[j] += normal; + m_normals[k] += normal; + + m_graph.m_vertices[f].m_surf = normal.GetNorm(); + m_graph.m_vertices[f].m_perimeter = u.GetNorm() + v.GetNorm() + w.GetNorm(); + + normal.Normalize(); + + m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(i,j)); + m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(j,k)); + m_graph.m_vertices[f].m_boudaryEdges.insert(GetEdgeIndex(k,i)); + if(m_addFacesPoints) + { + m_faceNormals[f] = normal; + m_facePoints[f] = (m_points[i] + m_points[j] + m_points[k]) / 3.0; + m_graph.m_vertices[f].m_distPoints[-static_cast(f)-1].m_distOnly = true; + } + if (m_addExtraDistPoints) + {// we need a kd-tree structure to accelerate this part! + long i1, j1, k1; + Vec3 u1, v1, normal1; + normal = -normal; + double distance = 0.0; + double distMin = 0.0; + size_t faceIndex = m_nTriangles; + Vec3 seedPoint((m_points[i] + m_points[j] + m_points[k]) / 3.0); + long nhit = 0; + for(size_t f1 = 0; f1 < m_nTriangles; f1++) + { + i1 = m_triangles[f1].X(); + j1 = m_triangles[f1].Y(); + k1 = m_triangles[f1].Z(); + u1 = m_points[j1] - m_points[i1]; + v1 = m_points[k1] - m_points[i1]; + normal1 = (u1 ^ v1); + if (normal * normal1 > 0.0) + { + nhit = IntersectRayTriangle(Vec3(seedPoint.X(), seedPoint.Y(), seedPoint.Z()), + Vec3(normal.X(), normal.Y(), normal.Z()), + Vec3(m_points[i1].X(), m_points[i1].Y(), m_points[i1].Z()), + Vec3(m_points[j1].X(), m_points[j1].Y(), m_points[j1].Z()), + Vec3(m_points[k1].X(), m_points[k1].Y(), m_points[k1].Z()), + distance); + if ((nhit==1) && ((distMin > distance) || (faceIndex == m_nTriangles))) + { + distMin = distance; + faceIndex = f1; + } + + } + } + if (faceIndex < m_nTriangles ) + { + i1 = m_triangles[faceIndex].X(); + j1 = m_triangles[faceIndex].Y(); + k1 = m_triangles[faceIndex].Z(); + m_graph.m_vertices[f].m_distPoints[i1].m_distOnly = true; + m_graph.m_vertices[f].m_distPoints[j1].m_distOnly = true; + m_graph.m_vertices[f].m_distPoints[k1].m_distOnly = true; + if (m_addFacesPoints) + { + m_graph.m_vertices[f].m_distPoints[-static_cast(faceIndex)-1].m_distOnly = true; + } + } + } + } + for (size_t v = 0; v < m_nPoints; v++) + { + m_normals[v].Normalize(); + } + } + + void HACD::NormalizeData() + { + if (m_nPoints == 0) + { + return; + } + m_barycenter = m_points[0]; + Vec3 min = m_points[0]; + Vec3 max = m_points[0]; + Real x, y, z; + for (size_t v = 1; v < m_nPoints ; v++) + { + m_barycenter += m_points[v]; + x = m_points[v].X(); + y = m_points[v].Y(); + z = m_points[v].Z(); + if ( x < min.X()) min.X() = x; + else if ( x > max.X()) max.X() = x; + if ( y < min.Y()) min.Y() = y; + else if ( y > max.Y()) max.Y() = y; + if ( z < min.Z()) min.Z() = z; + else if ( z > max.Z()) max.Z() = z; + } + m_barycenter /= static_cast(m_nPoints); + m_diag = (max-min).GetNorm(); + const Real invDiag = static_cast(2.0 * m_scale / m_diag); + if (m_diag != 0.0) + { + for (size_t v = 0; v < m_nPoints ; v++) + { + m_points[v] = (m_points[v] - m_barycenter) * invDiag; + } + } + } + void HACD::DenormalizeData() + { + if (m_nPoints == 0) + { + return; + } + if (m_diag != 0.0) + { + const Real diag = static_cast(m_diag / (2.0 * m_scale)); + for (size_t v = 0; v < m_nPoints ; v++) + { + m_points[v] = m_points[v] * diag + m_barycenter; + } + } + } + HACD::HACD(void) + { + m_convexHulls = 0; + m_triangles = 0; + m_points = 0; + m_normals = 0; + m_nTriangles = 0; + m_nPoints = 0; + m_nClusters = 0; + m_concavity = 0.0; + m_diag = 1.0; + m_barycenter = Vec3(0.0, 0.0,0.0); + m_alpha = 0.1; + m_beta = 0.1; + m_nVerticesPerCH = 30; + m_callBack = 0; + m_addExtraDistPoints = false; + m_addNeighboursDistPoints = false; + m_scale = 1000.0; + m_partition = 0; + m_nMinClusters = 3; + m_facePoints = 0; + m_faceNormals = 0; + } + HACD::~HACD(void) + { + delete [] m_normals; + delete [] m_convexHulls; + delete [] m_partition; + delete [] m_facePoints; + delete [] m_faceNormals; + } + int iteration = 0; + void HACD::ComputeEdgeCost(size_t e) + { + GraphEdge & gE = m_graph.m_edges[e]; + long v1 = gE.m_v1; + long v2 = gE.m_v2; + + if (m_graph.m_vertices[v2].m_distPoints.size()>m_graph.m_vertices[v1].m_distPoints.size()) + { + gE.m_v1 = v2; + gE.m_v2 = v1; + std::swap(v1, v2); + } + GraphVertex & gV1 = m_graph.m_vertices[v1]; + GraphVertex & gV2 = m_graph.m_vertices[v2]; + + // delete old convex-hull + delete gE.m_convexHull; + // create the edge's convex-hull + ICHull * ch = new ICHull; + gE.m_convexHull = ch; + (*ch) = (*gV1.m_convexHull); + + // update distPoints + gE.m_distPoints = gV1.m_distPoints; + std::map::iterator itDP(gV2.m_distPoints.begin()); + std::map::iterator itDPEnd(gV2.m_distPoints.end()); + std::map::iterator itDP1; + + for(; itDP != itDPEnd; ++itDP) + { + itDP1 = gE.m_distPoints.find(itDP->first); + if (itDP1 == gE.m_distPoints.end()) + { + gE.m_distPoints[itDP->first].m_distOnly = (itDP->second).m_distOnly; + if ( !(itDP->second).m_distOnly ) + { + ch->AddPoint(m_points[itDP->first], itDP->first); + } + } + else + { + if ( (itDP1->second).m_distOnly && !(itDP->second).m_distOnly) + { + gE.m_distPoints[itDP->first].m_distOnly = false; + ch->AddPoint(m_points[itDP->first], itDP->first); + } + } + } + + ch->SetDistPoints(&gE.m_distPoints); + // create the convex-hull + while (ch->Process() == ICHullErrorInconsistent) // if we face problems when constructing the visual-hull. really ugly!!!! + { +// if (m_callBack) (*m_callBack)("\t Problem with convex-hull construction [HACD::ComputeEdgeCost]\n", 0.0, 0.0, 0); + ch = new ICHull; + CircularList & verticesCH = (gE.m_convexHull)->GetMesh().m_vertices; + size_t nV = verticesCH.GetSize(); + long ptIndex = 0; + verticesCH.Next(); + for(size_t v = 1; v < nV; ++v) + { + ptIndex = verticesCH.GetHead()->GetData().m_name; + ch->AddPoint(m_points[ptIndex], ptIndex); + verticesCH.Next(); + } + delete gE.m_convexHull; + gE.m_convexHull = ch; + } + double volume = 0.0; + double concavity = 0.0; + if (ch->IsFlat()) + { + bool insideHull; + std::map::iterator itDP(gE.m_distPoints.begin()); + std::map::iterator itDPEnd(gE.m_distPoints.end()); + for(; itDP != itDPEnd; ++itDP) + { + if (itDP->first >= 0) + { + concavity = std::max(concavity, ch->ComputeDistance(itDP->first, m_points[itDP->first], m_normals[itDP->first], insideHull, false)); + } + } + } + else + { + if (m_addNeighboursDistPoints) + { // add distance points from adjacent clusters + std::set eEdges; + std::set_union(gV1.m_edges.begin(), + gV1.m_edges.end(), + gV2.m_edges.begin(), + gV2.m_edges.end(), + std::inserter( eEdges, eEdges.begin() ) ); + + std::set::const_iterator ed(eEdges.begin()); + std::set::const_iterator itEnd(eEdges.end()); + long a, b, c; + for(; ed != itEnd; ++ed) + { + a = m_graph.m_edges[*ed].m_v1; + b = m_graph.m_edges[*ed].m_v2; + if ( a != v2 && a != v1) + { + c = a; + } + else if ( b != v2 && b != v1) + { + c = b; + } + else + { + c = -1; + } + if ( c > 0) + { + GraphVertex & gVC = m_graph.m_vertices[c]; + std::map::iterator itDP(gVC.m_distPoints.begin()); + std::map::iterator itDPEnd(gVC.m_distPoints.end()); + std::map::iterator itDP1; + for(; itDP != itDPEnd; ++itDP) + { + itDP1 = gE.m_distPoints.find(itDP->first); + if (itDP1 == gE.m_distPoints.end()) + { + if (itDP->first >= 0 && itDP1 == gE.m_distPoints.end() && ch->IsInside(m_points[itDP->first])) + { + gE.m_distPoints[itDP->first].m_distOnly = true; + } + else if (itDP->first < 0 && ch->IsInside(m_facePoints[-itDP->first-1])) + { + gE.m_distPoints[itDP->first].m_distOnly = true; + } + } + } + } + } + } + concavity = Concavity(*ch, gE.m_distPoints); + } + + // compute boudary edges + double perimeter = 0.0; + double surf = 1.0; + if (m_alpha > 0.0) + { + gE.m_boudaryEdges.clear(); + std::set_symmetric_difference (gV1.m_boudaryEdges.begin(), + gV1.m_boudaryEdges.end(), + gV2.m_boudaryEdges.begin(), + gV2.m_boudaryEdges.end(), + std::inserter( gE.m_boudaryEdges, + gE.m_boudaryEdges.begin() ) ); + + std::set::const_iterator itBE(gE.m_boudaryEdges.begin()); + std::set::const_iterator itBEEnd(gE.m_boudaryEdges.end()); + for(; itBE != itBEEnd; ++itBE) + { + perimeter += (m_points[static_cast((*itBE) >> 32)] - + m_points[static_cast((*itBE) & 0xFFFFFFFFULL)]).GetNorm(); + } + surf = gV1.m_surf + gV2.m_surf; + } + double ratio = perimeter * perimeter / (4.0 * sc_pi * surf); + gE.m_volume = (m_beta == 0.0)?0.0:ch->ComputeVolume()/pow(m_scale, 3.0); // cluster's volume + gE.m_surf = surf; // cluster's area + gE.m_perimeter = perimeter; // cluster's perimeter + gE.m_concavity = concavity; // cluster's concavity + gE.m_error = static_cast(concavity + m_alpha * ratio + m_beta * volume); // cluster's priority + } + bool HACD::InitializePriorityQueue() + { + m_pqueue.reserve(m_graph.m_nE + 100); + for (size_t e=0; e < m_graph.m_nE; ++e) + { + ComputeEdgeCost(e); + m_pqueue.push(GraphEdgePriorityQueue(e, m_graph.m_edges[e].m_error)); + } + return true; + } + void HACD::Simplify() + { + long v1 = -1; + long v2 = -1; + double progressOld = -1.0; + double progress = 0.0; + double globalConcavity = 0.0; + char msg[1024]; + double ptgStep = 1.0; + while ( (globalConcavity < m_concavity) && + (m_graph.GetNVertices() > m_nMinClusters) && + (m_graph.GetNEdges()> 1)) + { + progress = 100.0-m_graph.GetNVertices() * 100.0 / m_nTriangles; + if (fabs(progress-progressOld) > ptgStep && m_callBack) + { + sprintf(msg, "%3.2f %% V = %lu \t C = %f \t \t \r", progress, static_cast(m_graph.GetNVertices()), globalConcavity); + (*m_callBack)(msg, progress, globalConcavity, m_graph.GetNVertices()); + progressOld = progress; + if (progress > 99.0) + { + ptgStep = 0.01; + } + else if (progress > 90.0) + { + ptgStep = 0.1; + } + } + + GraphEdgePriorityQueue currentEdge(0,0.0); + bool done = false; + do + { + done = false; + if (m_pqueue.size() == 0) + { + done = true; + break; + } + else + { + currentEdge = m_pqueue.top(); + m_pqueue.pop(); + } + } + while ( m_graph.m_edges[currentEdge.m_name].m_deleted || + m_graph.m_edges[currentEdge.m_name].m_error != currentEdge.m_priority); + + if (m_graph.m_edges[currentEdge.m_name].m_concavity < m_concavity && !done) + { + globalConcavity = std::max(globalConcavity ,m_graph.m_edges[currentEdge.m_name].m_concavity); + v1 = m_graph.m_edges[currentEdge.m_name].m_v1; + v2 = m_graph.m_edges[currentEdge.m_name].m_v2; + // update vertex info + m_graph.m_vertices[v1].m_error = m_graph.m_edges[currentEdge.m_name].m_error; + m_graph.m_vertices[v1].m_surf = m_graph.m_edges[currentEdge.m_name].m_surf; + m_graph.m_vertices[v1].m_volume = m_graph.m_edges[currentEdge.m_name].m_volume; + m_graph.m_vertices[v1].m_concavity = m_graph.m_edges[currentEdge.m_name].m_concavity; + m_graph.m_vertices[v1].m_perimeter = m_graph.m_edges[currentEdge.m_name].m_perimeter; + m_graph.m_vertices[v1].m_distPoints = m_graph.m_edges[currentEdge.m_name].m_distPoints; + (*m_graph.m_vertices[v1].m_convexHull) = (*m_graph.m_edges[currentEdge.m_name].m_convexHull); + (m_graph.m_vertices[v1].m_convexHull)->SetDistPoints(&(m_graph.m_vertices[v1].m_distPoints)); + m_graph.m_vertices[v1].m_boudaryEdges = m_graph.m_edges[currentEdge.m_name].m_boudaryEdges; + + // We apply the optimal ecol +// std::cout << "v1 " << v1 << " v2 " << v2 << std::endl; + m_graph.EdgeCollapse(v1, v2); + // recompute the adjacent edges costs + std::set::const_iterator itE(m_graph.m_vertices[v1].m_edges.begin()), + itEEnd(m_graph.m_vertices[v1].m_edges.end()); + for(; itE != itEEnd; ++itE) + { + size_t e = *itE; + ComputeEdgeCost(e); + m_pqueue.push(GraphEdgePriorityQueue(e, m_graph.m_edges[e].m_error)); + } + } + else + { + break; + } + } + while (!m_pqueue.empty()) + { + m_pqueue.pop(); + } + + m_cVertices.clear(); + m_cVertices.reserve(m_nClusters); + m_nClusters = m_graph.GetNVertices(); + for (size_t p=0, v = 0; v != m_graph.m_vertices.size(); ++v) + { + if (!m_graph.m_vertices[v].m_deleted) + { + if (m_callBack) + { + char msg[1024]; + sprintf(msg, "\t CH \t %lu \t %lf \t %lf\n", static_cast(p), m_graph.m_vertices[v].m_concavity, m_graph.m_vertices[v].m_error); + (*m_callBack)(msg, 0.0, 0.0, m_nClusters); + p++; + } + m_cVertices.push_back(v); + } + } + if (m_callBack) + { + sprintf(msg, "# clusters = %lu \t C = %f\n", static_cast(m_nClusters), globalConcavity); + (*m_callBack)(msg, progress, globalConcavity, m_graph.GetNVertices()); + } + + } + bool HACD::Compute(bool fullCH, bool exportDistPoints) + { + if ( !m_points || !m_triangles || !m_nPoints || !m_nTriangles) + { + return false; + } + size_t nV = m_nTriangles; + if (m_callBack) + { + std::ostringstream msg; + msg << "+ Mesh" << std::endl; + msg << "\t # vertices \t" << m_nPoints << std::endl; + msg << "\t # triangles \t" << m_nTriangles << std::endl; + msg << "+ Parameters" << std::endl; + msg << "\t min # of clusters \t" << m_nMinClusters << std::endl; + msg << "\t max concavity \t" << m_concavity << std::endl; + msg << "\t compacity weigth \t" << m_alpha << std::endl; + msg << "\t volume weigth \t" << m_beta << std::endl; + msg << "\t # vertices per convex-hull \t" << m_nVerticesPerCH << std::endl; + msg << "\t Scale \t" << m_scale << std::endl; + msg << "\t Add extra distance points \t" << m_addExtraDistPoints << std::endl; + msg << "\t Add neighbours distance points \t" << m_addNeighboursDistPoints << std::endl; + msg << "\t Add face distance points \t" << m_addFacesPoints << std::endl; + msg << "\t Produce full convex-hulls \t" << fullCH << std::endl; + (*m_callBack)(msg.str().c_str(), 0.0, 0.0, nV); + } + if (m_callBack) (*m_callBack)("+ Normalizing Data\n", 0.0, 0.0, nV); + NormalizeData(); + if (m_callBack) (*m_callBack)("+ Creating Graph\n", 0.0, 0.0, nV); + CreateGraph(); + // Compute the surfaces and perimeters of all the faces + if (m_callBack) (*m_callBack)("+ Initializing Dual Graph\n", 0.0, 0.0, nV); + InitializeDualGraph(); + if (m_callBack) (*m_callBack)("+ Initializing Priority Queue\n", 0.0, 0.0, nV); + InitializePriorityQueue(); + // we simplify the graph + if (m_callBack) (*m_callBack)("+ Simplification ...\n", 0.0, 0.0, m_nTriangles); + Simplify(); + if (m_callBack) (*m_callBack)("+ Denormalizing Data\n", 0.0, 0.0, m_nClusters); + DenormalizeData(); + if (m_callBack) (*m_callBack)("+ Computing final convex-hulls\n", 0.0, 0.0, m_nClusters); + delete [] m_convexHulls; + m_convexHulls = new ICHull[m_nClusters]; + delete [] m_partition; + m_partition = new long [m_nTriangles]; + for (size_t p = 0; p != m_cVertices.size(); ++p) + { + size_t v = m_cVertices[p]; + m_partition[v] = p; + for(size_t a = 0; a < m_graph.m_vertices[v].m_ancestors.size(); a++) + { + m_partition[m_graph.m_vertices[v].m_ancestors[a]] = p; + } + // compute the convex-hull + const std::map & pointsCH = m_graph.m_vertices[v].m_distPoints; + std::map::const_iterator itCH(pointsCH.begin()); + std::map::const_iterator itCHEnd(pointsCH.end()); + for(; itCH != itCHEnd; ++itCH) + { + if (!(itCH->second).m_distOnly) + { + m_convexHulls[p].AddPoint(m_points[itCH->first], itCH->first); + } + } + m_convexHulls[p].SetDistPoints(&m_graph.m_vertices[v].m_distPoints); + if (fullCH) + { + m_convexHulls[p].Process(); + } + else + { + m_convexHulls[p].Process(m_nVerticesPerCH); + } + if (exportDistPoints) + { + itCH = pointsCH.begin(); + for(; itCH != itCHEnd; ++itCH) + { + if ((itCH->second).m_distOnly) + { + if (itCH->first >= 0) + { + m_convexHulls[p].AddPoint(m_points[itCH->first], itCH->first); + } + else + { + m_convexHulls[p].AddPoint(m_facePoints[-itCH->first-1], itCH->first); + } + } + } + } + } + return true; + } + + size_t HACD::GetNTrianglesCH(size_t numCH) const + { + if (numCH >= m_nClusters) + { + return 0; + } + return m_convexHulls[numCH].GetMesh().GetNTriangles(); + } + size_t HACD::GetNPointsCH(size_t numCH) const + { + if (numCH >= m_nClusters) + { + return 0; + } + return m_convexHulls[numCH].GetMesh().GetNVertices(); + } + + bool HACD::GetCH(size_t numCH, Vec3 * const points, Vec3 * const triangles) + { + if (numCH >= m_nClusters) + { + return false; + } + m_convexHulls[numCH].GetMesh().GetIFS(points, triangles); + return true; + } + + bool HACD::Save(const char * fileName, bool uniColor, long numCluster) const + { + std::ofstream fout(fileName); + if (fout.is_open()) + { + if (m_callBack) + { + char msg[1024]; + sprintf(msg, "Saving %s\n", fileName); + (*m_callBack)(msg, 0.0, 0.0, m_graph.GetNVertices()); + } + Material mat; + if (numCluster < 0) + { + for (size_t p = 0; p != m_nClusters; ++p) + { + if (!uniColor) + { + mat.m_diffuseColor.X() = mat.m_diffuseColor.Y() = mat.m_diffuseColor.Z() = 0.0; + while (mat.m_diffuseColor.X() == mat.m_diffuseColor.Y() || + mat.m_diffuseColor.Z() == mat.m_diffuseColor.Y() || + mat.m_diffuseColor.Z() == mat.m_diffuseColor.X() ) + { + mat.m_diffuseColor.X() = (rand()%100) / 100.0; + mat.m_diffuseColor.Y() = (rand()%100) / 100.0; + mat.m_diffuseColor.Z() = (rand()%100) / 100.0; + } + } + m_convexHulls[p].GetMesh().SaveVRML2(fout, mat); + } + } + else if (numCluster < static_cast(m_cVertices.size())) + { + m_convexHulls[numCluster].GetMesh().SaveVRML2(fout, mat); + } + fout.close(); + return true; + } + else + { + if (m_callBack) + { + char msg[1024]; + sprintf(msg, "Error saving %s\n", fileName); + (*m_callBack)(msg, 0.0, 0.0, m_graph.GetNVertices()); + } + return false; + } + } +} \ No newline at end of file diff --git a/Extras/HACD/hacdHACD.h b/Extras/HACD/hacdHACD.h new file mode 100644 index 000000000..446719b9a --- /dev/null +++ b/Extras/HACD/hacdHACD.h @@ -0,0 +1,276 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_HACD_H +#define HACD_HACD_H +#include +#include +#include +#include +#include +#include +#include + +namespace HACD +{ + const double sc_pi = 3.14159265; + class HACD; + + // just to be able to set the capcity of the container + + template, class _Pr = std::less > + class reservable_priority_queue: public std::priority_queue<_Ty, _Container, _Pr> + { + typedef typename std::priority_queue<_Ty, _Container, _Pr>::size_type size_type; + public: + reservable_priority_queue(size_type capacity = 0) { reserve(capacity); }; + void reserve(size_type capacity) { this->c.reserve(capacity); } + size_type capacity() const { return this->c.capacity(); } + }; + + //! priority queque element + class GraphEdgePriorityQueue + { + public: + //! Constructor + //! @param name edge's id + //! @param priority edge's priority + GraphEdgePriorityQueue(long name, Real priority) + { + m_name = name; + m_priority = priority; + } + //! Destructor + ~GraphEdgePriorityQueue(void){} + private: + long m_name; //!< edge name + Real m_priority; //!< priority + //! Operator < for GraphEdgePQ + friend bool operator<(const GraphEdgePriorityQueue & lhs, const GraphEdgePriorityQueue & rhs); + //! Operator > for GraphEdgePQ + friend bool operator>(const GraphEdgePriorityQueue & lhs, const GraphEdgePriorityQueue & rhs); + friend class HACD; + }; + inline bool operator<(const GraphEdgePriorityQueue & lhs, const GraphEdgePriorityQueue & rhs) + { + if (lhs.m_priority!=rhs.m_priority) return lhs.m_priority(const GraphEdgePriorityQueue & lhs, const GraphEdgePriorityQueue & rhs) + { + if (lhs.m_priority!=rhs.m_priority) return lhs.m_priority>rhs.m_priority; + return lhs.m_name > rhs.m_name; + } + typedef void (*CallBackFunction)(const char *, double, double, size_t); + + //! Provides an implementation of the Hierarchical Approximate Convex Decomposition (HACD) technique described in "A Simple and Efficient Approach for 3D Mesh Approximate Convex Decomposition" Game Programming Gems 8 - Chapter 2.8, p.202. A short version of the chapter was published in ICIP09 and is available at ftp://ftp.elet.polimi.it/users/Stefano.Tubaro/ICIP_USB_Proceedings_v2/pdfs/0003501.pdf + class HACD + { + public: + + //! Gives the triangles partitionas an array of size m_nTriangles where the i-th element specifies the cluster to which belong the i-th triangle + //! @return triangles partition + const long * const GetPartition() const { return m_partition;} + //! Sets the scale factor + //! @param scale scale factor + void SetScaleFactor(double scale) { m_scale = scale;} + //! Gives the scale factor + //! @return scale factor + const double GetScaleFactor() const { return m_scale;} + //! Sets the call-back function + //! @param callBack pointer to the call-back function + void SetCallBack(CallBackFunction callBack) { m_callBack = callBack;} + //! Gives the call-back function + //! @return pointer to the call-back function + const CallBackFunction GetCallBack() const { return m_callBack;} + + //! Specifies whether faces points should be added when computing the concavity + //! @param addFacesPoints true = faces points should be added + void SetAddFacesPoints(bool addFacesPoints) { m_addFacesPoints = addFacesPoints;} + //! Specifies wheter faces points should be added when computing the concavity + //! @return true = faces points should be added + const bool GetAddFacesPoints() const { return m_addFacesPoints;} + //! Specifies whether extra points should be added when computing the concavity + //! @param addExteraDistPoints true = extra points should be added + void SetAddExtraDistPoints(bool addExtraDistPoints) { m_addExtraDistPoints = addExtraDistPoints;} + //! Specifies wheter extra points should be added when computing the concavity + //! @return true = extra points should be added + const bool GetAddExtraDistPoints() const { return m_addExtraDistPoints;} + //! Specifies whether extra points should be added when computing the concavity + //! @param addExteraDistPoints true = extra points should be added + void SetAddNeighboursDistPoints(bool addNeighboursDistPoints) { m_addNeighboursDistPoints = addNeighboursDistPoints;} + //! Specifies wheter extra points should be added when computing the concavity + //! @return true = extra points should be added + const bool GetAddNeighboursDistPoints() const { return m_addExtraDistPoints;} + //! Sets the points of the input mesh (Remark: the input points will be scaled and shifted. Use DenormalizeData() to invert those operations) + //! @param points pointer to the input points + void SetPoints(Vec3 * points) { m_points = points;} + //! Gives the points of the input mesh (Remark: the input points will be scaled and shifted. Use DenormalizeData() to invert those operations) + //! @return pointer to the input points + const Vec3 * GetPoints() const { return m_points;} + //! Sets the triangles of the input mesh. + //! @param triangles points pointer to the input points + void SetTriangles(Vec3 * triangles) { m_triangles = triangles;} + //! Gives the triangles in the input mesh + //! @return pointer to the input triangles + const Vec3 * GetTriangles() const { return m_triangles;} + //! Sets the number of points in the input mesh. + //! @param nPoints number of points the input mesh + void SetNPoints(size_t nPoints) { m_nPoints = nPoints;} + //! Gives the number of points in the input mesh. + //! @return number of points the input mesh + const size_t GetNPoints() const { return m_nPoints;} + //! Sets the number of triangles in the input mesh. + //! @param nTriangles number of triangles in the input mesh + void SetNTriangles(size_t nTriangles) { m_nTriangles = nTriangles;} + //! Gives the number of triangles in the input mesh. + //! @return number of triangles the input mesh + const size_t GetNTriangles() const { return m_nTriangles;} + //! Sets the minimum number of clusters to be generated. + //! @param nClusters minimum number of clusters + void SetNClusters(size_t nClusters) { m_nMinClusters = nClusters;} + //! Gives the number of generated clusters. + //! @return number of generated clusters + const size_t GetNClusters() const { return m_nClusters;} + //! Sets the maximum allowed concavity. + //! @param concavity maximum concavity + void SetConcavity(double concavity) { m_concavity = concavity;} + //! Gives the maximum allowed concavity. + //! @return maximum concavity + double GetConcavity() const { return m_concavity;} + //! Sets the volume weight. + //! @param beta volume weight + void SetVolumeWeight(double beta) { m_beta = beta;} + //! Gives the volume weight. + //! @return volume weight + double GetVolumeWeight() const { return m_beta;} + //! Sets the compacity weight (i.e. parameter alpha in ftp://ftp.elet.polimi.it/users/Stefano.Tubaro/ICIP_USB_Proceedings_v2/pdfs/0003501.pdf). + //! @param alpha compacity weight + void SetCompacityWeight(double alpha) { m_alpha = alpha;} + //! Gives the compacity weight (i.e. parameter alpha in ftp://ftp.elet.polimi.it/users/Stefano.Tubaro/ICIP_USB_Proceedings_v2/pdfs/0003501.pdf). + //! @return compacity weight + double GetCompacityWeight() const { return m_alpha;} + //! Sets the maximum number of vertices for each generated convex-hull. + //! @param nVerticesPerCH maximum # vertices per CH + void SetNVerticesPerCH(size_t nVerticesPerCH) { m_nVerticesPerCH = nVerticesPerCH;} + //! Gives the maximum number of vertices for each generated convex-hull. + //! @return maximum # vertices per CH + const size_t GetNVerticesPerCH() const { return m_nVerticesPerCH;} + //! Gives the number of vertices for the cluster number numCH. + //! @return number of vertices + size_t GetNPointsCH(size_t numCH) const; + //! Gives the number of triangles for the cluster number numCH. + //! @param numCH cluster's number + //! @return number of triangles + size_t GetNTrianglesCH(size_t numCH) const; + //! Gives the vertices and the triangles of the cluster number numCH. + //! @param numCH cluster's number + //! @param points pointer to the vector of points to be filled + //! @param triangles pointer to the vector of triangles to be filled + //! @return true if sucess + bool GetCH(size_t numCH, Vec3 * const points, Vec3 * const triangles); + //! Computes the HACD decomposition. + //! @param fullCH specifies whether to generate convex-hulls with a full or limited (i.e. < m_nVerticesPerCH) number of vertices + //! @param exportDistPoints specifies wheter distance points should ne exported or not (used only for debugging). + //! @return true if sucess + bool Compute(bool fullCH=false, bool exportDistPoints=false); + //! Saves the generated convex-hulls in a VRML 2.0 file. + //! @param fileName the output file name + //! @param uniColor specifies whether the different convex-hulls should have the same color or not + //! @param numCluster specifies the cluster to be saved, if numCluster < 0 export all clusters + //! @return true if sucess + bool Save(const char * fileName, bool uniColor, long numCluster=-1) const; + //! Shifts and scales to the data to have all the coordinates between 0.0 and 1000.0. + void NormalizeData(); + //! Inverse the operations applied by NormalizeData(). + void DenormalizeData(); + //! Constructor. + HACD(void); + //! Destructor. + ~HACD(void); + + private: + //! Gives the edge index. + //! @param a first vertex id + //! @param b second vertex id + //! @return edge's index + static unsigned long long GetEdgeIndex(unsigned long long a, unsigned long long b) + { + if (a > b) return (a << 32) + b; + else return (b << 32) + a; + } + //! Computes the concavity of a cluster. + //! @param ch the cluster's convex-hull + //! @param distPoints the cluster's points + //! @return cluster's concavity + double Concavity(ICHull & ch, std::map & distPoints); + //! Computes the perimeter of a cluster. + //! @param triIndices the cluster's triangles + //! @param distPoints the cluster's points + //! @return cluster's perimeter + double ComputePerimeter(const std::vector & triIndices) const; + //! Creates the Graph by associating to each mesh triangle a vertex in the graph and to each couple of adjacent triangles an edge in the graph. + void CreateGraph(); + //! Initializes the graph costs and computes the vertices normals + void InitializeDualGraph(); + //! Computes the cost of an edge + //! @param e edge's id + void ComputeEdgeCost(size_t e); + //! Initializes the priority queue + //! @param fast specifies whether fast mode is used + //! @return true if success + bool InitializePriorityQueue(); + //! Cleans the intersection between convex-hulls + void CleanClusters(); + //! Computes convex-hulls from partition information + //! @param fullCH specifies whether to generate convex-hulls with a full or limited (i.e. < m_nVerticesPerCH) number of vertices + void ComputeConvexHulls(bool fullCH); + //! Simplifies the graph + //! @param fast specifies whether fast mode is used + void Simplify(); + + private: + double m_scale; //>! scale factor used for NormalizeData() and DenormalizeData() + Vec3 * m_triangles; //>! pointer the triangles array + Vec3 * m_points; //>! pointer the points array + Vec3 * m_facePoints; //>! pointer to the faces points array + Vec3 * m_faceNormals; //>! pointer to the faces normals array + Vec3 * m_normals; //>! pointer the normals array + size_t m_nTriangles; //>! number of triangles in the original mesh + size_t m_nPoints; //>! number of vertices in the original mesh + size_t m_nClusters; //>! number of clusters + size_t m_nMinClusters; //>! minimum number of clusters + double m_concavity; //>! maximum concavity + double m_alpha; //>! compacity weigth + double m_beta; //>! volume weigth + double m_diag; //>! length of the BB diagonal + Vec3 m_barycenter; //>! barycenter of the mesh + std::vector< long > m_cVertices; //>! array of vertices each belonging to a different cluster + ICHull * m_convexHulls; //>! convex-hulls associated with the final HACD clusters + Graph m_graph; //>! simplification graph + size_t m_nVerticesPerCH; //>! maximum number of vertices per convex-hull + reservable_priority_queue, + std::greater::value_type> > m_pqueue; //!> priority queue + HACD(const HACD & rhs); + CallBackFunction m_callBack; //>! call-back function + long * m_partition; //>! array of size m_nTriangles where the i-th element specifies the cluster to which belong the i-th triangle + bool m_addFacesPoints; //>! specifies whether to add faces points or not + bool m_addExtraDistPoints; //>! specifies whether to add extra points for concave shapes or not + bool m_addNeighboursDistPoints; //>! specifies whether to add extra points from adjacent clusters or not + + }; +} +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdICHull.cpp b/Extras/HACD/hacdICHull.cpp new file mode 100644 index 000000000..1df733673 --- /dev/null +++ b/Extras/HACD/hacdICHull.cpp @@ -0,0 +1,1010 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#include +#include +namespace HACD +{ + const long ICHull::sc_dummyIndex = static_cast(-1); + ICHull::ICHull(void) + { + m_distPoints = 0; + m_isFlat = false; + m_dummyVertex = 0; + } + bool ICHull::AddPoints(const Vec3 * points, size_t nPoints) + { + if (!points) + { + return false; + } + CircularListElement * vertex = NULL; + for (size_t i = 0; i < nPoints; i++) + { + vertex = m_mesh.AddVertex(); + vertex->GetData().m_pos.X() = points[i].X(); + vertex->GetData().m_pos.Y() = points[i].Y(); + vertex->GetData().m_pos.Z() = points[i].Z(); + vertex->GetData().m_name = i; + } + return true; + } + bool ICHull::AddPoints(std::vector< Vec3 > points) + { + CircularListElement * vertex = NULL; + for (size_t i = 0; i < points.size(); i++) + { + vertex = m_mesh.AddVertex(); + vertex->GetData().m_pos.X() = points[i].X(); + vertex->GetData().m_pos.Y() = points[i].Y(); + vertex->GetData().m_pos.Z() = points[i].Z(); + } + return true; + } + + bool ICHull::AddPoint(const Vec3 & point, long id) + { + if (AddPoints(&point, 1)) + { + m_mesh.m_vertices.GetData().m_name = id; + return true; + } + return false; + } + + ICHullError ICHull::Process() + { + unsigned long addedPoints = 0; + if (m_mesh.GetNVertices() < 3) + { + return ICHullErrorNotEnoughPoints; + } + if (m_mesh.GetNVertices() == 3) + { + m_isFlat = true; + CircularListElement * t1 = m_mesh.AddTriangle(); + CircularListElement * t2 = m_mesh.AddTriangle(); + CircularListElement * v0 = m_mesh.m_vertices.GetHead(); + CircularListElement * v1 = v0->GetNext(); + CircularListElement * v2 = v1->GetNext(); + // Compute the normal to the plane + Vec3 p0 = v0->GetData().m_pos; + Vec3 p1 = v1->GetData().m_pos; + Vec3 p2 = v2->GetData().m_pos; + m_normal = (p1-p0) ^ (p2-p0); + m_normal.Normalize(); + t1->GetData().m_vertices[0] = v0; + t1->GetData().m_vertices[1] = v1; + t1->GetData().m_vertices[2] = v2; + t2->GetData().m_vertices[0] = v1; + t2->GetData().m_vertices[1] = v2; + t2->GetData().m_vertices[2] = v2; + return ICHullErrorOK; + } + if (m_isFlat) + { + m_mesh.m_edges.Clear(); + m_mesh.m_triangles.Clear(); + m_isFlat = false; + } + if (m_mesh.GetNTriangles() == 0) // we have to create the first polyhedron + { + ICHullError res = DoubleTriangle(); + if (res != ICHullErrorOK) + { + return res; + } + else + { + addedPoints += 3; + } + } + CircularList & vertices = m_mesh.GetVertices(); + // go to the first added and not processed vertex + while (!(vertices.GetHead()->GetPrev()->GetData().m_tag)) + { + vertices.Prev(); + } + while (!vertices.GetData().m_tag) // not processed + { + vertices.GetData().m_tag = true; + if (ProcessPoint()) + { + addedPoints++; + CleanUp(addedPoints); + vertices.Next(); + if (!GetMesh().CheckConsistancy()) + { + return ICHullErrorInconsistent; + } + } + } + if (m_isFlat) + { + std::vector< CircularListElement * > trianglesToDuplicate; + size_t nT = m_mesh.GetNTriangles(); + for(size_t f = 0; f < nT; f++) + { + TMMTriangle & currentTriangle = m_mesh.m_triangles.GetHead()->GetData(); + if( currentTriangle.m_vertices[0]->GetData().m_name == sc_dummyIndex || + currentTriangle.m_vertices[1]->GetData().m_name == sc_dummyIndex || + currentTriangle.m_vertices[2]->GetData().m_name == sc_dummyIndex ) + { + m_trianglesToDelete.push_back(m_mesh.m_triangles.GetHead()); + for(int k = 0; k < 3; k++) + { + for(int h = 0; h < 2; h++) + { + if (currentTriangle.m_edges[k]->GetData().m_triangles[h] == m_mesh.m_triangles.GetHead()) + { + currentTriangle.m_edges[k]->GetData().m_triangles[h] = 0; + break; + } + } + } + } + else + { + trianglesToDuplicate.push_back(m_mesh.m_triangles.GetHead()); + } + m_mesh.m_triangles.Next(); + } + size_t nE = m_mesh.GetNEdges(); + for(size_t e = 0; e < nE; e++) + { + TMMEdge & currentEdge = m_mesh.m_edges.GetHead()->GetData(); + if( currentEdge.m_triangles[0] == 0 && currentEdge.m_triangles[1] == 0) + { + m_edgesToDelete.push_back(m_mesh.m_edges.GetHead()); + } + m_mesh.m_edges.Next(); + } + m_mesh.m_vertices.Delete(m_dummyVertex); + m_dummyVertex = 0; + size_t nV = m_mesh.GetNVertices(); + CircularList & vertices = m_mesh.GetVertices(); + for(size_t v = 0; v < nV; ++v) + { + vertices.GetData().m_tag = false; + vertices.Next(); + } + CleanEdges(); + CleanTriangles(); + CircularListElement * newTriangle; + for(size_t t = 0; t < trianglesToDuplicate.size(); t++) + { + newTriangle = m_mesh.AddTriangle(); + newTriangle->GetData().m_vertices[0] = trianglesToDuplicate[t]->GetData().m_vertices[1]; + newTriangle->GetData().m_vertices[1] = trianglesToDuplicate[t]->GetData().m_vertices[0]; + newTriangle->GetData().m_vertices[2] = trianglesToDuplicate[t]->GetData().m_vertices[2]; + } + } + return ICHullErrorOK; + } + ICHullError ICHull::Process(unsigned long nPointsCH) + { + unsigned long addedPoints = 0; + if (nPointsCH < 3 || m_mesh.GetNVertices() < 3) + { + return ICHullErrorNotEnoughPoints; + } + if (m_mesh.GetNVertices() == 3) + { + m_isFlat = true; + CircularListElement * t1 = m_mesh.AddTriangle(); + CircularListElement * t2 = m_mesh.AddTriangle(); + CircularListElement * v0 = m_mesh.m_vertices.GetHead(); + CircularListElement * v1 = v0->GetNext(); + CircularListElement * v2 = v1->GetNext(); + // Compute the normal to the plane + Vec3 p0 = v0->GetData().m_pos; + Vec3 p1 = v1->GetData().m_pos; + Vec3 p2 = v2->GetData().m_pos; + m_normal = (p1-p0) ^ (p2-p0); + m_normal.Normalize(); + t1->GetData().m_vertices[0] = v0; + t1->GetData().m_vertices[1] = v1; + t1->GetData().m_vertices[2] = v2; + t2->GetData().m_vertices[0] = v1; + t2->GetData().m_vertices[1] = v0; + t2->GetData().m_vertices[2] = v2; + return ICHullErrorOK; + } + + if (m_isFlat) + { + m_mesh.m_triangles.Clear(); + m_mesh.m_edges.Clear(); + m_isFlat = false; + } + + if (m_mesh.GetNTriangles() == 0) // we have to create the first polyhedron + { + ICHullError res = DoubleTriangle(); + if (res != ICHullErrorOK) + { + return res; + } + else + { + addedPoints += 3; + } + } + CircularList & vertices = m_mesh.GetVertices(); + while (!vertices.GetData().m_tag && addedPoints < nPointsCH) // not processed + { + if (!FindMaxVolumePoint()) + { + break; + } + vertices.GetData().m_tag = true; + if (ProcessPoint()) + { + addedPoints++; + CleanUp(addedPoints); + if (!GetMesh().CheckConsistancy()) + { + return ICHullErrorInconsistent; + } + vertices.Next(); + } + } + // delete remaining points + while (!vertices.GetData().m_tag) + { + vertices.Delete(); + } + if (m_isFlat) + { + std::vector< CircularListElement * > trianglesToDuplicate; + size_t nT = m_mesh.GetNTriangles(); + for(size_t f = 0; f < nT; f++) + { + TMMTriangle & currentTriangle = m_mesh.m_triangles.GetHead()->GetData(); + if( currentTriangle.m_vertices[0]->GetData().m_name == sc_dummyIndex || + currentTriangle.m_vertices[1]->GetData().m_name == sc_dummyIndex || + currentTriangle.m_vertices[2]->GetData().m_name == sc_dummyIndex ) + { + m_trianglesToDelete.push_back(m_mesh.m_triangles.GetHead()); + for(int k = 0; k < 3; k++) + { + for(int h = 0; h < 2; h++) + { + if (currentTriangle.m_edges[k]->GetData().m_triangles[h] == m_mesh.m_triangles.GetHead()) + { + currentTriangle.m_edges[k]->GetData().m_triangles[h] = 0; + break; + } + } + } + } + else + { + trianglesToDuplicate.push_back(m_mesh.m_triangles.GetHead()); + } + m_mesh.m_triangles.Next(); + } + size_t nE = m_mesh.GetNEdges(); + for(size_t e = 0; e < nE; e++) + { + TMMEdge & currentEdge = m_mesh.m_edges.GetHead()->GetData(); + if( currentEdge.m_triangles[0] == 0 && currentEdge.m_triangles[1] == 0) + { + m_edgesToDelete.push_back(m_mesh.m_edges.GetHead()); + } + m_mesh.m_edges.Next(); + } + m_mesh.m_vertices.Delete(m_dummyVertex); + m_dummyVertex = 0; + size_t nV = m_mesh.GetNVertices(); + CircularList & vertices = m_mesh.GetVertices(); + for(size_t v = 0; v < nV; ++v) + { + vertices.GetData().m_tag = false; + vertices.Next(); + } + CleanEdges(); + CleanTriangles(); + CircularListElement * newTriangle; + for(size_t t = 0; t < trianglesToDuplicate.size(); t++) + { + newTriangle = m_mesh.AddTriangle(); + newTriangle->GetData().m_vertices[0] = trianglesToDuplicate[t]->GetData().m_vertices[1]; + newTriangle->GetData().m_vertices[1] = trianglesToDuplicate[t]->GetData().m_vertices[0]; + newTriangle->GetData().m_vertices[2] = trianglesToDuplicate[t]->GetData().m_vertices[2]; + } + } + return ICHullErrorOK; + } + bool ICHull::FindMaxVolumePoint() + { + CircularList & vertices = m_mesh.GetVertices(); + CircularListElement * vMaxVolume = 0; + CircularListElement * vHeadPrev = vertices.GetHead()->GetPrev(); + + double maxVolume = 0.0; + double volume = 0.0; + + while (!vertices.GetData().m_tag) // not processed + { + if (ComputePointVolume(volume, false)) + { + if ( maxVolume < volume) + { + maxVolume = volume; + vMaxVolume = vertices.GetHead(); + } + vertices.Next(); + } + } + CircularListElement * vHead = vHeadPrev->GetNext(); + vertices.GetHead() = vHead; + + if (!vMaxVolume) + { + return false; + } + + if (vMaxVolume != vHead) + { + Vec3 pos = vHead->GetData().m_pos; + long id = vHead->GetData().m_name; + vHead->GetData().m_pos = vMaxVolume->GetData().m_pos; + vHead->GetData().m_name = vMaxVolume->GetData().m_name; + vMaxVolume->GetData().m_pos = pos; + vHead->GetData().m_name = id; + } + + + return true; + } + ICHullError ICHull::DoubleTriangle() + { + // find three non colinear points + m_isFlat = false; + CircularList & vertices = m_mesh.GetVertices(); + CircularListElement * v0 = vertices.GetHead(); + while( Colinear(v0->GetData().m_pos, + v0->GetNext()->GetData().m_pos, + v0->GetNext()->GetNext()->GetData().m_pos)) + { + if ( (v0 = v0->GetNext()) == vertices.GetHead()) + { + return ICHullErrorCoplanarPoints; + } + } + CircularListElement * v1 = v0->GetNext(); + CircularListElement * v2 = v1->GetNext(); + // mark points as processed + v0->GetData().m_tag = v1->GetData().m_tag = v2->GetData().m_tag = true; + + // create two triangles + CircularListElement * f0 = MakeFace(v0, v1, v2, 0); + MakeFace(v2, v1, v0, f0); + + // find a fourth non-coplanar point to form tetrahedron + CircularListElement * v3 = v2->GetNext(); + vertices.GetHead() = v3; + + double vol = Volume(v0->GetData().m_pos, v1->GetData().m_pos, v2->GetData().m_pos, v3->GetData().m_pos); + while (vol == 0.0 && !v3->GetNext()->GetData().m_tag) + { + v3 = v3->GetNext(); + vol = Volume(v0->GetData().m_pos, v1->GetData().m_pos, v2->GetData().m_pos, v3->GetData().m_pos); + } + if (vol == 0.0) + { + // compute the barycenter + Vec3 bary(0.0,0.0,0.0); + CircularListElement * vBary = v0; + do + { + bary += vBary->GetData().m_pos; + } + while ( (vBary = vBary->GetNext()) != v0); + bary /= static_cast(vertices.GetSize()); + + // Compute the normal to the plane + Vec3 p0 = v0->GetData().m_pos; + Vec3 p1 = v1->GetData().m_pos; + Vec3 p2 = v2->GetData().m_pos; + m_normal = (p1-p0) ^ (p2-p0); + m_normal.Normalize(); + // add dummy vertex placed at (bary + normal) + vertices.GetHead() = v2; + Vec3 newPt = bary + m_normal; + AddPoint(newPt, sc_dummyIndex); + m_dummyVertex = vertices.GetHead(); + m_isFlat = true; + v3 = v2->GetNext(); + vol = Volume(v0->GetData().m_pos, v1->GetData().m_pos, v2->GetData().m_pos, v3->GetData().m_pos); + return ICHullErrorOK; + } + else if (v3 != vertices.GetHead()) + { + TMMVertex temp; + temp.m_name = v3->GetData().m_name; + temp.m_pos = v3->GetData().m_pos; + v3->GetData().m_name = vertices.GetHead()->GetData().m_name; + v3->GetData().m_pos = vertices.GetHead()->GetData().m_pos; + vertices.GetHead()->GetData().m_name = temp.m_name; + vertices.GetHead()->GetData().m_pos = temp.m_pos; + } + return ICHullErrorOK; + } + CircularListElement * ICHull::MakeFace(CircularListElement * v0, + CircularListElement * v1, + CircularListElement * v2, + CircularListElement * fold) + { + CircularListElement * e0; + CircularListElement * e1; + CircularListElement * e2; + long index = 0; + if (!fold) // if first face to be created + { + e0 = m_mesh.AddEdge(); // create the three edges + e1 = m_mesh.AddEdge(); + e2 = m_mesh.AddEdge(); + } + else // otherwise re-use existing edges (in reverse order) + { + e0 = fold->GetData().m_edges[2]; + e1 = fold->GetData().m_edges[1]; + e2 = fold->GetData().m_edges[0]; + index = 1; + } + e0->GetData().m_vertices[0] = v0; e0->GetData().m_vertices[1] = v1; + e1->GetData().m_vertices[0] = v1; e1->GetData().m_vertices[1] = v2; + e2->GetData().m_vertices[0] = v2; e2->GetData().m_vertices[1] = v0; + // create the new face + CircularListElement * f = m_mesh.AddTriangle(); + f->GetData().m_edges[0] = e0; f->GetData().m_edges[1] = e1; f->GetData().m_edges[2] = e2; + f->GetData().m_vertices[0] = v0; f->GetData().m_vertices[1] = v1; f->GetData().m_vertices[2] = v2; + // link edges to face f + e0->GetData().m_triangles[index] = e1->GetData().m_triangles[index] = e2->GetData().m_triangles[index] = f; + return f; + } + CircularListElement * ICHull::MakeConeFace(CircularListElement * e, CircularListElement * p) + { + // create two new edges if they don't already exist + CircularListElement * newEdges[2]; + for(int i = 0; i < 2; ++i) + { + if ( !( newEdges[i] = e->GetData().m_vertices[i]->GetData().m_duplicate ) ) + { // if the edge doesn't exits add it and mark the vertex as duplicated + newEdges[i] = m_mesh.AddEdge(); + newEdges[i]->GetData().m_vertices[0] = e->GetData().m_vertices[i]; + newEdges[i]->GetData().m_vertices[1] = p; + e->GetData().m_vertices[i]->GetData().m_duplicate = newEdges[i]; + } + } + // make the new face + CircularListElement * newFace = m_mesh.AddTriangle(); + newFace->GetData().m_edges[0] = e; + newFace->GetData().m_edges[1] = newEdges[0]; + newFace->GetData().m_edges[2] = newEdges[1]; + MakeCCW(newFace, e, p); + for(int i=0; i < 2; ++i) + { + for(int j=0; j < 2; ++j) + { + if ( ! newEdges[i]->GetData().m_triangles[j] ) + { + newEdges[i]->GetData().m_triangles[j] = newFace; + break; + } + } + } + return newFace; + } + bool ICHull::ComputePointVolume(double &totalVolume, bool markVisibleFaces) + { + // mark visible faces + CircularListElement * fHead = m_mesh.GetTriangles().GetHead(); + CircularListElement * f = fHead; + CircularList & vertices = m_mesh.GetVertices(); + CircularListElement * vertex0 = vertices.GetHead(); + bool visible = false; + Vec3 pos0 = Vec3(vertex0->GetData().m_pos.X(), + vertex0->GetData().m_pos.Y(), + vertex0->GetData().m_pos.Z()); + double vol = 0.0; + totalVolume = 0.0; + Vec3 ver0, ver1, ver2; + do + { + ver0.X() = f->GetData().m_vertices[0]->GetData().m_pos.X(); + ver0.Y() = f->GetData().m_vertices[0]->GetData().m_pos.Y(); + ver0.Z() = f->GetData().m_vertices[0]->GetData().m_pos.Z(); + ver1.X() = f->GetData().m_vertices[1]->GetData().m_pos.X(); + ver1.Y() = f->GetData().m_vertices[1]->GetData().m_pos.Y(); + ver1.Z() = f->GetData().m_vertices[1]->GetData().m_pos.Z(); + ver2.X() = f->GetData().m_vertices[2]->GetData().m_pos.X(); + ver2.Y() = f->GetData().m_vertices[2]->GetData().m_pos.Y(); + ver2.Z() = f->GetData().m_vertices[2]->GetData().m_pos.Z(); + vol = Volume(ver0, ver1, ver2, pos0); + if ( vol < 0.0 ) + { + vol = fabs(vol); + totalVolume += vol; + if (markVisibleFaces) + { + f->GetData().m_visible = true; + m_trianglesToDelete.push_back(f); + } + visible = true; + } + f = f->GetNext(); + } + while (f != fHead); + + if (m_trianglesToDelete.size() == m_mesh.m_triangles.GetSize()) + { + for(size_t i = 0; i < m_trianglesToDelete.size(); i++) + { + m_trianglesToDelete[i]->GetData().m_visible = false; + } + visible = false; + } + // if no faces visible from p then p is inside the hull + if (!visible && markVisibleFaces) + { + vertices.Delete(); + m_trianglesToDelete.clear(); + return false; + } + return true; + } + bool ICHull::ProcessPoint() + { + double totalVolume = 0.0; + if (!ComputePointVolume(totalVolume, true)) + { + return false; + } + // Mark edges in interior of visible region for deletion. + // Create a new face based on each border edge + CircularListElement * v0 = m_mesh.GetVertices().GetHead(); + CircularListElement * eHead = m_mesh.GetEdges().GetHead(); + CircularListElement * e = eHead; + CircularListElement * tmp = 0; + long nvisible = 0; + m_edgesToDelete.clear(); + m_edgesToUpdate.clear(); + do + { + tmp = e->GetNext(); + nvisible = 0; + for(int k = 0; k < 2; k++) + { + if ( e->GetData().m_triangles[k]->GetData().m_visible ) + { + nvisible++; + } + } + if ( nvisible == 2) + { + m_edgesToDelete.push_back(e); + } + else if ( nvisible == 1) + { + e->GetData().m_newFace = MakeConeFace(e, v0); + m_edgesToUpdate.push_back(e); + } + e = tmp; + } + while (e != eHead); + return true; + } + bool ICHull::MakeCCW(CircularListElement * f, + CircularListElement * e, + CircularListElement * v) + { + // the visible face adjacent to e + CircularListElement * fv; + if (e->GetData().m_triangles[0]->GetData().m_visible) + { + fv = e->GetData().m_triangles[0]; + } + else + { + fv = e->GetData().m_triangles[1]; + } + + // set vertex[0] and vertex[1] to have the same orientation as the corresponding vertices of fv. + long i; // index of e->m_vertices[0] in fv + CircularListElement * v0 = e->GetData().m_vertices[0]; + CircularListElement * v1 = e->GetData().m_vertices[1]; + for(i = 0; fv->GetData().m_vertices[i] != v0; i++); + + if ( fv->GetData().m_vertices[(i+1) % 3] != e->GetData().m_vertices[1] ) + { + f->GetData().m_vertices[0] = v1; + f->GetData().m_vertices[1] = v0; + } + else + { + f->GetData().m_vertices[0] = v0; + f->GetData().m_vertices[1] = v1; + // swap edges + CircularListElement * tmp = f->GetData().m_edges[0]; + f->GetData().m_edges[0] = f->GetData().m_edges[1]; + f->GetData().m_edges[1] = tmp; + } + f->GetData().m_vertices[2] = v; + return true; + } + bool ICHull::CleanUp(unsigned long & addedPoints) + { + bool r0 = CleanEdges(); + bool r1 = CleanTriangles(); + bool r2 = CleanVertices(addedPoints); + return r0 && r1 && r2; + } + bool ICHull::CleanEdges() + { + // integrate the new faces into the data structure + CircularListElement * e; + const std::vector *>::iterator itEndUpdate = m_edgesToUpdate.end(); + for(std::vector *>::iterator it = m_edgesToUpdate.begin(); it != itEndUpdate; ++it) + { + e = *it; + if ( e->GetData().m_newFace ) + { + if ( e->GetData().m_triangles[0]->GetData().m_visible) + { + e->GetData().m_triangles[0] = e->GetData().m_newFace; + } + else + { + e->GetData().m_triangles[1] = e->GetData().m_newFace; + } + e->GetData().m_newFace = 0; + } + } + // delete edges maked for deletion + CircularList & edges = m_mesh.GetEdges(); + const std::vector *>::iterator itEndDelete = m_edgesToDelete.end(); + for(std::vector *>::iterator it = m_edgesToDelete.begin(); it != itEndDelete; ++it) + { + edges.Delete(*it); + } + m_edgesToDelete.clear(); + m_edgesToUpdate.clear(); + return true; + } + bool ICHull::CleanTriangles() + { + CircularList & triangles = m_mesh.GetTriangles(); + const std::vector *>::iterator itEndDelete = m_trianglesToDelete.end(); + for(std::vector *>::iterator it = m_trianglesToDelete.begin(); it != itEndDelete; ++it) + { + if (m_distPoints) + { + if (m_isFlat) + { + // to be updated + } + else + { + std::set::const_iterator itPEnd((*it)->GetData().m_incidentPoints.end()); + std::set::const_iterator itP((*it)->GetData().m_incidentPoints.begin()); + std::map::iterator itPoint; + for(; itP != itPEnd; ++itP) + { + itPoint = m_distPoints->find(*itP); + if (itPoint != m_distPoints->end()) + { + itPoint->second.m_computed = false; + } + } + } + } + triangles.Delete(*it); + } + m_trianglesToDelete.clear(); + return true; + } + bool ICHull::CleanVertices(unsigned long & addedPoints) + { + // mark all vertices incident to some undeleted edge as on the hull + CircularList & edges = m_mesh.GetEdges(); + CircularListElement * e = edges.GetHead(); + size_t nE = edges.GetSize(); + for(size_t i = 0; i < nE; i++) + { + e->GetData().m_vertices[0]->GetData().m_onHull = true; + e->GetData().m_vertices[1]->GetData().m_onHull = true; + e = e->GetNext(); + } + // delete all the vertices that have been processed but are not on the hull + CircularList & vertices = m_mesh.GetVertices(); + CircularListElement * vHead = vertices.GetHead(); + CircularListElement * v = vHead; + v = v->GetPrev(); + do + { + if (v->GetData().m_tag && !v->GetData().m_onHull) + { + CircularListElement * tmp = v->GetPrev(); + vertices.Delete(v); + v = tmp; + addedPoints--; + } + else + { + v->GetData().m_duplicate = 0; + v->GetData().m_onHull = false; + v = v->GetPrev(); + } + } + while (v->GetData().m_tag && v != vHead); + return true; + } + void ICHull::Clear() + { + m_mesh.Clear(); + m_edgesToDelete = std::vector *>(); + m_edgesToUpdate = std::vector *>(); + m_trianglesToDelete= std::vector *>(); + m_isFlat = false; + } + const ICHull & ICHull::operator=(ICHull & rhs) + { + if (&rhs != this) + { + m_mesh.Copy(rhs.m_mesh); + m_edgesToDelete = rhs.m_edgesToDelete; + m_edgesToUpdate = rhs.m_edgesToUpdate; + m_trianglesToDelete = rhs.m_trianglesToDelete; + m_isFlat = rhs.m_isFlat; + } + return (*this); + } + double ICHull::ComputeVolume() + { + size_t nV = m_mesh.m_vertices.GetSize(); + if (nV == 0 || m_isFlat) + { + return 0.0; + } + Vec3 bary(0.0, 0.0, 0.0); + for(size_t v = 0; v < nV; v++) + { + bary.X() += m_mesh.m_vertices.GetHead()->GetData().m_pos.X(); + bary.Y() += m_mesh.m_vertices.GetHead()->GetData().m_pos.Y(); + bary.Z() += m_mesh.m_vertices.GetHead()->GetData().m_pos.Z(); + m_mesh.m_vertices.Next(); + } + bary /= nV; + + size_t nT = m_mesh.m_triangles.GetSize(); + Vec3 ver0, ver1, ver2; + double totalVolume = 0.0; + for(size_t t = 0; t < nT; t++) + { + ver0.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.X(); + ver0.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Y(); + ver0.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Z(); + ver1.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.X(); + ver1.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Y(); + ver1.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Z(); + ver2.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.X(); + ver2.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Y(); + ver2.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Z(); + totalVolume += Volume(ver0, ver1, ver2, bary); + m_mesh.m_triangles.Next(); + } + return totalVolume; + } + bool ICHull::IsInside(const Vec3 & pt0) + { + const Vec3 pt(pt0.X(), pt0.Y(), pt0.Z()); + if (m_isFlat) + { + size_t nT = m_mesh.m_triangles.GetSize(); + Vec3 ver0, ver1, ver2, a, b, c; + double u,v; + for(size_t t = 0; t < nT; t++) + { + ver0.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.X(); + ver0.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Y(); + ver0.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Z(); + ver1.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.X(); + ver1.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Y(); + ver1.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Z(); + ver2.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.X(); + ver2.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Y(); + ver2.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Z(); + a = ver1 - ver0; + b = ver2 - ver0; + c = pt - ver0; + u = c * a; + v = c * b; + if ( u >= 0.0 && u <= 1.0 && v >= 0.0 && u+v <= 1.0) + { + return true; + } + m_mesh.m_triangles.Next(); + } + return false; + } + else + { + size_t nT = m_mesh.m_triangles.GetSize(); + Vec3 ver0, ver1, ver2; + for(size_t t = 0; t < nT; t++) + { + ver0.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.X(); + ver0.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Y(); + ver0.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[0]->GetData().m_pos.Z(); + ver1.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.X(); + ver1.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Y(); + ver1.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[1]->GetData().m_pos.Z(); + ver2.X() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.X(); + ver2.Y() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Y(); + ver2.Z() = m_mesh.m_triangles.GetHead()->GetData().m_vertices[2]->GetData().m_pos.Z(); + if (Volume(ver0, ver1, ver2, pt) < 0.0) + { + return false; + } + m_mesh.m_triangles.Next(); + } + return true; + } + } + double ICHull::ComputeDistance(long name, const Vec3 & pt, const Vec3 & normal, bool & insideHull, bool updateIncidentPoints) + { + Vec3 ptNormal(static_cast(normal.X()), + static_cast(normal.Y()), + static_cast(normal.Z())); + Vec3 p0( static_cast(pt.X()), + static_cast(pt.Y()), + static_cast(pt.Z())); + + if (m_isFlat) + { + double distance = 0.0; + Vec3 chNormal(static_cast(m_normal.X()), + static_cast(m_normal.Y()), + static_cast(m_normal.Z())); + ptNormal -= (ptNormal * chNormal) * chNormal; + if (ptNormal.GetNorm() > 0.0) + { + ptNormal.Normalize(); + long nameVE1; + long nameVE2; + Vec3 pa, pb, d0, d1, d2, d3; + Vec3 p1 = p0 + ptNormal; + Vec3 p2, p3; + double mua, mub, s; + const double EPS = 0.00000000001; + size_t nE = m_mesh.GetNEdges(); + for(size_t e = 0; e < nE; e++) + { + TMMEdge & currentEdge = m_mesh.m_edges.GetHead()->GetData(); + nameVE1 = currentEdge.m_vertices[0]->GetData().m_name; + nameVE2 = currentEdge.m_vertices[1]->GetData().m_name; + if (currentEdge.m_triangles[0] == 0 || currentEdge.m_triangles[1] == 0) + { + if ( nameVE1==name || nameVE2==name ) + { + return 0.0; + } + /* + if (debug) std::cout << "V" << name + << " E " << nameVE1 << " " << nameVE2 << std::endl; + */ + + p2.X() = currentEdge.m_vertices[0]->GetData().m_pos.X(); + p2.Y() = currentEdge.m_vertices[0]->GetData().m_pos.Y(); + p2.Z() = currentEdge.m_vertices[0]->GetData().m_pos.Z(); + p3.X() = currentEdge.m_vertices[1]->GetData().m_pos.X(); + p3.Y() = currentEdge.m_vertices[1]->GetData().m_pos.Y(); + p3.Z() = currentEdge.m_vertices[1]->GetData().m_pos.Z(); + d0 = p3 - p2; + if (d0.GetNorm() > 0.0) + { + if (IntersectLineLine(p0, p1, p2, p3, pa, pb, mua, mub)) + { + d1 = pa - p2; + d2 = pa - pb; + d3 = pa - p0; + mua = d1.GetNorm()/d0.GetNorm(); + mub = d1*d0; + s = d3*ptNormal; + if (d2.GetNorm() < EPS && mua <= 1.0 && mub>=0.0 && s>0.0) + { + distance = std::max(distance, d3.GetNorm()); + } + } + } + } + m_mesh.m_edges.Next(); + } + } + return distance; + } + else + { + Vec3 ptNormal(static_cast(normal.X()), + static_cast(normal.Y()), + static_cast(normal.Z())); + + Vec3 impact; + long nhit; + double dist; + double distance = 0.0; + size_t nT = m_mesh.GetNTriangles(); + insideHull = false; + CircularListElement * face = 0; + Vec3 ver0, ver1, ver2; + for(size_t f = 0; f < nT; f++) + { + TMMTriangle & currentTriangle = m_mesh.m_triangles.GetHead()->GetData(); + /* + if (debug) std::cout << "T " << currentTriangle.m_vertices[0]->GetData().m_name << " " + << currentTriangle.m_vertices[1]->GetData().m_name << " " + << currentTriangle.m_vertices[2]->GetData().m_name << std::endl; + */ + if (currentTriangle.m_vertices[0]->GetData().m_name == name || + currentTriangle.m_vertices[1]->GetData().m_name == name || + currentTriangle.m_vertices[2]->GetData().m_name == name) + { + nhit = 1; + dist = 0.0; + } + else + { + ver0.X() = currentTriangle.m_vertices[0]->GetData().m_pos.X(); + ver0.Y() = currentTriangle.m_vertices[0]->GetData().m_pos.Y(); + ver0.Z() = currentTriangle.m_vertices[0]->GetData().m_pos.Z(); + ver1.X() = currentTriangle.m_vertices[1]->GetData().m_pos.X(); + ver1.Y() = currentTriangle.m_vertices[1]->GetData().m_pos.Y(); + ver1.Z() = currentTriangle.m_vertices[1]->GetData().m_pos.Z(); + ver2.X() = currentTriangle.m_vertices[2]->GetData().m_pos.X(); + ver2.Y() = currentTriangle.m_vertices[2]->GetData().m_pos.Y(); + ver2.Z() = currentTriangle.m_vertices[2]->GetData().m_pos.Z(); + nhit = IntersectRayTriangle(p0, ptNormal, ver0, ver1, ver2, dist); + } + + if (nhit == 1 && distance <= dist) + { + distance = dist; + insideHull = true; + face = m_mesh.m_triangles.GetHead(); +/* + std::cout << name << " -> T " << currentTriangle.m_vertices[0]->GetData().m_name << " " + << currentTriangle.m_vertices[1]->GetData().m_name << " " + << currentTriangle.m_vertices[2]->GetData().m_name << " Dist " + << dist << " P " << currentTriangle.m_normal * normal << std::endl; +*/ + if (dist > 0.1) + { + break; + } + } + m_mesh.m_triangles.Next(); + } + if (updateIncidentPoints && face && m_distPoints) + { + (*m_distPoints)[name].m_dist = static_cast(distance); + face->GetData().m_incidentPoints.insert(name); + } + return distance; + } + } +} + diff --git a/Extras/HACD/hacdICHull.h b/Extras/HACD/hacdICHull.h new file mode 100644 index 000000000..29ae6d53c --- /dev/null +++ b/Extras/HACD/hacdICHull.h @@ -0,0 +1,120 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_ICHULL_H +#define HACD_ICHULL_H +#include +#include +#include +#include +#include +namespace HACD +{ + class DPoint; + class HACD; + //! Incremental Convex Hull algorithm (cf. http://maven.smith.edu/~orourke/books/ftp.html ). + enum ICHullError + { + ICHullErrorOK = 0, + ICHullErrorCoplanarPoints, + ICHullErrorNoVolume, + ICHullErrorInconsistent, + ICHullErrorNotEnoughPoints + }; + class ICHull + { + public: + //! + bool IsFlat() { return m_isFlat;} + //! + std::map * GetDistPoints() const { return m_distPoints;} + //! + void SetDistPoints(std::map * distPoints) { m_distPoints = distPoints;} + //! Returns the computed mesh + TMMesh & GetMesh() { return m_mesh;} + //! Add one point to the convex-hull + bool AddPoint(const Vec3 & point) {return AddPoints(&point, 1);} + //! Add one point to the convex-hull + bool AddPoint(const Vec3 & point, long id); + //! Add points to the convex-hull + bool AddPoints(const Vec3 * points, size_t nPoints); + bool AddPoints(std::vector< Vec3 > points); + //! + ICHullError Process(); + //! + ICHullError Process(unsigned long nPointsCH); + //! + double ComputeVolume(); + //! + bool IsInside(const Vec3 & pt0); + //! + double ComputeDistance(long name, const Vec3 & pt, const Vec3 & normal, bool & insideHull, bool updateIncidentPoints); + //! + const ICHull & operator=(ICHull & rhs); + + //! Constructor + ICHull(void); + //! Destructor + virtual ~ICHull(void) {}; + + private: + //! DoubleTriangle builds the initial double triangle. It first finds 3 noncollinear points and makes two faces out of them, in opposite order. It then finds a fourth point that is not coplanar with that face. The vertices are stored in the face structure in counterclockwise order so that the volume between the face and the point is negative. Lastly, the 3 newfaces to the fourth point are constructed and the data structures are cleaned up. + ICHullError DoubleTriangle(); + //! MakeFace creates a new face structure from three vertices (in ccw order). It returns a pointer to the face. + CircularListElement * MakeFace(CircularListElement * v0, + CircularListElement * v1, + CircularListElement * v2, + CircularListElement * fold); + //! + CircularListElement * MakeConeFace(CircularListElement * e, CircularListElement * v); + //! + bool ProcessPoint(); + //! + bool ComputePointVolume(double &totalVolume, bool markVisibleFaces); + //! + bool FindMaxVolumePoint(); + //! + bool CleanEdges(); + //! + bool CleanVertices(unsigned long & addedPoints); + //! + bool CleanTriangles(); + //! + bool CleanUp(unsigned long & addedPoints); + //! + bool MakeCCW(CircularListElement * f, + CircularListElement * e, + CircularListElement * v); + void Clear(); + private: + static const long sc_dummyIndex; + static const double sc_distMin; + TMMesh m_mesh; + std::vector *> m_edgesToDelete; + std::vector *> m_edgesToUpdate; + std::vector *> m_trianglesToDelete; + std::map * m_distPoints; + CircularListElement * m_dummyVertex; + Vec3 m_normal; + bool m_isFlat; + + + ICHull(const ICHull & rhs); + + friend class HACD; + }; + +} +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdManifoldMesh.cpp b/Extras/HACD/hacdManifoldMesh.cpp new file mode 100644 index 000000000..e3bb55c0c --- /dev/null +++ b/Extras/HACD/hacdManifoldMesh.cpp @@ -0,0 +1,577 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#include +using namespace std; + + +namespace HACD +{ + Material::Material(void) + { + m_diffuseColor.X() = 0.5; + m_diffuseColor.Y() = 0.5; + m_diffuseColor.Z() = 0.5; + m_specularColor.X() = 0.5; + m_specularColor.Y() = 0.5; + m_specularColor.Z() = 0.5; + m_ambientIntensity = 0.4; + m_emissiveColor.X() = 0.0; + m_emissiveColor.Y() = 0.0; + m_emissiveColor.Z() = 0.0; + m_shininess = 0.4; + m_transparency = 0.0; + } + + TMMVertex::TMMVertex(void) + { + m_name = 0; + m_id = 0; + m_duplicate = 0; + m_onHull = false; + m_tag = false; + } + TMMVertex::~TMMVertex(void) + { + } + TMMEdge::TMMEdge(void) + { + m_id = 0; + m_triangles[0] = m_triangles[1] = m_newFace = 0; + m_vertices[0] = m_vertices[1] = 0; + } + TMMEdge::~TMMEdge(void) + { + } + TMMTriangle::TMMTriangle(void) + { + m_id = 0; + for(int i = 0; i < 3; i++) + { + m_edges[i] = 0; + m_vertices[0] = 0; + } + m_visible = false; + } + TMMTriangle::~TMMTriangle(void) + { + } + TMMesh::TMMesh(void) + { + m_barycenter = Vec3(0,0,0); + m_diag = 1; + } + TMMesh::~TMMesh(void) + { + } + + void TMMesh::Print() + { + size_t nV = m_vertices.GetSize(); + std::cout << "-----------------------------" << std::endl; + std::cout << "vertices (" << nV << ")" << std::endl; + for(size_t v = 0; v < nV; v++) + { + const TMMVertex & currentVertex = m_vertices.GetData(); + std::cout << currentVertex.m_id << ", " + << currentVertex.m_pos.X() << ", " + << currentVertex.m_pos.Y() << ", " + << currentVertex.m_pos.Z() << std::endl; + m_vertices.Next(); + } + + + size_t nE = m_edges.GetSize(); + std::cout << "edges (" << nE << ")" << std::endl; + for(size_t e = 0; e < nE; e++) + { + const TMMEdge & currentEdge = m_edges.GetData(); + const CircularListElement * v0 = currentEdge.m_vertices[0]; + const CircularListElement * v1 = currentEdge.m_vertices[1]; + const CircularListElement * f0 = currentEdge.m_triangles[0]; + const CircularListElement * f1 = currentEdge.m_triangles[1]; + + std::cout << "-> (" << v0->GetData().m_name << ", " << v1->GetData().m_name << ")" << std::endl; + std::cout << "-> F0 (" << f0->GetData().m_vertices[0]->GetData().m_name << ", " + << f0->GetData().m_vertices[1]->GetData().m_name << ", " + << f0->GetData().m_vertices[2]->GetData().m_name <<")" << std::endl; + std::cout << "-> F1 (" << f1->GetData().m_vertices[0]->GetData().m_name << ", " + << f1->GetData().m_vertices[1]->GetData().m_name << ", " + << f1->GetData().m_vertices[2]->GetData().m_name << ")" << std::endl; + m_edges.Next(); + } + size_t nT = m_triangles.GetSize(); + std::cout << "triangles (" << nT << ")" << std::endl; + for(size_t t = 0; t < nT; t++) + { + const TMMTriangle & currentTriangle = m_triangles.GetData(); + const CircularListElement * v0 = currentTriangle.m_vertices[0]; + const CircularListElement * v1 = currentTriangle.m_vertices[1]; + const CircularListElement * v2 = currentTriangle.m_vertices[2]; + const CircularListElement * e0 = currentTriangle.m_edges[0]; + const CircularListElement * e1 = currentTriangle.m_edges[1]; + const CircularListElement * e2 = currentTriangle.m_edges[2]; + + std::cout << "-> (" << v0->GetData().m_name << ", " << v1->GetData().m_name << ", "<< v2->GetData().m_name << ")" << std::endl; + + std::cout << "-> E0 (" << e0->GetData().m_vertices[0]->GetData().m_name << ", " + << e0->GetData().m_vertices[1]->GetData().m_name << ")" << std::endl; + std::cout << "-> E1 (" << e1->GetData().m_vertices[0]->GetData().m_name << ", " + << e1->GetData().m_vertices[1]->GetData().m_name << ")" << std::endl; + std::cout << "-> E2 (" << e2->GetData().m_vertices[0]->GetData().m_name << ", " + << e2->GetData().m_vertices[1]->GetData().m_name << ")" << std::endl; + m_triangles.Next(); + } + } + bool TMMesh::Save(const char *fileName) + { + std::ofstream fout(fileName); + std::cout << "Saving " << fileName << std::endl; + if (SaveVRML2(fout)) + { + fout.close(); + return true; + } + return false; + } + bool TMMesh::SaveVRML2(std::ofstream &fout) + { + return SaveVRML2(fout, Material()); + } + bool TMMesh::SaveVRML2(std::ofstream &fout, const Material & material) + { + if (fout.is_open()) + { + size_t nV = m_vertices.GetSize(); + size_t nT = m_triangles.GetSize(); + fout <<"#VRML V2.0 utf8" << std::endl; + fout <<"" << std::endl; + fout <<"# Vertices: " << nV << std::endl; + fout <<"# Triangles: " << nT << std::endl; + fout <<"" << std::endl; + fout <<"Group {" << std::endl; + fout <<" children [" << std::endl; + fout <<" Shape {" << std::endl; + fout <<" appearance Appearance {" << std::endl; + fout <<" material Material {" << std::endl; + fout <<" diffuseColor " << material.m_diffuseColor.X() << " " + << material.m_diffuseColor.Y() << " " + << material.m_diffuseColor.Z() << std::endl; + fout <<" ambientIntensity " << material.m_ambientIntensity << std::endl; + fout <<" specularColor " << material.m_specularColor.X() << " " + << material.m_specularColor.Y() << " " + << material.m_specularColor.Z() << std::endl; + fout <<" emissiveColor " << material.m_emissiveColor.X() << " " + << material.m_emissiveColor.Y() << " " + << material.m_emissiveColor.Z() << std::endl; + fout <<" shininess " << material.m_shininess << std::endl; + fout <<" transparency " << material.m_transparency << std::endl; + fout <<" }" << std::endl; + fout <<" }" << std::endl; + fout <<" geometry IndexedFaceSet {" << std::endl; + fout <<" ccw TRUE" << std::endl; + fout <<" solid TRUE" << std::endl; + fout <<" convex TRUE" << std::endl; + if (GetNVertices() > 0) { + fout <<" coord DEF co Coordinate {" << std::endl; + fout <<" point [" << std::endl; + for(size_t v = 0; v < nV; v++) + { + TMMVertex & currentVertex = m_vertices.GetData(); + fout <<" " << currentVertex.m_pos.X() << " " + << currentVertex.m_pos.Y() << " " + << currentVertex.m_pos.Z() << "," << std::endl; + currentVertex.m_id = v; + m_vertices.Next(); + } + fout <<" ]" << std::endl; + fout <<" }" << std::endl; + } + if (GetNTriangles() > 0) { + fout <<" coordIndex [ " << std::endl; + for(size_t f = 0; f < nT; f++) + { + TMMTriangle & currentTriangle = m_triangles.GetData(); + fout <<" " << currentTriangle.m_vertices[0]->GetData().m_id << ", " + << currentTriangle.m_vertices[1]->GetData().m_id << ", " + << currentTriangle.m_vertices[2]->GetData().m_id << ", -1," << std::endl; + m_triangles.Next(); + } + fout <<" ]" << std::endl; + } + fout <<" }" << std::endl; + fout <<" }" << std::endl; + fout <<" ]" << std::endl; + fout <<"}" << std::endl; + } + return true; + } + void TMMesh::GetIFS(Vec3 * const points, Vec3 * const triangles) + { + size_t nV = m_vertices.GetSize(); + size_t nT = m_triangles.GetSize(); + + for(size_t v = 0; v < nV; v++) + { + points[v] = m_vertices.GetData().m_pos; + m_vertices.GetData().m_id = v; + m_vertices.Next(); + } + for(size_t f = 0; f < nT; f++) + { + TMMTriangle & currentTriangle = m_triangles.GetData(); + triangles[f].X() = currentTriangle.m_vertices[0]->GetData().m_id; + triangles[f].Y() = currentTriangle.m_vertices[1]->GetData().m_id; + triangles[f].Z() = currentTriangle.m_vertices[2]->GetData().m_id; + m_triangles.Next(); + } + } + void TMMesh::Clear() + { + m_vertices.Clear(); + m_edges.Clear(); + m_triangles.Clear(); + } + void TMMesh::Copy(TMMesh & mesh) + { + Clear(); + // updating the id's + size_t nV = mesh.m_vertices.GetSize(); + size_t nE = mesh. m_edges.GetSize(); + size_t nT = mesh.m_triangles.GetSize(); + for(size_t v = 0; v < nV; v++) + { + mesh.m_vertices.GetData().m_id = v; + mesh.m_vertices.Next(); + } + for(size_t e = 0; e < nE; e++) + { + mesh.m_edges.GetData().m_id = e; + mesh.m_edges.Next(); + + } + for(size_t f = 0; f < nT; f++) + { + mesh.m_triangles.GetData().m_id = f; + mesh.m_triangles.Next(); + } + // copying data + m_vertices = mesh.m_vertices; + m_edges = mesh.m_edges; + m_triangles = mesh.m_triangles; + + // generating mapping + CircularListElement ** vertexMap = new CircularListElement * [nV]; + CircularListElement ** edgeMap = new CircularListElement * [nE]; + CircularListElement ** triangleMap = new CircularListElement * [nT]; + for(size_t v = 0; v < nV; v++) + { + vertexMap[v] = m_vertices.GetHead(); + m_vertices.Next(); + } + for(size_t e = 0; e < nE; e++) + { + edgeMap[e] = m_edges.GetHead(); + m_edges.Next(); + } + for(size_t f = 0; f < nT; f++) + { + triangleMap[f] = m_triangles.GetHead(); + m_triangles.Next(); + } + + // updating pointers + for(size_t v = 0; v < nV; v++) + { + if (vertexMap[v]->GetData().m_duplicate) + { + vertexMap[v]->GetData().m_duplicate = edgeMap[vertexMap[v]->GetData().m_duplicate->GetData().m_id]; + } + } + for(size_t e = 0; e < nE; e++) + { + if (edgeMap[e]->GetData().m_newFace) + { + edgeMap[e]->GetData().m_newFace = triangleMap[edgeMap[e]->GetData().m_newFace->GetData().m_id]; + } + if (nT > 0) + { + for(int f = 0; f < 2; f++) + { + if (edgeMap[e]->GetData().m_triangles[f]) + { + edgeMap[e]->GetData().m_triangles[f] = triangleMap[edgeMap[e]->GetData().m_triangles[f]->GetData().m_id]; + } + } + } + for(int v = 0; v < 2; v++) + { + if (edgeMap[e]->GetData().m_vertices[v]) + { + edgeMap[e]->GetData().m_vertices[v] = vertexMap[edgeMap[e]->GetData().m_vertices[v]->GetData().m_id]; + } + } + } + for(size_t f = 0; f < nT; f++) + { + if (nE > 0) + { + for(int e = 0; e < 3; e++) + { + if (triangleMap[f]->GetData().m_edges[e]) + { + triangleMap[f]->GetData().m_edges[e] = edgeMap[triangleMap[f]->GetData().m_edges[e]->GetData().m_id]; + } + } + } + for(int v = 0; v < 3; v++) + { + if (triangleMap[f]->GetData().m_vertices[v]) + { + triangleMap[f]->GetData().m_vertices[v] = vertexMap[triangleMap[f]->GetData().m_vertices[v]->GetData().m_id]; + } + } + } + delete [] vertexMap; + delete [] edgeMap; + delete [] triangleMap; + + } + long IntersectRayTriangle(const Vec3 & P0, const Vec3 & dir, + const Vec3 & V0, const Vec3 & V1, + const Vec3 & V2, double &t) + { + Vec3 edge1, edge2, edge3; + double det, invDet; + edge1 = V1 - V2; + edge2 = V2 - V0; + Vec3 pvec = dir ^ edge2; + det = edge1 * pvec; + if (det == 0.0) + return 0; + invDet = 1.0/det; + Vec3 tvec = P0 - V0; + Vec3 qvec = tvec ^ edge1; + t = (edge2 * qvec) * invDet; + if (t < 0.0) + { + return 0; + } + edge3 = V0 - V1; + Vec3 I(P0 + t * dir); + Vec3 s0 = (I-V0) ^ edge3; + Vec3 s1 = (I-V1) ^ edge1; + Vec3 s2 = (I-V2) ^ edge2; + if (s0*s1 >= 0.0 && s2*s1 >= 0.0) + { + return 1; + } + return 0; + } + + bool IntersectLineLine(const Vec3 & p1, const Vec3 & p2, + const Vec3 & p3, const Vec3 & p4, + Vec3 & pa, Vec3 & pb, + double & mua, double & mub) + { + Vec3 p13,p43,p21; + double d1343,d4321,d1321,d4343,d2121; + double numer,denom; + + p13.X() = p1.X() - p3.X(); + p13.Y() = p1.Y() - p3.Y(); + p13.Z() = p1.Z() - p3.Z(); + p43.X() = p4.X() - p3.X(); + p43.Y() = p4.Y() - p3.Y(); + p43.Z() = p4.Z() - p3.Z(); + if (p43.X()==0.0 && p43.Y()==0.0 && p43.Z()==0.0) + return false; + p21.X() = p2.X() - p1.X(); + p21.Y() = p2.Y() - p1.Y(); + p21.Z() = p2.Z() - p1.Z(); + if (p21.X()==0.0 && p21.Y()==0.0 && p21.Z()==0.0) + return false; + + d1343 = p13.X() * p43.X() + p13.Y() * p43.Y() + p13.Z() * p43.Z(); + d4321 = p43.X() * p21.X() + p43.Y() * p21.Y() + p43.Z() * p21.Z(); + d1321 = p13.X() * p21.X() + p13.Y() * p21.Y() + p13.Z() * p21.Z(); + d4343 = p43.X() * p43.X() + p43.Y() * p43.Y() + p43.Z() * p43.Z(); + d2121 = p21.X() * p21.X() + p21.Y() * p21.Y() + p21.Z() * p21.Z(); + + denom = d2121 * d4343 - d4321 * d4321; + if (denom==0.0) + return false; + numer = d1343 * d4321 - d1321 * d4343; + + mua = numer / denom; + mub = (d1343 + d4321 * (mua)) / d4343; + + pa.X() = p1.X() + mua * p21.X(); + pa.Y() = p1.Y() + mua * p21.Y(); + pa.Z() = p1.Z() + mua * p21.Z(); + pb.X() = p3.X() + mub * p43.X(); + pb.Y() = p3.Y() + mub * p43.Y(); + pb.Z() = p3.Z() + mub * p43.Z(); + + return true; + } + + long IntersectRayTriangle2(const Vec3 & P0, const Vec3 & dir, + const Vec3 & V0, const Vec3 & V1, + const Vec3 & V2, double &r) + { + Vec3 u, v, n; // triangle vectors + Vec3 w0, w; // ray vectors + double a, b; // params to calc ray-plane intersect + + // get triangle edge vectors and plane normal + u = V1 - V0; + v = V2 - V0; + n = u ^ v; // cross product + if (n.GetNorm() == 0.0) // triangle is degenerate + return -1; // do not deal with this case + + w0 = P0 - V0; + a = - n * w0; + b = n * dir; + if (fabs(b) <= 0.0) { // ray is parallel to triangle plane + if (a == 0.0) // ray lies in triangle plane + return 2; + else return 0; // ray disjoint from plane + } + + // get intersect point of ray with triangle plane + r = a / b; + if (r < 0.0) // ray goes away from triangle + return 0; // => no intersect + // for a segment, also test if (r > 1.0) => no intersect + + Vec3 I = P0 + r * dir; // intersect point of ray and plane + + // is I inside T? + double uu, uv, vv, wu, wv, D; + uu = u * u; + uv = u * v; + vv = v * v; + w = I - V0; + wu = w * u; + wv = w * v; + D = uv * uv - uu * vv; + + // get and test parametric coords + double s, t; + s = (uv * wv - vv * wu) / D; + if (s < 0.0 || s > 1.0) // I is outside T + return 0; + t = (uv * wu - uu * wv) / D; + if (t < 0.0 || (s + t) > 1.0) // I is outside T + return 0; + return 1; // I is in T + } + + + bool TMMesh::CheckConsistancy() + { + size_t nE = m_edges.GetSize(); + size_t nT = m_triangles.GetSize(); + for(size_t e = 0; e < nE; e++) + { + for(int f = 0; f < 2; f++) + { + if (!m_edges.GetHead()->GetData().m_triangles[f]) + { + return false; + } + } + m_edges.Next(); + } + + for(size_t f = 0; f < nT; f++) + { + for(int e = 0; e < 3; e++) + { + int found = 0; + for(int k = 0; k < 2; k++) + { + if (m_triangles.GetHead()->GetData().m_edges[e]->GetData().m_triangles[k] == m_triangles.GetHead()) + { + found++; + } + } + if (found != 1) + { + return false; + } + } + m_triangles.Next(); + } + + return true; + } + bool TMMesh::Normalize() + { + size_t nV = m_vertices.GetSize(); + if (nV == 0) + { + return false; + } + m_barycenter = m_vertices.GetHead()->GetData().m_pos; + Vec3 min = m_barycenter; + Vec3 max = m_barycenter; + Real x, y, z; + for(size_t v = 1; v < nV; v++) + { + m_barycenter += m_vertices.GetHead()->GetData().m_pos; + x = m_vertices.GetHead()->GetData().m_pos.X(); + y = m_vertices.GetHead()->GetData().m_pos.Y(); + z = m_vertices.GetHead()->GetData().m_pos.Z(); + if ( x < min.X()) min.X() = x; + else if ( x > max.X()) max.X() = x; + if ( y < min.Y()) min.Y() = y; + else if ( y > max.Y()) max.Y() = y; + if ( z < min.Z()) min.Z() = z; + else if ( z > max.Z()) max.Z() = z; + m_vertices.Next(); + } + m_barycenter /= static_cast(nV); + m_diag = static_cast(0.001 * (max-min).GetNorm()); + const Real invDiag = static_cast(1.0 / m_diag); + if (m_diag != 0.0) + { + for(size_t v = 0; v < nV; v++) + { + m_vertices.GetHead()->GetData().m_pos = (m_vertices.GetHead()->GetData().m_pos - m_barycenter) * invDiag; + m_vertices.Next(); + } + } + return true; + } + bool TMMesh::Denormalize() + { + size_t nV = m_vertices.GetSize(); + if (nV == 0) + { + return false; + } + if (m_diag != 0.0) + { + for(size_t v = 0; v < nV; v++) + { + m_vertices.GetHead()->GetData().m_pos = m_vertices.GetHead()->GetData().m_pos * m_diag + m_barycenter; + m_vertices.Next(); + } + } + return false; + } +} \ No newline at end of file diff --git a/Extras/HACD/hacdManifoldMesh.h b/Extras/HACD/hacdManifoldMesh.h new file mode 100644 index 000000000..3031480a9 --- /dev/null +++ b/Extras/HACD/hacdManifoldMesh.h @@ -0,0 +1,250 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) +All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. +*/ + +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_MANIFOLD_MESH_H +#define HACD_MANIFOLD_MESH_H +#include +#include +#include +#include +#include +#include +namespace HACD +{ + class TMMTriangle; + class TMMEdge; + class TMMesh; + class ICHull; + class HACD; + + class DPoint + { + public: + DPoint(Real dist=0, bool computed=false, bool distOnly=false) + :m_dist(dist), + m_computed(computed), + m_distOnly(distOnly){}; + ~DPoint(){}; + private: + Real m_dist; + bool m_computed; + bool m_distOnly; + friend class TMMTriangle; + friend class TMMesh; + friend class GraphVertex; + friend class GraphEdge; + friend class Graph; + friend class ICHull; + friend class HACD; + }; + + //! Vertex data structure used in a triangular manifold mesh (TMM). + class TMMVertex + { + public: + TMMVertex(void); + ~TMMVertex(void); + + private: + Vec3 m_pos; + long m_name; + size_t m_id; + CircularListElement * m_duplicate; // pointer to incident cone edge (or NULL) + bool m_onHull; + bool m_tag; + TMMVertex(const TMMVertex & rhs); + + friend class HACD; + friend class ICHull; + friend class TMMesh; + friend class TMMTriangle; + friend class TMMEdge; + }; + + //! Edge data structure used in a triangular manifold mesh (TMM). + class TMMEdge + { + public: + TMMEdge(void); + ~TMMEdge(void); + private: + size_t m_id; + CircularListElement * m_triangles[2]; + CircularListElement * m_vertices[2]; + CircularListElement * m_newFace; + + + TMMEdge(const TMMEdge & rhs); + + friend class HACD; + friend class ICHull; + friend class TMMTriangle; + friend class TMMVertex; + friend class TMMesh; + }; + + //! Triangle data structure used in a triangular manifold mesh (TMM). + class TMMTriangle + { + public: + TMMTriangle(void); + ~TMMTriangle(void); + private: + size_t m_id; + CircularListElement * m_edges[3]; + CircularListElement * m_vertices[3]; + std::set m_incidentPoints; + bool m_visible; + + TMMTriangle(const TMMTriangle & rhs); + + friend class HACD; + friend class ICHull; + friend class TMMesh; + friend class TMMVertex; + friend class TMMEdge; + }; + + class Material + { + public: + Material(void); + ~Material(void){} +// private: + Vec3 m_diffuseColor; + double m_ambientIntensity; + Vec3 m_specularColor; + Vec3 m_emissiveColor; + double m_shininess; + double m_transparency; + + friend class TMMesh; + friend class HACD; + }; + + //! triangular manifold mesh data structure. + class TMMesh + { + public: + + //! Returns the number of vertices> + inline size_t GetNVertices() const { return m_vertices.GetSize();} + //! Returns the number of edges + inline size_t GetNEdges() const { return m_edges.GetSize();} + //! Returns the number of triangles + inline size_t GetNTriangles() const { return m_triangles.GetSize();} + //! Returns the vertices circular list + inline const CircularList & GetVertices() const { return m_vertices;} + //! Returns the edges circular list + inline const CircularList & GetEdges() const { return m_edges;} + //! Returns the triangles circular list + inline const CircularList & GetTriangles() const { return m_triangles;} + //! Returns the vertices circular list + inline CircularList & GetVertices() { return m_vertices;} + //! Returns the edges circular list + inline CircularList & GetEdges() { return m_edges;} + //! Returns the triangles circular list + inline CircularList & GetTriangles() { return m_triangles;} + //! Add vertex to the mesh + CircularListElement * AddVertex() {return m_vertices.Add();} + //! Add vertex to the mesh + CircularListElement * AddEdge() {return m_edges.Add();} + //! Add vertex to the mesh + CircularListElement * AddTriangle() {return m_triangles.Add();} + //! Print mesh information + void Print(); + //! + void GetIFS(Vec3 * const points, Vec3 * const triangles); + //! Save mesh + bool Save(const char *fileName); + //! Save mesh to VRML 2.0 format + bool SaveVRML2(std::ofstream &fout); + //! Save mesh to VRML 2.0 format + bool SaveVRML2(std::ofstream &fout, const Material & material); + //! + void Clear(); + //! + void Copy(TMMesh & mesh); + //! + bool CheckConsistancy(); + //! + bool Normalize(); + //! + bool Denormalize(); + //! Constructor + TMMesh(void); + //! Destructor + virtual ~TMMesh(void); + + private: + CircularList m_vertices; + CircularList m_edges; + CircularList m_triangles; + Real m_diag; //>! length of the BB diagonal + Vec3 m_barycenter; //>! barycenter of the mesh + + // not defined + TMMesh(const TMMesh & rhs); + friend class ICHull; + friend class HACD; + }; + //! IntersectRayTriangle(): intersect a ray with a 3D triangle + //! Input: a ray R, and a triangle T + //! Output: *I = intersection point (when it exists) + //! 0 = disjoint (no intersect) + //! 1 = intersect in unique point I1 + long IntersectRayTriangle( const Vec3 & P0, const Vec3 & dir, + const Vec3 & V0, const Vec3 & V1, + const Vec3 & V2, double &t); + + // intersect_RayTriangle(): intersect a ray with a 3D triangle + // Input: a ray R, and a triangle T + // Output: *I = intersection point (when it exists) + // Return: -1 = triangle is degenerate (a segment or point) + // 0 = disjoint (no intersect) + // 1 = intersect in unique point I1 + // 2 = are in the same plane + long IntersectRayTriangle2(const Vec3 & P0, const Vec3 & dir, + const Vec3 & V0, const Vec3 & V1, + const Vec3 & V2, double &r); + + /* + Calculate the line segment PaPb that is the shortest route between + two lines P1P2 and P3P4. Calculate also the values of mua and mub where + Pa = P1 + mua (P2 - P1) + Pb = P3 + mub (P4 - P3) + Return FALSE if no solution exists. + */ + bool IntersectLineLine(const Vec3 & p1, const Vec3 & p2, + const Vec3 & p3, const Vec3 & p4, + Vec3 & pa, Vec3 & pb, + double & mua, double &mub); +} +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdVector.h b/Extras/HACD/hacdVector.h new file mode 100644 index 000000000..4f4d3ed79 --- /dev/null +++ b/Extras/HACD/hacdVector.h @@ -0,0 +1,67 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_VECTOR_H +#define HACD_VECTOR_H +#include +#include +#include + +namespace HACD +{ + typedef float Real; + //! Vector dim 3. + template < typename T > class Vec3 + { + public: + T & X(); + T & Y(); + T & Z(); + const T & X() const; + const T & Y() const; + const T & Z() const; + void Normalize(); + T GetNorm() const; + void operator= (const Vec3 & rhs); + void operator+=(const Vec3 & rhs); + void operator-=(const Vec3 & rhs); + void operator-=(T a); + void operator+=(T a); + void operator/=(T a); + void operator*=(T a); + Vec3 operator^ (const Vec3 & rhs) const; + T operator* (const Vec3 & rhs) const; + Vec3 operator+ (const Vec3 & rhs) const; + Vec3 operator- (const Vec3 & rhs) const; + Vec3 operator- () const; + Vec3 operator* (T rhs) const; + Vec3 operator/ (T rhs) const; + Vec3(); + Vec3(T a); + Vec3(T x, T y, T z); + Vec3(const Vec3 & rhs); + /*virtual*/ ~Vec3(void); + + private: + T m_data[3]; + }; + template + const bool Colinear(const Vec3 & a, const Vec3 & b, const Vec3 & c); + template + const T Volume(const Vec3 & a, const Vec3 & b, const Vec3 & c, const Vec3 & d); + +} +#include "hacdVector.inl" // template implementation +#endif \ No newline at end of file diff --git a/Extras/HACD/hacdVector.inl b/Extras/HACD/hacdVector.inl new file mode 100644 index 000000000..3532741c3 --- /dev/null +++ b/Extras/HACD/hacdVector.inl @@ -0,0 +1,178 @@ +#pragma once +#ifndef HACD_VECTOR_INL +#define HACD_VECTOR_INL +namespace HACD +{ + template + inline Vec3 operator*(T lhs, const Vec3 & rhs) + { + return Vec3(lhs * rhs.X(), lhs * rhs.Y(), lhs * rhs.Z()); + } + template + inline T & Vec3::X() + { + return m_data[0]; + } + template + inline T & Vec3::Y() + { + return m_data[1]; + } + template + inline T & Vec3::Z() + { + return m_data[2]; + } + template + inline const T & Vec3::X() const + { + return m_data[0]; + } + template + inline const T & Vec3::Y() const + { + return m_data[1]; + } + template + inline const T & Vec3::Z() const + { + return m_data[2]; + } + template + inline void Vec3::Normalize() + { + T n = sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]+m_data[2]*m_data[2]); + if (n != 0.0) (*this) /= n; + } + template + inline T Vec3::GetNorm() const + { + return sqrt(m_data[0]*m_data[0]+m_data[1]*m_data[1]+m_data[2]*m_data[2]); + } + template + inline void Vec3::operator= (const Vec3 & rhs) + { + this->m_data[0] = rhs.m_data[0]; + this->m_data[1] = rhs.m_data[1]; + this->m_data[2] = rhs.m_data[2]; + } + template + inline void Vec3::operator+=(const Vec3 & rhs) + { + this->m_data[0] += rhs.m_data[0]; + this->m_data[1] += rhs.m_data[1]; + this->m_data[2] += rhs.m_data[2]; + } + template + inline void Vec3::operator-=(const Vec3 & rhs) + { + this->m_data[0] -= rhs.m_data[0]; + this->m_data[1] -= rhs.m_data[1]; + this->m_data[2] -= rhs.m_data[2]; + } + template + inline void Vec3::operator-=(T a) + { + this->m_data[0] -= a; + this->m_data[1] -= a; + this->m_data[2] -= a; + } + template + inline void Vec3::operator+=(T a) + { + this->m_data[0] += a; + this->m_data[1] += a; + this->m_data[2] += a; + } + template + inline void Vec3::operator/=(T a) + { + this->m_data[0] /= a; + this->m_data[1] /= a; + this->m_data[2] /= a; + } + template + inline void Vec3::operator*=(T a) + { + this->m_data[0] *= a; + this->m_data[1] *= a; + this->m_data[2] *= a; + } + template + inline Vec3 Vec3::operator^ (const Vec3 & rhs) const + { + return Vec3(m_data[1] * rhs.m_data[2] - m_data[2] * rhs.m_data[1], + m_data[2] * rhs.m_data[0] - m_data[0] * rhs.m_data[2], + m_data[0] * rhs.m_data[1] - m_data[1] * rhs.m_data[0]); + } + template + inline T Vec3::operator*(const Vec3 & rhs) const + { + return (m_data[0] * rhs.m_data[0] + m_data[1] * rhs.m_data[1] + m_data[2] * rhs.m_data[2]); + } + template + inline Vec3 Vec3::operator+(const Vec3 & rhs) const + { + return Vec3(m_data[0] + rhs.m_data[0],m_data[1] + rhs.m_data[1],m_data[2] + rhs.m_data[2]); + } + template + inline Vec3 Vec3::operator-(const Vec3 & rhs) const + { + return Vec3(m_data[0] - rhs.m_data[0],m_data[1] - rhs.m_data[1],m_data[2] - rhs.m_data[2]) ; + } + template + inline Vec3 Vec3::operator-() const + { + return Vec3(-m_data[0],-m_data[1],-m_data[2]) ; + } + + template + inline Vec3 Vec3::operator*(T rhs) const + { + return Vec3(rhs * this->m_data[0], rhs * this->m_data[1], rhs * this->m_data[2]); + } + template + inline Vec3 Vec3::operator/ (T rhs) const + { + return Vec3(m_data[0] / rhs, m_data[1] / rhs, m_data[2] / rhs); + } + template + inline Vec3::Vec3(T a) + { + m_data[0] = m_data[1] = m_data[2] = a; + } + template + inline Vec3::Vec3(T x, T y, T z) + { + m_data[0] = x; + m_data[1] = y; + m_data[2] = z; + } + template + inline Vec3::Vec3(const Vec3 & rhs) + { + m_data[0] = rhs.m_data[0]; + m_data[1] = rhs.m_data[1]; + m_data[2] = rhs.m_data[2]; + } + template + inline Vec3::~Vec3(void){}; + + template + inline Vec3::Vec3() {} + + template + inline const bool Colinear(const Vec3 & a, const Vec3 & b, const Vec3 & c) + { + return ((c.Z() - a.Z()) * (b.Y() - a.Y()) - (b.Z() - a.Z()) * (c.Y() - a.Y()) == 0.0 /*EPS*/) && + ((b.Z() - a.Z()) * (c.X() - a.X()) - (b.X() - a.X()) * (c.Z() - a.Z()) == 0.0 /*EPS*/) && + ((b.X() - a.X()) * (c.Y() - a.Y()) - (b.Y() - a.Y()) * (c.X() - a.X()) == 0.0 /*EPS*/); + } + + template + inline const T Volume(const Vec3 & a, const Vec3 & b, const Vec3 & c, const Vec3 & d) + { + return (a-d) * ((b-d) ^ (c-d)); + } +} +#endif //HACD_VECTOR_INL \ No newline at end of file diff --git a/Extras/HACD/hacdVersion.h b/Extras/HACD/hacdVersion.h new file mode 100644 index 000000000..3ce69b6f5 --- /dev/null +++ b/Extras/HACD/hacdVersion.h @@ -0,0 +1,20 @@ +/* Copyright (c) 2011 Khaled Mamou (kmamou at gmail dot com) + All rights reserved. + + + Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + + 2. 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. + + 3. The names of the contributors may not 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 HOLDER 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. + */ +#pragma once +#ifndef HACD_VERSION_H +#define HACD_VERSION_H +#define HACD_VERSION_MAJOR 0 +#define HACD_VERSION_MINOR 0 +#endif \ No newline at end of file