User Manual

Author

François Protais

Volume Mesh Smoothing and Geometric Fitting

This package implements an optimization algorithm for improving the quality of volumetric meshes and fitting them to geometric targets, based on the method described in [1]. It can be used to improve an already valid mesh, untangle an invalid mesh, or deform a mesh so that its boundary follows a prescribed geometry.

The algorithm only modifies vertex coordinates. It does not insert or remove vertices, change cell connectivity, or alter the combinatorial structure of the input mesh. Consequently, cell indices, material labels, adjacency relations, and other data attached to the mesh are preserved throughout the optimization. When changes to the mesh connectivity or element sizing are required, the package Tetrahedral Remeshing should be used instead.

Smoothing and Fitting Algorithm

The algorithm optimizes vertex positions according to an element-quality energy. It currently uses a conformal energy (MIPS3D) that improves the dihedral angles of the cells. Interior vertices are moved to improve the volume mesh, while boundary vertices can additionally be attracted towards a target geometry.

The energy incorporates a barrier term preventing element inversion. For an initially valid mesh, all accepted optimization steps preserve element orientation, providing a validity guarantee while the element quality is improved. The use of a penalization approach also allows the optimizer to start from an invalid mesh and recover a valid configuration by untangling inverted elements.

Geometric targets can be specified at several dimensions:

• surface targets constrain boundary polygons to locally estimated tangent planes;
• curve targets constrain selected mesh edges to target tangent directions;
• point targets attract individual vertices to prescribed positions.

Surface patches and curves may be assigned different identifiers, allowing different parts of the mesh to use different target geometries. Constraints are soft and weighted by default, so geometric fidelity can be balanced against element quality. Vertices, or individual coordinate dimensions, can also be locked when hard constraints are required. Combining surface, curve, and point targets allows smooth regions, sharp curves, corners, and user handles to be treated within the same optimization.

By recovering the tangent planes of the target geometry, the smoother can recover curvature discontinuities, enabling automatic feature recovery and preservation.

API

The main function of the package is CGAL::boundary_aware_mesh_smoothing(), which takes a model of CGAL::MeshComplex_3InTriangulation_3 as input mesh and a model of ConstructTangentSpace for re-projections. The vertex coordinates are then updated to improve element quality and to fit geometric targets.

CGAL::Mesh_smoothing_3::C3t3_mesh_projector provides a model of ConstructTangentSpace to re-project the input mesh into another mesh represented by a model of CGAL::MeshComplex_3InTriangulation_3.

Warning

The algorithm relocates vertices without updating the connectivity of the triangulation. Therefore, if the input triangulation satisfies a Delaunay or regular triangulation property, this property is not maintained and must be considered lost after smoothing.

Untangling an invalid mesh without any boundary constraints or fitting terms will result in a collapsed mesh. It is a current limitation of the optimization approach.

Examples

Direct Smoothing of a C3t3

The following example demonstrates the direct use of the smoother on a CGAL::Mesh_complex_3_in_triangulation_3. The tetrahedral mesh is read from a Medit file with its surface patches, and CGAL::boundary_aware_mesh_smoothing() is then called to improve the mesh quality while preserving the input surface patches.

Note that the CGAL reader does not currently read feature edges from Medit files, so the current example only preserves input patches.

Figure 69.1 Running the smoother on mambo_m3.mesh slightly improves the dihedral angles, as the initial mesh is already of good quality for its current geometric target. The elements can slide along their respective patches while preserving sharpness and key features. The bunny.mesh file only contains volume elements; consequently, the smoother does not attempt to preserve any boundary surface. The resulting boundary is therefore less geometrically regular, but the quality of the volume elements is significantly improved.

Example: Mesh_smoothing_3/c3t3_smooth.cpp

Show / Hide

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
#include <CGAL/Triangulation_3.h>
#include <CGAL/Triangulation_data_structure_3.h>
#include <CGAL/Simplicial_mesh_cell_base_3.h>
#include <CGAL/Simplicial_mesh_vertex_base_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
 
#include <CGAL/tags.h>
 
