Matrix singularity in cMCF computation in Step 1

[dkobayas-cyber]

Hi, thank you so much for sharing your great work!
I'm trying to run your example Jupyter notebook on my Linux machine with Intel chip. I didn't encounter any issues in the first notebook (Step0_binary_cell_segmentation_surface_meshing_and_curvature_measurement.ipynb) but got an error in the second notebook (Step1_conformalized_mean_curvature_flow_and_S_ref_creation.ipynb). The error happens in the second cell where the conformalized mean curvature flow (cMCF) is computed:

  0%|          | 0/50 [00:00<?, ?it/s]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 2
      1 # run cMCF on the mesh for 50 iterations, delta=5e-4 using the standard cotangent Laplacian 
----> 2 Usteps_MCF, mesh_F, MCF_measures_dict = meshtools.conformalized_mean_curvature_flow(mesh, 
      3                                                                                     max_iter=50, 
      4                                                                                     delta=5e-4, 
      5                                                                                     conformalize=True, # set this flag to run cMCF, else it will run MCF
      6                                                                                     robust_L =False, # if set, runs MCF/cMCF using the robust Laplacian instead of cotangent Laplacian
      7                                                                                     mollify_factor=0,
      8                                                                                     solver='pardiso')  # this is a parameter used in the robust Laplacian
     10 # cMCF is run as a free-form deformation, if we want to have the mesh represent a gradual smoothing in accordance with the physical size of the cell and the image coordinates we need to put back the correct scaling factors.
     11 Usteps_MCF_img = meshtools.recover_img_coordinate_conformal_mean_flow(Usteps_MCF, MCF_measures_dict) # recovers the proper associated image coordinates 

File ~/LC-shadow/u-unwrap3D/unwrap3D/Mesh/meshtools.py:5305, in conformalized_mean_curvature_flow(mesh, max_iter, delta, rescale_output, min_diff, conformalize, robust_L, mollify_factor, solver)
   5301 b = M.dot(U)
   5302 # b = U.copy()
   5303 
   5304 # U = spsparse.linalg.spsolve(S,b)
-> 5305 U = spsolver.spsolve(S,b)
   5307 if np.sum(np.isnan(U[:,0])) > 0:
   5308     # if we detect nan, no good. 
   5309     break 

File /depot/cfrueh/apps/env_shadow/lib/python3.11/site-packages/pypardiso/scipy_aliases.py:44, in spsolve(A, b, factorize, squeeze, solver, *args, **kwargs)
     41 if sp.issparse(A) and A.format == "csc":
     42     A = A.tocsr()  # fixes issue with brightway2 technosphere matrix
---> 44 solver._check_A(A)
     45 if factorize and not solver._is_already_factorized(A):
     46     solver.factorize(A)

File /depot/cfrueh/apps/env_shadow/lib/python3.11/site-packages/pypardiso/pardiso_wrapper.py:238, in PyPardisoSolver._check_A(self, A)
    236 if not np.diff(A.indptr).all():
    237     row_col = 'column' if self._solve_transposed else 'row'
--> 238     raise ValueError('Matrix A is singular, because it contains empty {}(s)'.format(row_col))
    240 if A.dtype != np.float64:
    241     raise TypeError('PyPardiso currently only supports float64, but matrix A has dtype: {}'.format(A.dtype))

ValueError: Matrix A is singular, because it contains empty row(s)

Strangely, the matrix S computed in this line becomes singular for some reason. I'm wondering what a potential cause for this error would be. I'm using the example data (bleb_example.tif) as input.

[fyz11]

Heya,

that's a funky error, for the example data. I just encountered this with one of our interns recently, I've not seen that before... In this case, there was actually no mesh being input into the function, and was a python environment issue (since it was working on his local machine). Did you get a proper binary mask and mesh in step 0 and the input paths are all good?

