mesh_geometry.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_MESH_MESH_GEOMETRY
#define GEOGRAM_MESH_MESH_GEOMETRY
 
#include <geogram/basic/common.h>
#include <geogram/mesh/mesh.h>
#include <geogram/basic/geometry.h>
#include <geogram/basic/geometry_nd.h>
 
namespace GEO {
 
    namespace Geom {
 
        inline const vec3& mesh_vertex(const Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 3);
            return *(const vec3*) (M.vertices.point_ptr(v));
        }
 
        inline const vec3& mesh_vertex_ref(const Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 3);
            return *(vec3 const *) (M.vertices.point_ptr(v));
        }
 
        inline vec3& mesh_vertex_ref(Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 3);
            return *(vec3*) (M.vertices.point_ptr(v));
        }
 
        inline const vec3& mesh_corner_vertex(const Mesh& M, index_t c) {
            return mesh_vertex(M, M.facet_corners.vertex(c));
        }
 
        inline vec3& mesh_corner_vertex_ref(Mesh& M, index_t c) {
            return mesh_vertex_ref(M, M.facet_corners.vertex(c));
        }
 
        inline const vec3& mesh_vertex_normal(const Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 6);
            return *(const vec3*) (M.vertices.point_ptr(v) + 3);
        }
 
        inline vec3& mesh_vertex_normal_ref(Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 6);
            return *(vec3*) (M.vertices.point_ptr(v) + 3);
        }
 
        inline const vec3& mesh_vertex_normal_ref(const Mesh& M, index_t v) {
            geo_debug_assert(M.vertices.dimension() >= 6);
            return *(vec3 const *) (M.vertices.point_ptr(v) + 3);
        }
 
        inline double mesh_facet_area(const Mesh& M, index_t f, index_t dim=0) {
            geo_debug_assert(dim <= M.vertices.dimension());
            if(dim == 0) {
                dim = M.vertices.dimension();
            }
            double result = 0.0;
            // Check for empty facet, should not happen.
            if(M.facets.corners_end(f) == M.facets.corners_begin(f)) {
                return result;
            }
            const double* p0 = M.vertices.point_ptr(
                M.facet_corners.vertex(M.facets.corners_begin(f))
            );
            for(
                index_t i = M.facets.corners_begin(f) + 1;
                i + 1 < M.facets.corners_end(f); i++
            ) {
                result += GEO::Geom::triangle_area(
                    p0,
                    M.vertices.point_ptr(M.facet_corners.vertex(i)),
                    M.vertices.point_ptr(M.facet_corners.vertex(i + 1)),
                    coord_index_t(dim)
                );
            }
            return result;
        }
 
        vec3 GEOGRAM_API  mesh_facet_normal(const Mesh& M, index_t f);
 
        inline vec3 mesh_facet_center(const Mesh& M, index_t f) {
            vec3 result(0.0, 0.0, 0.0);
            double count = 0.0;
            for(index_t c = M.facets.corners_begin(f);
                c < M.facets.corners_end(f); ++c) {
                result += Geom::mesh_corner_vertex(M, c);
                count += 1.0;
            }
            return (1.0 / count) * result;
        }
 
        inline vec3 mesh_cell_center(const Mesh& M, index_t c) {
            vec3 result(0.0, 0.0, 0.0);
            for(index_t lv=0; lv<M.cells.nb_vertices(c); ++lv) {
                index_t v = M.cells.vertex(c,lv);
                result += vec3(M.vertices.point_ptr(v));
            }
            return (1.0 / double(M.cells.nb_vertices(c))) * result;
        }
 
 
        inline vec3 mesh_tet_center(const Mesh& M, index_t t) {
            index_t iv1 = M.cells.vertex(t, 0);
            index_t iv2 = M.cells.vertex(t, 1);
            index_t iv3 = M.cells.vertex(t, 2);
            index_t iv4 = M.cells.vertex(t, 3);
            const vec3& v1 = Geom::mesh_vertex(M, iv1);
            const vec3& v2 = Geom::mesh_vertex(M, iv2);
            const vec3& v3 = Geom::mesh_vertex(M, iv3);
            const vec3& v4 = Geom::mesh_vertex(M, iv4);
            return 0.25 * (v1 + v2 + v3 + v4);
        }
 
        inline vec3 mesh_corner_vector(const Mesh& M, index_t c1) {
            geo_debug_assert(M.facets.are_simplices());
            index_t c2 = M.facets.next_corner_around_facet(c1/3, c1);
            index_t v1 = M.facet_corners.vertex(c1);
            index_t v2 = M.facet_corners.vertex(c2);
            return mesh_vertex(M,v2) - mesh_vertex(M,v1);
        }
 
        double GEOGRAM_API mesh_normal_angle(const Mesh& M, index_t c);
 
        double GEOGRAM_API mesh_unsigned_normal_angle(
            const Mesh& M, index_t f1, index_t f2
        );
 
        double GEOGRAM_API mesh_area(const Mesh& M, index_t dim);
 
        inline double mesh_area(const Mesh& M) {
            return mesh_area(M, M.vertices.dimension());
        }
 
        double GEOGRAM_API mesh_enclosed_volume(const Mesh& M);
    }
 
    void GEOGRAM_API compute_normals(Mesh& M);
 
    void GEOGRAM_API simple_Laplacian_smooth(
        Mesh& M, index_t nb_iter, bool normals_only
    );
 
    void GEOGRAM_API get_bbox(const Mesh& M, double* xyzmin, double* xyzmax);
 
    double GEOGRAM_API bbox_diagonal(const Mesh& M);
 
    void GEOGRAM_API set_anisotropy(Mesh& M, double s);
 
    void GEOGRAM_API unset_anisotropy(Mesh& M);
 
    void GEOGRAM_API compute_sizing_field(
        Mesh& M, double gradation = 1.0, index_t nb_lfs_samples = 0
    );
 
    void GEOGRAM_API normalize_embedding_area(Mesh& M);
 
    double GEOGRAM_API mesh_cell_volume(
        const Mesh& M, index_t c
    );
 
    double GEOGRAM_API mesh_cells_volume(const Mesh& M);
 
 
    vec3 GEOGRAM_API mesh_cell_facet_normal(
        const Mesh& M, index_t c, index_t lf
    );
 
    double GEOGRAM_API surface_average_edge_length(
        const Mesh& M
    );
}
 
#endif