#include <CGAL/IO/File_medit.h>
#include <fstream>
 
#include <CGAL/Mesh_smoothing_3/boundary_aware_mesh_smoothing.h>
#include <CGAL/Mesh_smoothing_3/projectors.h>
 
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
 
using Subdomain_index = int;
using Surface_patch_index = unsigned;
using Curve_index = unsigned;
using Corner_index = unsigned;
 
using Cb = CGAL::Simplicial_mesh_cell_base_3<K, Subdomain_index, Surface_patch_index>;
using Vb = CGAL::Simplicial_mesh_vertex_base_3<K, Subdomain_index, Surface_patch_index,
                                                  Curve_index, Corner_index>;
 
using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb, CGAL::Sequential_tag>;
using Triangulation = CGAL::Triangulation_3<K, Tds>;
 
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Triangulation>;
 
 
int main(int argc, char* argv[])
{
    std::string filename = (argc > 1) ? std::string(argv[1])
                                      : "data/mambo_m3.mesh";
 
    C3t3 c3t3;
 
    std::ifstream is(filename, std::ios_base::in);
    if(!CGAL::IO::read_MEDIT(is,c3t3.triangulation()))
    {
      std::cerr << "Failed to read" << std::endl;
      return EXIT_FAILURE;
    }
    c3t3.rescan_after_load_of_triangulation();
 
    std::ofstream os("c3t3_initial.mesh");
    CGAL::IO::write_MEDIT(os, c3t3.triangulation(), CGAL::parameters::all_vertices(true));
    os.close();
 
    auto result =
    CGAL::boundary_aware_mesh_smoothing(
        c3t3,
        CGAL::Mesh_smoothing_3::C3t3_mesh_projector(c3t3),
        CGAL::parameters::verbose(true)
    );
 
    std::cout << "Number of inverted elements: " << result.nb_invalid_elements << std::endl;
    std::cout << "Number of vertex updates: " << result.nb_vertex_updates << std::endl;
    std::cout << "Number of metric evaluations: " << result.nb_metric_evaluations << std::endl;
    std::cout << "Pre-processing time: " << result.pre_processing_time << std::endl;
    std::cout << "Smoothing time: " << result.optimization_time << std::endl;
 
    std::ofstream os2("c3t3_smoothed.mesh");
    CGAL::IO::write_MEDIT(os2, c3t3.triangulation(), CGAL::parameters::all_vertices(true));
    os2.close();
 
    std::cout << "Done" << std::endl;
    return EXIT_SUCCESS;
 
}
 
 
 
 

Feature Recovery on a Mesh Generated from an Implicit Domain

This example builds on the example presented in Section Construction from a Vector of Implicit Functions and a Vector of Strings of the 3D Mesh Generation package. After generating a mesh of the implicit domain using CGAL::make_mesh_3, we smooth the mesh using CGAL::Mesh_smoothing_3::Signed_distance_function_projector to project its boundary back onto the implicit domain. This recovers the sharp features formed by the intersection of the implicit surfaces.

Figure 69.2 Mesh generated by CGAL::make_mesh_3 on an implicit domain (left) and result of the smoothing algorithm (right).

Example: Mesh_smoothing_3/implicit_domain_feature_recovery.cpp

Show / Hide

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Implicit_to_labeling_function_wrapper.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Mesh_3/Dump_c3t3.h>
 
#include <CGAL/Mesh_smoothing_3/boundary_aware_mesh_smoothing.h>
#include <CGAL/Mesh_smoothing_3/projectors.h>
 
#include <utility>
#include <vector>
#include <string>
 
// Kernel
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
using FT = K::FT;
using Point_3 = K::Point_3;
using Vector_3 = K::Vector_3;
 
using Value_and_gradient = std::pair<FT, Vector_3>;
 
// -----------------------------------------------------------------------------
// Implicit functions
// -----------------------------------------------------------------------------
 
