From ac4efc84ef30c25c206fc957c656624cdfeef43e Mon Sep 17 00:00:00 2001 From: cnr-isti-vclab Date: Thu, 21 Sep 2006 18:10:05 +0000 Subject: [PATCH] initial commit --- vcg/complex/vertexmesh/update/normal.h | 60 +++++ vcg/math/disjoint_set.h | 123 ++++++++++ vcg/space/normal_extrapolation.h | 328 +++++++++++++++++++++++++ 3 files changed, 511 insertions(+) create mode 100644 vcg/complex/vertexmesh/update/normal.h create mode 100644 vcg/math/disjoint_set.h create mode 100644 vcg/space/normal_extrapolation.h diff --git a/vcg/complex/vertexmesh/update/normal.h b/vcg/complex/vertexmesh/update/normal.h new file mode 100644 index 00000000..2916ab3d --- /dev/null +++ b/vcg/complex/vertexmesh/update/normal.h @@ -0,0 +1,60 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef __VCG_VERTEXMESH_UPDATE_NORMAL +#define __VCG_VERTEXMESH_UPDATE_NORMAL + +#include + +namespace vcg { + namespace vertex { + + /** \addtogroup vertexmesh */ + /* @{ */ + + /*! + * This class is used to update the normals of a Vertex mesh + */ + template < class VERTEX_CONTAINER > + class UpdateNormal + { + public: + typedef VERTEX_CONTAINER VertexContainer; + typedef typename VERTEX_CONTAINER::value_type VertexType; + typedef typename VertexType::ScalarType ScalarType; + typedef typename VERTEX_CONTAINER::iterator VertexIterator; + + /*! + */ + static void UpdateNormals(const VertexIterator &begin, const VertexIterator &end, int k) + { + vcg::NormalExtrapolation< VertexContainer >::ExtrapolateNormlas(begin, end, k); + }; + + }; //end of class UpdateNormal + + /*! @} */ + }; //end of namespace vertex +}; //end of namespace vcg + +#endif //__VCG_VERTEXMESH_UPDATE_NORMAL diff --git a/vcg/math/disjoint_set.h b/vcg/math/disjoint_set.h new file mode 100644 index 00000000..e271f4a9 --- /dev/null +++ b/vcg/math/disjoint_set.h @@ -0,0 +1,123 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef VCG_MATH_UNIONSET_H +#define VCG_MATH_UNIONSET_H + +#include +#include +#include + +namespace vcg +{ + /*! + * Given a set of elements, it is often useful to break them up or partition them into a number of separate, nonoverlapping groups. + * A disjoint-set data structure is a data structure that keeps track of such a partitioning. See + * Diskoint-set data structure on Wikipedia for more details. + */ + template + class DisjointSet + { + /************************************************* + * Inner class definitions + **************************************************/ + struct DisjointSetNode + { + DisjointSetNode(OBJECT_TYPE *x) {obj=x; parent=obj; rank=0;} + OBJECT_TYPE *obj; + OBJECT_TYPE *parent; + int rank; + }; + + typedef OBJECT_TYPE* ObjectPointer; + typedef std::pair< ObjectPointer, int > hPair; + typedef typename stdext::hash_map< ObjectPointer, int >::iterator hIterator; + typedef std::pair< hIterator, bool > hInsertResult; + + public: + /*! + * Default constructor + */ + DisjointSet() {} + + /*! + * Makes a group containing only a given element (a singleton). + */ + void MakeSet(OBJECT_TYPE *x) + { + int object_count = int(inserted_objects.size()); + assert(inserted_objects.find(x)==inserted_objects.end()); //the map mustn't already contain the object x + nodes.push_back(DisjointSetNode(x)); + inserted_objects.insert( hPair(x,object_count) ); + } + + /*! + * Combine or merge two groups into a single group. + */ + void Union(OBJECT_TYPE *x, OBJECT_TYPE *y) + { + OBJECT_TYPE *s0 = FindSet(x); + OBJECT_TYPE *s1 = FindSet(y); + Link(s0, s1); + } + + /*! + * Determine which group a particular element is in. + */ + OBJECT_TYPE* FindSet(OBJECT_TYPE *x) + { + hIterator pos = inserted_objects.find(x); + assert(pos!=inserted_objects.end()); + DisjointSetNode *node = &nodes[pos->second]; + if (node->parent!=x) + node->parent = FindSet(node->parent); + return node->parent; + } + + private: + /* + */ + void Link(OBJECT_TYPE *x, OBJECT_TYPE *y) + { + hIterator xPos = inserted_objects.find(x); + hIterator yPos = inserted_objects.find(y); + assert(xPos!=inserted_objects.end() && yPos!=inserted_objects.end()); + DisjointSetNode *xNode = &nodes[xPos->second]; + DisjointSetNode *yNode = &nodes[yPos->second]; + if (xNode->rank>yNode->rank) + xNode->parent = y; + else + { + yNode->parent = x; + if (xNode->rank==yNode->rank) + yNode->rank++; + } + } + + protected: + stdext::hash_map< OBJECT_TYPE*, int > inserted_objects; + std::vector< DisjointSetNode > nodes; + }; +};// end of namespace vcg + +#endif //VCG_MATH_UNIONSET_H \ No newline at end of file diff --git a/vcg/space/normal_extrapolation.h b/vcg/space/normal_extrapolation.h new file mode 100644 index 00000000..c59df052 --- /dev/null +++ b/vcg/space/normal_extrapolation.h @@ -0,0 +1,328 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +#ifndef VCG_SPACE_NORMAL_EXTRAPOLATION_H +#define VCG_SPACE_NORMAL_EXTRAPOLATION_H + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace vcg +{ + /*! + */ + template < class VERTEX_CONTAINER > + class NormalExtrapolation + { + public: + typedef typename VERTEX_CONTAINER::value_type VertexType; + typedef typename VertexType *VertexPointer; + typedef typename VERTEX_CONTAINER::iterator VertexIterator; + typedef typename VertexType::CoordType CoordType; + typedef typename VertexType::NormalType NormalType; + typedef typename VertexType::ScalarType ScalarType; + typedef typename vcg::Box3< ScalarType > BoundingBoxType; + typedef typename vcg::Matrix33 MatrixType; + + enum NormalOrientation {IsCorrect=0, MustBeFlipped=1}; + + public: + /*! + */ + static void ExtrapolateNormlas(const VertexIterator &begin, const VertexIterator &end, int k, const int root_index=-1, NormalOrientation orientation=IsCorrect) + { + /************************************************* + * Inner class definitions + **************************************************/ + // Dummy class: no object marker is needed + class DummyObjectMarker {}; + // Object functor: return the bounding-box enclosing a given vertex + struct BoundingBoxForVertexFunctor + { + inline BoundingBoxType operator()( const VertexType &vertex ) const + { BoundingBoxType bb; bb.Set(vertex.P()); return bb; } + }; + + // Object functor: compute the distance between a vertex and a point + struct VertPointDistanceFunctor + { + inline bool operator()(const VertexType &v, const CoordType &p, ScalarType &d, CoordType &q) const + { + float distance = vcg::Distance(p, v.P()); + if (distance>d) + return false; + + d = distance; + q = v.P(); + return true; + } + }; + // Plane structure: identify a plain as a pair + struct Plane + { + Plane() { center.Zero(); normal.Zero();}; + CoordType center; + NormalType normal; + int index; + }; + + // Object functor: compute the distance between a point and the plane + struct PlanePointDistanceFunctor + { + inline bool operator()(const Plane &plane, const vcg::Point3f &p, float &d, vcg::Point3f &q) const + { + float distance = vcg::Distance(p, plane.center); + if (distance>d) + return false; + + d = distance; + q = plane.center; + return true; + } + }; + + // Object functor: return the bounding-box enclosing a given plane + struct BoundingBoxForPlaneFunctor + { + inline BoundingBoxType operator()( const Plane &plane ) const + { BoundingBoxType bb; bb.Set(plane.center); return bb; } + }; + + // Represent an edge in the Riemannian graph + struct RiemannianEdge + { + RiemannianEdge(Plane *p=NULL, ScalarType w=std::numeric_limits::max()) {plane=p; weight=w; } + + Plane *plane; + ScalarType weight; + }; + // Represent an edge in the MST tree + struct MSTEdge + { + MSTEdge(Plane *p0=NULL, Plane *p1=NULL, ScalarType w=std::numeric_limits::max()) {u=p0; v=p1; weight=w;}; + inline bool operator<(const MSTEdge &e) const {return weight sons; + }; + + /************************************************* + * The Algorithm + **************************************************/ + BoundingBoxType dataset_bb; + for (VertexIterator iter=begin; iter!=end; iter++) + dataset_bb.Add(iter->P()); + float max_distance = dataset_bb.Diag(); + + // Step 1: identify the tangent planes used to locally approximate the surface + int vertex_count = int( std::distance(begin, end) ); + std::vector< Plane > tangent_planes(vertex_count); + vcg::Octree< VertexType, ScalarType > octree_for_planes; + octree_for_planes.Set< VertexIterator , BoundingBoxForVertexFunctor >(begin, end, dataset_bb, BoundingBoxForVertexFunctor()); + + std::vector< VertexPointer > nearest_vertices; + std::vector< CoordType > nearest_points; + std::vector< ScalarType > distances; + for (VertexIterator iter=begin; iter!=end; iter++) + { + octree_for_planes.GetKClosest, std::vector, std::vector > + (VertPointDistanceFunctor(), DummyObjectMarker(), k, iter->P(), max_distance, nearest_vertices, distances, nearest_points); + + // for each vertex *iter, compute the centroid as avarege of the k-nearest vertices of *iter + Plane *plane = &tangent_planes[ std::distance(begin, iter) ]; + for (int n=0; ncenter += nearest_points[n]; + plane->center /= float(k); + + // then, identity the normal associated to the centroid + MatrixType covariance_matrix; + CoordType diff; + covariance_matrix.SetZero(); + for (int n=0; ncenter; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + covariance_matrix[i][j]+=diff[i]*diff[j]; + } + + CoordType eigenvalues; + MatrixType eigenvectors; + int required_rotations; + vcg::Jacobi< MatrixType, CoordType >(covariance_matrix, eigenvalues, eigenvectors, required_rotations); + vcg::SortEigenvaluesAndEigenvectors< MatrixType, CoordType >(eigenvalues, eigenvectors); + for (int d=0; d<3; d++) + plane->normal[d] = eigenvectors[d][2]; + plane->normal.Normalize(); + + plane->index = int( std::distance(begin, iter) ); + } + + // Step 2: build the Riemannian graph, i.e. the graph where each point is connected to the k-nearest neigbours. + dataset_bb.SetNull(); + std::vector< Plane >::iterator ePlane = tangent_planes.end(); + for (std::vector< Plane >::iterator iPlane=tangent_planes.begin(); iPlane!=ePlane; iPlane++) + dataset_bb.Add(iPlane->center); + max_distance = dataset_bb.Diag(); + + vcg::Octree< Plane, ScalarType > octree_for_plane; + octree_for_plane.Set< std::vector::iterator, BoundingBoxForPlaneFunctor >(tangent_planes.begin(), tangent_planes.end(), dataset_bb, BoundingBoxForPlaneFunctor()); + std::vector< Plane* > nearest_planes(distances.size()); + std::vector< std::vector< RiemannianEdge > > riemannian_graph(vertex_count); //it's probably that we are wasting the last position... + for (std::vector< Plane >::iterator iPlane=tangent_planes.begin(); iPlane!=ePlane; iPlane++) + { + octree_for_plane.GetKClosest< PlanePointDistanceFunctor, DummyObjectMarker, std::vector< Plane* >, std::vector< ScalarType >, std::vector< CoordType > > + (PlanePointDistanceFunctor(), DummyObjectMarker(), k, iPlane->center, max_distance, nearest_planes, distances, nearest_points, true, false); + + for (int n=0; nindexindex) + riemannian_graph[iPlane->index].push_back( RiemannianEdge( nearest_planes[n], 1.0f - fabs(iPlane->normal * nearest_planes[n]->normal)) ); + } + + // Step 3: compute the minimum spanning tree (MST) over the Riemannian graph (we use the Kruskal algorithm) + std::vector< MSTEdge > E; + std::vector< std::vector< RiemannianEdge > >::iterator iRiemannian = riemannian_graph.begin(); + std::vector< RiemannianEdge >::iterator iRiemannianEdge, eRiemannianEdge; + for (int i=0; ibegin(), eRiemannianEdge=iRiemannian->end(); iRiemannianEdge!=eRiemannianEdge; iRiemannianEdge++) + E.push_back(MSTEdge(&tangent_planes[i], iRiemannianEdge->plane, iRiemannianEdge->weight)); + + std::sort( E.begin(), E.end() ); + vcg::DisjointSet set; + + for (std::vector< Plane >::iterator iPlane=tangent_planes.begin(); iPlane!=ePlane; iPlane++) + set.MakeSet( &*iPlane ); + + std::vector< MSTEdge >::iterator iMSTEdge = E.begin(); + std::vector< MSTEdge >::iterator eMSTEdge = E.end(); + std::vector< MSTEdge > unoriented_tree; + Plane *u, *v; + for ( ; iMSTEdge!=eMSTEdge; iMSTEdge++) + if ((u=set.FindSet(iMSTEdge->u))!=(v=set.FindSet(iMSTEdge->v))) + unoriented_tree.push_back( *iMSTEdge ), set.Union(u, v); + E.clear(); + + // compute for each plane the list of sorting edges + std::vector< std::vector< int > > incident_edges(vertex_count); + iMSTEdge = unoriented_tree.begin(); + eMSTEdge = unoriented_tree.end(); + for ( ; iMSTEdge!=eMSTEdge; iMSTEdge++) + { + int u_index = int(iMSTEdge->u->index); + int v_index = int(iMSTEdge->v->index); + incident_edges[ u_index ].push_back( v_index ), + incident_edges[ v_index ].push_back( u_index ); + } + + // Traverse the incident_edges vector and build the MST + VertexIterator iCurrentVertex, iSonVertex; + std::vector< MSTNode > MST(vertex_count); + + std::vector< Plane >::iterator iFirstPlane = tangent_planes.begin(); + std::vector< Plane >::iterator iCurrentPlane, iSonPlane; + + MSTNode *mst_root; + int r_index = (root_index!=-1)? root_index : rand()*vertex_count/RAND_MAX; + mst_root = &MST[ r_index ]; + mst_root->parent = mst_root; //the parent of the root is the root itself + + if (orientation==MustBeFlipped) + { + iCurrentVertex = begin; + std::advance(iCurrentVertex, r_index); + iCurrentVertex->N() = iCurrentVertex->N()*ScalarType(-1.0f); + } + + { // just to limit the scope of the variable border + std::queue< int > border; + border.push(r_index); + while (!border.empty()) + { + int current_node_index = border.front(); border.pop(); + + MSTNode *current_node = &MST[current_node_index]; //retrieve the pointer to the current MST node + std::advance((iCurrentVertex=begin), current_node_index); //retrieve the pointer to the correspective vertex + current_node->vertex = &*iCurrentVertex; //and associate it to the MST node + + std::vector< int >::iterator iSon = incident_edges[ current_node_index ].begin(); + std::vector< int >::iterator eSon = incident_edges[ current_node_index ].end(); + for ( ; iSon!=eSon; iSon++) + { + MSTNode *son = &MST[ *iSon ]; + if (son->parent==NULL) // the node hasn't been visited + { + son->parent = current_node; // Update the MST nodes + current_node->sons.push_back(son); + //std::advance((iSonVertex=begin), *iSon);//retrieve the pointer to the Vertex associated to son + border.push( *iSon ); + } + } + } + } + + // and finally visit the MST tree in order to propagate the normals + { + std::queue< MSTNode* > border; + border.push(mst_root); + while (!border.empty()) + { + MSTNode *current_node = border.front(); border.pop(); + //std::vector< MSTNode* >::iterator iMSTSon = current_node->sons.begin(); + //std::vector< MSTNode* >::iterator eMSTSon = current_node->sons.end(); + for (int s=0; ssons.size()); s++) + { + if (current_node->vertex->N()*current_node->sons[s]->vertex->N()sons[s]->vertex->N() *= ScalarType(-1.0f); + border.push( current_node->sons[s] ); + } + } + } + }; + }; + +};//end of namespace vcg + +#endif //end of VCG_SPACE_NORMAL_EXTRAPOLATION_H \ No newline at end of file