lindstrom_turk_core_impl.h

很多二维 三维几何计算算法 C++ 类库

// 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 Surface_mesh_simplification license may use this file in
// accordance with the Surface_mesh_simplification 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Surface_mesh_simplification/include/CGAL/Surface_mesh_simplification/Policies/Edge_collapse/Detail/Lindstrom_Turk_core_impl.h $
// $Id: Lindstrom_Turk_core_impl.h 37469 2007-03-26 08:11:41Z afabri $
//
// Author(s)     : Fernando Cacciola <[e-mail]>
//
#ifndef CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROM_TURK_CORE_IMPL_H
#define CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_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 Surface_mesh_simplification
{

template<class ECM>
LindstromTurkCore<ECM>::LindstromTurkCore( Params const& aParams, Profile const& aProfile )
  :
   mParams(aParams)
  ,mProfile(aProfile)
{
  Extract_triangle_data();
  Extract_boundary_data();
}

template<class ECM>
void LindstromTurkCore<ECM>::Extract_boundary_data()
{
  for ( const_border_edge_iterator it = mProfile.border_edges().begin(), eit = mProfile.border_edges().end() ; it != eit ; ++ it )
  {
    edge_descriptor border_edge = *it ;
    
    edge_descriptor face_edge = opposite_edge(border_edge,surface()) ;
        
    vertex_descriptor sv = source(face_edge,surface());
    vertex_descriptor tv = target(face_edge,surface());
    
    Point const& sp = get_point(sv);
    Point const& tp = get_point(tv);
    
    Vector v = tp - sp ;
    Vector n = Point_cross_product(tp,sp) ;
    
    CGAL_ECMS_LT_TRACE(1,"Boundary edge. S:" << xyz_to_string(sp) << " T:" << xyz_to_string(tp)
                      << " V:" << xyz_to_string(v) << " N:" << xyz_to_string(n) 
                      ) ;
    
    mBdry_data.push_back( Boundary_data(sp,tp,v,n) ) ;
  }  
}

template<class ECM>
void LindstromTurkCore<ECM>::Extract_triangle_data()
{
  for ( const_triangle_iterator it = mProfile.triangles().begin(), eit = mProfile.triangles().end() ; it != eit ; ++ it )
  {
    Triangle const& tri = *it ;
    
    Point const& p0 = get_point(tri.v0);
    Point const& p1 = get_point(tri.v1);
    Point const& p2 = get_point(tri.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_ECMS_LT_TRACE(1,"Extracting triangle v" << tri.v0->id() << "->v" << tri.v1->id() << "->v" << tri.v2->id()
                      << " N:" << xyz_to_string(lNormalV) << " L:" << lNormalL
                      );
    
    mTriangle_data.push_back(Triangle_data(lNormalV,lNormalL));
  }
}

template<class ECM>
typename LindstromTurkCore<ECM>::Optional_point LindstromTurkCore<ECM>::compute_placement()
{
  Optional_point  rPlacementP ;
  Optional_vector lPlacementV ;
  
  CGAL_ECMS_LT_TRACE(0,"Computing LT data for E" << mProfile.v0v1()->id() << " (V" << mProfile.v0()->id() << "->V" << mProfile.v1()->id() << ")" );
  
  //
  // 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 ( mBdry_data.size() > 0 )
    Add_boundary_preservation_constrians(mBdry_data);
  
  if ( mConstrians.n < 3 )
    Add_volume_preservation_constrians(mTriangle_data);
    
  if ( mConstrians.n < 3 )
    Add_boundary_and_volume_optimization_constrians(mBdry_data,mTriangle_data); 
    
  if ( mConstrians.n < 3 )
    Add_shape_optimization_constrians(mProfile.link());
  
  // 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 ;
      
      CGAL_ECMS_LT_TRACE(2,"     b: " << xyz_to_string(mConstrians.b) );
      CGAL_ECMS_LT_TRACE(2,"inv(A): " << matrix_to_string(lAi) );
      
      lPlacementV = mConstrians.b * lAi ;
      
      CGAL_ECMS_LT_TRACE(0,"New vertex point: " << xyz_to_string(*lPlacementV) );
    }
    else
      CGAL_ECMS_LT_TRACE(1,"Can't solve optimization, singular system.");
  }
  else
    CGAL_ECMS_LT_TRACE(1,"Can't solve optimization, not enough alpha-compatible constrians.");
  
  if ( lPlacementV )
    rPlacementP = filter_infinity( ORIGIN + (*lPlacementV) );
    
  return rPlacementP;
}