Value_and_gradient torus_distance_and_gradient(const Point_3& p)
{
    const FT major_radius = FT(1.5);
    const FT minor_radius = FT(0.5);
 
    const FT x = p.x();
    const FT y = p.y();
    const FT z = p.z();
 
    const FT radial = CGAL::sqrt(x * x + z * z);
    const FT dr = radial - major_radius;
    const FT tube_distance = CGAL::sqrt(dr * dr + y * y);
 
    CGAL_precondition(radial != FT(0));
    CGAL_precondition(tube_distance != FT(0));
 
    const FT value = tube_distance - minor_radius;
 
    const Vector_3 gradient(
        dr * x / (tube_distance * radial),
        y / tube_distance,
        dr * z / (tube_distance * radial));
 
    return {value, gradient};
}
 
Value_and_gradient sphere_distance_and_gradient(const Point_3& p)
{
    const FT x = p.x();
    const FT y = p.y();
    const FT z = p.z();
 
    const FT distance_to_center =
        CGAL::sqrt(x * x + y * y + z * z);
 
    CGAL_precondition(distance_to_center != FT(0));
 
    const FT radius = CGAL::sqrt(FT(3));
 
    return {
        distance_to_center - radius,
        Vector_3(
            x / distance_to_center,
            y / distance_to_center,
            z / distance_to_center)
    };
}
 
// Scalar versions used by Mesh_3.
FT torus_function(const Point_3& p)
{
    return torus_distance_and_gradient(p).first;
}
 
FT sphere_function(const Point_3& p)
{
    return sphere_distance_and_gradient(p).first;
}
 
// -----------------------------------------------------------------------------
// Mesh domain
// -----------------------------------------------------------------------------
 
using Function = FT (*)(const Point_3&);
 
using Function_wrapper =
    CGAL::Implicit_multi_domain_to_labeling_function_wrapper<Function>;
 
using Function_vector = Function_wrapper::Function_vector;
using Mesh_domain = CGAL::Labeled_mesh_domain_3<K>;
 
// Triangulation
using Tr = CGAL::Mesh_triangulation_3<Mesh_domain>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
 
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
using Facet_criteria = Mesh_criteria::Facet_criteria;
using Cell_criteria = Mesh_criteria::Cell_criteria;
 
int main()
{
    namespace params = CGAL::parameters;
 
    // The same domain as mesh_implicit_domains_2.cpp:
    //
    //     torus_function > 0
    //     sphere_function < 0
    //
    Function_vector functions;
    functions.push_back(&torus_function);
    functions.push_back(&sphere_function);
 
    std::vector<std::string> signs;
    signs.push_back("+-");
 
    Mesh_domain domain(
        Function_wrapper(functions, signs),
        K::Sphere_3(
            CGAL::ORIGIN,
            CGAL::square(FT(5))),
        params::relative_error_bound(1e-6));
 
    // Same mesh criteria as the original example.
    Facet_criteria facet_criteria(
        30,     // angle
        0.2,    // size
        0.02);  // approximation
 
    Cell_criteria cell_criteria(
        2.,     // radius-edge ratio
        0.4);   // size
 
    Mesh_criteria criteria(
        facet_criteria,
        cell_criteria);
 
    C3t3 c3t3 =
        CGAL::make_mesh_3<C3t3>(
            domain,
            criteria,
            params::no_exude().no_perturb());
 
    CGAL::dump_c3t3(c3t3, "implicit_initial");
 
    // -------------------------------------------------------------------------
    // Projection target
    // -------------------------------------------------------------------------
    //
    // The domain corresponds to:
    //
    //     torus > 0  &&  sphere < 0
    //
    // Reorient both functions so that the desired domain is negative:
    //
    //     -torus < 0 && sphere < 0.
    //
    // max(-torus, sphere) therefore represents the entire boundary as one
    // implicit function. No patch identifier is used by the projector.
    //
    auto projection_function = [](const Point_3& p) -> Value_and_gradient {
        const auto [torus_distance, torus_gradient] =
            torus_distance_and_gradient(p);
 
        const auto [sphere_distance, sphere_gradient] =
            sphere_distance_and_gradient(p);
 
        // Desired domain:
        //
        //   outside torus  -> -torus_distance < 0
        //   inside sphere  ->  sphere_distance < 0
        //
        const FT outside_torus_distance = -torus_distance;
 
        if(outside_torus_distance > sphere_distance)
        {
            return {
                outside_torus_distance,
                -torus_gradient
            };
        }
 
        return {
            sphere_distance,
            sphere_gradient
        };
    };
 
    using Projector =
        CGAL::Mesh_smoothing_3::Signed_distance_function_projector<
            K,
            decltype(projection_function)>;
 
    // tolerance is expressed in the units of the input geometry.
    Projector projector(
        projection_function,
        10,     // maximum projection iterations
        1e-8);  // positional tolerance
 
    const auto result =
        CGAL::boundary_aware_mesh_smoothing(
            c3t3,
            projector,
            CGAL::parameters::verbose(true));
 
    std::cout << "Number of inverted elements: "
              << result.nb_invalid_elements << '\n';
    std::cout << "Number of vertex updates: "
              << result.nb_vertex_updates << '\n';
    std::cout << "Number of metric evaluations: "
              << result.nb_metric_evaluations << '\n';
    std::cout << "Smoothing time: "
              << result.total_time << " s." << '\n';
 
    CGAL::dump_c3t3(c3t3, "implicit_smoothed");
 
    return EXIT_SUCCESS;
}