The S matrix is constructed based on the mesh laplacian matrix, L. This matrix is known to be singular if the mesh has poor triangulation which affects matrix inversion. There is likely a triangle whose internal angle is very close to 0. You can check this using the measure_triangle_props function to compute statistics of the triangle. If this is indeed the case, proceed to solution 1 below i.e. remesh with a lower remesh_samples fraction.

unwrap3D.Mesh.meshtools.measure_triangle_props (c.f. https://github.com/DanuserLab/u-unwrap3D/blob/master/unwrap3D/Mesh/meshtools.py)

There are two solutions then:

1. in the tutorial i.e. step 0, please use a lower remesh_samples in the marching cubes step. I think currently it is 0.9 -> use 0.5 instead. This makes the output mesh use less triangles. The output mesh when viewed in meshlab should have equilateral triangles. There may have been a change in the pyacvd dependency which meant this parameter setting no longer is giving equilateral meshing.....
   

2. You can often still make the cMCF step (step1) work by setting the keyword argument (robust_L=True) in meshtools.conformalized_mean_curvature_flow.

Let me know if this then works.

Felix

[dkobayas-cyber]

Thank you so much for your help! Sadly, I still have no luck figuring out the error so far... Below is a list of things I have tried so far:

Mesh generation

In step 0, I was using remesh_samples=0.9 . As you pointed out, this parameter caused a poor triangulation in the mesh, which gave me the following statistics:

'min_angle': np.float64(0.0012567934950694312), 
'avg_angle': np.float64(60.0),
 'max_angle': np.float64(179.994805757994), 
'std_dev_angle': np.float64(26.108316772153273), 
'min_quality': np.float64(3.0148344741940624e-09), 
'avg_quality': np.float64(0.7568483073357946), 
'max_quality': np.float64(0.9999972296523654),

After I change it to remesh_samples=0.5, the smallest internal angle as well as the quality has improved:

'min_angle': np.float64(20.52528587052864),
 'avg_angle': np.float64(60.00000000000001), 
'max_angle': np.float64(131.59532255198266), 
'std_dev_angle': np.float64(12.524908818513676), 
'min_quality': np.float64(0.31317252309741744), 
'avg_quality': np.float64(0.9321358765361495), 
'max_quality': np.float64(0.9999994116434059)

However, even with the improved mesh, cMCF calculation by meshtools.conformalized_mean_curvature_flow() still failed with the same error message.

Dependency on the solver?

When mapping the surface proximal surface signals in Step 0, I actually got the same error message (i.e., singular matrix error) using the example notebook. However, when I changed the solver from pardiso to 'scipy.sparse`, the mapping calculation was successful:

n_samples = 1./ .104 # total number of steps
stepsize = 0.5 # voxels
    
# flip the mesh vertex coordinates so that it aligns with the volume size 
img_binary_surf_mesh.vertices = img_binary_surf_mesh.vertices[:,::-1].copy()

# run the active contour cMCF to get the coordinates at different depths into the cell according to the external image gradient given by the gradient of the signed distance function.
v_depth = meshtools.parametric_mesh_constant_img_flow(img_binary_surf_mesh, 
                    external_img_gradient = H_sdf_vol_normal.transpose(1,2,3,0), 
                    niters=int(n_samples/stepsize), 
                    deltaL=5e-5, # delta which controls the stiffness of the mesh
                    step_size=stepsize, 
                    method='implicit', # this specifies the cMCF solver.
                    conformalize=True, # ensure we use the cMCF Laplacian
                    robust_L=False,
                    solver='other') # <--Used to be pardiso, but now using scipy.sparse

Output:

0%|          | 0/19 [00:00<?, ?it/s]/home/dkobayas/LC-shadow/u-unwrap3D/unwrap3D/Mesh/meshtools.py:4033: MatrixRankWarning: Matrix is exactly singular
  U = spsolver.spsolve(S,b)
100%|██████████| 19/19 [00:02<00:00,  7.83it/s]

This motivated me to also use the scipy.sparse in meshtools.conformalized_mean_curvature_flow() in Step 1, but I ended up getting a new error message...

 0%|          | 0/50 [00:00<?, ?it/s]/home/dkobayas/LC-shadow/u-unwrap3D/unwrap3D/Mesh/meshtools.py:5305: MatrixRankWarning: Matrix is exactly singular
  U = spsolver.spsolve(S,b)
  0%|          | 0/50 [00:00<?, ?it/s]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 9
      1 # run cMCF on the mesh for 50 iterations, delta=5e-4 using the standard cotangent Laplacian 
      2 # Usteps_MCF, mesh_F, MCF_measures_dict = meshtools.conformalized_mean_curvature_flow(mesh, 
      3 #                                                 max_iter=50, 
   (...)
      7 #                                                 mollify_factor=1e-6,
      8 #                                                 solver='pardiso')  # this is a parameter used in the robust Laplacian
----> 9 Usteps_MCF, mesh_F, MCF_measures_dict = meshtools.conformalized_mean_curvature_flow(mesh, 
     10                                                 max_iter=50, 
     11                                                 delta=5e-4, 
     12                                                 conformalize=True, # set this flag to run cMCF, else it will run MCF
     13                                                 robust_L =False, # if set, runs MCF/cMCF using the robust Laplacian instead of cotangent Laplacian
     14                                                 mollify_factor=1e-6,
     15                                                 solver='other')  # this is a parameter used in the robust Laplacian
     17 # cMCF is run as a free-form deformation, if we want to have the mesh represent a gradual smoothing in accordance with the physical size of the cell and the image coordinates we need to put back the correct scaling factors.
     18 Usteps_MCF_img = meshtools.recover_img_coordinate_conformal_mean_flow(Usteps_MCF, MCF_measures_dict) # recovers the proper associated image coordinates 

File ~/LC-shadow/u-unwrap3D/unwrap3D/Mesh/meshtools.py:5355, in conformalized_mean_curvature_flow(mesh, max_iter, delta, rescale_output, min_diff, conformalize, robust_L, mollify_factor, solver)
   5352         else:
   5353             Usteps[...,ii+1] = U.copy() # copy the current iteration into it. 
-> 5355 c_all = np.vstack(c_all)
   5356 area_all = np.hstack(area_all)
   5357 Usteps = Usteps[..., :len(area_all)+1]

File /depot/cfrueh/apps/env_shadow/lib/python3.11/site-packages/numpy/_core/shape_base.py:291, in vstack(tup, dtype, casting)
    289 if not isinstance(arrs, tuple):
    290     arrs = (arrs,)
--> 291 return _nx.concatenate(arrs, 0, dtype=dtype, casting=casting)

ValueError: need at least one array to concatenate

Robust Laplacian (L_robust=True)

Lastly, I also tried the cMCF step (step1) with the keyword argument, L_robust=True. However, unfortunately, it caused a segmentaiton fault and crashed my JupyterLab notebook! In the first place, the robust-Laplacian package doesn't seem to be working in my environment. I tested a toy example using the Stanford bunny model but it also caused the segmentation fault. The error was reproduced for both MacOS and Linux, so the problem may be in the package itself. I created a new issue in their repository (nmwsharp/robust-laplacians-py#14 (comment)). Hopefully, I will get a response and find a fix in this problem...

[dkobayas-cyber]

The problem has been resolved. A root cause turned out to be in the libigl Python Bindings (igl). The function igl.cotmatrix() returned an empty Laplacian matrix, regardless of the input mesh model (e.g., Stanford Bunny model).

For some reason, the igl package didn't work as expected with Numpy version 2.1.0. After downgrading the Numpy to 1.26.4, the igl started to work and the cMCF calculation became successful in Step 1, even with remesh_samples=0.9.

[fyz11]

thanks for letting us know ! - that's a pretty wild error coming from igl and numpy. I wonder if this numpy also causes robust_laplacian to fail too..... I will need to test these and perhaps change the install requirements. Cheers!