from pybug.io import auto_import import numpy as np mesh = auto_import('/vol/atlas/databases/frgc/spring2003/04201d302.abs')[0] tris = mesh.trilist points = mesh.points print mesh def normalise(vec): # Avoid divisions by almost 0 numbers # np.spacing(1) is equivalent to Matlab's eps d = np.sqrt(np.sum(vec ** 2, axis=1)) d[d < np.spacing(1)] = 1.0 return vec / d[..., None] def py_compute_normal(vertex, face): nface = face.shape[0] nvert = vertex.shape[0] # Calculate the cross product (per-face normal) normalf = np.cross(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) normalf = normalise(normalf) # Calculate per-vertex normal normal = np.zeros([nvert, 3]) for i in xrange(nface): f = face[i, :] for j in xrange(3): normal[f[j], :] += normalf[i, :] # Normalize normal = normalise(normal) # Enforce that the normal are outward v = vertex - np.mean(vertex)[..., None] s = np.sum(v * normal, axis=1) if np.sum(np.greater(s, 0)) < np.sum(np.less(s, 0)): # flip normal = -normal normalf = -normalf return normal, normalf %timeit py_compute_normal(points, tris) %load_ext cythonmagic %%cython import numpy as np cimport cython cimport numpy as np %%cython import numpy as np cimport cython cimport numpy as np def my_pow(double x): return np.power(x, 2) print my_pow(2.0) %%cython import numpy as np cimport numpy as np cimport cython cdef np.ndarray[np.float64_t, ndim=2] cy_normalise_naive(np.ndarray[np.float64_t, ndim=2] vec): # Avoid divisions by almost 0 numbers cdef np.ndarray[np.float64_t, ndim=1] d = np.sqrt(np.sum(vec ** 2, axis=1)) d[d < np.spacing(1)] = 1.0 return vec / d[..., None] cpdef cy_compute_normal_naive(np.ndarray[np.float64_t, ndim=2] vertex, np.ndarray[int, ndim=2] face): cdef int nface = face.shape[0] cdef int nvert = vertex.shape[0] # Calculate the cross product (per-face normal) cdef np.ndarray[np.float64_t, ndim=2] normalf = np.cross(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) normalf = cy_normalise_naive(normalf) # Calculate per-vertex normal cdef np.ndarray[np.float64_t, ndim=2] normal = np.zeros([nvert, 3]) cdef np.ndarray[int, ndim=1] f for i in xrange(nface): f = face[i, :] for j in xrange(3): normal[f[j], :] += normalf[i, :] # Normalize normal = cy_normalise_naive(normal) # Enforce that the normal are outward cdef np.ndarray[np.float64_t, ndim=2] v = vertex - np.mean(vertex)[..., None] cdef np.ndarray[np.float64_t, ndim=1] s = np.sum(v * normal, axis=1) if np.sum(np.greater(s, 0)) < np.sum(np.less(s, 0)): # flip normal = -normal normalf = -normalf return normal, normalf %timeit cy_compute_normal_naive(points, tris) %%cython -a import numpy as np cimport numpy as np cimport cython cdef np.ndarray[np.float64_t, ndim=2] cy_normalise_naive(np.ndarray[np.float64_t, ndim=2] vec): # Avoid divisions by almost 0 numbers cdef np.ndarray[np.float64_t, ndim=1] d = np.sqrt(np.sum(vec ** 2, axis=1)) d[d < np.spacing(1)] = 1.0 return vec / d[..., None] cpdef cy_compute_normal_naive(np.ndarray[np.float64_t, ndim=2] vertex, np.ndarray[int, ndim=2] face): cdef int nface = face.shape[0] cdef int nvert = vertex.shape[0] # unit normals to the faces cdef np.ndarray[np.float64_t, ndim=2] normalf = np.cross(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) normalf = cy_normalise_naive(normalf) # unit normal to the vertex cdef np.ndarray[np.float64_t, ndim=2] normal = np.zeros([nvert, 3]) cdef double[:] f for i in xrange(nface): f = face[i, :] for j in xrange(3): normal[f[j], :] += normalf[i, :] # normalize normal = cy_normalise_naive(normal) # enforce that the normal are outward cdef np.ndarray[np.float64_t, ndim=2] v = vertex - np.mean(vertex)[..., None] cdef np.ndarray[np.float64_t, ndim=1] s = np.sum(v * normal, axis=1) if np.sum(np.greater(s, 0)) < np.sum(np.less(s, 0)): # flip normal = -normal normalf = -normalf return normal, normalf %%cython -a import numpy as np cimport numpy as np cimport cython cdef np.ndarray[np.float64_t, ndim=2] cy_normalise_naive(np.ndarray[np.float64_t, ndim=2] vec): # Avoid divisions by almost 0 numbers cdef np.ndarray[np.float64_t, ndim=1] d = np.sqrt(np.sum(vec ** 2, axis=1)) d[d < np.spacing(1)] = 1.0 return vec / d[..., None] cpdef cy_compute_normal_better(np.ndarray[np.float64_t, ndim=2] vertex, np.ndarray[int, ndim=2] face): cdef int nface = face.shape[0] cdef int nvert = vertex.shape[0] # unit normals to the faces cdef np.ndarray[np.float64_t, ndim=2] normalf = np.cross(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) normalf = cy_normalise_naive(normalf) # unit normal to the vertex cdef np.ndarray[np.float64_t, ndim=2] normal = np.zeros([nvert, 3]) cdef int f0, f1, f2 for i in range(nface): f0 = face[i, 0] f1 = face[i, 1] f2 = face[i, 2] for j in range(3): normal[f0, j] += normalf[i, j] normal[f1, j] += normalf[i, j] normal[f2, j] += normalf[i, j] # normalize normal = cy_normalise_naive(normal) # enforce that the normal are outward cdef np.ndarray[np.float64_t, ndim=2] v = vertex - np.mean(vertex)[..., None] cdef np.ndarray[np.float64_t, ndim=1] s = np.sum(v * normal, axis=1) if np.sum(np.greater(s, 0)) < np.sum(np.less(s, 0)): # flip normal = -normal normalf = -normalf return normal, normalf %timeit cy_compute_normal_better(points, tris) p = %prun -r cy_compute_normal_better(points, tris) p.print_stats() %%cython -a import numpy as np cimport numpy as np cimport cython cdef np.ndarray[np.float64_t, ndim=2] normalise(np.ndarray[np.float64_t, ndim=2] vec): # Avoid divisions by almost 0 numbers cdef np.ndarray[np.float64_t, ndim=1] d = np.sqrt(np.sum(vec ** 2, axis=1)) d[d < np.spacing(1)] = 1.0 return vec / d[..., None] cdef inline np.ndarray[np.float64_t, ndim=2] cross(double[:, :] x, double[:, :] y): cdef np.ndarray[np.float64_t, ndim=2] z = np.empty_like(x) cdef int n = x.shape[0] for i in range(n): z[i, 0] = x[i, 1] * y[i, 2] - x[i, 2] * y[i, 1] z[i, 1] = x[i, 2] * y[i, 0] - x[i, 0] * y[i, 2] z[i, 2] = x[i, 0] * y[i, 1] - x[i, 1] * y[i, 0] return z cpdef cy_compute_normal(np.ndarray[np.float64_t, ndim=2] vertex, np.ndarray[int, ndim=2] face): cdef int nface = face.shape[0] cdef int nvert = vertex.shape[0] # unit normals to the faces cdef np.ndarray[np.float64_t, ndim=2] normalf = cross(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) normalf = normalise(normalf) # unit normal to the vertex cdef np.ndarray[np.float64_t, ndim=2] normal = np.zeros([nvert, 3]) cdef int f0, f1, f2 for i in range(nface): f0 = face[i, 0] f1 = face[i, 1] f2 = face[i, 2] for j in range(3): normal[f0, j] += normalf[i, j] normal[f1, j] += normalf[i, j] normal[f2, j] += normalf[i, j] # normalize normal = normalise(normal) # enforce that the normal are outward cdef np.ndarray[np.float64_t, ndim=2] v = vertex - np.mean(vertex)[..., None] cdef np.ndarray[np.float64_t, ndim=1] s = np.sum(v * normal, axis=1) if np.sum(np.greater(s, 0)) < np.sum(np.less(s, 0)): # flip normal = -normal normalf = -normalf return normal, normalf %timeit cy_compute_normal(points, tris) print "" p = %prun -r cy_compute_normal(points, tris) p.print_stats() cy_normal, cy_normalf = cy_compute_normal(points, tris) py_normal, py_normalf = py_compute_normal(points, tris) print np.allclose(cy_normal, py_normal) print np.allclose(cy_normalf, py_normalf) import numba from numba.decorators import autojit, jit from math import sqrt # Had to unroll this otherwise it complained about some strange python object # coercion error @autojit def numba_normalise(vec): # Avoid divisions by almost 0 numbers # np.spacing(1) is equivalent to Matlab's eps n = vec.shape[0] for i in range(n): d = sqrt(vec[i, 0] * vec[i, 0] + vec[i, 1] * vec[i, 1] + vec[i, 2] * vec[i, 2]) if d < np.spacing(1): d = 1.0 vec[i, 0] /= d vec[i, 1] /= d vec[i, 2] /= d # If I didn't roll my own cross product then computing # the normals actually takes LONGER than the pure Python implementation @jit(argtypes=(numba.double[:, :], numba.double[:, :])) def cross_numba(x, y): output = np.empty_like(x) n = x.shape[0] for i in range(n): output[i, 0] = x[i, 1] * y[i, 2] - x[i, 2] * y[i, 1] output[i, 1] = x[i, 2] * y[i, 0] - x[i, 0] * y[i, 2] output[i, 2] = x[i, 0] * y[i, 1] - x[i, 1] * y[i, 0] return output @jit(argtypes=(numba.double[:, :], numba.int32[:, :])) def numba_compute_normal(vertex, face): nface = face.shape[0] # Calculate the cross product (per-face normal) normalf = cross_numba(vertex[face[:, 1], :] - vertex[face[:, 0], :], vertex[face[:, 2], :] - vertex[face[:, 0], :]) numba_normalise(normalf) # Calculate per-vertex normal normal = np.zeros_like(vertex) for i in range(nface): f = face[i, :] for j in range(3): normal[f[j], :] += normalf[i, :] # Normalize numba_normalise(normal) # Enforce that the normal are outward v = vertex - np.mean(vertex)[..., None] s = np.sum(v * normal, axis=1) s_gt_sum = 0 s_lt_sum = 0 # Had to expand this loop otherwise numba complained: # 'only length-1 arrays can be converted to Python scalars' # On the calls to np.greater and np.less for i in range(s.shape[0]): if s[i] > 0: s_gt_sum += 1 elif s[i] < 0: s_lt_sum += 1 if s_gt_sum < s_lt_sum: # flip normal = -normal normalf = -normalf return normal, normalf %timeit numba_compute_normal(points, tris) numba_normal, numba_normalf = numba_compute_normal(points, tris) py_normal, py_normalf = py_compute_normal(points, tris) print np.allclose(numba_normal, py_normal) print np.allclose(numba_normalf, py_normalf)