diff --git a/.gitattributes b/.gitattributes index 8eeb2c30046..394bc9618d6 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1900,7 +1900,7 @@ Surface_mesh_parameterization/test/Surface_mesh_parameterization/data/cube.off - Surface_mesh_parameterization/test/Surface_mesh_parameterization/data/high_genus.off -text svneol=unset#application/octet-stream Surface_mesh_parameterization/test/Surface_mesh_parameterization/data/knot2.off -text svneol=unset#application/octet-stream Surface_mesh_parameterization/test/Surface_mesh_parameterization/data/oni.off -text svneol=unset#application/octet-stream -Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse_functor_base.h -text +Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_collapse_functor_base.h -text Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Polyhedron_edge_cached_pointer_map.h -text Surface_mesh_simplification/test/Surface_mesh_simplification/Midpoint_edge_collapse_test.cpp -text Surface_mesh_simplification/test/Surface_mesh_simplification/edge_collapse_test.cpp -text diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Count_ratio_stop_pred.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_pred.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Count_ratio_stop_pred.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_ratio_stop_pred.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Count_stop_pred.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_pred.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Count_stop_pred.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Count_stop_pred.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h new file mode 100644 index 00000000000..bd7be108078 --- /dev/null +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core.h @@ -0,0 +1,231 @@ +// Copyright (c) 2005, 2006 Fernando Luis Cacciola Carballal. All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Fernando Cacciola +// +#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H +#define CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H 1 + +#include + +#include +#include + +CGAL_BEGIN_NAMESPACE + +// +// This should be in +// +// Implementation of the collapsing cost and placement strategy from: +// +// "Fast and Memory Efficient Polygonal Symplification" +// Peter Lindstrom, Greg Turk +// + +namespace Triangulated_surface_mesh { namespace Simplification +{ + +template +class LindstromTurkCore +{ +public: + + typedef Collapse_data_ Collapse_data ; + + typedef typename Collapse_data::TSM TSM ; + typedef typename Collapse_data::vertex_descriptor vertex_descriptor ; + typedef typename Collapse_data::edge_descriptor edge_descriptor ; + typedef typename Collapse_data::Params Params ; + + typedef boost::graph_traits GraphTraits ; + + typedef typename GraphTraits::in_edge_iterator in_edge_iterator ; + + typedef typename Surface_geometric_traits::Kernel Kernel ; + + typedef typename Kernel::Point_3 Point ; + typedef typename Kernel::Vector_3 Vector ; + typedef typename Kernel::FT FT ; + + typedef optional Optional_FT ; + typedef optional Optional_point ; + typedef optional Optional_vector ; + + typedef MatrixC33 Matrix ; + +public: + + LindstromTurkCore( Params const& aParams + , vertex_descriptor const& aP + , vertex_descriptor const& aQ + , bool aIsPFixed + , bool aIsQFixed + , edge_descriptor const& aP_Q + , edge_descriptor const& aQ_P + , TSM& aSurface + ) ; + + void compute( Collapse_data& aData ) ; + +private : + + struct Triangle + { + Triangle() {} + + Triangle( Vector const& aNormalV, FT const& aNormalL ) : NormalV(aNormalV), NormalL(aNormalL) {} + + Vector NormalV ; + FT NormalL ; + } ; + + typedef std::vector Triangles ; + typedef std::vector Link ; + typedef std::vector edge_descriptor_vector ; + + struct BoundaryEdge + { + BoundaryEdge ( Point s_, Point t_, Vector const& v_, Vector const& n_ ) : s(s_), t(t_), v(v_), n(n_) {} + + Point s, t ; + Vector v, n ; + } ; + typedef std::vector BoundaryEdges ; + + class Constrians + { + public: + + Constrians() : n(0), A(NULL_MATRIX), b(NULL_VECTOR) {} + + void Add_if_alpha_compatible( Vector const& Ai, FT const& bi ) ; + + void Add_from_gradient ( Matrix const& H, Vector const& c ) ; + + int n ; + Matrix A ; + Vector b ; + + private: + + // alpha = 1 degree + static double squared_cos_alpha() { return 0.999695413509 ; } + static double squared_sin_alpha() { return 3.04586490453e-4; } + } ; + +private : + + void Add_boundary_preservation_constrians( BoundaryEdges const& aBdry ) ; + void Add_volume_preservation_constrians( Triangles const& aTriangles ); + void Add_boundary_and_volume_optimization_constrians( BoundaryEdges const& aBdry, Triangles const& aTriangles ) ; + void Add_shape_optimization_constrians( Link const& aLink ) ; + + FT Compute_boundary_cost( Vector const& v, BoundaryEdges const& aBdry ) ; + FT Compute_volume_cost ( Vector const& v, Triangles const& aTriangles ) ; + FT Compute_shape_cost ( Point const& p, Link const& aLink ) ; + + bool is_border ( edge_descriptor const& edge ) const + { + edge_is_border_t is_border_property ; + return get(is_border_property,mSurface,edge) ; + } + bool is_undirected_edge_a_border ( edge_descriptor const& edge ) const + { + return is_border(edge) || is_border(opposite_edge(edge,mSurface)) ; + } + + Point get_point( vertex_descriptor const& v ) const + { + vertex_point_t vertex_point_property ; + return get(vertex_point_property,mSurface,v) ; + } + + + static Vector Point_cross_product ( Point const& a, Point const& b ) + { + return cross_product(a-ORIGIN,b-ORIGIN); + } + + // This is the (uX)(Xu) product described in the Lindstrom-Turk paper + static Matrix LT_product( Vector const& u ) + { + FT a00 = ( u.y()*u.y() ) + ( u.z()*u.z() ) ; + FT a01 = -u.x()*u.y(); + FT a02 = -u.x()*u.z(); + + FT a10 = a01 ; + FT a11 = ( u.x()*u.x() ) + ( u.z()*u.z() ) ; + FT a12 = - u.y() * u.z(); + + FT a20 = a02 ; + FT a21 = a12 ; + FT a22 = ( u.x()*u.x() ) + ( u.y()*u.y() ) ; + + return Matrix(a00,a01,a02 + ,a10,a11,a12 + ,a20,a21,a22 + ); + } + + Triangle Get_triangle ( vertex_descriptor const& v0 + , vertex_descriptor const& v1 + , vertex_descriptor const& v2 + ) ; + + void Extract_triangle( vertex_descriptor const& v0 + , vertex_descriptor const& v1 + , vertex_descriptor const& v2 + , edge_descriptor const& e02 + , Triangles& rTriangles + ) ; + + void Extract_triangles_and_link( Triangles& rTriangles, Link& rLink ); + + void Extract_boundary_edge ( edge_descriptor edge, BoundaryEdges& rBdry ) ; + void Extract_boundary_edges( vertex_descriptor const& v + , edge_descriptor_vector& rCollected + , BoundaryEdges& rBdry + ) ; + void Extract_boundary_edges( BoundaryEdges& rBdry ) ; + + +private: + + Params const& mParams ; + vertex_descriptor const& mP ; + vertex_descriptor const& mQ ; + bool mIsPFixed ; + bool mIsQFixed ; + edge_descriptor const& mP_Q ; + edge_descriptor const& mQ_P ; + TSM& mSurface ; + +private: + + Constrians mConstrians ; + + Point mV ; + +}; + +} } // namespace Triangulated_surface_mesh::Simplification + +CGAL_END_NAMESPACE + +#include + +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H // +// EOF // + diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h new file mode 100644 index 00000000000..00c279832c1 --- /dev/null +++ b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h @@ -0,0 +1,700 @@ +// Copyright (c) 2005, 2006 Fernando Luis Cacciola Carballal. All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Fernando Cacciola +// +#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_IMPL_H +#define CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_IMPL_H 1 + +CGAL_BEGIN_NAMESPACE + +// +// Implementation of the strategy from: +// +// "Fast and Memory Efficient Polygonal Symplification" +// Peter Lindstrom, Greg Turk +// + +namespace Triangulated_surface_mesh { namespace Simplification +{ + +template +LindstromTurkCore::LindstromTurkCore( Params const& aParams + , vertex_descriptor const& aP + , vertex_descriptor const& aQ + , bool aIsPFixed + , bool aIsQFixed + , edge_descriptor const& aP_Q + , edge_descriptor const& aQ_P + , TSM& aSurface + ) + : + mParams(aParams) + ,mP(aP) + ,mQ(aQ) + ,mIsPFixed(aIsPFixed) + ,mIsQFixed(aIsQFixed) + ,mP_Q(aP_Q) + ,mQ_P(aQ_P) + ,mSurface(aSurface) +{ +} + +template +void LindstromTurkCore::compute( Collapse_data& rData ) +{ + + Optional_FT lOptionalCost ; + Optional_point lOptionalP ; + + if ( !mIsPFixed || !mIsQFixed ) + { + CGAL_TSMS_LT_TRACE(2,"Computing LT data for E" << mP_Q->ID << " (V" << mP->ID << "->V" << mQ->ID << ")" ); + + // Volume preservation and optimization constrians are based on the normals to the triangles in the star of the collapsing egde + // Triangle shape optimization constrians are based on the link of the collapsing edge (the cycle of vertices around the edge) + Triangles lTriangles; + Link lLink; + + lTriangles.reserve(16); + lLink .reserve(16); + + Extract_triangles_and_link(lTriangles,lLink); + +#ifdef CGAL_SURFACE_SIMPLIFICATION_ENABLE_LT_TRACE + std::ostringstream ss ; + for( typename Link::const_iterator it = lLink.begin(), eit = lLink.end() ; it != eit ; ++it ) + ss << "v" << (*it)->ID << " " ; + std::string s = ss.str(); + CGAL_TSMS_LT_TRACE(3,"Link: " << s ); +#endif + + BoundaryEdges lBdry ; + + Extract_boundary_edges(lBdry); + + Point const& lP = get_point(mP) ; + Point const& lQ = get_point(mQ) ; + + Optional_vector lOptionalV ; + + if ( mIsPFixed ) + { + lOptionalV = Optional_vector(lP - ORIGIN); + } + else if ( mIsQFixed ) + { + lOptionalV = Optional_vector(lQ - ORIGIN); + } + else + { + // + // Each vertex constrian is an equation of the form: Ai * v = bi + // Where 'v' is a vector representing the vertex, + // 'Ai' is a (row) vector and 'bi' a scalar. + // + // The vertex is completely determined with 3 such constrians, + // so is the solution to the folloing system: + // + // A.r0(). * v = b0 + // A1 * v = b1 + // A2 * v = b2 + // + // Which in matrix form is : A * v = b + // + // (with 'A' a 3x3 matrix and 'b' a vector) + // + // The member variable mConstrinas contains A and b. Indidivual constrians (Ai,bi) can be added to it. + // Once 3 such constrians have been added 'v' is directly solved a: + // + // v = b*inverse(A) + // + // A constrian (Ai,bi) must be alpha-compatible with the previously added constrians (see Paper); if it's not, is discarded. + // + if ( lBdry.size() > 0 ) + Add_boundary_preservation_constrians(lBdry); + + if ( mConstrians.n < 3 ) + Add_volume_preservation_constrians(lTriangles); + + if ( mConstrians.n < 3 ) + Add_boundary_and_volume_optimization_constrians(lBdry,lTriangles); + + if ( mConstrians.n < 3 ) + Add_shape_optimization_constrians(lLink); + + // It might happen that there were not enough alpha-compatible constrians. + // In that case there is simply no good vertex placement + if ( mConstrians.n == 3 ) + { + // If the matrix is singular it's inverse cannot be computed so an 'absent' value is returned. + optional lOptional_Ai = inverse_matrix(mConstrians.A); + if ( lOptional_Ai ) + { + Matrix const& lAi = *lOptional_Ai ; + + lOptionalV = mConstrians.b * lAi ; + + } + else + CGAL_TSMS_LT_TRACE(1,"Can't solve optimization, singular system."); + } + else + CGAL_TSMS_LT_TRACE(1,"Can't solve optimization, not enough alpha-compatible constrians."); + } + + if ( lOptionalV ) + { + // + // lP is the optimized new vertex position. Now we compute the collapsing cost: + // + lOptionalP = Optional_point( ORIGIN + (*lOptionalV) ); + + FT lSquaredLength = squared_distance(lP,lQ); + + CGAL_TSMS_LT_TRACE(1,"Squared edge length: " << lSquaredLength ) ; + + FT lBdryCost = Compute_boundary_cost(*lOptionalV,lBdry); + FT lVolumeCost = Compute_volume_cost (*lOptionalV,lTriangles); + FT lShapeCost = Compute_shape_cost (*lOptionalP,lLink); + + FT lTotalCost = mParams.VolumeWeight * lVolumeCost + + mParams.BoundaryWeight * lBdryCost * lSquaredLength + + mParams.ShapeWeight * lShapeCost * lSquaredLength * lSquaredLength ; + + lOptionalCost = Optional_FT(lTotalCost); + + CGAL_TSMS_LT_TRACE(1,"New vertex point: " << xyz_to_string(*lOptionalP) << "\nTotal cost: " << lTotalCost ); + + CGAL_TSMS_LT_TRACE(2,"\nSquared edge length: " << lSquaredLength + << "\nBoundary cost: " << lBdryCost + << "\nVolume cost: " << lVolumeCost + << "\nShape cost: " << lShapeCost + << "\nTOTAL COST: " << lTotalCost + ); + } + + } + else + CGAL_TSMS_LT_TRACE(1,"The edge is a fixed edge."); + + + rData = Collapse_data(mP,mQ,mIsPFixed,mIsQFixed,mP_Q,mSurface,lOptionalCost,lOptionalP) ; +} + + +// +// Caches the "local boundary", that is, the sequence of 3 border edges: o->p, p->q, q->e +// +template +void LindstromTurkCore::Extract_boundary_edge( edge_descriptor edge, BoundaryEdges& rBdry ) +{ + edge_descriptor face_edge = is_border(edge) ? opposite_edge(edge,mSurface) : edge ; + + vertex_descriptor sv = source(face_edge,mSurface); + vertex_descriptor tv = target(face_edge,mSurface); + + Point const& sp = get_point(sv); + Point const& tp = get_point(tv); + + Vector v = tp - sp ; + Vector n = Point_cross_product(tp,sp) ; + + CGAL_TSMS_LT_TRACE(3,"Boundary edge. S:" << xyz_to_string(sp) << " T:" << xyz_to_string(tp) + << " V:" << xyz_to_string(v) << " N:" << xyz_to_string(n) + ) ; + + rBdry.push_back( BoundaryEdge(sp,tp,v,n) ) ; + +} + +template +void LindstromTurkCore::Extract_boundary_edges( vertex_descriptor const& v + , edge_descriptor_vector& rCollected + , BoundaryEdges& rBdry + ) +{ + in_edge_iterator eb, ee ; + for ( tie(eb,ee) = in_edges(v,mSurface) ; eb != ee ; ++ eb ) + { + edge_descriptor edge = *eb ; + + if ( is_undirected_edge_a_border(edge) && std::find(rCollected.begin(),rCollected.end(),edge) == rCollected.end() ) + { + Extract_boundary_edge(edge,rBdry); + rCollected.push_back(edge); + rCollected.push_back(opposite_edge(edge,mSurface)); + } + } +} + +template +void LindstromTurkCore::Extract_boundary_edges( BoundaryEdges& rBdry ) +{ + edge_descriptor_vector lCollected ; + Extract_boundary_edges(mP,lCollected,rBdry); + Extract_boundary_edges(mQ,lCollected,rBdry); +} + +// +// Calculates the normal of the triangle (v0,v1,v2) (both vector and its length as (v0xv1).v2) +// +template +typename LindstromTurkCore::Triangle LindstromTurkCore::Get_triangle( vertex_descriptor const& v0 + , vertex_descriptor const& v1 + , vertex_descriptor const& v2 + ) +{ + Point const& p0 = get_point(v0); + Point const& p1 = get_point(v1); + Point const& p2 = get_point(v2); + + Vector v01 = p1 - p0 ; + Vector v02 = p2 - p0 ; + + Vector lNormalV = cross_product(v01,v02); + + FT lNormalL = Point_cross_product(p0,p1) * (p2-ORIGIN); + + CGAL_TSMS_LT_TRACE(3,"Extracting triangle v" << v0->ID << "->v" << v1->ID << "->v" << v2->ID + << " N:" << xyz_to_string(lNormalV) << " L:" << lNormalL + ); + + return Triangle(lNormalV,lNormalL); +} + +// +// If (v0,v1,v2) is a finite triangular facet of the mesh, that is, NONE of these vertices are boundary vertices, +// the triangle (properly oriented) is added to rTriangles. +// The triangle is encoded as its normal, calculated using the actual facet orientation [(v0,v1,v2) or (v0,v2,v1)] +// +template +void LindstromTurkCore::Extract_triangle( vertex_descriptor const& v0 + , vertex_descriptor const& v1 + , vertex_descriptor const& v2 + , edge_descriptor const& e02 + , Triangles& rTriangles + ) +{ + // The 3 vertices are obtained by circulating ccw around v0, that is, e02 = next_ccw(e01). + // Since these vertices are NOT obtained by circulating the face, the actual triangle orientation is unspecified. + + // The triangle is oriented v0->v2->v1 if the next edge that follows e02 (which is the edge v0->v2) is v2->v1. + if ( target(next_edge(e02,mSurface),mSurface) == v1 ) + { + // The triangle is oriented v0->v2->v1. + // In this case e02 is an edge of the facet. + // If this facet edge is a border edge then this triangle is not in the mesh . + if ( !is_border(e02) ) + rTriangles.push_back(Get_triangle(v0,v2,v1) ) ; + } + else + { + // The triangle is oriented v0->v1->v2. + // In this case, e20 and not e02, is an edge of the facet. + // If this facet edge is a border edge then this triangle is not in the mesh . + if ( !is_border(opposite_edge(e02,mSurface)) ) + rTriangles.push_back(Get_triangle(v0,v1,v2) ) ; + } +} + +// +// Extract all triangles (its normals) and vertices (the link) around the collpasing edge p_q +// +template +void LindstromTurkCore::Extract_triangles_and_link( Triangles& rTriangles, Link& rLink ) +{ + // + // Extract around mP, CCW + // + vertex_descriptor v0 = mP; + vertex_descriptor v1 = mQ; + + edge_descriptor e02 = mP_Q; + + do + { + e02 = next_edge_ccw(e02,mSurface); + + vertex_descriptor v2 = target(e02,mSurface); + + if ( v2 != mQ ) + { + CGAL_expensive_assertion ( std::find(rLink.begin(),rLink.end(),v2) == rLink.end() ) ; + rLink.push_back(v2) ; + } + + Extract_triangle(v0,v1,v2,e02,rTriangles); + + v1 = v2 ; + } + while ( e02 != mP_Q ) ; + + // + // Extract around mQ, CCW + // + + v0 = mQ; + + e02 = next_edge_ccw(mQ_P,mSurface); + + v1 = target(e02,mSurface); + + // This could have been added to the link while circulating around mP + if ( v1 != mP && std::find(rLink.begin(),rLink.end(),v1) == rLink.end() ) + rLink.push_back(v1) ; + + e02 = next_edge_ccw(e02,mSurface); + + do + { + vertex_descriptor v2 = target(e02,mSurface); + + // Any of the vertices found around mP can be reached again around mQ, but we can't duplicate them here. + if ( v2 != mP && std::find(rLink.begin(),rLink.end(),v2) == rLink.end() ) + rLink.push_back(v2) ; + + Extract_triangle(v0,v1,v2,e02,rTriangles); + + v1 = v2 ; + + e02 = next_edge_ccw(e02,mSurface); + + } + while ( e02 != mQ_P ) ; +} + +template +void LindstromTurkCore::Add_boundary_preservation_constrians( BoundaryEdges const& aBdry ) +{ + + if ( aBdry.size() > 0 ) + { + Vector e1 = NULL_VECTOR ; + Vector e2 = NULL_VECTOR ; + + for ( typename BoundaryEdges::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it ) + { + e1 = e1 + it->v ; + e2 = e2 + it->n ; + } + + CGAL_TSMS_LT_TRACE(2,"Adding boundary preservation constrians. SumV=" << xyz_to_string(e1) << " SumN=" << xyz_to_string(e2) ); + + Matrix H = LT_product(e1); + + Vector c = cross_product(e1,e2); + + mConstrians.Add_from_gradient(H,c); + } +} + +template +void LindstromTurkCore::Add_volume_preservation_constrians( Triangles const& aTriangles ) +{ + CGAL_TSMS_LT_TRACE(2,"Adding volume preservation constrians. " << aTriangles.size() << " triangles."); + + Vector lSumV = NULL_VECTOR ; + FT lSumL(0) ; + + for( typename Triangles::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it ) + { + lSumV = lSumV + it->NormalV ; + lSumL = lSumL + it->NormalL ; + } + + mConstrians.Add_if_alpha_compatible(lSumV,lSumL); + +} + +template +void LindstromTurkCore::Add_boundary_and_volume_optimization_constrians( BoundaryEdges const& aBdry, Triangles const& aTriangles ) +{ + CGAL_TSMS_LT_TRACE(2,"Adding boundary and volume optimization constrians. "); + + Matrix H = NULL_MATRIX ; + Vector c = NULL_VECTOR ; + + // + // Volume optimization + // + for( typename Triangles::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it ) + { + Triangle const& lTri = *it ; + + H += direct_product(lTri.NormalV,lTri.NormalV) ; + + c = c - ( lTri.NormalL * lTri.NormalV ) ; + } + + CGAL_TSMS_LT_TRACE(3,"Hv:" << matrix_to_string(H) << "\n cv:" << xyz_to_string(c) ) ; + + + if ( aBdry.size() > 0 ) + { + // + // Boundary optimization + // + Matrix Hb = NULL_MATRIX ; + Vector cb = NULL_VECTOR ; + + for ( typename BoundaryEdges::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it ) + { + Matrix H = LT_product(it->v); + Vector c = cross_product(it->v,it->n); + + Hb += H ; + cb = cb + c ; + } + + CGAL_TSMS_LT_TRACE(3,"Hb:" << matrix_to_string(Hb) << "\n cb:" << xyz_to_string(cb) ) ; + + // + // Weighted average + // + FT lScaledBoundaryWeight = FT(9) * mParams.BoundaryWeight * squared_distance ( get_point(mP), get_point(mQ) ) ; + + H *= mParams.VolumeWeight ; + c = c * mParams.VolumeWeight ; + + H += lScaledBoundaryWeight * Hb ; + c = c + ( lScaledBoundaryWeight * cb ) ; + + CGAL_TSMS_LT_TRACE(3,"VolW=" << mParams.VolumeWeight << " BdryW=" << mParams.BoundaryWeight << " ScaledBdryW=" << lScaledBoundaryWeight ) ; + + } + + mConstrians.Add_from_gradient(H,c); +} + +template +void LindstromTurkCore::Add_shape_optimization_constrians( Link const& aLink ) +{ + FT s(aLink.size()); + + Matrix H (s,0,0 + ,0,s,0 + ,0,0,s + ); + + Vector c = NULL_VECTOR ; + + for( typename Link::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it ) + c = c + (ORIGIN - get_point(*it)) ; + + CGAL_TSMS_LT_TRACE(2,"Adding shape optimization constrians: Shape vector: " << xyz_to_string(c) ); + + mConstrians.Add_from_gradient(H,c); +} + +template +typename LindstromTurkCore::FT +LindstromTurkCore::Compute_boundary_cost( Vector const& v, BoundaryEdges const& aBdry ) +{ + FT rCost(0); + for ( typename BoundaryEdges::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it ) + { + Vector u = (it->t - ORIGIN ) - v ; + Vector c = cross_product(it->v,u); + rCost += c*c; + } + return rCost / FT(4) ; +} + +template +typename LindstromTurkCore::FT +LindstromTurkCore::Compute_volume_cost( Vector const& v, Triangles const& aTriangles ) +{ + FT rCost(0); + + for( typename Triangles::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it ) + { + Triangle const& lTri = *it ; + + FT lF = lTri.NormalV * v - lTri.NormalL ; + + rCost += ( lF * lF ) ; + + } + + return rCost / FT(36) ; +} + +template +typename LindstromTurkCore::FT +LindstromTurkCore::Compute_shape_cost( Point const& p, Link const& aLink ) +{ + FT rCost(0); + + for( typename Link::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it ) + rCost += squared_distance(p,get_point(*it)) ; + + return rCost ; +} + +template +void LindstromTurkCore::Constrians::Add_if_alpha_compatible( Vector const& Ai, FT const& bi ) +{ + double slai = to_double(Ai*Ai) ; + if ( slai > 0.0 ) + { + double l = CGAL_NTS sqrt(slai) ; + + Vector Ain = Ai / l ; + FT bin = bi / l ; + + bool lAddIt = true ; + + if ( n == 1 ) + { + FT d01 = A.r0() * Ai ; + + double sla0 = to_double(A.r0() * A.r0()) ; + double sd01 = to_double(d01 * d01) ; + + if ( sd01 > ( sla0 * slai * squared_cos_alpha() ) ) + lAddIt = false ; + } + else if ( n == 2 ) + { + Vector N = cross_product(A.r0(),A.r1()); + + FT dc012 = N * Ai ; + + double slc01 = to_double(N * N) ; + double sdc012 = to_double(dc012 * dc012); + + if ( sdc012 <= slc01 * slai * squared_sin_alpha() ) + lAddIt = false ; + } + + if ( lAddIt ) + { + switch ( n ) + { + case 0 : + A.r0() = Ain ; + b = Vector(bin,b.y(),b.z()); + break ; + case 1 : + A.r1() = Ain ; + b = Vector(b.x(),bin,b.z()); + break ; + case 2 : + A.r2() = Ain ; + b = Vector(b.x(),b.y(),bin); + break ; + } + ++ n ; + + CGAL_TSMS_LT_TRACE(1,"Constrains.A:" << matrix_to_string(A) << "\nConstrains.b:" << xyz_to_string(b) ) ; + } + } +} + +template +int index_of_max_component ( V const& v ) +{ + typedef typename Kernel_traits::Kernel::FT FT ; + + int i = 0 ; + FT max = v.x(); + if ( max < v.y() ) + { + max = v.y(); + i = 1 ; + } + if ( max < v.z() ) + { + max = v.z(); + i = 2 ; + } + return i ; +} + +template +void LindstromTurkCore::Constrians::Add_from_gradient ( Matrix const& H, Vector const& c ) +{ + CGAL_precondition(n >= 0 && n<=2 ); + + switch(n) + { + case 0 : + + Add_if_alpha_compatible(H.r0(),-c.x()); + Add_if_alpha_compatible(H.r1(),-c.y()); + Add_if_alpha_compatible(H.r2(),-c.z()); + + break; + + case 1 : + { + Vector const& A0 = A.r0(); + + Vector A02( A0.x()*A0.x() + , A0.y()*A0.y() + , A0.z()*A0.z() + ); + + Vector Q0; + switch ( index_of_max_component(A02) ) + { + case 0: Q0 = Vector(- A0.z()/A0.x(),0 ,1 ); break; + case 1: Q0 = Vector(0 ,- A0.z()/A0.y(),1 ); break; + case 2: Q0 = Vector(1 ,0 ,- A0.x()/A0.z()); break; + + default : Q0 = NULL_VECTOR ; // This should never happen! + } + + CGAL_assertion( Q0 != NULL_VECTOR ) ; + + Vector Q1 = cross_product(A0,Q0); + + Vector A1 = H * Q0 ; + Vector A2 = H * Q1 ; + FT b1 = - ( Q0 * c ) ; + FT b2 = - ( Q1 * c ) ; + + Add_if_alpha_compatible(A1,b1); + Add_if_alpha_compatible(A2,b2); + + } + break ; + + case 2: + { + + Vector Q = cross_product(A.r0(),A.r1()); + + Vector A2 = H * Q ; + + FT b2 = - ( Q * c ) ; + + Add_if_alpha_compatible(A2,b2); + + } + break ; + + } +} + +} } // namespace Triangulated_surface_mesh::Simplification + +CGAL_END_NAMESPACE + +#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROMTURK_CORE_IMPL_H // +// EOF // + + diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse_functor_base.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_collapse_functor_base.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse_functor_base.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_collapse_functor_base.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_length_cost.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_length_cost.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Edge_length_cost.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_collapse_data.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_collapse_data.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_collapse_data.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_collapse_data.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_cost.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_cost.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_cost.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_set_collapse_data.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_set_collapse_data.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_set_collapse_data.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_set_collapse_data.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_vertex_placement.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_vertex_placement.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/LindstromTurk_vertex_placement.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/LindstromTurk_vertex_placement.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Midpoint_vertex_placement.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_vertex_placement.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Midpoint_vertex_placement.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Midpoint_vertex_placement.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Minimal_collapse_data.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Minimal_collapse_data.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Minimal_collapse_data.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Minimal_collapse_data.h diff --git a/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Set_minimal_collapse_data.h b/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Set_minimal_collapse_data.h similarity index 100% rename from Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Set_minimal_collapse_data.h rename to Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Set_minimal_collapse_data.h