Feature Recovery on a Mesh Generated from a Hybrid Domain

This example builds on the example presented in Section Construction of a Hybrid Domain : From an Implicit and a Polyhedral Domain of the 3D Mesh Generation package, but does not explicitly provide the 1-dimensional features of the domain. We extend the hybrid domain so that it becomes a model of ConstructTangentSpace, using CGAL::Mesh_smoothing_3::Polyhedral_mesh_domain_projector for its polyhedral component. The smoothing algorithm can then recover the sharp features from the surface constraints alone.

Figure 69.3 Mesh generated by CGAL::make_mesh_3 on a hybrid domain (left) and result of the smoothing algorithm (center). The close-up on the right shows that, although the features are recovered without creating inverted elements, the local topology of the surface patches can still lead to highly distorted surface facets.

Example: Mesh_smoothing_3/hybrid_domain_feature_recovery.cpp

Show / Hide

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Polyhedral_mesh_domain_3.h>
#include <CGAL/make_mesh_3.h>
 
#include <CGAL/Mesh_smoothing_3/boundary_aware_mesh_smoothing.h>
#include <CGAL/Mesh_smoothing_3/projectors.h>
 
#include <CGAL/centroid.h>
#include <CGAL/Kernel/global_functions_3.h>
 
#include <CGAL/IO/File_medit.h>
 
#include <array>
#include <cmath>
#include <fstream>
#include <iostream>
#include <tuple>
#include <vector>
 
// Kernel
using K = CGAL::Exact_predicates_inexact_constructions_kernel;
 
// Implicit domain
using Implicit_domain = CGAL::Labeled_mesh_domain_3<K>;
 
// Polyhedral domain
using Polyhedron = CGAL::Polyhedron_3<K>;
using Polyhedron_domain = CGAL::Polyhedral_mesh_domain_3<Polyhedron, K>;
using Polyhedron_projector = CGAL::Mesh_smoothing_3::Polyhedral_mesh_domain_projector<Polyhedron_domain>;
 
class Hybrid_domain
{
    const Implicit_domain& implicit_domain;
    const Polyhedron_domain& polyhedron_domain;
    Polyhedron_projector polyhedron_projector;
 
public:
    Hybrid_domain(const Implicit_domain& implicit_domain,
                  const Polyhedron_domain& polyhedron_domain)
        : implicit_domain(implicit_domain)
        , polyhedron_domain(polyhedron_domain)
        , polyhedron_projector(polyhedron_domain)
    {}
 
    // Types required by MeshDomain_3
    using Surface_patch_index = int;
    using Subdomain_index = int;
    using Index = int;
 
    using R = K;
    using Point_3 = K::Point_3;
    using Vector_3 = K::Vector_3;
    using FT = K::FT;
    using Intersection = std::tuple<Point_3, Index, int>;
 
