OpenMesh version of Nathan Reed and lemon:
import numpy as np
import openmesh as om
from vedo import *
def array_angle(array1, array2):
"""
INPUT: N x 2/3/...
"""
res = np.sum(array1 * array2, axis=1)
res /= np.linalg.norm(array1, axis=1)
res /= np.linalg.norm(array2, axis=1)
res = np.clip(res, -1.0, 1.0)
res = np.arccos(res)
return res
def run(path_mesh):
mesh_om = om.read_trimesh(path_mesh, vertex_normal=True)
v = mesh_om.points()
f = mesh_om.face_vertex_indices()
n = mesh_om.vertex_normals()
if True:
"""
https://computergraphics.stackexchange.com/a/1719/17639
"""
e = mesh_om.ev_indices()
ev = v[e[:, 1]] - v[e[:, 0]]
en = n[e[:, 1]] - n[e[:, 0]]
e_curv = np.sum(ev * en, axis=1) / np.linalg.norm(ev, axis=1)
v_curv = np.zeros(v.shape[0])
for idx_e in range(e.shape[0]):
v_curv[e[idx_e, 0]] += e_curv[idx_e]
v_curv[e[idx_e, 1]] += e_curv[idx_e]
else:
"""
https://blender.stackexchange.com/a/147371/82691
"""
vf_indices = mesh_om.vf_indices()
ev_indices = mesh_om.ev_indices()
fe_indices = mesh_om.fe_indices()
ev_01 = v[ev_indices[:, 1]] - v[ev_indices[:, 0]]
en_01 = n[ev_indices[:, 1]] - n[ev_indices[:, 0]]
e_curv = np.sum(ev_01 * en_01, axis=1) / np.linalg.norm(ev_01, axis=1)
e_01 = v[f[:, 1]] - v[f[:, 0]]
e_12 = v[f[:, 2]] - v[f[:, 1]]
e_20 = v[f[:, 0]] - v[f[:, 2]]
f_angles = np.zeros_like(f).astype(np.float32)
f_angles[:, 0] = array_angle(e_01, -e_20)
f_angles[:, 1] = array_angle(-e_01, e_12)
f_angles[:, 2] = array_angle(e_20, -e_12)
v_curv = np.zeros(v.shape[0])
for idx_v in range(v.shape[0]):
va = 0
for i in range(vf_indices.shape[1]):
idx_f = vf_indices[idx_v, i]
if idx_f == -1:
break
idx_v_in_f = np.where(f[idx_f] == idx_v)[0][0]
a = f_angles[idx_f, idx_v_in_f]
va += a
idx_e0 = fe_indices[idx_f, (0+idx_v_in_f) % 3]
idx_e2 = fe_indices[idx_f, (2+idx_v_in_f) % 3]
v_curv[idx_v] += a * (e_curv[idx_e0] + e_curv[idx_e2])
v_curv[idx_v] /= (2 * va)
if True:
#! vis
mean_v_curv = np.mean(v_curv)
std_v_curv = np.std(v_curv)
k_min = mean_v_curv - 3 * std_v_curv
k_max = mean_v_curv + 3 * std_v_curv
v_curv = np.clip(v_curv, k_min, k_max)
M = Mesh([v, f])
M.cmap('jet', v_curv, vmin=k_min, vmax=k_max, on="points").addScalarBar()
M.show(resetcam=True)
if name == "main":
path_mesh = "/home/lab9/Mesh/bunny.obj"
run(path_mesh)
When reading with Openmesh, the following message appears, indicating that the mesh may be flawed. At this time, the values of vertices and faces read by Openmesh and Blender may be different, and preprocessing may be required.
PolyMeshT::add_face: complex vertex
PolyMeshT::add_face: complex edge
PolyMeshT::add_face: patch re-linking failed
PolyMeshT::add_face: complex vertex is usually caused by the "wrong" faces, so we just need to delete some of the faces.
One possible solution is to use igl to read the mesh(keep the content of the file) with texcoord and normal(excluded here), and use openmesh to form the mesh, which will automatically delete the bad faces, we can calculate the mask and select valid f/ft/fn from the original array.
import numpy as np
import igl # conda install -c conda-forge igl
import openmesh as om
def write_mesh(path_mesh, v, vc=None, vt=None, vn=None, f=None, ft=None, fn=None, str_mtllib="", str_usemtl=""):
assert v.ndim == 2 and v.shape[1] == 3, ("\033[1;31mv.shape=%s\n\033[0m") % str(v.shape)
if vc is not None and len(vc) == 0:
vc = None #! vc = []
if vt is not None and len(vt) == 0:
vt = None
if vn is not None and len(vn) == 0:
vn = None
if ft is not None and len(ft) == 0:
ft = None
if fn is not None and len(fn) == 0:
fn = None
with open(path_mesh, "w") as fp:
if str_mtllib != "":
fp.write("mtllib " + str_mtllib + "\n")
if vc is not None:
v = np.concatenate([v, vc / 255.0], 1).reshape(-1, 6)
fp.write(("v {:f} {:f} {:f} {:.3f} {:.3f} {:.3f}\n" * v.shape[0]).format(*v.reshape(-1)))
else:
fp.write(("v {:f} {:f} {:f}\n" * v.shape[0]).format(*v.reshape(-1)))
if vt is not None:
fp.write(("vt {:f} {:f}\n" * vt.shape[0]).format(*vt.reshape(-1)))
if vn is not None:
fp.write(("vn {:f} {:f} {:f}\n" * vn.shape[0]).format(*vn.reshape(-1)))
if str_usemtl != "":
fp.write("usemtl " + str_usemtl + "\n")
if (ft is not None) and (fn is not None) and (f is not None):
num_f = f.shape[0]
# f v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...
for idx_f in range(num_f):
fp.write(("f {:d}/{:d}/{:d} {:d}/{:d}/{:d} {:d}/{:d}/{:d}\n").format(
f[idx_f, 0] + 1, ft[idx_f, 0] + 1, fn[idx_f, 0] + 1,
f[idx_f, 1] + 1, ft[idx_f, 1] + 1, fn[idx_f, 1] + 1,
f[idx_f, 2] + 1, ft[idx_f, 2] + 1, fn[idx_f, 2] + 1,
))
elif (ft is None) and (fn is not None) and (f is not None):
num_f = f.shape[0]
# f v1//vn1 v2//vn2 v3//vn3 ...
for idx_f in range(num_f):
fp.write(("f {:d}//{:d} {:d}//{:d} {:d}//{:d}\n").format(
f[idx_f, 0] + 1, fn[idx_f, 0] + 1,
f[idx_f, 1] + 1, fn[idx_f, 1] + 1,
f[idx_f, 2] + 1, fn[idx_f, 2] + 1,
))
elif (ft is not None) and (fn is None) and (f is not None):
num_f = f.shape[0]
# f v1/vt1 v2/vt2 v3/vt3 ...
for idx_f in range(num_f):
fp.write(("f {:d}/{:d} {:d}/{:d} {:d}/{:d}\n").format(
f[idx_f, 0] + 1, ft[idx_f, 0] + 1,
f[idx_f, 1] + 1, ft[idx_f, 1] + 1,
f[idx_f, 2] + 1, ft[idx_f, 2] + 1,
))
elif (ft is None) and (fn is None) and (f is not None):
# f v1 v2 v3 ....
fp.write(("f {:d} {:d} {:d}\n" * f.shape[0]).format(*(f.reshape(-1) + 1)))
if name == "main":
path_mesh = "XXX.obj"
path_mesh_iglom = path_mesh[:-4] + "_IGLOM.obj"
v_igl, vt_igl, _, f_igl, ft_igl, _ = igl.read_obj(path_mesh)
#! check
if True:
dbl_area = igl.doublearea(v_igl, f_igl)
f_mask0 = (dbl_area > 2e-9).astype(np.bool_) # ! 2e-9
f_igl = f_igl[f_mask0]
if len(vt_igl) != 0:
ft_igl = ft_igl[f_mask0]
f_mask = np.zeros(f_igl.shape[0]).astype(np.bool_)
mesh_iglom = om.PolyMesh(v_igl, f_igl)
v_iglom = mesh_iglom.points()
f_iglom = mesh_iglom.face_vertex_indices()
offset_f_igl = 0
for idx_f_iglom in range(f_iglom.shape[0]):
is_true = True
while is_true:
is_equ = True
is_equ &= f_igl[offset_f_igl, 0] == f_iglom[idx_f_iglom, 0]
is_equ &= f_igl[offset_f_igl, 1] == f_iglom[idx_f_iglom, 1]
is_equ &= f_igl[offset_f_igl, 2] == f_iglom[idx_f_iglom, 2]
if is_equ:
f_mask[offset_f_igl] = True
offset_f_igl += 1
is_true = False
else:
offset_f_igl += 1
if len(vt_igl) != 0:
ft_iglom = ft_igl[f_mask]
if len(vt_igl) != 0:
write_mesh(path_mesh_iglom, v=v_iglom, vt=vt_igl, f=f_iglom, ft=ft_iglom)
else:
write_mesh(path_mesh_iglom, v=v_iglom, f=f_iglom)
After:
bpy.opscommands used is scary. You might attract a more proficient audience if you rewrote the codebase using thebmeshmodule. (However, I don't know if all commands [bisecte.g.] are directly translatable.) – Leander Aug 01 '19 at 13:02bpy.opsis actually not defining or really relevant to the question, since this is a question which could be dealt with in a manner which is entirely independent from the surface shape. – Robert Roth Aug 01 '19 at 13:13discrete_mean_curvature_measure(mesh, points, radius)https://github.com/mikedh/trimesh/blob/master/trimesh/curvature.py – Robert Roth Aug 02 '19 at 07:07bpy.opsis bad? – Mike Aug 02 '19 at 07:24bpy.ops. Most coders write blender scripts without bpy.ops, since it is slower, context-sensitive and more difficult to read. This answer provides some generalized explanation. Every bpy.ops triggers a scene update. – Leander Aug 02 '19 at 07:29