EdgeWeightEvaluators.hh 6.63 KB
Newer Older
Valentin Nigolian's avatar
Valentin Nigolian committed
1 2 3 4 5 6 7 8 9 10 11 12 13
#pragma once

#include <OpenVolumeMesh/Core/GeometryKernel.hh>
#include <OpenVolumeMesh/Mesh/TetrahedralMesh.hh>
#include <cfloat>


namespace OpenVolumeMesh{

namespace Laplacian{



14

Valentin Nigolian's avatar
Valentin Nigolian committed
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
/** \brief Gives all-1 halfedge weights*/
template<class _polyhedral_mesh>
class UniformEdgeWeightEvaluator{

public:

    using Scalar = typename _polyhedral_mesh::PointT::value_type;

    static Scalar halfedge_weight(_polyhedral_mesh& mesh,
                                  const HalfEdgeHandle& edge){
        return 1.;
    }

};



/** \brief Laplacian based on the "Dual Face" of edges.
 * See [Alexa 2020] Properties of Laplace Operators for Tetrahedral Meshes. */
template<class _tetrahedral_mesh>
class DualEdgeWeightEvaluator{


public:

    using VecT   = typename _tetrahedral_mesh::PointT;
    using Scalar = typename VecT::value_type;


    /** See paper mentioned above for computations and naming detail */
    static Scalar halfedge_weight(_tetrahedral_mesh& mesh,
                                  const HalfEdgeHandle& edge){

        Scalar weight(0);

Valentin Nigolian's avatar
Valentin Nigolian committed
50
        double alpha(0.f), beta(0.f), theta(0.f);
Valentin Nigolian's avatar
Valentin Nigolian committed
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

        auto hehf_iter = mesh.hehf_iter(edge);

        /*the idea is to iterate through the halffaces around the halfedge to
         * compute the successive beta angles and then use them as alpha angles
         * for the next iteration. */
        VertexHandle xi = mesh.halfedge_vertices(edge)[0];
        VertexHandle xj = mesh.halfedge_vertices(edge)[1];

        VertexHandle xk;
        for(const auto& v_it: mesh.halfface_vertices(*hehf_iter)){
            if(v_it != xi && v_it != xj){
                xk = v_it;
            }
        }

        VecT xk_xi = (mesh.vertex(xi) - mesh.vertex(xk)).normalized();
        VecT xk_xj = (mesh.vertex(xj) - mesh.vertex(xk)).normalized();

        VecT nijk = xk_xi.cross(xk_xj).normalized();

        alpha = acos(xk_xi.dot(xk_xj));

        double first_alpha = alpha;
        VecT first_nijk = nijk;

        hehf_iter++;

79 80
        int tet_count(0);

Valentin Nigolian's avatar
Valentin Nigolian committed
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
        while(hehf_iter.valid()){

            VertexHandle xl;
            for(const auto& v_it: mesh.halfface_vertices(*hehf_iter)){
                if(v_it != xi && v_it != xj){
                    xl = v_it;
                }
            }

            VecT xl_xi = (mesh.vertex(xi) - mesh.vertex(xl)).normalized();
            VecT xl_xj = (mesh.vertex(xj) - mesh.vertex(xl)).normalized();

            VecT nijl = xl_xi.cross(xl_xj).normalized();

            beta = acos(xl_xi.dot(xl_xj));
            theta = acos(nijk.dot(nijl));

            weight += tet_weight(alpha, beta, theta);

Valentin Nigolian's avatar
Valentin Nigolian committed
100
#if 0
101 102 103 104 105 106
            std::cout<<" ===== TET "<<tet_count<<std::endl;
            std::cout<<"  xl : "<<xl<<std::endl;
            std::cout<<"  xk : "<<xk<<std::endl;
            std::cout<<"  xl-xi : "<<xl_xi<<std::endl;
            std::cout<<"  xl-xj : "<<xl_xj<<std::endl;
            std::cout<<" -------"<<std::endl;
Valentin Nigolian's avatar
digging  
Valentin Nigolian committed
107
            std::cout<<"  alpha : "<<std::setprecision(20)<<alpha<<std::endl;
108 109 110 111 112
            std::cout<<"   beta : "<<beta<<std::endl;
            std::cout<<"  theta : "<<theta<<std::endl;
            std::cout<<"  cot alpha : "<<cot(alpha)<<std::endl;
            std::cout<<"  cot  beta : "<<cot(beta)<<std::endl;
            std::cout<<"  cot theta : "<<cot(theta)<<std::endl;
Valentin Nigolian's avatar
digging  
Valentin Nigolian committed
113 114
            double tet_w = tet_weight(alpha, beta, theta);
            std::cout<<"  tet weight : "<<tet_w<<std::endl;
115
            std::cout<<"  --> cummulated weight : "<<weight<<std::endl;
Valentin Nigolian's avatar
Valentin Nigolian committed
116
#endif
117

Valentin Nigolian's avatar
Valentin Nigolian committed
118 119 120 121 122 123
            //setting-up next iteration
            xk = xl;
            alpha = beta;
            nijk = nijl;

            hehf_iter++;
124
            tet_count++;
Valentin Nigolian's avatar
Valentin Nigolian committed
125 126 127 128 129 130 131 132 133
        }

        //compute the weight contribution of the last cell, which consists of the
        //last halfface and the first one
        beta = first_alpha;
        theta = acos(nijk.dot(first_nijk));

        weight += tet_weight(alpha, beta, theta);

Valentin Nigolian's avatar
Valentin Nigolian committed
134
#if 0
135 136 137 138 139 140 141
        std::cout<<" ===== TET "<<tet_count<<std::endl;
        std::cout<<"  alpha : "<<alpha<<std::endl;
        std::cout<<"   beta : "<<beta<<std::endl;
        std::cout<<"  theta : "<<theta<<std::endl;
        std::cout<<"  cot alpha : "<<cot(alpha)<<std::endl;
        std::cout<<"  cot  beta : "<<cot(beta)<<std::endl;
        std::cout<<"  cot theta : "<<cot(theta)<<std::endl;
Valentin Nigolian's avatar
digging  
Valentin Nigolian committed
142 143
        double tet_w = tet_weight(alpha, beta, theta);
        std::cout<<"  tet weight : "<<tet_w<<std::endl;
144
        std::cout<<"  --> cummulated weight : "<<weight<<std::endl;
Valentin Nigolian's avatar
Valentin Nigolian committed
145
#endif
146

Valentin Nigolian's avatar
Valentin Nigolian committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160
        //multiply by Vol(i,j) (i.e. edge length)
        weight *= (mesh.vertex(xi) - mesh.vertex(xj)).norm() / 8.;

        return weight;
    }


private:

