[yt-svn] commit/yt: xarthisius: Merged in xarthisius/yt (pull request #2028)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Mar 16 20:57:36 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/9d3023efc9e3/
Changeset: 9d3023efc9e3
Branch: yt
User: xarthisius
Date: 2016-03-17 03:57:30+00:00
Summary: Merged in xarthisius/yt (pull request #2028)
Pre-calculate the size of projection of rays onto image plane for perspective camera
Affected #: 7 files
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba setup.py
--- a/setup.py
+++ b/setup.py
@@ -127,7 +127,9 @@
["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"]),
+ libraries=["m"],
+ depends=["yt/utilities/lib/bounding_volume_hierarchy.pxd",
+ "yt/utilities/lib/vec3_ops.pxd"]),
Extension("yt.utilities.lib.contour_finding",
["yt/utilities/lib/contour_finding.pyx"],
include_dirs=["yt/utilities/lib/",
@@ -177,7 +179,8 @@
"yt/utilities/lib/kdtree.h",
"yt/utilities/lib/fixed_interpolator.h",
"yt/utilities/lib/fixed_interpolator.pxd",
- "yt/utilities/lib/field_interpolation_tables.pxd"]),
+ "yt/utilities/lib/field_interpolation_tables.pxd",
+ "yt/utilities/lib/vec3_ops.pxd"]),
Extension("yt.utilities.lib.element_mappings",
["yt/utilities/lib/element_mappings.pyx"],
libraries=["m"], depends=["yt/utilities/lib/element_mappings.pxd"]),
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/utilities/lib/bounding_volume_hierarchy.pyx
--- a/yt/utilities/lib/bounding_volume_hierarchy.pyx
+++ b/yt/utilities/lib/bounding_volume_hierarchy.pyx
@@ -4,6 +4,7 @@
from libc.math cimport fabs, fmax, fmin
from libc.stdlib cimport malloc, free
from cython.parallel import parallel, prange
+from vec3_ops cimport dot, subtract, cross
cdef extern from "mesh_construction.h":
enum:
@@ -15,35 +16,6 @@
cdef np.float64_t INF = np.inf
cdef np.int64_t LEAF_SIZE = 16
- 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]
-
@cython.boundscheck(False)
@cython.wraparound(False)
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/utilities/lib/grid_traversal.pxd
--- a/yt/utilities/lib/grid_traversal.pxd
+++ b/yt/utilities/lib/grid_traversal.pxd
@@ -27,6 +27,7 @@
np.float64_t[:,:] zbuffer
np.float64_t pdx, pdy
np.float64_t bounds[4]
+ np.float64_t[:,:] camera_data # position, width, unit_vec[0,2]
int nv[2]
np.float64_t *x_vec
np.float64_t *y_vec
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -17,9 +17,11 @@
cimport numpy as np
cimport cython
#cimport healpix_interface
+cdef extern from "limits.h":
+ cdef int SHRT_MAX
from libc.stdlib cimport malloc, calloc, free, abs
from libc.math cimport exp, floor, log2, \
- lrint, fabs, atan, atan2, asin, cos, sin, sqrt
+ lrint, fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, i64clip
from field_interpolation_tables cimport \
FieldInterpolationTable, FIT_initialize_table, FIT_eval_transfer,\
@@ -27,6 +29,7 @@
from fixed_interpolator cimport *
from cython.parallel import prange, parallel, threadid
+from vec3_ops cimport dot, subtract, L2_norm, fma
DEF Nch = 4
@@ -210,6 +213,108 @@
rv[2] = lrint((extrema[2] - cy - image.bounds[2])/image.pdy)
rv[3] = rv[2] + lrint((extrema[3] - extrema[2])/image.pdy)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef void calculate_extent_perspective(ImageContainer *image,
+ VolumeContainer *vc, np.int64_t rv[4]) nogil:
+
+ cdef np.float64_t cam_pos[3]
+ cdef np.float64_t cam_width[3]
+ cdef np.float64_t north_vector[3]
+ cdef np.float64_t east_vector[3]
+ cdef np.float64_t normal_vector[3]
+ cdef np.float64_t vertex[3]
+ cdef np.float64_t pos1[3]
+ cdef np.float64_t sight_vector[3]
+ cdef np.float64_t sight_center[3]
+ cdef np.float64_t corners[3][8]
+ cdef float sight_vector_norm, sight_angle_cos, sight_length, dx, dy
+ cdef int i, iv, px, py
+ cdef int min_px, min_py, max_px, max_py
+
+ min_px = SHRT_MAX
+ min_py = SHRT_MAX
+ max_px = -SHRT_MAX
+ max_py = -SHRT_MAX
+
+ # calculate vertices for 8 corners of vc
+ corners[0][0] = vc.left_edge[0]
+ corners[0][1] = vc.right_edge[0]
+ corners[0][2] = vc.right_edge[0]
+ corners[0][3] = vc.left_edge[0]
+ corners[0][4] = vc.left_edge[0]
+ corners[0][5] = vc.right_edge[0]
+ corners[0][6] = vc.right_edge[0]
+ corners[0][7] = vc.left_edge[0]
+
+ corners[1][0] = vc.left_edge[1]
+ corners[1][1] = vc.left_edge[1]
+ corners[1][2] = vc.right_edge[1]
+ corners[1][3] = vc.right_edge[1]
+ corners[1][4] = vc.left_edge[1]
+ corners[1][5] = vc.left_edge[1]
+ corners[1][6] = vc.right_edge[1]
+ corners[1][7] = vc.right_edge[1]
+
+ corners[2][0] = vc.left_edge[2]
+ corners[2][1] = vc.left_edge[2]
+ corners[2][2] = vc.left_edge[2]
+ corners[2][3] = vc.left_edge[2]
+ corners[2][4] = vc.right_edge[2]
+ corners[2][5] = vc.right_edge[2]
+ corners[2][6] = vc.right_edge[2]
+ corners[2][7] = vc.right_edge[2]
+
+ # This code was ported from
+ # yt.visualization.volume_rendering.lens.PerspectiveLens.project_to_plane()
+ for i in range(3):
+ cam_pos[i] = image.camera_data[0, i]
+ cam_width[i] = image.camera_data[1, i]
+ east_vector[i] = image.camera_data[2, i]
+ north_vector[i] = image.camera_data[3, i]
+ normal_vector[i] = image.camera_data[4, i]
+
+ for iv in range(8):
+ vertex[0] = corners[0][iv]
+ vertex[1] = corners[1][iv]
+ vertex[2] = corners[2][iv]
+
+ subtract(vertex, cam_pos, sight_vector)
+ fma(cam_width[2], normal_vector, cam_pos, sight_center)
+
+ sight_vector_norm = L2_norm(sight_vector)
+
+ if sight_vector_norm != 0:
+ for i in range(3):
+ sight_vector[i] /= sight_vector_norm
+
+ sight_angle_cos = dot(sight_vector, normal_vector)
+ sight_angle_cos = fclip(sight_angle_cos, -1.0, 1.0)
+
+ if acos(sight_angle_cos) < 0.5 * M_PI and sight_angle_cos != 0.0:
+ sight_length = cam_width[2] / sight_angle_cos
+ else:
+ sight_length = sqrt(cam_width[0]**2 + cam_width[1]**2)
+ sight_length = sight_length / sqrt(1.0 - sight_angle_cos**2)
+
+ fma(sight_length, sight_vector, cam_pos, pos1)
+ subtract(pos1, sight_center, pos1)
+ dx = dot(pos1, east_vector)
+ dy = dot(pos1, north_vector)
+
+ px = int(image.nv[0] * 0.5 + image.nv[0] / cam_width[0] * dx)
+ py = int(image.nv[1] * 0.5 + image.nv[1] / cam_width[1] * dy)
+ min_px = min(min_px, px)
+ max_px = max(max_px, px)
+ min_py = min(min_py, py)
+ max_py = max(max_py, py)
+
+ rv[0] = max(min_px, 0)
+ rv[1] = min(max_px, image.nv[0])
+ rv[2] = max(min_py, 0)
+ rv[3] = min(max_py, image.nv[1])
+
+
# We do this for a bunch of lenses. Fallback is to grab them from the vector
# info supplied.
@@ -271,6 +376,13 @@
*args, **kwargs):
self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
cdef np.float64_t[:,:] zbuffer
+ cdef np.float64_t[:,:] camera_data
+ cdef int i
+
+ camera_data = kwargs.pop("camera_data", None)
+ if camera_data is not None:
+ self.image.camera_data = camera_data
+
zbuffer = kwargs.pop("zbuffer", None)
if zbuffer is None:
zbuffer = np.ones((image.shape[0], image.shape[1]), "float64")
@@ -281,15 +393,19 @@
else:
if not (vp_pos.shape[0] == vp_dir.shape[0] == image.shape[0]) or \
not (vp_pos.shape[1] == vp_dir.shape[1] == image.shape[1]):
- msg = "Bad lens shape / direction for %s\n" % (self.lens_type)
+ msg = "Bad lense shape / direction for %s\n" % (self.lens_type)
msg += "Shapes: (%s - %s - %s) and (%s - %s - %s)" % (
vp_pos.shape[0], vp_dir.shape[0], image.shape[0],
vp_pos.shape[1], vp_dir.shape[1], image.shape[1])
raise RuntimeError(msg)
- self.extent_function = calculate_extent_null
+
+ if camera_data is not None and "perspective" in self.lens_type:
+ self.extent_function = calculate_extent_perspective
+ else:
+ self.extent_function = calculate_extent_null
self.vector_function = generate_vector_info_null
+
self.sampler = NULL
- cdef int i
# These assignments are so we can track the objects and prevent their
# de-allocation from reference counts. Note that we do this to the
# "atleast_3d" versions. Also, note that we re-assign the input
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/utilities/lib/vec3_ops.pxd
--- /dev/null
+++ b/yt/utilities/lib/vec3_ops.pxd
@@ -0,0 +1,50 @@
+cimport cython
+import numpy as np
+cimport numpy as np
+from libc.math cimport sqrt
+
+ 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 inline void fma(const np.float64_t f,
+ const np.float64_t a[3],
+ const np.float64_t b[3],
+ np.float64_t c[3]) nogil:
+ c[0] = f * a[0] + b[0]
+ c[1] = f * a[1] + b[1]
+ c[2] = f * a[2] + b[2]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline np.float64_t L2_norm(const np.float64_t a[3]) nogil:
+ return sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2])
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/visualization/volume_rendering/camera.py
--- a/yt/visualization/volume_rendering/camera.py
+++ b/yt/visualization/volume_rendering/camera.py
@@ -280,6 +280,8 @@
def _get_sampler_params(self, render_source):
lens_params = self.lens._get_sampler_params(self, render_source)
lens_params.update(width=self.width)
+ lens_params.update(camera_data=np.vstack(
+ (self.position.d, self.width.d, self.unit_vectors.d)))
return lens_params
def set_lens(self, lens_type):
diff -r 1cc2beb515d168c04a39df000d33bd4f06f987d4 -r 9d3023efc9e320ecc1cb048819162f701c9487ba yt/visualization/volume_rendering/utils.py
--- a/yt/visualization/volume_rendering/utils.py
+++ b/yt/visualization/volume_rendering/utils.py
@@ -49,6 +49,8 @@
params['num_samples'],
)
kwargs = {'lens_type': params['lens_type']}
+ if "camera_data" in params:
+ kwargs['camera_data'] = params['camera_data']
if render_source.zbuffer is not None:
kwargs['zbuffer'] = render_source.zbuffer.z
args[4][:] = np.reshape(render_source.zbuffer.rgba[:], \
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