lindstrom_turk_core_impl.h

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

    for ( typename Boundary_data_vector::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_ECMS_LT_TRACE(2,"Hb:" << matrix_to_string(Hb) << "\ncb:" << xyz_to_string(cb) ) ;
    
    //
    // Weighted average
    //
    FT lScaledBoundaryWeight = FT(9) * FT(mParams.BoundaryWeight) * squared_distance(mProfile.p0(),mProfile.p1())  ;
    
    H *= mParams.VolumeWeight ;
    c = c * mParams.VolumeWeight ;
    
    H += lScaledBoundaryWeight * Hb ;
    c = c + ( lScaledBoundaryWeight * cb ) ;
    
    CGAL_ECMS_LT_TRACE(2,"H:" << matrix_to_string(H) << "\nc:" << xyz_to_string(c) ) ;
    CGAL_ECMS_LT_TRACE(2,"VolW:" << mParams.VolumeWeight << " BdryW:" << mParams.BoundaryWeight << " ScaledBdryW:" << lScaledBoundaryWeight ) ;
    
  }
  
  mConstrians.Add_from_gradient(H,c);
}

template<class ECM>
void LindstromTurkCore<ECM>::Add_shape_optimization_constrians( vertex_descriptor_vector const& aLink )
{
  FT s((double)aLink.size());
  
  Matrix H (  s,0.0,0.0
           ,0.0,  s,0.0
           ,0.0,0.0,  s
           );
           
  Vector c = NULL_VECTOR ;
  
  for( typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it )
    c = c + (ORIGIN - get_point(*it)) ;  
           
  CGAL_ECMS_LT_TRACE(1,"Adding shape optimization constrians. Shape vector: " << xyz_to_string(c) );
  
  mConstrians.Add_from_gradient(H,c);
}

template<class ECM>
typename LindstromTurkCore<ECM>::FT
LindstromTurkCore<ECM>::Compute_boundary_cost( Vector const& v, Boundary_data_vector const& aBdry )
{
  FT rCost(0);
  for ( typename Boundary_data_vector::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 ECM>
typename LindstromTurkCore<ECM>::FT
LindstromTurkCore<ECM>::Compute_volume_cost( Vector const& v, Triangle_data_vector const& aTriangles )
{
  FT rCost(0);
  
  for( typename Triangle_data_vector::const_iterator it = aTriangles.begin(), eit = aTriangles.end() ; it != eit ; ++it )
  {
    Triangle_data const& lTri = *it ;
    
    FT lF = lTri.NormalV * v - lTri.NormalL ;
    
    rCost += ( lF * lF ) ;
    
  }
    
  return rCost / FT(36) ;
}

template<class ECM>
typename LindstromTurkCore<ECM>::FT
LindstromTurkCore<ECM>::Compute_shape_cost( Point const& p, vertex_descriptor_vector const& aLink )
{
  FT rCost(0);
  
  for( typename vertex_descriptor_vector::const_iterator it = aLink.begin(), eit = aLink.end() ; it != eit ; ++it )
    rCost += squared_distance(p,get_point(*it)) ;  
    
  return rCost ;
}

template<class ECM>
void LindstromTurkCore<ECM>::Constrians::Add_if_alpha_compatible( Vector const& Ai, FT const& bi )
{
  FT slai = Ai*Ai ;
  
  CGAL_ECMS_LT_TRACE(3,"[constrians] Adding new if alpha-compatble.\nslai: " << slai );
  
  FT l = CGAL_NTS sqrt(slai) ;
  
  CGAL_ECMS_LT_TRACE(3,"[constrians] Adding new if alpha-compatble.\nslai: " << slai );
  
  Vector Ain ;
  FT     bin ;
  
  if ( !CGAL_NTS is_zero(l) )
  {
    Ain = Ai / l ;
    bin = bi / l ;
  }
  else
  {
    CGAL_ECMS_LT_TRACE(3,"[constrians] l is ZERO." );
    
    Ain = Vector( big_value(), big_value(), big_value() ) ;
    bin = big_value() ;
  }
  
  CGAL_ECMS_LT_TRACE(3,"[constrians] Ain: " << xyz_to_string(Ain) << " bin:" << bin );
  
  bool lAddIt = true ;
  
  if ( n == 1 )
  {
    FT d01 = A.r0() * Ai  ;
    
    FT sla0 = A.r0() * A.r0() ;
    FT sd01 = d01 * d01 ;

    FT max = sla0 * slai * squared_cos_alpha() ;     
    
    CGAL_ECMS_LT_TRACE(3,"[constrians] Second constrain. d01: " << d01 << " sla0:" << sla0 << " sd01:" << sd01 << " max:" << max );
    
    if ( sd01 > max )
      lAddIt = false ;
  }
  else if ( n == 2 )
  {
    Vector N = cross_product(A.r0(),A.r1());
    
    FT dc012 = N * Ai ;
    
    FT slc01  = N * N ;
    FT sdc012 = dc012 * dc012;       

    FT min = slc01 * slai * squared_sin_alpha() ;
    
    CGAL_ECMS_LT_TRACE(3,"[constrians] Thirds constrain. N: " << xyz_to_string(N) << " dc012:" << dc012 << " slc01:" << slc01 
                       << " sdc012:" << sdc012 << " min:" << min );
    
    if ( sdc012 <= min )
      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 ;
    }
    
    CGAL_ECMS_LT_TRACE(3,"[constrains] Accepting # " << n << " A:" << matrix_to_string(A) << " b:" << xyz_to_string(b) ) ;
    
    ++ n ;
    
  }
  else
  {
    CGAL_ECMS_LT_TRACE(3,"[constrains] INCOMPATIBLE. Discarded" ) ;
  }
}

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 ECM>
void LindstromTurkCore<ECM>::Constrians::Add_from_gradient ( Matrix const& H, Vector const& c )
{
  CGAL_ECMS_LT_TRACE(3,"[constrains] Adding from gradient. Current n=" << n ) ;
  
  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();
        
        CGAL_assertion( A0 != NULL_VECTOR ) ;
        
        Vector AbsA0( CGAL_NTS abs(A0.x())
                    , CGAL_NTS abs(A0.y())
                    , CGAL_NTS abs(A0.z())
                    );
           
        Vector Q0;      
        switch ( index_of_max_component(AbsA0) ) 
        {
          // Since A0 is guaranteed to be non-zero, the denominators here are known to be non-zero too.
          
          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_ECMS_LT_TRACE(3,"[constrains] Q0:" << xyz_to_string(Q0) ) ;
        
        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 ) ;
        
        CGAL_ECMS_LT_TRACE(3,"[constrains] Q1:" << xyz_to_string(Q1) << " A1: " << xyz_to_string(A1) << " A2:" << xyz_to_string(A2) << "\nb1:" << b1 << " b2:" << b2 ) ;
        
        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 ) ;
        
        CGAL_ECMS_LT_TRACE(3,"[constrains] Q:" << xyz_to_string(Q) << " A2: " << xyz_to_string(A2) << " b2:" << b2 ) ;
        Add_if_alpha_compatible(A2,b2);
        
      }
      break ;
      
  }
}

} // namespace Surface_mesh_simplification

CGAL_END_NAMESPACE

#endif // CGAL_SURFACE_MESH_SIMPLIFICATION_POLICIES_EDGE_COLLAPSE_DETAIL_LINDSTROMTURK_CORE_IMPL_H //
// EOF //