    static Scalar cot(Scalar x){
        return cos(x)/sin(x);
    }


Valentin Nigolian's avatar
Valentin Nigolian committed
161

Valentin Nigolian's avatar
Valentin Nigolian committed
162 163
    static Scalar tet_weight(Scalar alpha, Scalar beta, Scalar theta) {

Valentin Nigolian's avatar
Valentin Nigolian committed
164 165 166 167
        if(fabs(sin(alpha)) > std::numeric_limits<Scalar>::epsilon()  &&
                fabs(sin(beta)) > std::numeric_limits<Scalar>::epsilon() &&
                fabs(sin(theta)) > std::numeric_limits<Scalar>::epsilon() &&
                fabs(cos(theta)) > std::numeric_limits<Scalar>::epsilon()){
Valentin Nigolian's avatar
Valentin Nigolian committed
168 169 170 171 172 173 174

            double cotan_alpha = cot(alpha);
            double cotan_beta  = cot(beta);

            //weight computation following the paper's formula
            return cot(theta) * (2. * (cotan_alpha * cotan_beta)/cos(theta) - cotan_alpha*cotan_alpha - cotan_beta*cotan_beta);;
        }else{
Valentin Nigolian's avatar
Valentin Nigolian committed
175
#if 0
Valentin Nigolian's avatar
digging  
Valentin Nigolian committed
176
            std::cout<<" can't compute tet weight : "<<std::endl;
Valentin Nigolian's avatar
Valentin Nigolian committed
177
            std::cout<<"  epsilon   : "<<std::numeric_limits<Scalar>::epsilon() <<std::endl;
Valentin Nigolian's avatar
digging  
Valentin Nigolian committed
178 179 180 181
            std::cout<<"  sin alpha : "<<std::setprecision(20)<<sin(alpha)<<std::endl;
            std::cout<<"  sin  beta : "<<std::setprecision(20)<<sin(beta)<<std::endl;
            std::cout<<"  sin theta : "<<std::setprecision(20)<<sin(theta)<<std::endl;
            std::cout<<"  cos theta : "<<std::setprecision(20)<<cos(theta)<<std::endl;
Valentin Nigolian's avatar
Valentin Nigolian committed
182 183 184 185 186
            std::cout<<"  abs sin alpha : "<<std::setprecision(20)<<abs(sin(alpha))<<std::endl;
            std::cout<<"  abs sin  beta : "<<std::setprecision(20)<<abs(sin(beta))<<std::endl;
            std::cout<<"  abs sin theta : "<<std::setprecision(20)<<abs(sin(theta))<<std::endl;
            std::cout<<"  abs cos theta : "<<std::setprecision(20)<<abs(cos(theta))<<std::endl;
            std::cout<<" comparisons : "<<
Valentin Nigolian's avatar
Valentin Nigolian committed
187 188 189 190
                       (fabs(sin(alpha)) > std::numeric_limits<Scalar>::epsilon()) << " " <<
                       (fabs(sin(beta))  > std::numeric_limits<Scalar>::epsilon()) << " "<<
                       (fabs(sin(theta)) > std::numeric_limits<Scalar>::epsilon()) << " "<<
                       (fabs(cos(theta)) > std::numeric_limits<Scalar>::epsilon()) << std::endl;
Valentin Nigolian's avatar
Valentin Nigolian committed
191
#endif
Valentin Nigolian's avatar
Valentin Nigolian committed
192 193 194 195 196 197 198 199 200 201
            return 0.f;
        }
    }

};


}
}