template<class ECM>
typename LindstromTurkCore<ECM>::Optional_FT LindstromTurkCore<ECM>::compute_cost( Optional_point const& aP )
{
  Optional_FT rCost ;

  if ( aP )
  {
    Vector lV = (*aP) - ORIGIN ;

    FT lSquaredLength = squared_distance(mProfile.p0(),mProfile.p1());
   
    FT lBdryCost   = Compute_boundary_cost(lV ,mBdry_data);
    FT lVolumeCost = Compute_volume_cost  (lV ,mTriangle_data);
    FT lShapeCost  = Compute_shape_cost   (*aP,mProfile.link());
    
    FT lTotalCost =   FT(mParams.VolumeWeight)   * lVolumeCost
                    + FT(mParams.BoundaryWeight) * lBdryCost   * lSquaredLength
                    + FT(mParams.ShapeWeight)    * lShapeCost  * lSquaredLength * lSquaredLength ;
    
    rCost = filter_infinity(lTotalCost);
    
    CGAL_ECMS_LT_TRACE(0,"\nSquared edge length: " << lSquaredLength 
                      << "\nBoundary cost: " << lBdryCost << " weight: " << mParams.BoundaryWeight
                      << "\nVolume cost: " << lVolumeCost << " weight: " << mParams.VolumeWeight
                      << "\nShape cost: " << lShapeCost  << " weight: " << mParams.ShapeWeight
                      << "\nTOTAL COST: " << lTotalCost 
                      );
                      
  }
    
  return rCost;
}


template<class ECM>
void LindstromTurkCore<ECM>::Add_boundary_preservation_constrians( Boundary_data_vector const& aBdry )
{
 
  if ( aBdry.size() > 0 )
  {
    Vector e1 = NULL_VECTOR ; 
    Vector e2 = NULL_VECTOR ;

    FT e1x = FT(0.0), e1y = FT(0.0), e1z = FT(0.0) ;
    
    for ( typename Boundary_data_vector::const_iterator it = aBdry.begin() ; it != aBdry.end() ; ++ it )
    {
      e1 = e1 + it->v ;
      e2 = e2 + it->n ;
      
      FT vx = it->v.x() ;
      FT vy = it->v.y() ;
      FT vz = it->v.z() ;
      
      e1x = e1x + vx;
      e1y = e1y + vy;
      e1z = e1z + vz;
      
      CGAL_ECMS_LT_TRACE(1,"vx:" << vx << " vy:" << vy << " vz:" << vz << " e1x:" << e1x << " e1y:" << e1y << " e1z:" << e1z );
    }     
    
    Matrix H = LT_product(e1);
    
    Vector c = cross_product(e1,e2);
    
    CGAL_ECMS_LT_TRACE(1
                      ,"Adding boundary preservation constrians. SumV:" << xyz_to_string(e1) 
                       << " SumN:" << xyz_to_string(e2) 
                       << "\nH:" << matrix_to_string(H)
                       << "\nc:" << xyz_to_string(c)                       
                      );
    
    mConstrians.Add_from_gradient(H,c);
  }
}

template<class ECM>
void LindstromTurkCore<ECM>::Add_volume_preservation_constrians( Triangle_data_vector const& aTriangles )
{
  CGAL_ECMS_LT_TRACE(1,"Adding volume preservation constrians. " << aTriangles.size() << " triangles.");
  
  Vector lSumV = NULL_VECTOR ;
  FT     lSumL(0) ;
  
  for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
  {
    lSumV = lSumV + it->NormalV ;
    lSumL = lSumL + it->NormalL ;  
  }   
  
  CGAL_ECMS_LT_TRACE(1, " SumV:" << xyz_to_string(lSumV) << " SumL:" << lSumL );
  
  mConstrians.Add_if_alpha_compatible(lSumV,lSumL);   

}

template<class ECM>
void LindstromTurkCore<ECM>::Add_boundary_and_volume_optimization_constrians( Boundary_data_vector const& aBdry
                                                                            , Triangle_data_vector const& aTriangles
                                                                            )
{
  CGAL_ECMS_LT_TRACE(1,"Adding boundary and volume optimization constrians. ");
  
  Matrix H = NULL_MATRIX ;
  Vector c = NULL_VECTOR ;

  //
  // Volume optimization  
  //
  for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
  {
    Triangle_data const& lTri = *it ;
    
    H += direct_product(lTri.NormalV,lTri.NormalV) ;
    
    c = c - ( lTri.NormalL * lTri.NormalV ) ;
  }   
  
  CGAL_ECMS_LT_TRACE(2,"Hv:" << matrix_to_string(H) << "\ncv:" << xyz_to_string(c) ) ;
  
  
  if ( aBdry.size() > 0 )
  {
    //
    // Boundary optimization
    //
    Matrix Hb = NULL_MATRIX ;
    Vector cb = NULL_VECTOR ;