Gobal changes to optimize memory consumption -- UNFINISHED WORK

This commit is contained in:
Fernando Cacciola 2006-08-24 13:30:51 +00:00
parent e2cc288e30
commit 42766fef0a
15 changed files with 932 additions and 1 deletions

2
.gitattributes vendored
View File

@ -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

View File

@ -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 <fernando_cacciola@ciudad.com.ar>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H 1
#include <vector>
#include <CGAL/Surface_mesh_simplification/TSMS_common.h>
#include <CGAL/Surface_mesh_simplification/Policies/LindstromTurk_collapse_data.h>
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 Collapse_data_>
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<TSM> GraphTraits ;
typedef typename GraphTraits::in_edge_iterator in_edge_iterator ;
typedef typename Surface_geometric_traits<TSM>::Kernel Kernel ;
typedef typename Kernel::Point_3 Point ;
typedef typename Kernel::Vector_3 Vector ;
typedef typename Kernel::FT FT ;
typedef optional<FT> Optional_FT ;
typedef optional<Point> Optional_point ;
typedef optional<Vector> Optional_vector ;
typedef MatrixC33<Kernel> 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<Triangle> Triangles ;
typedef std::vector<vertex_descriptor> Link ;
typedef std::vector<edge_descriptor> 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<BoundaryEdge> 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 <CGAL/Surface_mesh_simplification/Policies/Detail/Lindstrom_Turk_core_impl.h>
#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_LINDSTROM_TURK_CORE_H //
// EOF //

View File

@ -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 <fernando_cacciola@ciudad.com.ar>
//
#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<class CD>
LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<Matrix> 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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
typename LindstromTurkCore<CD>::Triangle LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class CD>
typename LindstromTurkCore<CD>::FT
LindstromTurkCore<CD>::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<class CD>
typename LindstromTurkCore<CD>::FT
LindstromTurkCore<CD>::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<class CD>
typename LindstromTurkCore<CD>::FT
LindstromTurkCore<CD>::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<class CD>
void LindstromTurkCore<CD>::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<class V>
int index_of_max_component ( V const& v )
{
typedef typename Kernel_traits<V>::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<class CD>
void LindstromTurkCore<CD>::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 //