Contour lines generation

The generation of contour lines is a crucial step in accurately representing the terrain features within a mesh. In this project, the ContourConstraint class is responsible for implementing the contour lines computation. This class is defined in contour_constraint.hpp and implemented in contour_constraint.cpp.

1. Class definition

The ContourConstraint class is designed to work with an existing mesh, analyzing its vertices and edges to generate contour lines at specified intervals of elevation (z values). Below is the class definition:


class ContourConstraint
{
public:
    using mesh_type = Mesh;
    using edge_type = mesh_type::edge_type;
    using poly_type = Polygon;
    using value_type = mesh_type::point_type::value_type;
    using cache_type = std::unordered_map<std::string, value_type>;
    ContourConstraint(mesh_type const &mesh);
    std::unique_ptr<mesh_type> apply(value_type precision);
    std::pair<value_type, value_type> getMinMaxZ() const;
    poly_type getBorder() const;

private:
    cache_type M_cache;
    mesh_type const &M_meshref;
    poly_type M_border;

    std::vector<std::unique_ptr<mesh_type>> computeAllContoursMeshes(value_type interval);
    std::unique_ptr<mesh_type> assembleContourMeshes(std::vector<std::unique_ptr<mesh_type>> const &meshes);
    std::unique_ptr<mesh_type> triangulateAssembledMesh(mesh_type const &nonTriangulated);
    std::unique_ptr<mesh_type> computeContourMesh(value_type z);
    void registerZs(mesh_type const &mesh);
    void postprocess(mesh_type &mesh);
    void restoreHeight(mesh_type &target);
    void reset_cache();
    value_type findClosestZ(mesh_type::point_type const &vertex);
    void registerVertex(mesh_type::point_type const &vertex);
    std::string hash_xy(value_type const &x, value_type const &y) const;
    decltype(M_cache)::const_iterator getZFromCache(value_type const &x, value_type const &y) const;
};

2. Implementation details

The ContourConstraint class works by iterating through the mesh’s triangles and computing where each triangle intersects with the horizontal planes at the specified contour levels. These intersections form the contour lines. This process primarily involves the following steps:

2.1. computeAllContoursMeshes

This function first determines the minimum and maximum elevation (z values) present in the mesh using the getMinMaxZ() method. It then computes contour lines for each level within this range based on the specified precision. The main work of generating contour lines is done by computeContourMesh(), which finds intersections of mesh triangles with horizontal planes at the specified z levels.


std::vector<std::unique_ptr<ContourConstraint::mesh_type>> ContourConstraint::computeAllContoursMeshes(value_type interval)
{
    auto [min, max] = getMinMaxZ();
    std::vector<std::unique_ptr<mesh_type>> meshes;
    for (auto z = min; z <= max; z += interval)
    {
        auto msh = computeContourMesh(z);
        meshes.push_back(std::move(msh));
    }
    return meshes;
}

2.2. getMinMaxZ

This function scans through all the nodes of the provided mesh to find the lowest and highest z values. These values are then used to determine the range of elevations at which contour lines will be generated.


std::pair<ContourConstraint::value_type, ContourConstraint::value_type> ContourConstraint::getMinMaxZ() const
{
    auto const &points = M_meshref.points();
    std::vector<value_type> zs;
    zs.reserve(points.size());
    std::transform(points.begin(), points.end(), std::back_insert_iterator(zs), [](auto const &it)
                   { return it.second->z(); });
    auto mimas = std::minmax_element(zs.cbegin(), zs.cend());
    value_type min = *(mimas.first);
    value_type max = *(mimas.second);
    return {min, max};
}

2.3. computeContourMesh

This method is responsable for finding the intersection points between the meshs triangles and a horizontal plane at a specific z level iterating over all the triangles of the mesh. While going through one triangle the method checks for the following cases:

  • The line coincides with the edge of the triangle

  • The line intersects with the triangle

  • One of the vertices of the triangle is at the contour level

The case where all 3 vertices of the triangle are at the contour level is not considered as it represents a flat area.

If the intersection points are found, they are stored as segments.


