Step 2: Quasi-conformal spherical parameterization of the genus-0 reference surface, \(S_{ref}(x,y,z)\)#
Theoretically, any genus-0 3D surface can be mapped to a sphere. Practically, there is no unique way to do this. Various ways can be found to do this in the literature. What is good depends on the downstream application and whether we care about speed and size of mesh (i.e. the number of vertices and faces).
u-Unwrap3D implements a fast iterative relaxation algorithm for quasi-conformal spherical parameterization. This does not give the lowest conformal error however it results for better triangle quality to enable area distortion relaxation in the next step to get a quasi-equiareal parameterization.
Note: u-Unwrap3D also provides implementation of the spherical quasiconformal map of Choi et al. which guarantees a bijective mapping with no foldovers, and has lower conformal error but we find yields poorer triangle quality (smaller area triangles and smaller minimum angle), and is very slow, scaling poorly for larger meshes (>100k vertices).
1. Load genus-0 \(S_{ref}(x,y,z)\) mesh#
We assume the user has worked through step 1 which generates and saves a genus-0 reference surface, \(S_{ref}(x,y,z)\) for an input cell surface mesh, \(S(x,y,z)\) to the folder example_results/bleb_example/step1_cMCF_reference. Moreover this surface has been colored by the mean curvature of \(S(x,y,z)\)
import unwrap3D.Utility_Functions.file_io as fio
import unwrap3D.Mesh.meshtools as meshtools
import numpy as np
import os
import scipy.io as spio
import igl
# example cell used
imgfolder = '../../data/img'
imgfile = os.path.join(imgfolder, 'bleb_example.tif')
basefname = os.path.split(imgfile)[-1].split('.tif')[0] # get the filename with extension
# create the analysis save folder for this step
savefolder = os.path.join('example_results',
basefname,
'step2_conformal_sphere_param')
fio.mkdir(savefolder) # auto generates the specified folder structure if doesn't currently exist.
# load the pre-computed genus-0 reference surface mesh for the cell
S_ref_folder = 'example_results/%s/step1_cMCF_reference' %(basefname)
S_ref_file = os.path.join(S_ref_folder,
'unwrap_cMCF_Sref_mesh_H_color.obj')
mesh = meshtools.read_mesh(S_ref_file)
2. Conformal mapping of \(S_{ref}(x,y,z)\) onto the 3D unit sphere to get \(S_{\text{ref}}^{\mathcal{Q}}(x,y,z)\)#
We will perform the spherical mapping using the i) recommended iterative approach and ii) the quasiconformal mapping of Choi et al. to show the difference in timing and in the deformation errors.
i) Iterative conformal spherical map (preferred if performing equiareal parameterization, and large meshes)
This method first directly projects the mesh onto the unit sphere by normalizing the vertex displacement relative to the mesh centroid \((v-c)\), i.e. \(v\leftarrow \frac{v-c}{|v-c|}\). We then apply cMCF iteratively to ‘untangle’ the mesh directly on the sphere using the uniform Laplacian, and the mass matrix from the robust Laplacian construction. The result is always a spherical map. A valid bijective spherical map is one with no foldovers i.e no inverted faces. For most shapes, we get fast convergence within 5 iterations.
ii) Direct quasiconformal spherical map of Choi et al.
This method leverages stereographic projection and attempts to find a quasiconformal map by solving the Beltrami equation. Not robust for triangle quality. Will error for non genus-0 input, is memory-intensive and slow. Result can be more challenging to apply area-distortion relaxation to (more parameter tuning). Use for simpler or more downsamples meshes, or when the objective is to get a lower quasiconformal error.
# explicitly expose the vertex and faces.
v = mesh.vertices.copy()
f = mesh.faces.copy()
# import time module for timing
import time
"""
1) Iterative conformal spherical map (*recommended default*, faster, more robust, better triangle quality)
"""
t1_iterative = time.time()
sphere_xyz_iterative, n_inverted_faces = meshtools.iterative_tutte_spherical_map(v,f,
deltaL = 5e-3, # like the cMCF, controls speed of convergence.
min_iter=5, # return as soon as iteration number > min_iter, and number of inverted faces is 0
max_iter=25, # run this many iterations to see if number of inverted reduces to 0
mollify_factor=1e-5,
scale=1.) # gives just the vertices.
assert(n_inverted_faces[-1]==0) # the output takes the last iteration.
t2_iterative = time.time()
print('printed number of inverted faces per iteration')
print(n_inverted_faces)
"""
2) Direct quasiconformal spherical map of Choi et al. (only for obtaining lower conformal error)
"""
t1_direct = time.time()
sphere_xyz_direct = meshtools.direct_spherical_conformal_map(v,f) # gives just the vertices.
t2_direct = time.time()
print('recommended iterative conformal: ', t2_iterative-t1_iterative)
print('Choi et al. direct conformal: ', t2_direct-t1_direct)
printed number of inverted faces per iteration
[0, 0, 0, 0, 0, 0]
recommended iterative conformal: 4.810135841369629
Choi et al. direct conformal: 14.325770616531372
# create mesh for iterative
sphere_mesh_iterative = mesh.copy() # this helps to copy the existing colormap
sphere_mesh_iterative.vertices = sphere_xyz_iterative.copy()
tmp = sphere_mesh_iterative.export(os.path.join(savefolder, 'iterative-conformal_sphere_param_Sref_%s.obj' %(basefname)))
triangle_props_iterative = meshtools.measure_triangle_props(sphere_mesh_iterative, clean=True)
print('triangle props iterative:')
print(triangle_props_iterative)
print('=========================')
# create mesh for Choi direct method
sphere_mesh_direct = mesh.copy() # this helps to copy the existing colormap
sphere_mesh_direct.vertices = sphere_xyz_direct.copy()
tmp = sphere_mesh_direct.export(os.path.join(savefolder, 'direct-conformal_sphere_param_Sref_%s.obj' %(basefname)))
triangle_props_direct = meshtools.measure_triangle_props(sphere_mesh_direct, clean=True)
print('triangle props direct:')
print(triangle_props_direct)
print('=========================')
triangle props iterative:
{'min_angle': 37.72755825012703, 'avg_angle': 60.00000000000001, 'max_angle': 96.77951462608578, 'std_dev_angle': 4.602416979613976, 'min_quality': 0.7544342656414823, 'avg_quality': 0.9904537832406153, 'max_quality': 0.9999998443829548, 'quality': array([0.99022527, 0.98870827, 0.98488669, ..., 0.98583723, 0.9777652 ,
0.99995536]), 'angles': array([[0.99094327, 0.98729905, 1.16335033],
[1.16315121, 0.9516156 , 1.02682584],
[0.9635418 , 0.98645905, 1.1915918 ],
...,
[1.18472648, 1.00075562, 0.95611055],
[0.92452656, 0.99967113, 1.21739496],
[1.03953202, 1.05172123, 1.0503394 ]])}
=========================
triangle props direct:
{'min_angle': 30.930843003120742, 'avg_angle': 60.00000000000001, 'max_angle': 98.11151126528594, 'std_dev_angle': 7.829824117694427, 'min_quality': 0.6931658075816856, 'avg_quality': 0.9728591140223615, 'max_quality': 0.9999983516674269, 'quality': array([0.99567969, 0.99565083, 0.99566051, ..., 0.96767302, 0.91652935,
0.90920887]), 'angles': array([[1.09158677, 0.97249596, 1.07750992],
[1.07778518, 1.09157683, 0.97223065],
[1.09195766, 0.97236561, 1.07726938],
...,
[1.01109782, 0.88512555, 1.24536928],
[0.74086023, 1.09982487, 1.30090755],
[1.08785054, 1.32211514, 0.73162698]])}
=========================
"""
Plotting using matplotlib for comparison. The final version should have vertices more evenly distributed over the sphere.
"""
import unwrap3D.Visualisation.plotting as plotting
import pylab as plt
sampling = 1 # plot every just so its not fully dense in the plot so we don't see anything!
# plotting the before
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(projection='3d')
plt.title('Iterative spherical parameterization')
ax.set_box_aspect(aspect = (1,1,1))
ax.scatter(sphere_mesh_iterative.vertices[::sampling,2],
sphere_mesh_iterative.vertices[::sampling,1],
sphere_mesh_iterative.vertices[::sampling,0],
s=0.5,
c=sphere_mesh_iterative.visual.vertex_colors[::sampling,:3]/255.)
ax.view_init(0, 0)
plotting.set_axes_equal(ax)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
# plotting the after
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(projection='3d')
plt.title('Direct spherical parameterization')
ax.set_box_aspect(aspect = (1,1,1))
ax.scatter(sphere_mesh_direct.vertices[::sampling,2],
sphere_mesh_direct.vertices[::sampling,1],
sphere_mesh_direct.vertices[::sampling,0],
s=0.5,
c=sphere_mesh_direct.visual.vertex_colors[::sampling,:3]/255.)
ax.view_init(0, 0)
plotting.set_axes_equal(ax)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
3. Quantifying and colormapping the conformal and equiareal distortion error of the spherical parameterization#
We compute the conformal and equiareal errors for the iterative spherical parameterization to be used in the next step. We will also show how to map this onto the mesh.
import unwrap3D.Visualisation.colors as vol_colors
from matplotlib import cm
"""
Conformal error
"""
conformal_error = meshtools.quasi_conformal_error(v,
sphere_xyz_iterative,
f)
# color the error
conformal_error_colors = vol_colors.get_colors(conformal_error[1], colormap=cm.jet, vmin=1., vmax=1.5)
sphere_conformal_error = meshtools.create_mesh(vertices=sphere_xyz_iterative,
faces=f,
face_colors=np.uint8(255*conformal_error_colors[...,:3]))
tmp = sphere_conformal_error.export(os.path.join(savefolder,
'conformal_sphere_param_Sref_%s_conformal-error.obj' %(basefname)))
"""
Area distortion error
"""
area_distortion_error = meshtools.area_distortion_measure(v,
sphere_xyz_iterative,
f)
area_distortion_error_colors = vol_colors.get_colors(np.log(area_distortion_error),
colormap=cm.coolwarm,
vmin=-3,
vmax=3)
sphere_areal_error = sphere_conformal_error.copy()
sphere_areal_error.visual.face_colors = np.uint8(255*area_distortion_error_colors[...,:3]).copy()
tmp = sphere_areal_error.export(os.path.join(savefolder,
'conformal_sphere_param_Sref_%s_log_equiareal-error.obj' %(basefname)))
# save the statistics
spio.savemat(os.path.join(savefolder,
'conformal_sphere-param_conformal_area_errors.mat'),
{'conformal_Jac_eigvals': conformal_error[0],
'conformal_stretch_factors': conformal_error[1],
'mean_conformal_stretch_factors': conformal_error[2],
'conformal_tri_area_pairs': conformal_error[3],
'conformal_area_distortion_factor': area_distortion_error})
We now also compute the errors for the Choi direct quasiconformal map for comparison.
mean_conformal_error_iterative = conformal_error[2] # ideal is 1
mean_equiareal_error_iterative = np.nanmean(area_distortion_error) # ideal is 1
print('quasiconformal error of iterative method:', mean_conformal_error_iterative)
print('mean equiareal error of iterative method:', mean_equiareal_error_iterative)
print('-----------------------------------------')
# compute mean conformal and equiareal for direct method
conformal_error_direct = meshtools.quasi_conformal_error(v,
sphere_xyz_direct,
f)
area_distortion_error_direct = meshtools.area_distortion_measure(v,
sphere_xyz_direct,
f)
mean_conformal_error_direct = conformal_error_direct[2] # ideal is 1
mean_equiareal_error_direct = np.nanmean(area_distortion_error_direct) # ideal is 1
print('quasiconformal error of direct method:', mean_conformal_error_direct)
print('mean equiareal error of direct method:', mean_equiareal_error_direct)
print('-----------------------------------------')
quasiconformal error of iterative method: 1.1873036453396435
mean equiareal error of iterative method: 1.4135101985988752
-----------------------------------------
quasiconformal error of direct method: 1.0192368494880721
mean equiareal error of direct method: 1.5193297616244414
-----------------------------------------
The direct quasiconformal map achieves lower conformal error but at the expense of greater area distortion.