Problems with mesh contraction, NaN values gets generated

[mesabg]

I'm having issues with this implementation. I already read the paper and have everything in place. Sometimes my linear algebra equation solver rans into NaN values, so basically some vertices became unstable. How is the proper way to handle this ?

Note: I also have a basic mesh degeneration detection which takes the one ring edges length and mark the vertices so they can preserve their current position, but even though there's a moment where the linear algebra solver returns NaN in some vertices

Also another question is related to the mesh characteristics for this to work. Are there any special constraints ? Cause I have notice, that scaling up or down the model gives different results.

Environment: Also C++ with Eigen library. Tested with all Sparse Solvers that the library provides

[harrisonyei]

It yields different result with scaled input models may be caused by the weights of constraints, you will have to apply proper weights for different models.
Our solution doesn’t handle different models either, we manually pitch those values for each case.
As the mesh contracts, it will inevitably facing numerical errors, your proposed method above for detecting edge length might help avoiding this error.
Although, how did you implement your one-ring length detection as stop criteria for vertices, does it calculate every iteration for all vertices and how you form your linear system, I might need more information.

[mesabg]

That's right, the weighted constraints affects the different models in different ways. At the moment I'm trying to get a set of models which I can use for testing purposes with the same sL and constraints values. Can you provide me with the models that you used for testing ? And also the values used for that generation, it would be definitely helpful.

So, my approach to handle that is scaling the model into a range, let's say [-1, 1] or [-100, 100] to minimize that calculation error and avoid the calculation issues.

About the one ring edges lengths calculation, I recalculate for all of them after getting the Vi solution. Pretty much someting like:

for e in edges:
length = calculate edge length
vertices_one_ring_edge_lengths[u] += length
vertices_one_ring_edge_lengths[v] += length

And then I just evaluate if the length is less than a threshold so, I can block the position of that vertex, causing not moving on the following iterations

But as I said, even with all this apporaches, I still not able to see a fully functional skeleton

[mesabg]

About the linear system I use the following code, is Eigen, same library you used:

typedef SparseMatrix<double> SparseMatrixXd;
typedef Eigen::SimplicialLDLT<SparseMatrixXd> CholeskyLDLTSolver;

MatrixXd choleskyLdltSolver(SparseMatrixXd &A, MatrixXd &b) {
    SparseMatrixXd At = A.transpose();
    SparseMatrixXd AtA = At * A;

    // Solve System
    CholeskyLDLTSolver solver;
    solver.compute(AtA);

    return solver.solve(At * b);
}

[harrisonyei]

By "block the position of the vertex", you mean that the matrices treat those vertices' positions as constants?

Irrelevant to avoiding numerical errors the paper propose using volume as an indicator,
but in our case instead we just run many times and select an acceptable result mesh from these iterations to enforce the second step - edge collapsing.

[mesabg]

Yes, that's right making them fixed

Do you use 1e-6 as threshold (ratio volume) to stop the process ?

And yes, I haven't done yet the edge collapsing, but I will

[harrisonyei]

No we didnt implement the volume stop criteria, we run iterations until numerical error occurs and cannot contract anymore, then select from iteration results.
This project is an easy implementation for mesh skeleton extraction.