std::unique_ptr<ContourConstraint::mesh_type> ContourConstraint::computeContourMesh(value_type contourLevel)
{
    auto medges = std::make_unique<mesh_type>();

    for (const auto &triangle : M_meshref.triangles())
    {
        std::vector<mesh_type::point_type> localIntersections;
        for (int i = 0; i < 3; i++)
        {
            int j = (i + 1) % 3;
            const auto &pt1 = triangle->point(i);
            const auto &pt2 = triangle->point(j);

            value_type z1 = pt1.z();
            value_type z2 = pt2.z();

            // CASE 1: Line coincides with edge
            if (std::fabs(z1 - contourLevel) < TOLERANCE && std::fabs(z2 - contourLevel) < TOLERANCE)
            {
                localIntersections.emplace_back(pt1.x(), pt1.y(), pt1.z());
                localIntersections.emplace_back(pt2.x(), pt2.y(), pt2.z());
                continue;
            }

            // CASE 2: Intersection with triangle
            if ((z1 < contourLevel && z2 > contourLevel) || (z1 > contourLevel && z2 < contourLevel))
            {
                value_type t = (contourLevel - z1) / (z2 - z1);
                value_type x = pt1.x() + t * (pt2.x() - pt1.x());
                value_type y = pt1.y() + t * (pt2.y() - pt1.y());
                localIntersections.emplace_back(x, y, contourLevel);
            }

            // CASE 3: One vertex at contour level
            if (std::fabs(z1 - contourLevel) < TOLERANCE && std::fabs(z2 - contourLevel) >= TOLERANCE)
            {
                localIntersections.emplace_back(pt1.x(), pt1.y(), pt1.z());
            }
            if (std::fabs(z2 - contourLevel) < TOLERANCE && std::fabs(z1 - contourLevel) >= TOLERANCE)
            {
                localIntersections.emplace_back(pt2.x(), pt2.y(), pt2.z());
            }
        }

        // Store the intersection as a segment if there are exactly two points
        if (localIntersections.size() == 2)
        {
            auto pt1 = std::make_unique<mesh_type::point_type>(localIntersections[0].x(), localIntersections[0].y(), localIntersections[0].z());
            auto pt2 = std::make_unique<mesh_type::point_type>(localIntersections[1].x(), localIntersections[1].y(), localIntersections[1].z());
            auto pt1b = medges->addPoint(std::move(pt1));
            auto pt2b = medges->addPoint(std::move(pt2));

            medges->addEdge(std::make_unique<mesh_type::edge_type>(pt1b, pt2b));
        }
    }

    algorithms::mergeDuplicatePoints(*medges);

    return medges;
}

The methods assembleContourMeshes and triangulateAssembledMesh are used for combining the computed contour meshes and preparing them for further processing, but these details will be discussed in the next subsection.

In summary, the ContourConstraint class effectively generates contour lines from a given mesh by finding the intersections between the mesh’s triangles and horizontal planes at specified intervals, creating a series of contour meshes that represent these intersections.

3. Results

Here are visualizations of the contour lines for the mesh generated with the LambdaGenerator using the \(z(x, y) = \sqrt{r^2 - (x - x_0)^2 - (y - y_0)^2}\) function:


contour lambda 1 points
Figure 1. Lambda 1 contour lines points

contour lambda 1 edges
Figure 2. Lambda 1 contour lines edges

Here are visualizations of the contour lines for the mesh generated with the LambdaGenerator using the \(z(x, y) = sin(x) + sin(y)\) function:


contour lambda 2 points
Figure 3. Lambda 2 contour lines points

contour lambda 2 edges
Figure 4. Lambda 2 contour lines edges

Here are visualizations of contour lines for the mesh generated with the GpsGenerator class for the tiles containing the tile with coordinates (x, y) = (33809, 23527) at zoom level 16, which includes the location at longitude 5.7232 and latitude 45.1835, located in Grenoble, France:


contour strasbourg 16 1
Figure 5. GPS contour lines points

contour strasbourg 16 2
Figure 6. GPS contour lines edges

All the meshes visualized here were exported in the MSH format, using the ExporterMeshGmsh class, and were visualized using the Gmsh software.


Next: Re-triangulation


References