From 90d868e36a304dba2e162e4581a5779ce7a67337 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20M=C3=B6bius?= Date: Tue, 20 Feb 2024 15:06:17 +0100 Subject: [PATCH] Initial Holefiller files from OpenFlipper --- src/OpenMesh/Tools/HoleFiller/HoleFillerT.hh | 200 +++++ .../Tools/HoleFiller/HoleFillerT_impl.hh | 784 ++++++++++++++++++ 2 files changed, 984 insertions(+) create mode 100644 src/OpenMesh/Tools/HoleFiller/HoleFillerT.hh create mode 100644 src/OpenMesh/Tools/HoleFiller/HoleFillerT_impl.hh diff --git a/src/OpenMesh/Tools/HoleFiller/HoleFillerT.hh b/src/OpenMesh/Tools/HoleFiller/HoleFillerT.hh new file mode 100644 index 00000000..c9cbcb78 --- /dev/null +++ b/src/OpenMesh/Tools/HoleFiller/HoleFillerT.hh @@ -0,0 +1,200 @@ +/*===========================================================================*\ +* * +* OpenFlipper * + * Copyright (c) 2001-2015, RWTH-Aachen University * + * Department of Computer Graphics and Multimedia * + * All rights reserved. * + * www.openflipper.org * + * * + *---------------------------------------------------------------------------* + * This file is part of OpenFlipper. * + *---------------------------------------------------------------------------* + * * + * 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. Neither the name of the copyright holder nor the names of its * + * contributors may 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 + +#include +#include +#include +#include +#include "OpenMeshUtils.hh" +#include +#include + +//============================================================================= + +template< class TheMesh > +class HoleFiller +{ +public: + typedef TheMesh Mesh; + + import_om_abbreviations( typename Mesh ); + + // A weight is a tuple of area and maximum dihedral angle + // + + class Weight { + public: + + Weight() : angle_( 180 ), area_( FLT_MAX ) {} + Weight( Scalar _angle, Scalar _area ) : angle_( _angle ), area_( _area ) {} + ~Weight() {} + + Scalar angle() const { return angle_; } + Scalar area() const { return area_; } + + Weight operator+( const Weight & _other ) const { + return Weight( std::max( angle(), _other.angle() ), + area() + _other.area() ); + } + + bool operator<( const Weight & _rhs ) const { + return ( angle() < _rhs.angle() || + ( angle() == _rhs.angle() && area() < _rhs.area() ) ); + } + + private: + Scalar angle_; + Scalar area_; + }; + + + // Ctors + explicit HoleFiller( Mesh & _mesh ); + ~HoleFiller(); + + // Identify and fill all holes of the mesh. + void fill_all_holes( int _stages = 3 ); + + + // Fill a hole which is identified by one of its boundary edges. + void fill_hole( EH _eh, int _stages = 3 ); + + // Fair a filling + //void fairing( std::vector< FH >& _faceHandles ); + void fairing( std::vector< OpenMesh::SmartFaceHandle >& _faceHandles ); + + // Remove degenerated faces from the filling + void removeDegeneratedFaces( std::vector< FH >& _faceHandles ); + +private: + + // Refine a face + bool refine( FH _fh ); + + // Relax an edge + bool relax_edge( OpenMesh::SmartEdgeHandle _eh ); + + // Test whether a point _x lies in the circumsphere of _a,_b,_c. + bool in_circumsphere( const Point & _x, + const Point & _a, + const Point & _b, + const Point & _c ) const; + + // Create the triangulation for polygon (_i,...,_j). + bool fill( int _i, int _j ); + + // Compute the weight of the triangle (_i,_j,_k). + Weight weight( int _i, int _j, int _k ); + + // Does edge (_u,_v) already exist? + bool exists_edge( OpenMesh::SmartVertexHandle _u, VH _w ); + + // Compute the area of the triangle (_a,_b,_c). + Scalar area( VH _a, VH _b, VH _c ); + + // Compute the dihedral angle (in degrees) between triangle + // (_u,_v,_a) and triangle (_v,_u,_b). + Scalar dihedral_angle( VH _u, VH _v, VH _a, VH _b ); + + + // The mesh, with each vertex we associate a scale factor that is + // needed for remeshing + + Mesh & mesh_; + OpenMesh::VPropHandleT< Scalar > scale_; + + /* + HOLE + + boundary_vertex_ + | + V + ==*=======*=======*== BOUNDARY + / \ / \ / \ + / \ / \ / \ + \ / \ / + * * <- opposite_vertex_ + */ + + + typedef std::vector< VH > VHVec; + typedef typename std::vector< VH >::iterator VHVecIter; + typedef typename std::vector< VH >::const_iterator CVHVecIter; + + typedef std::vector< FH > FHVec; + typedef typename std::vector< FH >::iterator FHVecIter; + typedef typename std::vector< FH >::const_iterator CFHVecIter; + + + // This vector contains all vertices of the hole (in order) + std::vector< OpenMesh::SmartVertexHandle > boundary_vertex_; + + // This vector contains all vertices that are opposite to an edge of the hole + VHVec opposite_vertex_; + + // This vector contains all edges of the hole (in order) + std::vector< OpenMesh::SmartEdgeHandle > hole_edge_; + + // This vector stores handles to all triangles of the current hole + std::vector< OpenMesh::SmartFaceHandle > hole_triangle_; + + // These are the two central arrays that are needed for the dynamic + // programming approach to hole filling. + // w_[i][j] : stores the minimal weight that can be achieved + // for a triangulation of the polygon + // boundary_vertex_[i],...,boundary_vertex_[j] + // l_[i][j] : stores the third index of the triangle + // + // that is needed for reconstructing the minimal triangulation + + std::vector< std::vector< Weight > > w_; + std::vector< std::vector< int > > l_; +}; + +//============================================================================= +#ifndef HOLEFILLER_CC + #include "HoleFillerT_impl.hh" +#endif +//============================================================================= + diff --git a/src/OpenMesh/Tools/HoleFiller/HoleFillerT_impl.hh b/src/OpenMesh/Tools/HoleFiller/HoleFillerT_impl.hh new file mode 100644 index 00000000..ead575ff --- /dev/null +++ b/src/OpenMesh/Tools/HoleFiller/HoleFillerT_impl.hh @@ -0,0 +1,784 @@ +/*===========================================================================*\ +* * +* OpenFlipper * + * Copyright (c) 2001-2015, RWTH-Aachen University * + * Department of Computer Graphics and Multimedia * + * All rights reserved. * + * www.openflipper.org * + * * + *---------------------------------------------------------------------------* + * This file is part of OpenFlipper. * + *---------------------------------------------------------------------------* + * * + * 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. Neither the name of the copyright holder nor the names of its * + * contributors may 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 "HoleFillerT.hh" +//============================================================================= + +template< class MeshT > +HoleFiller< MeshT >::HoleFiller( Mesh & _mesh ) + : mesh_( _mesh ) +{ + mesh_.request_vertex_status(); + mesh_.request_edge_status(); + + if (! mesh_.get_property_handle(scale_,"scale") ) + mesh_.add_property( scale_ , "scale" ); +} + + +//============================================================================= + + + +template< class MeshT > +HoleFiller< MeshT >::~HoleFiller() +{ + mesh_.release_vertex_status(); + mesh_.release_edge_status(); + + if ( mesh_.get_property_handle(scale_,"scale") ) + mesh_.remove_property( scale_ ); +} + + +//============================================================================= +// +// Identify and fill all holes of the mesh. +// +//============================================================================= + + +template< class MeshT > +void +HoleFiller< MeshT >::fill_all_holes( int _stages ) +{ + // Collect all boundary edges + + std::vector< EH > bdry_edge; + + for (auto ei : mesh_.edges()) + if ( ei.is_boundary() ) + bdry_edge.push_back( ei ); + + // Fill holes + + int cnt = 0; + for (auto i : bdry_edge) + if ( mesh_.is_boundary( i ) ) + { + ++cnt; + std::cerr << "Filling hole " << cnt << "\n"; + fill_hole( i, _stages ); + } + + + // Smooth fillings + + if ( _stages <= 2 ) + return; + + std::cerr << "Stage 3 : Smoothing the hole fillings ... "; + + OpenMesh::Smoother::JacobiLaplaceSmootherT< Mesh > smoother( mesh_ ); + smoother.initialize( OpenMesh::Smoother::SmootherT< Mesh >:: + Tangential_and_Normal, + OpenMesh::Smoother::SmootherT< Mesh >::C1 ); + + smoother.smooth( 500 ); + std::cerr << "ok\n"; +} + + +//============================================================================= +// +// Fill a hole which is identified by one of its boundary edges. +// +//============================================================================= + + +template< class MeshT > +void +HoleFiller< MeshT >::fill_hole( EH _eh, int _stages ) +{ + std::cerr << " Stage 1 : Computing a minimal triangulation ... "; + + //remember last vertex for selection of new ones + typename MeshT::VertexHandle old_last_handle = *(--mesh_.vertices_end()); + + // No boundary edge, no hole + + if ( ! mesh_.is_boundary( _eh ) ) + return; + + + // Get boundary halfedge + + OpenMesh::SmartHalfedgeHandle hh = make_smart(_eh, mesh_).h0(); + + if ( ! hh.is_boundary() ) + hh = hh.opp(); + + + // Collect boundary vertices + + boundary_vertex_.clear(); + opposite_vertex_.clear(); + + OpenMesh::SmartHalfedgeHandle ch = hh; + + do { + boundary_vertex_.push_back( ch.from() ); + opposite_vertex_.push_back( ch.opp().next().to() ); + //check number of outgoing boundary HEH's at Vertex + int c = 0; + OpenMesh::SmartVertexHandle vh = ch.to(); + + for (auto voh_it : vh.outgoing_halfedges()) + if ( voh_it.is_boundary() ) + c++; + + if ( c >= 2){ + OpenMesh::SmartHalfedgeHandle op = ch.opp(); + typename MeshT::VertexOHalfedgeIter voh_it(mesh_,op); + + ch = *(++voh_it); + }else + ch = ch.next(); + + } while ( ch != hh ); + + + int nv = boundary_vertex_.size(); + + // Compute an initial triangulation + + w_.clear(); + w_.resize( nv, std::vector( nv, Weight() ) ); + + l_.clear(); + l_.resize( nv, std::vector( nv, 0 ) ); + + + for ( int i = 0; i < nv - 1; ++i ) + w_[i][i+1] = Weight( 0, 0 ); + + for ( int j = 2; j < nv; ++j ) + { + #pragma omp parallel for shared(j, nv) + for(int i = 0; i < nv-j; ++i) + { + Weight valmin; + int argmin = -1; + for ( int m = i + 1; m < i + j; ++m ) + { + Weight newval = w_[i][m] + w_[m][i+j] + weight( i, m, i+j ); + if ( newval < valmin ) + { + valmin = newval; + argmin = m; + } + } + + w_[i][i+j] = valmin; + l_[i][i+j] = argmin; + } + } + + + // Actually fill the hole. We collect all triangles and edges of + // this filling for further processing. + + hole_edge_.clear(); + hole_triangle_.clear(); + if ( fill( 0, nv - 1 ) ){ + + std::cerr << "ok\n"; + + if ( _stages <= 1 ) + return; + + std::cerr << " Stage 2 : Fairing the filling ... "; + + std::vector< OpenMesh::SmartFaceHandle > handles = hole_triangle_; + + fairing(handles); + + //select all new vertices + typename MeshT::VertexIter old_end = ++typename MeshT::VertexIter(mesh_,old_last_handle); + typename MeshT::VertexIter v_end = mesh_.vertices_end(); + + for(; old_end != v_end; ++old_end) + if ( !mesh_.status(*old_end).deleted() ) + mesh_.status(*old_end).set_selected( true ); + + std::cerr << "ok\n"; + }else + std::cerr << "Could not create triangulation" << std::endl; +} + + +/// path fairing +template< class MeshT > +void +HoleFiller< MeshT >::fairing( std::vector< OpenMesh::SmartFaceHandle >& _faceHandles ){ + + //generate vector of all edges + hole_edge_.clear(); + + hole_triangle_ = _faceHandles; + + OpenMesh::EPropHandleT< bool > edgeProp; + OpenMesh::VPropHandleT< bool > vertexProp; + OpenMesh::FPropHandleT< bool > faceProp; + OpenMesh::FPropHandleT< bool > orderProp; + + if (! mesh_.get_property_handle(edgeProp,"edgeProp") ) + mesh_.add_property( edgeProp, "edgeProp" ); + if (! mesh_.get_property_handle(vertexProp,"vertexProp") ) + mesh_.add_property( vertexProp, "vertexProp" ); + if (! mesh_.get_property_handle(faceProp,"faceProp") ) + mesh_.add_property( faceProp, "faceProp" ); + if (! mesh_.get_property_handle(orderProp,"orderProp") ) + mesh_.add_property( orderProp, "orderProp" ); + + //init properties by setting all of them to false + for (auto fIt : mesh_.faces()) { + mesh_.property( orderProp, fIt ) = false; + mesh_.property( faceProp, fIt ) = false; + } + + for (auto eIt : mesh_.edges()) + mesh_.property( edgeProp, eIt ) = false; + + for (auto vIt : mesh_.vertices()) { + mesh_.property( vertexProp, vIt ) = false; + } + + //set face property + for (uint i = 0; i < hole_triangle_.size(); i++){ + mesh_.property( faceProp, hole_triangle_[i] ) = true; + } + + //set properties + for (unsigned int i = 0; i < hole_triangle_.size(); i++){ + for (auto fei : hole_triangle_[i].edges()) { + mesh_.status( fei ).set_locked(true); + //set edge property for all edges inside the hole (eg not on the hole boundary) + if (mesh_.property( faceProp, fei.h0().face() ) && + mesh_.property( faceProp, fei.h1().face() ) ){ + + mesh_.property( edgeProp, fei ) = true; + hole_edge_.push_back( fei ); + mesh_.status( fei ).set_locked(false); + } + } + + /// @TODO, strange iterator at property! + for (auto fvi : hole_triangle_[i].vertices()){ + //set vertex property for all vertices of the hole + for ( auto vfi : fvi.faces() ) + mesh_.property( vertexProp, fvi ) = true; + } + } + + //calculate scaling weights for vertices + for (auto vIt : mesh_.vertices()) + if (mesh_.property(vertexProp, vIt)){ + + Scalar cnt = 0; + Scalar scale = 0; + + for ( auto voh_it : vIt.outgoing_halfedges()) + { + + if (voh_it.face().is_valid() && + voh_it.opp().face().is_valid() && + mesh_.property(faceProp, voh_it.face() ) && + mesh_.property(faceProp, voh_it.opp().face() )) + continue; + + cnt += 1.0f; + Point p0 = mesh_.point( vIt ); + Point p1 = mesh_.point( voh_it.to() ); + scale += ( p1 - p0 ).norm(); + + } + + scale /= cnt; + + mesh_.property( scale_, vIt ) = scale; + } + + mesh_.remove_property(edgeProp); + mesh_.remove_property(vertexProp); + mesh_.remove_property(faceProp); + mesh_.remove_property(orderProp); + + // Do the patch fairing + + bool did_refine = true; + + for ( int k = 0; k < 40 && did_refine; ++k ) + { + uint end = hole_triangle_.size(); + + did_refine = false; + for ( unsigned int j = 0; j < end; ++j ) + did_refine |= refine( hole_triangle_[j] ); + + for ( int i = 0; i < 10; ++i ) + for ( unsigned int j = 0; j < hole_edge_.size(); ++j ) + relax_edge( hole_edge_[j] ); + } + + // unlock everything + for ( auto ei : mesh_.edges()) + mesh_.status( ei ).set_locked( false ); +} + +//============================================================================= +// +// Refine a face: Apply a 1-to-3 split if the edge lengths of the +// face do not match the interpolated edge lengths of the hole +// boundary. +// +//============================================================================= + + +template< class MeshT > +bool +HoleFiller< MeshT >::refine( FH _fh ) +{ + + // Collect the three edges of the face into e0, e1, e2 + + FEI fei = mesh_.fe_iter( _fh ); + OpenMesh::SmartEdgeHandle e0 = *fei; ++fei; + OpenMesh::SmartEdgeHandle e1 = *fei; ++fei; + OpenMesh::SmartEdgeHandle e2 = *fei; ++fei; + + + // Collect the vertices, vertex positions and scale factors of the face. + + FVI fvi = mesh_.fv_iter( _fh ); + + VH v0 = *fvi; ++fvi; + VH v1 = *fvi; ++fvi; + VH v2 = *fvi; ++fvi; + + Point p0 = mesh_.point( v0 ); + Point p1 = mesh_.point( v1 ); + Point p2 = mesh_.point( v2 ); + + Scalar scale0 = mesh_.property( scale_, v0 ); + Scalar scale1 = mesh_.property( scale_, v1 ); + Scalar scale2 = mesh_.property( scale_, v2 ); + + // Interpolate the scale factor. + + Scalar scale = ( scale0 + scale1 + scale2 ) / 3.0f; + Point center = typename MeshT::Scalar(1.0/3.0) * ( p0 + p1 + p2 ); + + Scalar d0 = 1.0f * ( p0 - center ).norm(); + Scalar d1 = 1.0f * ( p1 - center ).norm(); + Scalar d2 = 1.0f * ( p2 - center ).norm(); + + + //dont split triangles which tend to degenerate + if ( (d0 + d1 + d2) / 3.0f < scale) return false; + + + // If the edge lengths differ too much from the scale, perform a + // triangle split. + + if ( ( d0 > scale && d0 > scale0 ) || + ( d1 > scale && d1 > scale1 ) || + ( d2 > scale && d2 > scale2 ) ) + { + // Split the face ... + OpenMesh::SmartVertexHandle ch = mesh_.add_vertex( center ); + mesh_.split( _fh, ch ); + + // ... put the new triangles into the global triangle list ... + + for ( auto vfi : ch.faces() ) + if ( vfi != _fh ) + hole_triangle_.push_back( vfi ); + + // ... put the new edges into the global edge list ... + + for ( auto vei : ch.edges() ) + hole_edge_.push_back( vei ); + + // ... and set the appropriate scale factor for the new vertex. + + mesh_.property( scale_, ch ) = scale; + + relax_edge( e0 ); + relax_edge( e1 ); + relax_edge( e2 ); + + return true; + } + + return false; +} + + +//============================================================================= +// +// Relax an edge: Flip it if one of its opposing vertices lies in +// the circumsphere of the other three vertices. +// +//============================================================================= + + +template< class MeshT > +bool +HoleFiller< MeshT >::relax_edge( OpenMesh::SmartEdgeHandle _eh ) +{ + if ( mesh_.status( _eh ).locked() ) + return false; + + // Abbreviations for the two halfedges of _eh + + OpenMesh::SmartHalfedgeHandle h0 = _eh.h0(); + OpenMesh::SmartHalfedgeHandle h1 = _eh.h1(); + + // Get the two end-vertices u and v of the edge + + Point u( mesh_.point( h0.to() ) ); + Point v( mesh_.point( h1.to() ) ); + + // Get the two opposing vertices a and b + + Point a( mesh_.point( h0.next().to() ) ); + Point b( mesh_.point( h1.next().to() ) ); + + // If the circumsphere criterion is fullfilled AND if the flip is + // topologically admissible, we do it. + + if ( in_circumsphere( a, u, v, b ) || in_circumsphere( b, u, v, a ) ){ + if ( mesh_.is_flip_ok( _eh ) ) + { + + mesh_.flip( _eh ); + + return true; + }else + mesh_.status(_eh).set_selected( true ); +} + return false; +} + + +//============================================================================= +// +// Test whether a point _x lies in the circumsphere of _a,_b,_c. +// +//============================================================================= + + +template< class MeshT > +bool +HoleFiller< MeshT >::in_circumsphere( const Point & _x, + const Point & _a, + const Point & _b, + const Point & _c ) const +{ + Point ab = _b - _a; + Point ac = _c - _a; + + Scalar a00 = -2.0f * ( ab | _a ); + Scalar a01 = -2.0f * ( ab | _b ); + Scalar a02 = -2.0f * ( ab | _c ); + Scalar b0 = _a.sqrnorm() - _b.sqrnorm(); + + Scalar a10 = -2.0f * ( ac | _a ); + Scalar a11 = -2.0f * ( ac | _b ); + Scalar a12 = -2.0f * ( ac | _c ); + Scalar b1 = _a.sqrnorm() - _c.sqrnorm(); + + typename MeshT::Scalar alpha = -(-a11*a02+a01*a12-a12*b0+b1*a02+a11*b0-a01*b1) + / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12); + typename MeshT::Scalar beta = (a10*b0-a10*a02-a12*b0+a00*a12+b1*a02-a00*b1) + / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12); + typename MeshT::Scalar gamma = (-a11*a00-a10*b0+a00*b1+a11*b0+a01*a10-a01*b1) + / (-a11*a00+a11*a02-a10*a02+a00*a12+a01*a10-a01*a12); + + Point center = alpha * _a + beta * _b + gamma * _c; + + return ( _x - center ).sqrnorm() < ( _a - center ).sqrnorm(); +} + + +//============================================================================= +// +// Create the triangulation +// +// Recursively creates the triangulation for polygon (_i,...,_j). +// +//============================================================================= + + +template< class MeshT > +bool +HoleFiller< MeshT >::fill( int _i, int _j ) +{ + // If the two vertices _i and _j are adjacent, there is nothing to do. + + if ( _i + 1 == _j ) + return true; + + + // Create and store the middle triangle, store its edges. + + OpenMesh::SmartFaceHandle fh = mesh_.add_face( boundary_vertex_[_i], + boundary_vertex_[ l_[_i][_j] ], + boundary_vertex_[_j] ); + hole_triangle_.push_back( fh ); + + if (!fh.is_valid()) + return false; + + hole_edge_.push_back( mesh_.edge_handle + ( mesh_.find_halfedge( boundary_vertex_[_i], + boundary_vertex_[ l_[_i][_j] ] ) ) ); + hole_edge_.push_back( mesh_.edge_handle + ( mesh_.find_halfedge( boundary_vertex_[ l_[_i][_j] ], + boundary_vertex_[_j] ) ) ); + + + // Recursively create the left and right side of the + // triangulation. + + if (!fill( _i, l_[_i][_j] ) || !fill( l_[_i][_j], _j )) + return false; + else + return true; +} + + + +//============================================================================= +// +// Compute the weight of the triangle (_i,_j,_k). +// +//============================================================================= + + +template< class MeshT > +typename HoleFiller< MeshT >::Weight +HoleFiller< MeshT >::weight( int _i, int _j, int _k ) +{ + // Return an infinite weight if the insertion of this triangle + // would create complex edges. + + if ( exists_edge( boundary_vertex_[_i], boundary_vertex_[_j] ) || + exists_edge( boundary_vertex_[_j], boundary_vertex_[_k] ) || + exists_edge( boundary_vertex_[_k], boundary_vertex_[_i] ) ) + return Weight(); + + + // Return an infinite weight, if one of the neighboring patches + // could not be created. + + + if ( l_[_i][_j] == -1 ) return Weight(); + if ( l_[_j][_k] == -1 ) return Weight(); + + + // Compute the maxmimum dihedral angles to the adjacent triangles + // (if they exist) + + Scalar angle = 0.0f; + + if ( _i + 1 == _j ) + angle = std::max( angle, dihedral_angle( boundary_vertex_[_i], + boundary_vertex_[_j], + boundary_vertex_[_k], + opposite_vertex_[_i] ) ); + else + angle = std::max( angle, dihedral_angle( boundary_vertex_[_i], + boundary_vertex_[_j], + boundary_vertex_[_k], + boundary_vertex_[l_[_i][_j]] ) ); + + if ( _j + 1 == _k ) + angle = std::max( angle, dihedral_angle( boundary_vertex_[_j], + boundary_vertex_[_k], + boundary_vertex_[_i], + opposite_vertex_[_j] ) ); + else + angle = std::max( angle, dihedral_angle( boundary_vertex_[_j], + boundary_vertex_[_k], + boundary_vertex_[_i], + boundary_vertex_[l_[_j][_k]] ) ); + + if ( _i == 0 && _k == (int) boundary_vertex_.size() - 1 ) + angle = std::max( angle, dihedral_angle( boundary_vertex_[_k], + boundary_vertex_[_i], + boundary_vertex_[_j], + opposite_vertex_[_k] ) ); + + + return Weight( angle, area( boundary_vertex_[_i], + boundary_vertex_[_j], + boundary_vertex_[_k] ) ); +} + + + +//============================================================================= +// +// Does an edge from vertex _u to _v exist? +// +//============================================================================= + + +template< class MeshT > +bool +HoleFiller< MeshT >::exists_edge( OpenMesh::SmartVertexHandle _u, VH _w ) +{ + for ( auto vohi : _u.outgoing_halfedges() ) + if ( ! vohi.edge().is_boundary() ) + if ( vohi.to() == _w ) + return true; + + return false; +} + + + +//============================================================================= +// +// Compute the area of the triangle (_a,_b,_c). +// +//============================================================================= + + +template< class MeshT > +typename MeshT::Scalar +HoleFiller< MeshT >::area( VH _a, VH _b, VH _c ) +{ + Point a( mesh_.point( _a ) ); + Point b( mesh_.point( _b ) ); + Point c( mesh_.point( _c ) ); + + Point n( ( b - a ) % ( c - b ) ); + + return 0.5 * n.norm(); +} + + +//============================================================================= +// +// Compute a dihedral angle +// +// Computes the dihedral angle (in degrees) between triangle +// (_u,_v,_a) and triangle (_v,_u,_b), no matter whether these +// triangles actually exist in the current mesh or not). +// +//============================================================================= + + +template< class MeshT > +typename MeshT::Scalar +HoleFiller< MeshT >::dihedral_angle( VH _u, VH _v, VH _a, VH _b ) +{ + Point u( mesh_.point( _u ) ); + Point v( mesh_.point( _v ) ); + Point a( mesh_.point( _a ) ); + Point b( mesh_.point( _b ) ); + + Point n0( ( v - u ) % ( a - v ) ); + Point n1( ( u - v ) % ( b - u ) ); + + n0.normalize(); + n1.normalize(); + + return acos( n0 | n1 ) * 180.0 / M_PI; +} + + +/// remove degenerated faces +template< class MeshT > +void +HoleFiller< MeshT >::removeDegeneratedFaces( std::vector< FH >& _faceHandles ){ + + for (int i = _faceHandles.size()-1; i >= 0 ; i--){ + + if ( mesh_.status( _faceHandles[i] ).deleted() ){ + // face might be deleted because of a previous edge collapse + // erase the face from the vector + _faceHandles.erase( _faceHandles.begin() + i ); + continue; + } + + //get the vertices (works only on triMeshes) + FVI fvi = mesh_.fv_iter( _faceHandles[i] ); + Point v0 = mesh_.point( *fvi); + ++fvi; + Point v1 = mesh_.point( *fvi ); + ++fvi; + Point v2 = mesh_.point( *fvi ); + + //check if its degenerated + Point v0v1 = v1 - v0; + Point v0v2 = v2 - v0; + Point n = v0v1 % v0v2; // not normalized ! + double d = n.sqrnorm(); + + if (d < FLT_MIN && d > -FLT_MIN) { + // degenerated face found + FHI hIt = mesh_.fh_iter( _faceHandles[i] ); + + //try to collapse one of the edges + while (hIt.is_valid()){ + if ( mesh_.is_collapse_ok( *hIt ) ){ + // collapse the edge to remove the triangle + mesh_.collapse( *hIt ); + // and erase the corresponding face from the vector + _faceHandles.erase( _faceHandles.begin() + i ); + break; + } else { + ++hIt; + } + } + } + } +} + +//============================================================================= +