[yt-svn] commit/yt: MatthewTurk: Merged in atmyers/yt (pull request #1994)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Feb 24 09:37:20 PST 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/d18dc7748e19/
Changeset: d18dc7748e19
Branch: yt
User: MatthewTurk
Date: 2016-02-24 17:37:11+00:00
Summary: Merged in atmyers/yt (pull request #1994)
Pure Cython Ray Tracer for Triangle Meshes
Affected #: 5 files
diff -r aaa54f6b1a1b1f98ea52b203fccce7c51a1ccc91 -r d18dc7748e19ec77816ee9b44147e5c0597c3185 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -30,6 +30,7 @@
yt/utilities/lib/amr_kdtools.c
yt/utilities/lib/basic_octree.c
yt/utilities/lib/bitarray.c
+yt/utilities/lib/bounding_volume_hierarchy.c
yt/utilities/lib/contour_finding.c
yt/utilities/lib/depth_first_octree.c
yt/utilities/lib/element_mappings.c
diff -r aaa54f6b1a1b1f98ea52b203fccce7c51a1ccc91 -r d18dc7748e19ec77816ee9b44147e5c0597c3185 setup.py
--- a/setup.py
+++ b/setup.py
@@ -122,6 +122,11 @@
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=["m"], depends=["yt/utilities/lib/bitarray.pxd"]),
+ Extension("yt.utilities.lib.bounding_volume_hierarchy",
+ ["yt/utilities/lib/bounding_volume_hierarchy.pyx"],
+ extra_compile_args=omp_args,
+ extra_link_args=omp_args,
+ libraries=["m"], depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd"]),
Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],
include_dirs=["yt/utilities/lib/",
diff -r aaa54f6b1a1b1f98ea52b203fccce7c51a1ccc91 -r d18dc7748e19ec77816ee9b44147e5c0597c3185 yt/utilities/lib/bounding_volume_hierarchy.pxd
--- /dev/null
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pxd
@@ -0,0 +1,50 @@
+cimport cython
+import numpy as np
+cimport numpy as np
+
+# ray data structure
+cdef struct Ray:
+ np.float64_t origin[3]
+ np.float64_t direction[3]
+ np.float64_t inv_dir[3]
+ np.float64_t data_val
+ np.float64_t t_near
+ np.float64_t t_far
+ np.int64_t elem_id
+
+# axis-aligned bounding box
+cdef struct BBox:
+ np.float64_t left_edge[3]
+ np.float64_t right_edge[3]
+
+# node for the bounding volume hierarchy
+cdef struct BVHNode:
+ np.int64_t begin
+ np.int64_t end
+ BVHNode* left
+ BVHNode* right
+ BBox bbox
+
+# triangle data structure
+cdef struct Triangle:
+ np.float64_t p0[3]
+ np.float64_t p1[3]
+ np.float64_t p2[3]
+ np.float64_t d0, d1, d2
+ np.int64_t elem_id
+ np.float64_t centroid[3]
+ BBox bbox
+
+cdef class BVH:
+ cdef BVHNode* root
+ cdef Triangle* triangles
+ cdef np.int64_t leaf_size
+ cdef np.float64_t[:, ::1] vertices
+ cdef np.int64_t _partition(self, np.int64_t begin, np.int64_t end,
+ np.int64_t ax, np.float64_t split) nogil
+ cdef void intersect(self, Ray* ray) nogil
+ cdef void _get_node_bbox(self, BVHNode* node,
+ np.int64_t begin, np.int64_t end) nogil
+ cdef void _recursive_intersect(self, Ray* ray, BVHNode* node) nogil
+ cdef BVHNode* _recursive_build(self, np.int64_t begin, np.int64_t end) nogil
+ cdef void _recursive_free(self, BVHNode* node) nogil
diff -r aaa54f6b1a1b1f98ea52b203fccce7c51a1ccc91 -r d18dc7748e19ec77816ee9b44147e5c0597c3185 yt/utilities/lib/bounding_volume_hierarchy.pyx
--- /dev/null
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -0,0 +1,320 @@
+cimport cython
+import numpy as np
+cimport numpy as np
+from libc.math cimport fabs, fmax, fmin
+from libc.stdlib cimport malloc, free
+from cython.parallel import parallel, prange
+
+cdef extern from "mesh_construction.h":
+ enum:
+ MAX_NUM_TRI
+ int triangulate_hex[MAX_NUM_TRI][3]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline np.float64_t dot(const np.float64_t a[3],
+ const np.float64_t b[3]) nogil:
+ return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void cross(const np.float64_t a[3],
+ const np.float64_t b[3],
+ np.float64_t c[3]) nogil:
+ c[0] = a[1]*b[2] - a[2]*b[1]
+ c[1] = a[2]*b[0] - a[0]*b[2]
+ c[2] = a[0]*b[1] - a[1]*b[0]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void subtract(const np.float64_t a[3],
+ const np.float64_t b[3],
+ np.float64_t c[3]) nogil:
+ c[0] = a[0] - b[0]
+ c[1] = a[1] - b[1]
+ c[2] = a[2] - b[2]
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_triangle_intersect(Ray* ray, const Triangle* tri) nogil:
+# https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+
+ # edge vectors
+ cdef np.float64_t e1[3]
+ cdef np.float64_t e2[3]
+ subtract(tri.p1, tri.p0, e1)
+ subtract(tri.p2, tri.p0, e2)
+
+ cdef np.float64_t P[3]
+ cross(ray.direction, e2, P)
+
+ cdef np.float64_t det, inv_det
+ det = dot(e1, P)
+ if(det > -1.0e-10 and det < 1.0e-10):
+ return False
+ inv_det = 1.0 / det
+
+ cdef np.float64_t T[3]
+ subtract(ray.origin, tri.p0, T)
+
+ cdef np.float64_t u = dot(T, P) * inv_det
+ if(u < 0.0 or u > 1.0):
+ return False
+
+ cdef np.float64_t Q[3]
+ cross(T, e1, Q)
+
+ cdef np.float64_t v = dot(ray.direction, Q) * inv_det
+ if(v < 0.0 or u + v > 1.0):
+ return False
+
+ cdef np.float64_t t = dot(e2, Q) * inv_det
+
+ if(t > 1.0e-10 and t < ray.t_far):
+ ray.t_far = t
+ ray.data_val = (1.0 - u - v)*tri.d0 + u*tri.d1 + v*tri.d2
+ ray.elem_id = tri.elem_id
+ return True
+
+ return False
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef np.int64_t ray_bbox_intersect(Ray* ray, const BBox bbox) nogil:
+# https://tavianator.com/fast-branchless-raybounding-box-intersections/
+
+ cdef np.float64_t tmin = -1.0e300
+ cdef np.float64_t tmax = 1.0e300
+
+ cdef np.int64_t i
+ cdef np.float64_t t1, t2
+ for i in range(3):
+ t1 = (bbox.left_edge[i] - ray.origin[i])*ray.inv_dir[i]
+ t2 = (bbox.right_edge[i] - ray.origin[i])*ray.inv_dir[i]
+ tmin = fmax(tmin, fmin(t1, t2))
+ tmax = fmin(tmax, fmax(t1, t2))
+
+ return tmax >= fmax(tmin, 0.0)
+
+
+cdef class BVH:
+ '''
+
+ This class implements a bounding volume hierarchy (BVH), a spatial acceleration
+ structure for fast ray-tracing. A BVH is like a kd-tree, except that instead of
+ partitioning the *volume* of the parent to create the children, we partition the
+ triangles themselves into 'left' or 'right' sub-trees. The bounding volume for a
+ node is then determined by computing the bounding volume of the triangles that
+ belong to it. This allows us to quickly discard triangles that are not close
+ to intersecting a given ray.
+
+ '''
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __cinit__(self,
+ np.float64_t[:, ::1] vertices,
+ np.int64_t[:, ::1] indices,
+ np.float64_t[:, ::1] field_data):
+
+ self.leaf_size = 16
+ self.vertices = vertices
+ cdef np.int64_t num_elem = indices.shape[0]
+ cdef np.int64_t num_tri = 12*num_elem
+
+ # fill our array of triangles
+ cdef np.int64_t i, j, k
+ cdef np.int64_t offset, tri_index
+ cdef np.int64_t v0, v1, v2
+ cdef Triangle* tri
+ self.triangles = <Triangle*> malloc(num_tri * sizeof(Triangle))
+ for i in range(num_elem):
+ offset = 12*i
+ for j in range(12):
+ tri_index = offset + j
+ tri = &(self.triangles[tri_index])
+ tri.elem_id = i
+ v0 = indices[i][triangulate_hex[j][0]]
+ v1 = indices[i][triangulate_hex[j][1]]
+ v2 = indices[i][triangulate_hex[j][2]]
+ tri.d0 = field_data[i][triangulate_hex[j][0]]
+ tri.d1 = field_data[i][triangulate_hex[j][1]]
+ tri.d2 = field_data[i][triangulate_hex[j][2]]
+ for k in range(3):
+ tri.p0[k] = vertices[v0][k]
+ tri.p1[k] = vertices[v1][k]
+ tri.p2[k] = vertices[v2][k]
+ tri.centroid[k] = (tri.p0[k] + tri.p1[k] + tri.p2[k]) / 3.0
+ tri.bbox.left_edge[k] = fmin(fmin(tri.p0[k], tri.p1[k]), tri.p2[k])
+ tri.bbox.right_edge[k] = fmax(fmax(tri.p0[k], tri.p1[k]), tri.p2[k])
+
+ self.root = self._recursive_build(0, num_tri)
+
+ cdef void _recursive_free(self, BVHNode* node) nogil:
+ if node.end - node.begin > self.leaf_size:
+ self._recursive_free(node.left)
+ self._recursive_free(node.right)
+ free(node)
+
+ def __dealloc__(self):
+ self._recursive_free(self.root)
+ free(self.triangles)
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef np.int64_t _partition(self, np.int64_t begin, np.int64_t end,
+ np.int64_t ax, np.float64_t split) nogil:
+ # this re-orders the triangle array so that all of the triangles
+ # to the left of mid have centroids less than or equal to "split"
+ # along the direction "ax". All the triangles to the right of mid
+ # will have centroids *greater* than "split" along "ax".
+ cdef np.int64_t mid = begin
+ while (begin != end):
+ if self.triangles[mid].centroid[ax] > split:
+ mid += 1
+ elif self.triangles[begin].centroid[ax] > split:
+ self.triangles[mid], self.triangles[begin] = \
+ self.triangles[begin], self.triangles[mid]
+ mid += 1
+ begin += 1
+ return mid
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void _get_node_bbox(self, BVHNode* node,
+ np.int64_t begin, np.int64_t end) nogil:
+ cdef np.int64_t i, j
+ cdef BBox box = self.triangles[begin].bbox
+ for i in range(begin+1, end):
+ for j in range(3):
+ box.left_edge[j] = fmin(box.left_edge[j],
+ self.triangles[i].bbox.left_edge[j])
+ box.right_edge[j] = fmax(box.right_edge[j],
+ self.triangles[i].bbox.right_edge[j])
+ node.bbox = box
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void intersect(self, Ray* ray) nogil:
+ self._recursive_intersect(ray, self.root)
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef void _recursive_intersect(self, Ray* ray, BVHNode* node) nogil:
+
+ # check for bbox intersection:
+ if not ray_bbox_intersect(ray, node.bbox):
+ return
+
+ # check for leaf
+ cdef np.int64_t i, hit
+ cdef Triangle* tri
+ if (node.end - node.begin) <= self.leaf_size:
+ for i in range(node.begin, node.end):
+ tri = &(self.triangles[i])
+ hit = ray_triangle_intersect(ray, tri)
+ return
+
+ # if not leaf, intersect with left and right children
+ self._recursive_intersect(ray, node.left)
+ self._recursive_intersect(ray, node.right)
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef BVHNode* _recursive_build(self, np.int64_t begin, np.int64_t end) nogil:
+ cdef BVHNode *node = <BVHNode* > malloc(sizeof(BVHNode))
+ node.begin = begin
+ node.end = end
+
+ self._get_node_bbox(node, begin, end)
+
+ # check for leaf
+ if (end - begin) <= self.leaf_size:
+ return node
+
+ # we use the "split in the middle of the longest axis approach"
+ # see: http://www.vadimkravcenko.com/bvh-tree-building/
+
+ # compute longest dimension
+ cdef np.int64_t ax = 0
+ cdef np.float64_t d = fabs(node.bbox.right_edge[0] -
+ node.bbox.left_edge[0])
+ if fabs(node.bbox.right_edge[1] - node.bbox.left_edge[1]) > d:
+ ax = 1
+ if fabs(node.bbox.right_edge[2] - node.bbox.left_edge[2]) > d:
+ ax = 2
+
+ # split in half along that dimension
+ cdef np.float64_t split = 0.5*(node.bbox.right_edge[ax] +
+ node.bbox.left_edge[ax])
+
+ # sort triangle list
+ cdef np.int64_t mid = self._partition(begin, end, ax, split)
+
+ if(mid == begin or mid == end):
+ mid = begin + (end-begin)/2
+
+ # recursively build sub-trees
+ node.left = self._recursive_build(begin, mid)
+ node.right = self._recursive_build(mid, end)
+
+ return node
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void cast_rays(np.float64_t* image,
+ const np.float64_t* origins,
+ const np.float64_t* direction,
+ const int N,
+ BVH bvh) nogil:
+
+ cdef Ray* ray
+ cdef int i, j, k
+
+ with nogil, parallel():
+
+ ray = <Ray *> malloc(sizeof(Ray))
+
+ for k in range(3):
+ ray.direction[k] = direction[k]
+ ray.inv_dir[k] = 1.0 / direction[k]
+
+ for i in prange(N):
+ for j in range(3):
+ ray.origin[j] = origins[N*j + i]
+ ray.t_far = 1e30
+ ray.t_near = 0.0
+ ray.data_val = 0
+ bvh.intersect(ray)
+ image[i] = ray.data_val
+
+ free(ray)
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def test_ray_trace(np.ndarray[np.float64_t, ndim=1] image,
+ np.ndarray[np.float64_t, ndim=2] origins,
+ np.ndarray[np.float64_t, ndim=1] direction,
+ BVH bvh):
+
+ cdef int N = origins.shape[0]
+ cast_rays(&image[0], &origins[0, 0], &direction[0], N, bvh)
diff -r aaa54f6b1a1b1f98ea52b203fccce7c51a1ccc91 -r d18dc7748e19ec77816ee9b44147e5c0597c3185 yt/utilities/lib/tests/test_bounding_volume_hierarchy.py
--- /dev/null
+++ b/yt/utilities/lib/tests/test_bounding_volume_hierarchy.py
@@ -0,0 +1,47 @@
+import yt
+import numpy as np
+from yt.utilities.lib.bounding_volume_hierarchy import BVH, \
+ test_ray_trace
+from yt.visualization.volume_rendering.api import Camera
+from yt.testing import requires_file
+
+
+def get_rays(camera):
+
+ normal_vector = camera.unit_vectors[2].d
+ W = np.array([8.0, 8.0])
+ N = np.array([800, 800])
+ dx = W / N
+
+ x_points = np.linspace((-N[0]/2 + 0.5)*dx[0], (N[0]/2 - 0.5)*dx[0], N[0])
+ y_points = np.linspace((-N[1]/2 + 0.5)*dx[1], (N[1]/2 - 0.5)*dx[0], N[1])
+
+ X, Y = np.meshgrid(x_points, y_points)
+ result = np.dot(camera.unit_vectors[0:2].T, [X.ravel(), Y.ravel()])
+ vec_origins = result.T + camera.position
+ return np.array(vec_origins), np.array(normal_vector)
+
+
+fn = "MOOSE_sample_data/out.e-s010"
+
+
+ at requires_file(fn)
+def test_bounding_volume_hierarchy():
+ ds = yt.load(fn)
+ vertices = ds.index.meshes[0].connectivity_coords
+ indices = ds.index.meshes[0].connectivity_indices - 1
+
+ ad = ds.all_data()
+ field_data = ad['connect1', 'diffused']
+
+ bvh = BVH(vertices, indices, field_data)
+
+ cam = Camera()
+ cam.set_position(np.array([8.0, 8.0, 8.0]))
+ cam.focus = np.array([0.0, 0.0, 0.0])
+ origins, direction = get_rays(cam)
+
+ image = np.empty(800*800, np.float64)
+ test_ray_trace(image, origins, direction, bvh)
+ image = image.reshape((800, 800))
+ return image
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list