nmwsharp / geometry-central

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

Exception when traceBack geodesic curve. #133

Closed Cindy-xdZhang closed 1 year ago

Cindy-xdZhang commented 1 year ago

Hello there, thank you for providing this framework, I have encountered a crash without knowing what's the cause:

Code for building a gaussian-looklike mesh:

std::unique_ptr<geometrycentral::surface::SimplePolygonMesh> makeGaussianBumpMesh(const Vector2& extent, const Vector2& bumpCenter, float sigma)
{

    std::unique_ptr<geometrycentral::surface::SimplePolygonMesh> simpleMesh = std::make_unique<geometrycentral::surface::SimplePolygonMesh>();
    simpleMesh->clear();

    float xExtent = extent.x;
    float yExtent = extent.y;

    Vector3 startPos;
    startPos.x = xExtent * -0.5f;
    startPos.y = yExtent * -0.5f;
    startPos.z = 0.f;

    Vector3 offsetX = Vector3::zero();
    offsetX.x = xExtent / g_numOfVertCols;
    Vector3 offsetY = Vector3::zero();
    offsetY.y = yExtent / g_numOfVertRows;

    Vector3 bumpCenter3D;
    bumpCenter3D.x = bumpCenter.x;
    bumpCenter3D.y = bumpCenter.y;
    bumpCenter3D.z = 0.f;

    for (int j = 0; j < g_numOfVertRows; j++) {
        for (int i = 0; i < g_numOfVertCols; i++) {

            Vector3 position = startPos + offsetX * i + offsetY * j;
            Vector3 d = position - bumpCenter3D;
            float dist_squared = d.norm2();
            float val = 1.0 / (sigma * sqrt(2.0 * geometrycentral::PI)) * exp(-0.5 * dist_squared / (sigma * sigma));
            position.z = val;
            simpleMesh->vertexCoordinates.push_back(position);
        }
    }

    for (int j = 0; j < g_numOfVertRows - 1; j++) {
        for (int i = 0; i < g_numOfVertCols - 1; i++) {
            int linindx = j * g_numOfVertCols + i;
            int indx00 = linindx;
            int indx01 = linindx + 1;
            int indx10 = linindx + g_numOfVertCols;
            int indx11 = linindx + g_numOfVertCols + 1;
            std::vector<size_t> face1;
            face1.push_back(indx00);
            face1.push_back(indx01);
            face1.push_back(indx11);

            std::vector<size_t> face2;
            face2.push_back(indx00);
            face2.push_back(indx11);
            face2.push_back(indx10);
            simpleMesh->polygons.push_back(face1);
            simpleMesh->polygons.push_back(face2);
        }
    }

    // simpleMesh->stripUnusedVertices();

    return std::move(simpleMesh);
}

code for the geodesic curve:

 std::unique_ptr<geometrycentral::surface::SimplePolygonMesh> simpleMeshBump = makeGaussianBumpMesh(g_bumpMeshParameters.domain_extent, g_bumpMeshParameters.bump_center, g_bumpMeshParameters.bump_sigma);
    std::unique_ptr<geometrycentral::surface::SimplePolygonMesh> simpleMeshBump = makeGaussianBumpMesh(g_bumpMeshParameters.domain_extent, g_bumpMeshParameters.bump_center, g_bumpMeshParameters.bump_sigma);
    auto MeshBump = geometrycentral::surface::makeManifoldSurfaceMeshAndGeometry(simpleMeshBump->polygons, simpleMeshBump->vertexCoordinates);

std::unique_ptr<geometrycentral::surface::GeodesicAlgorithmExact> g_mmp= std::make_unique<geometrycentral::surface::GeodesicAlgorithmExact>(*g_mesh_bump, *g_geometry_bump);
//initialize a source point  via {faceid+barycentric coordinate}
 sourcePoints.push_back(geometrycentral::surface::SurfacePoint(face, fBary));
g_mmp->propagate(sourcePoints);
        for (int i = 0; i < g_mesh_bump->nVertices(); i += 1) {
            geometrycentral::surface::Vertex& v_bump = g_mesh_bump->vertex(i);
            geometrycentral::surface::SurfacePoint queryPoint(v_bump);
            double distance = g_mmp->getDistance(queryPoint); //geodesic distance  works
            geoDistance[v_bump] = distance;
            std::vector<SurfacePoint> geodesicCurve = g_mmp->traceBack(queryPoint); // trace geodesic curve will crash
        }

The trackBack here will throw:

GC_SAFETY_ASSERT FAILURE from \usr\sources\gc-polyscope-project-template\deps\geometry-central\src\surface\exact_geodesic_helpers.cpp:332 - compute_local_coordinates() err: Point [SurfacePoint: type=Vertex, vertex= v_0] not adjacent to e_512

Even when I change the source points or the query point, neither works for the traceBack. Do you have any clue on what might cause this?

Cindy-xdZhang commented 1 year ago

I've fixed this(not sure fully correct just add a filter to possible edges), but I am not sure whether this tricky fix might give you some hints on what's the problem.

void possibleEdgeSafetyfilter(SurfacePoint& point, std::vector<Edge>& storage) {
  std::vector<Edge> temp;
  temp.reserve(storage.size());

  for (const auto& e : storage) {
    SurfacePoint eDummyPt(e, 0.5);
    Face fShared = sharedFace(point, eDummyPt);
    if (fShared != Face()) {
      temp.emplace_back(e);
    }
  }
  storage.swap(temp);
}

then in 
std::vector<SurfacePoint> GeodesicAlgorithmExact::traceBack(const SurfacePoint& point) {
...........
 while (visible_from_source(path.back()) < 0) {
      SurfacePoint& q = path.back();

      possible_traceback_edges(q, possible_edges);
      possibleEdgeSafetyfilter(q, possible_edges);
...........
}
nmwsharp commented 1 year ago

Uh oh, that sounds like it might be a bug. Could you please do two things two help us debug?

Also paging @MarkGillespie in case he immediately knows what is up here.

Cindy-xdZhang commented 1 year ago

Sorry I didin't check github recently, I will give you what you want shortly.

Cindy-xdZhang commented 1 year ago

You can reproduce this exception in this repo branch "reproduce_exception": https://github.com/Cindy-xdZhang/geometry-central/tree/reproduce_exception

I am currently busy on some other projects until next year January 30 , I can work on this bug if it hasn't been fix at that time. again, thanks for the library.

nervoussystem commented 1 year ago

I also encountered this bug. I don't have a stack trace but it seems to have to do with interacting with a mesh boundary. In this image all the white vertices are where it fails

Screenshot 2023-03-17 21 41 52

MarkGillespie commented 1 year ago

Thanks for pointing this out! I made some changes which I think fix the issue, but feel free to reach out again if you find more weird behavior