1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# Copyright (c) 2014, Nicolas P. Rougier
# Distributed under the (new) BSD License. See LICENSE.txt for more info.
# -----------------------------------------------------------------------------
import numpy as np
def compact(vertices, indices, tolerance=1e-3):
""" Compact vertices and indices within given tolerance """
# Transform vertices into a structured array for np.unique to work
n = len(vertices)
V = np.zeros(n, dtype=[("pos", np.float32, 3)])
V["pos"][:, 0] = vertices[:, 0]
V["pos"][:, 1] = vertices[:, 1]
V["pos"][:, 2] = vertices[:, 2]
epsilon = 1e-3
decimals = int(np.log(epsilon)/np.log(1/10.))
# Round all vertices within given decimals
V_ = np.zeros_like(V)
X = V["pos"][:, 0].round(decimals=decimals)
X[np.where(abs(X) < epsilon)] = 0
V_["pos"][:, 0] = X
Y = V["pos"][:, 1].round(decimals=decimals)
Y[np.where(abs(Y) < epsilon)] = 0
V_["pos"][:, 1] = Y
Z = V["pos"][:, 2].round(decimals=decimals)
Z[np.where(abs(Z) < epsilon)] = 0
V_["pos"][:, 2] = Z
# Find the unique vertices AND the mapping
U, RI = np.unique(V_, return_inverse=True)
# Translate indices from original vertices into the reduced set (U)
indices = indices.ravel()
I_ = indices.copy().ravel()
for i in range(len(indices)):
I_[i] = RI[indices[i]]
I_ = I_.reshape(len(indices)/3, 3)
# Return reduced vertices set, transalted indices and mapping that allows
# to go from U to V
return U.view(np.float32).reshape(len(U), 3), I_, RI
def normals(vertices, indices):
"""
Compute normals over a triangulated surface
Parameters
----------
vertices : ndarray (n,3)
triangles vertices
indices : ndarray (p,3)
triangles indices
"""
# Compact similar vertices
vertices, indices, mapping = compact(vertices, indices)
T = vertices[indices]
N = np.cross(T[:, 1] - T[:, 0], T[:, 2]-T[:, 0])
L = np.sqrt(np.sum(N * N, axis=1))
L[L == 0] = 1.0 # prevent divide-by-zero
N /= L[:, np.newaxis]
normals = np.zeros_like(vertices)
normals[indices[:, 0]] += N
normals[indices[:, 1]] += N
normals[indices[:, 2]] += N
L = np.sqrt(np.sum(normals*normals, axis=1))
L[L == 0] = 1.0
normals /= L[:, np.newaxis]
return normals[mapping]
|
