eltopo3d/dynamicsurface.cpp
// ---------------------------------------------------------
///
/// Compute rank of the quadric metric tensor at a vertex
///
// ---------------------------------------------------------
Compute rank of the quadric metric tensor at a vertex
Compute rank of the quadric metric tensor at a vertex
unsigned int DynamicSurface::classify_vertex( size_t v ) const
{     
    if ( m_mesh.m_vertex_to_triangle_map[v].empty() )     { return 0; }
    const std::vector<size_t>& incident_triangles = m_mesh.m_vertex_to_triangle_map[v];
    Mat33d A(0,0,0,0,0,0,0,0,0);
    for ( size_t i = 0; i < incident_triangles.size(); ++i )
    {
        size_t triangle_index = incident_triangles[i];
        Vec3d normal = get_triangle_normal(triangle_index);
        double w = get_triangle_area(triangle_index);
        A(0,0) += normal[0] * w * normal[0];
        A(1,0) += normal[1] * w * normal[0];
        A(2,0) += normal[2] * w * normal[0];
        A(0,1) += normal[0] * w * normal[1];
        A(1,1) += normal[1] * w * normal[1];
        A(2,1) += normal[2] * w * normal[1];
        A(0,2) += normal[0] * w * normal[2];
        A(1,2) += normal[1] * w * normal[2];
        A(2,2) += normal[2] * w * normal[2];
    }
    double eigenvalues[3];
    double work[9];
    int info = ~0, n = 3, lwork = 9;
    LAPACK::get_eigen_decomposition( &n, A.a, &n, eigenvalues, work, &lwork, &info );
    if ( info != 0 )
    {
        if ( m_verbose )
        {
            std::cout << "Eigen decomp failed.  Incident triangles: " << std::endl;
            for ( size_t i = 0; i < incident_triangles.size(); ++i )
            {
                size_t triangle_index = incident_triangles[i];
                Vec3d normal = get_triangle_normal(triangle_index);
                double w = get_triangle_area(triangle_index);
                std::cout << "normal: ( " << normal << " )    ";  
                std::cout << "area: " << w << std::endl;
            }
        }
        return 4;
    }
    unsigned int rank = 0;
    for ( unsigned int i = 0; i < 3; ++i )
    {
        if ( eigenvalues[i] > G_EIGENVALUE_RANK_RATIO * eigenvalues[2] )
        {
            ++rank;
        }
    }
    return rank;
}
unsigned int DynamicSurface :: classify_vertex ( size_t v ) const { if ( m_mesh . m_vertex_to_triangle_map [ v ] . empty ( ) ) { return 0 ; } const std :: vector < size_t > & incident_triangles = m_mesh . m_vertex_to_triangle_map [ v ] ; Mat33d A ( 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) ; for ( size_t i = 0 ; i < incident_triangles . size ( ) ; ++ i ) { size_t triangle_index = incident_triangles [ i ] ; Vec3d normal = get_triangle_normal ( triangle_index ) ; double w = get_triangle_area ( triangle_index ) ; A ( 0 , 0 ) += normal [ 0 ] * w * normal [ 0 ] ; A ( 1 , 0 ) += normal [ 1 ] * w * normal [ 0 ] ; A ( 2 , 0 ) += normal [ 2 ] * w * normal [ 0 ] ; A ( 0 , 1 ) += normal [ 0 ] * w * normal [ 1 ] ; A ( 1 , 1 ) += normal [ 1 ] * w * normal [ 1 ] ; A ( 2 , 1 ) += normal [ 2 ] * w * normal [ 1 ] ; A ( 0 , 2 ) += normal [ 0 ] * w * normal [ 2 ] ; A ( 1 , 2 ) += normal [ 1 ] * w * normal [ 2 ] ; A ( 2 , 2 ) += normal [ 2 ] * w * normal [ 2 ] ; } double eigenvalues [ 3 ] ; double work [ 9 ] ; int info = ~ 0 , n = 3 , lwork = 9 ; LAPACK :: get_eigen_decomposition ( & n , A . a , & n , eigenvalues , work , & lwork , & info ) ; if ( info != 0 ) { if ( m_verbose ) { std :: cout << " " << std :: endl ; for ( size_t i = 0 ; i < incident_triangles . size ( ) ; ++ i ) { size_t triangle_index = incident_triangles [ i ] ; Vec3d normal = get_triangle_normal ( triangle_index ) ; double w = get_triangle_area ( triangle_index ) ; std :: cout << " " << normal << " " ; std :: cout << " " << w << std :: endl ; } } return 4 ; } unsigned int rank = 0 ; for ( unsigned int i = 0 ; i < 3 ; ++ i ) { if ( eigenvalues [ i ] > G_EIGENVALUE_RANK_RATIO * eigenvalues [ 2 ] ) { ++ rank ; } } return rank ; }
Compute rank of the quadric metric tensor at a vertex
Compute rank of the quadric metric tensor at a vertex
// get eigen decomposition // compute rank of primary space