nmwsharp / geometry-central

Applied 3D geometry in C++, with a focus on surface meshes.
https://geometry-central.net
MIT License
1.07k stars 149 forks source link

Calling FlipEdgeNetwork::delaunayRefine() doesn't seem to preserve the shape of the original mesh #162

Open dpar39 opened 1 year ago

dpar39 commented 1 year ago

I'm trying to compute a new mesh from the intrinsic triangulation, but results are not quite as expected. Any clue why some triangles near right edges are being modified like this? Below is the code I'm using to generate a new mesh on which the geodesic path are all surface points on vertices.

PathNetworkSPtr GeodesicPathExtractor::getPathNetworkOnDelaunayMesh(double areaThresh,
                                                                    size_t maxInsertions,
                                                                    double angleBound) const
{
    using namespace geometrycentral::surface;

    //_edgeNetwork->makeDelaunay();
    _edgeNetwork->delaunayRefine(areaThresh, maxInsertions, angleBound);

    auto & tri = _edgeNetwork->tri;
    auto newMeshCopy = tri->intrinsicMesh->copy();

    const auto & vertexData = _mesh->getGeometry().vertexPositions;
    auto * intrinsicMeshPtr = tri->intrinsicMesh.get();
    auto * meshCopyPtr = newMeshCopy.get();

    VertexData<Vector3> newPositions(*newMeshCopy);
    for (auto v : newMeshCopy->vertices()) {
        Vertex vi(intrinsicMeshPtr, v.getIndex());
        auto spOnOrig = tri->vertexLocations[vi];
        newPositions[v] = spOnOrig.interpolate(vertexData);
    }

    const auto vertexIndices = tri->intrinsicMesh->getVertexIndices();

    std::vector<std::vector<SurfacePoint>> newPathNetwork;
    newPathNetwork.reserve(_edgeNetwork->paths.size());
    for (auto & path : _edgeNetwork->paths) {
        const auto heList = path->getHalfedgeList();
        std::vector<SurfacePoint> newSubpath;
        newSubpath.reserve(heList.size() + 1);
        for (auto he : heList) {
            auto vA = Vertex(meshCopyPtr, vertexIndices[he.tailVertex()]);
            auto vB = Vertex(meshCopyPtr, vertexIndices[he.tipVertex()]);
            if (newSubpath.empty() || newSubpath.back() != vA) {
                newSubpath.emplace_back(vA);
            }
            newSubpath.emplace_back(vB);
        }
        newPathNetwork.emplace_back(std::move(newSubpath));
    }

    newMeshCopy->compress();

    auto newGeometry = std::make_unique<VertexPositionGeometry>(*newMeshCopy, newPositions);
    Mesh3DBaseSPtr newMesh = std::make_shared<Mesh3DBase>(std::move(newMeshCopy), std::move(newGeometry));
    return std::make_shared<PathNetwork>(std::move(newPathNetwork), newMesh);
}

image

nmwsharp commented 10 months ago

Hi, I think this is a misunderstanding of how intrinsic triangulations work. The intrinsic geometry of the original surface is indeed exactly preserved, in the sense that length/angles/areas/etc are unchanged. In this sense you are working with an abstract mesh along the surface of the input geometry. However, if you take the vertex positions and draw straight lines between them, you get a different surface back---as you've noticed here, the corners get cut off.

You could use the common subdivision to extract a mesh which preserves the geometry exactly, but that is generally not what you want.

Usually the whole point with intrinsic triangulations is to avoid realizing the actual intrinsic mesh sitting in space. You generally can't (without approximations like you see here), and if you do you are just back to doing ordinary mesh operations.

patrick-laurin commented 6 months ago

I believe it's about replicating section 4.9 Comparison to Traditional Remeshing in Geometry Processing with Intrinsic Triangulations

I was also wondering if there was an easy way to get the remeshed surface but it looks like it must be done by hand.