    // Type required by ConstructTangentSpace
    using Tangent_space = CGAL::Mesh_smoothing_3::Tangent_space<K>;
 
    CGAL::Bbox_3 bbox() const
    {
        return implicit_domain.bbox() + polyhedron_domain.bbox();
    }
 
    struct Construct_initial_points
    {
        Construct_initial_points(const Hybrid_domain& domain)
            : r_domain_(domain)
        {}
 
        template<class OutputIterator>
        OutputIterator operator()(OutputIterator pts, const int n = 20) const
        {
            using Implicit_Index = Implicit_domain::Index;
            std::vector<std::pair<Point_3, Implicit_Index> > implicit_points_vector;
 
            Implicit_domain::Construct_initial_points cstr_implicit_initial_points =
                r_domain_.implicit_domain.construct_initial_points_object();
 
            cstr_implicit_initial_points(
                std::back_inserter(implicit_points_vector), n / 2);
 
            for(const auto& p : implicit_points_vector)
                *pts++ = std::make_pair(p.first, 2);
 
            using Polyhedron_Index = Polyhedron_domain::Index;
            std::vector<std::pair<Point_3, Polyhedron_Index> > polyhedron_points_vector;
 
            Polyhedron_domain::Construct_initial_points cstr_polyhedron_initial_points =
                r_domain_.polyhedron_domain.construct_initial_points_object();
 
            cstr_polyhedron_initial_points(
                std::back_inserter(polyhedron_points_vector), n / 2);
 
            for(const auto& p : polyhedron_points_vector)
                *pts++ = std::make_pair(p.first, 1);
 
            return pts;
        }
 
    private:
        const Hybrid_domain& r_domain_;
    };
 
    Construct_initial_points construct_initial_points_object() const
    {
        return Construct_initial_points(*this);
    }
 
    struct Is_in_domain
    {
        Is_in_domain(const Hybrid_domain& domain)
            : r_domain_(domain)
        {}
 
        std::optional<Subdomain_index> operator()(const Point_3& p) const
        {
            const std::optional<Subdomain_index> subdomain_index =
                r_domain_.implicit_domain.is_in_domain_object()(p);
 
            if(subdomain_index)
                return 2;
 
            return r_domain_.polyhedron_domain.is_in_domain_object()(p);
        }
 
    private:
        const Hybrid_domain& r_domain_;
    };
 
    Is_in_domain is_in_domain_object() const
    {
        return Is_in_domain(*this);
    }
 
    struct Construct_intersection
    {
        Construct_intersection(const Hybrid_domain& domain)
            : r_domain_(domain)
        {}
 
        template<typename Query>
        Intersection operator()(const Query& query) const
        {
            using boost::get;
 
            const Implicit_domain::Intersection implicit_inter =
                r_domain_.implicit_domain.construct_intersection_object()(query);
 
            if(get<2>(implicit_inter) != 0)
                return Intersection(get<0>(implicit_inter), 2, get<2>(implicit_inter));
 
            const Polyhedron_domain::Intersection polyhedron_inter =
                r_domain_.polyhedron_domain.construct_intersection_object()(query);
 
            if(get<2>(polyhedron_inter) != 0)
            {
                const Point_3 inter_point = get<0>(polyhedron_inter);
 
                if(!r_domain_.implicit_domain.is_in_domain_object()(inter_point))
                    return Intersection(inter_point, 1, get<2>(polyhedron_inter));
            }
 
            return Intersection();
        }
 
    private:
        const Hybrid_domain& r_domain_;
    };
 
    Construct_intersection construct_intersection_object() const
    {
        return Construct_intersection(*this);
    }
 
    Index index_from_surface_patch_index(const Surface_patch_index& index) const
    {
        return index;
    }
 
    Index index_from_subdomain_index(const Subdomain_index& index) const
    {
        return index;
    }
 
    Surface_patch_index surface_patch_index(const Index& index) const
    {
        return index;
    }
 
    Subdomain_index subdomain_index(const Index& index) const
    {
        return index;
    }
 
