Mesh refining run forever
sonukiller opened this issue · 6 comments
When I am trying to implement 6_delaunay_refinement/, the step of refining before coarsening takes a lot of time (also refining after coarsening) and runs for forever. I am getting results quickly for the sample mesh boxpart.obj but for some big meshes (>30K), this problem is coming. Any solution?
Also, I am getting this error for dragon_fat.obj
malloc(): invalid next size (unsorted)
Aborted (core dumped)
I tried mesh simplification for dragon_fat.obj using Blender and reduced the vertices to 5.5K and it was working with --refinement_time=BEFORE and I got same error as above with --refinement_time=AFTER! Can you please help with this?
The infinite loop is due to the or condition. What do you think?
It might help to decrease the angle_threshold argument to the delaunay_refine function, which controls how much Delaunay refinement tries to improve triangle quality. The default value of 25 degrees has usually worked pretty well for me, but can cause problems on meshes with very pointy vertices
I tried to reduce the default value of 25 degrees to 5 degrees and it was working but not as expected.
I implemented Dealuney refinement using this repo and it is working for 25 degrees also.
What is the difference between the Delaunay refinement using this repo and here? And one more thing I want to ask, does your coarsening algorithm have implementation using geometry central library?
Ah, interesting! No, this code uses a reimplementation of Delaunay refinement, and in this implementation, I left out some code to handle complicated edge cases (specifically this part about parent faces) which I thought were pretty rare. But these edge cases could well come up more often than I thought. I'm a bit busy at the moment, but I can probably take a closer look at the code in a week or two.
In the meantime, here's some code we used to do Delaunay refinement using geometry-central while working on this project. It sounds like switching to using geometry-central's implementation will hopefully solve the problems you're having
Code
#include <Eigen/Core>
#include <Eigen/Dense>
#include "geometrycentral/surface/manifold_surface_mesh.h"
#include "geometrycentral/surface/vertex_position_geometry.h"
#include "geometrycentral/surface/signpost_intrinsic_triangulation.h"
#include "global_variables.h"
#include "glue_face_sides.h"
/*
doing delaunay refinement and transfer barycentric coordinates to the refined mesh
Inputs:
VO: input vertex list
FO: input face list
bary_coords: |BC|x3 barycentric coordinates
bary_faces: |BC|x1 index to FO
Outputs:
F: delaunay refined face list
l: delaunay refined edge lengths
bary_coords: |BC|x3 barycentric coordinates
bary_faces: |BC|x1 index to F
*/
void delaunay_refine(
const Eigen::MatrixXd & VO,
const Eigen::MatrixXi & FO,
Eigen::MatrixXi & F,
Eigen::MatrixXi & G,
Eigen::MatrixXd & l)
{
using namespace std;
using namespace Eigen;
using namespace geometrycentral;
using namespace surface;
using namespace global_variables;
ManifoldSurfaceMesh mesh(FO);
VertexPositionGeometry geometry(mesh, VO);
SignpostIntrinsicTriangulation tri(mesh, geometry);
tri.delaunayRefine();
tri.intrinsicMesh->validateConnectivity();
int nF = tri.intrinsicMesh->nFaces();
F.resize(nF,3);
l.resize(nF,3);
F = tri.intrinsicMesh->getFaceVertexMatrix<int>();
tri.requireEdgeLengths();
for(Face f : tri.intrinsicMesh->faces()) {
size_t iHe = 0;
for(Halfedge he : f.adjacentHalfedges()) {
if (f.getIndex() >= nF)
throw std::runtime_error("[Error] delaunayRefine deleted vertices");
l(f.getIndex(),iHe) = tri.edgeLengths[he.edge()];
iHe++;
}
}
// build glue map
G.resize(nF, 3*2);
G.fill(0);
for(Face f : tri.intrinsicMesh->faces()) {
size_t s = 0; // side
for(Halfedge he : f.adjacentHalfedges()) {
// the face side of "he"
Vector2i fs = {he.face().getIndex(), s};
// get the face side of "he.twin()"
Vector2i fs_twin;
if (he.twin().isInterior())
{
// if he.twin() is an inteior half-edge, then find its face-side
int s_twin = 0; // side of he_tmp
for (Halfedge he_tmp : he.twin().face().adjacentHalfedges())
{
if (he_tmp == he.twin()) // if found it, then break
break;
s_twin++;
}
fs_twin << he.twin().face().getIndex(), s_twin;
}
else
{
// if he.twin() is an exterior half-edge, then set it to ghost
fs_twin << GHOST_INDEX, GHOST_INDEX;
}
glue_face_sides(fs, fs_twin, G);
s++;
}
}
}
void delaunay_refine(
const Eigen::MatrixXd & VO,
const Eigen::MatrixXi & FO,
Eigen::MatrixXi & F,
Eigen::MatrixXi & G,
Eigen::MatrixXd & l,
Eigen::MatrixXd & bary_coords,
Eigen::VectorXi & bary_faces)
{
using namespace std;
using namespace Eigen;
using namespace geometrycentral;
using namespace surface;
using namespace global_variables;
ManifoldSurfaceMesh mesh(FO);
VertexPositionGeometry geometry(mesh, VO);
SignpostIntrinsicTriangulation tri(mesh, geometry);
tri.delaunayRefine();
tri.intrinsicMesh->validateConnectivity();
for (int ii=0; ii<bary_coords.rows(); ii++)
{
if (bary_faces(ii) != GHOST_INDEX) // if this bary is valid
{
int f = bary_faces(ii);
geometrycentral::Vector3 bc{bary_coords(ii,0), bary_coords(ii,1), bary_coords(ii,2)};
SurfacePoint bary_point(tri.inputMesh.face(f), bc);
SurfacePoint bary_point_intrinsic = tri.equivalentPointOnIntrinsic(bary_point);
bary_faces(ii) = bary_point_intrinsic.face.getIndex();
bary_coords.row(ii) << bary_point_intrinsic.faceCoords[0], bary_point_intrinsic.faceCoords[1], bary_point_intrinsic.faceCoords[2];
}
}
int nF = tri.intrinsicMesh->nFaces();
F.resize(nF,3);
F = tri.intrinsicMesh->getFaceVertexMatrix<int>();
l.resize(nF, 3);
tri.requireEdgeLengths();
for(Face f : tri.intrinsicMesh->faces()) {
size_t iHe = 0;
for(Halfedge he : f.adjacentHalfedges()) {
l(f.getIndex(),iHe) = tri.edgeLengths[he.edge()];
iHe++;
}
}
// build glue map
G.resize(nF, 3*2);
G.fill(0);
for(Face f : tri.intrinsicMesh->faces()) {
size_t s = 0; // side
for(Halfedge he : f.adjacentHalfedges()) {
// the face side of "he"
Vector2i fs = {he.face().getIndex(), s};
// get the face side of "he.twin()"
Vector2i fs_twin;
if (he.twin().isInterior())
{
// if he.twin() is an inteior half-edge, then find its face-side
int s_twin = 0; // side of he_tmp
for (Halfedge he_tmp : he.twin().face().adjacentHalfedges())
{
if (he_tmp == he.twin()) // if found it, then break
break;
s_twin++;
}
fs_twin << he.twin().face().getIndex(), s_twin;
}
else
{
// if he.twin() is an exterior half-edge, then set it to ghost
fs_twin << GHOST_INDEX, GHOST_INDEX;
}
glue_face_sides(fs, fs_twin, G);
s++;
}
}
}
void delaunay_refine(
const Eigen::MatrixXi & FO,
const Eigen::MatrixXd & lO,
Eigen::MatrixXi & F,
Eigen::MatrixXi & G,
Eigen::MatrixXd & l)
{
using namespace std;
using namespace Eigen;
using namespace geometrycentral;
using namespace surface;
using namespace global_variables;
ManifoldSurfaceMesh mesh(FO);
EdgeData<double> edgeLengths(mesh);
for(Face f : mesh.faces()) {
size_t iHe = 0;
for(Halfedge he : f.adjacentHalfedges()) {
edgeLengths[he.edge()] = lO(f.getIndex(),iHe);
iHe++;
}
}
EdgeLengthGeometry geometry(mesh, edgeLengths);
SignpostIntrinsicTriangulation tri(mesh, geometry);
tri.delaunayRefine();
tri.intrinsicMesh->validateConnectivity();
int nF = tri.intrinsicMesh->nFaces();
F.resize(nF,3);
l.resize(nF,3);
F = tri.intrinsicMesh->getFaceVertexMatrix<int>();
tri.requireEdgeLengths();
for(Face f : tri.intrinsicMesh->faces()) {
size_t iHe = 0;
for(Halfedge he : f.adjacentHalfedges()) {
if (f.getIndex() >= nF)
throw std::runtime_error("[Error] delaunayRefine deleted vertices");
l(f.getIndex(),iHe) = tri.edgeLengths[he.edge()];
iHe++;
}
}
// build glue map
G.resize(nF, 3*2);
G.fill(0);
for(Face f : tri.intrinsicMesh->faces()) {
size_t s = 0; // side
for(Halfedge he : f.adjacentHalfedges()) {
// the face side of "he"
Vector2i fs = {he.face().getIndex(), s};
// get the face side of "he.twin()"
Vector2i fs_twin;
if (he.twin().isInterior())
{
// if he.twin() is an inteior half-edge, then find its face-side
int s_twin = 0; // side of he_tmp
for (Halfedge he_tmp : he.twin().face().adjacentHalfedges())
{
if (he_tmp == he.twin()) // if found it, then break
break;
s_twin++;
}
fs_twin << he.twin().face().getIndex(), s_twin;
}
else
{
// if he.twin() is an exterior half-edge, then set it to ghost
fs_twin << GHOST_INDEX, GHOST_INDEX;
}
glue_face_sides(fs, fs_twin, G);
s++;
}
}
}Thanks for the code!
The coarsening algorithm is implemented before and after Delaunay refinement and it requires F, G, l, A, v2fs, BC, F2V. The above three delaunay_refine functions not fully output these values.
Can you please provide a delaunay_refine function using geometry-central which will output all these values required for the coarsening algorithm or if possible can you provide the implementation of the coarsening algorithm using geometry-central.
A possibility is to use the delaunay refinement from geometry central to obtain an intrinsic mesh, and you can then compute the information using these two functions
build_glue_map(F,G);
build_angular_coordinates(F,G,l,A,v2fs);