    template<typename Patch_face>
    Tangent_space
    patch_face_projection_plane(const Patch_face& patch_face,
                                const std::vector<Point_3>& face_points) const
    {
        const Point_3 face_center =
            CGAL::centroid(face_points.begin(), face_points.end());
 
        // Patch 1 is the polyhedral domain.
        if(patch_face.first == 1)
            return polyhedron_projector.patch_face_projection_plane(patch_face, face_points);
 
        // Patch 2 is the implicit sphere:
        // ||p - (1,1,1)||^2 - 1 = 0.
        CGAL_assertion(patch_face.first == 2);
 
        const Point_3 sphere_center(1., 1., 1.);
        const Vector_3 radial = face_center - sphere_center;
 
        CGAL_assertion(radial.squared_length() != FT(0));
 
        const FT length = CGAL::sqrt(radial.squared_length());
        const Vector_3 normal = radial / length;
        const Point_3 projected = sphere_center + normal;
 
        return Tangent_space{projected, normal};
    }
 
    template<typename Curve_edge>
    Tangent_space
    curve_edge_projection_line(const Curve_edge& curve_edge,
                               const std::array<Point_3, 2>& edge_points) const
    {
        return polyhedron_projector.curve_edge_projection_line(curve_edge, edge_points);
    }
};
 
using Domain = Hybrid_domain;
 
// Triangulation
using Tr = CGAL::Mesh_triangulation_3<Domain, K>::type;
using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
 
// Criteria
using Mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
using Facet_criteria = Mesh_criteria::Facet_criteria;
using Cell_criteria = Mesh_criteria::Cell_criteria;
 
using Point = K::Point_3;
using FT = K::FT;
 
FT sphere_centered_at_111(const Point& p)
{
    const FT dx = p.x() - 1;
    const FT dy = p.y() - 1;
    const FT dz = p.z() - 1;
 
    return dx * dx + dy * dy + dz * dz - 1;
}
 
namespace params = CGAL::parameters;
 
int main()
{
    const std::string fname = CGAL::data_file_path("meshes/cube.off");
 
    Polyhedron polyhedron;
    std::ifstream input(fname);
    input >> polyhedron;
 
    if(input.bad())
    {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }
 
    Polyhedron_domain polyhedron_domain(polyhedron);
 
    Implicit_domain sphere_domain =
        Implicit_domain::create_implicit_mesh_domain(
            sphere_centered_at_111,
            K::Sphere_3(K::Point_3(1, 1, 1), K::FT(2)));
 
    Domain domain(sphere_domain, polyhedron_domain);
 
    Facet_criteria facet_criteria(30, 0.08, 0.025);
    Cell_criteria cell_criteria(2, 0.1);
    Mesh_criteria criteria(facet_criteria, cell_criteria);
 
    C3t3 c3t3 =
        CGAL::make_mesh_3<C3t3>(
            domain,
            criteria,
            params::no_perturb().no_exude());
 
    dump_c3t3(c3t3, "hybrid_initial");
 
 
    const auto result =
        CGAL::boundary_aware_mesh_smoothing(
            c3t3,
            domain,
            CGAL::parameters::verbose(true));
 
    std::cout << "Number of inverted elements: "
              << result.nb_invalid_elements << '\n';
    std::cout << "Number of vertex updates: "
              << result.nb_vertex_updates << '\n';
    std::cout << "Number of metric evaluations: "
              << result.nb_metric_evaluations << '\n';
    std::cout << "Smoothing time: "
              << result.total_time << " s." << '\n';
 
    dump_c3t3(c3t3, "hybrid_smoothed");
 
    return EXIT_SUCCESS;
}

Implementation History

The optimization method implemented in this package was introduced by François Protais, Gianmarco Cherchi, and Marco Livesu in [1].

The original implementation was developed as part of the Mesh_optimization project. The CGAL implementation was subsequently adapted to operate directly on CGAL::Mesh_complex_3_in_triangulation_3 objects and to expose the geometric fitting mechanism through the concepts ConstructTangentSpace and TangentSpace.

It was initially published in CGAL-6.3.