[yt-svn] commit/yt: 6 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sat Jul 16 10:42:02 PDT 2016
6 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/caf2aa0abab3/
Changeset: caf2aa0abab3
Branch: yt
User: MatthewTurk
Date: 2016-07-07 19:57:29+00:00
Summary: Beginning to refactor volume rendering Cython code.
Affected #: 15 files
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -62,6 +62,11 @@
yt/utilities/lib/marching_cubes.c
yt/utilities/lib/png_writer.h
yt/utilities/lib/write_array.c
+yt/utilities/lib/perftools_wrap.c
+yt/utilities/lib/partitioned_grid.c
+yt/utilities/lib/volume_container.c
+yt/utilities/lib/lenses.c
+yt/utilities/lib/image_samplers.c
syntax: glob
*.pyc
*.pyd
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b setup.py
--- a/setup.py
+++ b/setup.py
@@ -186,7 +186,7 @@
"particle_mesh_operations", "depth_first_octree", "fortran_reader",
"interpolators", "misc_utilities", "basic_octree", "image_utilities",
"points_in_volume", "quad_tree", "ray_integrators", "mesh_utilities",
- "amr_kdtools"
+ "amr_kdtools", "image_samplers", "lenses"
]
for ext_name in lib_exts:
cython_extensions.append(
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -22,8 +22,10 @@
from .oct_container cimport OctreeContainer, OctAllocationContainer, Oct
cimport oct_visitors
from .oct_visitors cimport cind
+from yt.utilities.lib.volume_container cimport \
+ VolumeContainer
from yt.utilities.lib.grid_traversal cimport \
- VolumeContainer, sample_function, walk_volume
+ sampler_function, walk_volume
from yt.utilities.lib.bitarray cimport ba_get_value, ba_set_value
cdef extern from "math.h":
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/contour_finding.pyx
--- a/yt/utilities/lib/contour_finding.pyx
+++ b/yt/utilities/lib/contour_finding.pyx
@@ -26,8 +26,10 @@
from yt.geometry.oct_visitors cimport \
Oct
from .amr_kdtools cimport Node
-from .grid_traversal cimport VolumeContainer, PartitionedGrid, \
- vc_index, vc_pos_index
+from .partitioned_grid cimport \
+ PartitionedGrid
+from .volume_container cimport \
+ VolumeContainer, vc_index, vc_pos_index
import sys
cdef inline ContourID *contour_create(np.int64_t contour_id,
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/field_interpolation_tables.pxd
--- a/yt/utilities/lib/field_interpolation_tables.pxd
+++ b/yt/utilities/lib/field_interpolation_tables.pxd
@@ -18,8 +18,6 @@
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, fabs
from libc.stdlib cimport malloc
-DEF Nch = 4
-
cdef struct FieldInterpolationTable:
# Note that we make an assumption about retaining a reference to values
# externally.
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/grid_traversal.pxd
--- a/yt/utilities/lib/grid_traversal.pxd
+++ b/yt/utilities/lib/grid_traversal.pxd
@@ -17,22 +17,8 @@
import numpy as np
cimport numpy as np
cimport cython
-cimport kdtree_utils
-
-cdef struct ImageContainer:
- np.float64_t[:,:,:] vp_pos
- np.float64_t[:,:,:] vp_dir
- np.float64_t *center
- np.float64_t[:,:,:] image
- np.float64_t[:,:] zbuffer
- np.int64_t[:,:] image_used
- np.int64_t[:,:] mesh_lines
- 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
+from .image_samplers cimport ImageContainer, ImageSampler
+from .volume_container cimport VolumeContainer, vc_index, vc_pos_index
ctypedef void sampler_function(
VolumeContainer *vc,
@@ -43,82 +29,11 @@
int index[3],
void *data) nogil
-ctypedef void calculate_extent_function(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil
-
-cdef calculate_extent_function calculate_extent_plane_parallel
-
-ctypedef void generate_vector_info_function(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil
-
-cdef generate_vector_info_function generate_vector_info_plane_parallel
-cdef generate_vector_info_function generate_vector_info_null
-
-cdef class ImageSampler:
- cdef ImageContainer *image
- cdef sampler_function *sampler
- cdef public object acenter, aimage, ax_vec, ay_vec
- cdef public object azbuffer
- cdef public object aimage_used
- cdef public object amesh_lines
- cdef void *supp_data
- cdef np.float64_t width[3]
- cdef public object lens_type
- cdef calculate_extent_function *extent_function
- cdef generate_vector_info_function *vector_function
-
- cdef void setup(self, PartitionedGrid pg)
-
-cdef struct VolumeContainer:
- int n_fields
- np.float64_t **data
- # The mask has dimensions one fewer in each direction than data
- np.uint8_t *mask
- np.float64_t left_edge[3]
- np.float64_t right_edge[3]
- np.float64_t dds[3]
- np.float64_t idds[3]
- int dims[3]
-
-cdef class PartitionedGrid:
- cdef public object my_data
- cdef public object source_mask
- cdef public object LeftEdge
- cdef public object RightEdge
- cdef public int parent_grid_id
- cdef VolumeContainer *container
- cdef kdtree_utils.kdtree *star_list
- cdef np.float64_t star_er
- cdef np.float64_t star_sigma_num
- cdef np.float64_t star_coeff
- cdef void get_vector_field(self, np.float64_t pos[3],
- np.float64_t *vel, np.float64_t *vel_mag)
-
-ctypedef void sample_function(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil
-
cdef int walk_volume(VolumeContainer *vc,
np.float64_t v_pos[3],
np.float64_t v_dir[3],
- sample_function *sampler,
+ sampler_function *sampler,
void *data,
np.float64_t *return_t = *,
np.float64_t max_t = *) nogil
-cdef inline int vc_index(VolumeContainer *vc, int i, int j, int k):
- return (i*vc.dims[1]+j)*vc.dims[2]+k
-
-cdef inline int vc_pos_index(VolumeContainer *vc, np.float64_t *spos):
- cdef int index[3]
- cdef int i
- for i in range(3):
- index[i] = <int> ((spos[i] - vc.left_edge[i]) * vc.idds[i])
- return vc_index(vc, index[0], index[1], index[2])
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -17,8 +17,6 @@
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, \
fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
@@ -28,913 +26,23 @@
FIT_eval_transfer_with_light
from fixed_interpolator cimport *
-cdef extern from "platform_dep.h":
- long int lrint(double x) nogil
-
from cython.parallel import prange, parallel, threadid
-from vec3_ops cimport dot, subtract, L2_norm, fma
-
from cpython.exc cimport PyErr_CheckSignals
+from .image_samplers cimport \
+ ImageSampler, \
+ ImageContainer, \
+ VolumeRenderAccumulator
+
DEF Nch = 4
-cdef class PartitionedGrid:
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def __cinit__(self,
- int parent_grid_id, data,
- mask,
- np.ndarray[np.float64_t, ndim=1] left_edge,
- np.ndarray[np.float64_t, ndim=1] right_edge,
- np.ndarray[np.int64_t, ndim=1] dims,
- star_kdtree_container star_tree = None):
- # The data is likely brought in via a slice, so we copy it
- cdef np.ndarray[np.float64_t, ndim=3] tdata
- cdef np.ndarray[np.uint8_t, ndim=3] mask_data
- self.container = NULL
- self.parent_grid_id = parent_grid_id
- self.LeftEdge = left_edge
- self.RightEdge = right_edge
- self.container = <VolumeContainer *> \
- malloc(sizeof(VolumeContainer))
- cdef VolumeContainer *c = self.container # convenience
- cdef int n_fields = len(data)
- c.n_fields = n_fields
- for i in range(3):
- c.left_edge[i] = left_edge[i]
- c.right_edge[i] = right_edge[i]
- c.dims[i] = dims[i]
- c.dds[i] = (c.right_edge[i] - c.left_edge[i])/dims[i]
- c.idds[i] = 1.0/c.dds[i]
- self.my_data = data
- self.source_mask = mask
- mask_data = mask
- c.data = <np.float64_t **> malloc(sizeof(np.float64_t*) * n_fields)
- for i in range(n_fields):
- tdata = data[i]
- c.data[i] = <np.float64_t *> tdata.data
- c.mask = <np.uint8_t *> mask_data.data
- if star_tree is None:
- self.star_list = NULL
- else:
- self.set_star_tree(star_tree)
-
- def set_star_tree(self, star_kdtree_container star_tree):
- self.star_list = star_tree.tree
- self.star_sigma_num = 2.0*star_tree.sigma**2.0
- self.star_er = 2.326 * star_tree.sigma
- self.star_coeff = star_tree.coeff
-
- def __dealloc__(self):
- # The data fields are not owned by the container, they are owned by us!
- # So we don't need to deallocate them.
- if self.container == NULL: return
- if self.container.data != NULL: free(self.container.data)
- free(self.container)
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def integrate_streamline(self, pos, np.float64_t h, mag):
- cdef np.float64_t cmag[1]
- cdef np.float64_t k1[3]
- cdef np.float64_t k2[3]
- cdef np.float64_t k3[3]
- cdef np.float64_t k4[3]
- cdef np.float64_t newpos[3]
- cdef np.float64_t oldpos[3]
- for i in range(3):
- newpos[i] = oldpos[i] = pos[i]
- self.get_vector_field(newpos, k1, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + 0.5*k1[i]*h
-
- if not (self.LeftEdge[0] < newpos[0] and newpos[0] < self.RightEdge[0] and \
- self.LeftEdge[1] < newpos[1] and newpos[1] < self.RightEdge[1] and \
- self.LeftEdge[2] < newpos[2] and newpos[2] < self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k2, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + 0.5*k2[i]*h
-
- if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
- self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
- self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k3, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + k3[i]*h
-
- if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
- self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
- self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k4, cmag)
-
- for i in range(3):
- pos[i] = oldpos[i] + h*(k1[i]/6.0 + k2[i]/3.0 + k3[i]/3.0 + k4[i]/6.0)
-
- if mag is not None:
- for i in range(3):
- newpos[i] = pos[i]
- self.get_vector_field(newpos, k4, cmag)
- mag[0] = cmag[0]
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef void get_vector_field(self, np.float64_t pos[3],
- np.float64_t *vel, np.float64_t *vel_mag):
- cdef np.float64_t dp[3]
- cdef int ci[3]
- cdef VolumeContainer *c = self.container # convenience
-
- for i in range(3):
- ci[i] = (int)((pos[i]-self.LeftEdge[i])/c.dds[i])
- dp[i] = (pos[i] - ci[i]*c.dds[i] - self.LeftEdge[i])/c.dds[i]
-
- cdef int offset = ci[0] * (c.dims[1] + 1) * (c.dims[2] + 1) \
- + ci[1] * (c.dims[2] + 1) + ci[2]
-
- vel_mag[0] = 0.0
- for i in range(3):
- vel[i] = offset_interpolate(c.dims, dp, c.data[i] + offset)
- vel_mag[0] += vel[i]*vel[i]
- vel_mag[0] = np.sqrt(vel_mag[0])
- if vel_mag[0] != 0.0:
- for i in range(3):
- vel[i] /= vel_mag[0]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void calculate_extent_plane_parallel(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil:
- # We do this for all eight corners
- cdef np.float64_t temp
- cdef np.float64_t *edges[2]
- cdef np.float64_t cx, cy
- cdef np.float64_t extrema[4]
- cdef int i, j, k
- edges[0] = vc.left_edge
- edges[1] = vc.right_edge
- extrema[0] = extrema[2] = 1e300; extrema[1] = extrema[3] = -1e300
- for i in range(2):
- for j in range(2):
- for k in range(2):
- # This should rotate it into the vector plane
- temp = edges[i][0] * image.x_vec[0]
- temp += edges[j][1] * image.x_vec[1]
- temp += edges[k][2] * image.x_vec[2]
- if temp < extrema[0]: extrema[0] = temp
- if temp > extrema[1]: extrema[1] = temp
- temp = edges[i][0] * image.y_vec[0]
- temp += edges[j][1] * image.y_vec[1]
- temp += edges[k][2] * image.y_vec[2]
- if temp < extrema[2]: extrema[2] = temp
- if temp > extrema[3]: extrema[3] = temp
- cx = cy = 0.0
- for i in range(3):
- cx += image.center[i] * image.x_vec[i]
- cy += image.center[i] * image.y_vec[i]
- rv[0] = lrint((extrema[0] - cx - image.bounds[0])/image.pdx)
- rv[1] = rv[0] + lrint((extrema[1] - extrema[0])/image.pdx)
- 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]
-
- cam_width[1] = cam_width[0] * image.nv[1] / image.nv[0]
-
- 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.
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void calculate_extent_null(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil:
- rv[0] = 0
- rv[1] = image.nv[0]
- rv[2] = 0
- rv[3] = image.nv[1]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void generate_vector_info_plane_parallel(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- # Now outbound
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
- cdef int i
- cdef np.float64_t px, py
- px = width[0] * (<np.float64_t>vi)/(<np.float64_t>im.nv[0]-1) - width[0]/2.0
- py = width[1] * (<np.float64_t>vj)/(<np.float64_t>im.nv[1]-1) - width[1]/2.0
- # atleast_3d will add to beginning and end
- v_pos[0] = im.vp_pos[0,0,0]*px + im.vp_pos[0,3,0]*py + im.vp_pos[0,9,0]
- v_pos[1] = im.vp_pos[0,1,0]*px + im.vp_pos[0,4,0]*py + im.vp_pos[0,10,0]
- v_pos[2] = im.vp_pos[0,2,0]*px + im.vp_pos[0,5,0]*py + im.vp_pos[0,11,0]
- for i in range(3): v_dir[i] = im.vp_dir[0,i,0]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void generate_vector_info_null(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- # Now outbound
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
- cdef int i
- for i in range(3):
- # Here's a funny thing: we use vi here because our *image* will be
- # flattened. That means that im.nv will be a better one-d offset,
- # since vp_pos has funny strides.
- v_pos[i] = im.vp_pos[vi, vj, i]
- v_dir[i] = im.vp_dir[vi, vj, i]
-
-cdef struct ImageAccumulator:
- np.float64_t rgba[Nch]
- void *supp_data
-
-cdef class ImageSampler:
- def __init__(self,
- np.float64_t[:,:,:] vp_pos,
- np.float64_t[:,:,:] vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- *args, **kwargs):
- self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
- cdef np.float64_t[:,:] zbuffer
- cdef np.int64_t[:,:] image_used
- cdef np.int64_t[:,:] mesh_lines
- 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")
-
- image_used = np.zeros((image.shape[0], image.shape[1]), "int64")
- mesh_lines = np.zeros((image.shape[0], image.shape[1]), "int64")
-
- self.lens_type = kwargs.pop("lens_type", None)
- if self.lens_type == "plane-parallel":
- self.extent_function = calculate_extent_plane_parallel
- self.vector_function = generate_vector_info_plane_parallel
- 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 += "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)
-
- if camera_data is not None and self.lens_type == 'perspective':
- self.extent_function = calculate_extent_perspective
- else:
- self.extent_function = calculate_extent_null
- self.vector_function = generate_vector_info_null
-
- self.sampler = NULL
- # 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
- # arguments.
- self.image.vp_pos = vp_pos
- self.image.vp_dir = vp_dir
- self.image.image = self.aimage = image
- self.acenter = center
- self.image.center = <np.float64_t *> center.data
- self.ax_vec = x_vec
- self.image.x_vec = <np.float64_t *> x_vec.data
- self.ay_vec = y_vec
- self.image.y_vec = <np.float64_t *> y_vec.data
- self.image.zbuffer = self.azbuffer = zbuffer
- self.image.image_used = self.aimage_used = image_used
- self.image.mesh_lines = self.amesh_lines = mesh_lines
- self.image.nv[0] = image.shape[0]
- self.image.nv[1] = image.shape[1]
- for i in range(4): self.image.bounds[i] = bounds[i]
- self.image.pdx = (bounds[1] - bounds[0])/self.image.nv[0]
- self.image.pdy = (bounds[3] - bounds[2])/self.image.nv[1]
- for i in range(3):
- self.width[i] = width[i]
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def __call__(self, PartitionedGrid pg, int num_threads = 0):
- # This routine will iterate over all of the vectors and cast each in
- # turn. Might benefit from a more sophisticated intersection check,
- # like http://courses.csusm.edu/cs697exz/ray_box.htm
- cdef int vi, vj, hit, i, j
- cdef np.int64_t iter[4]
- cdef VolumeContainer *vc = pg.container
- cdef ImageContainer *im = self.image
- self.setup(pg)
- if self.sampler == NULL: raise RuntimeError
- cdef np.float64_t *v_pos
- cdef np.float64_t *v_dir
- cdef np.float64_t max_t
- hit = 0
- cdef np.int64_t nx, ny, size
- self.extent_function(self.image, vc, iter)
- iter[0] = i64clip(iter[0]-1, 0, im.nv[0])
- iter[1] = i64clip(iter[1]+1, 0, im.nv[0])
- iter[2] = i64clip(iter[2]-1, 0, im.nv[1])
- iter[3] = i64clip(iter[3]+1, 0, im.nv[1])
- nx = (iter[1] - iter[0])
- ny = (iter[3] - iter[2])
- size = nx * ny
- cdef ImageAccumulator *idata
- cdef np.float64_t width[3]
- cdef int chunksize = 100
- for i in range(3):
- width[i] = self.width[i]
- with nogil, parallel(num_threads = num_threads):
- idata = <ImageAccumulator *> malloc(sizeof(ImageAccumulator))
- idata.supp_data = self.supp_data
- v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- for j in prange(size, schedule="static", chunksize=chunksize):
- vj = j % ny
- vi = (j - vj) / ny + iter[0]
- vj = vj + iter[2]
- # Dynamically calculate the position
- self.vector_function(im, vi, vj, width, v_dir, v_pos)
- for i in range(Nch):
- idata.rgba[i] = im.image[vi, vj, i]
- max_t = fclip(im.zbuffer[vi, vj], 0.0, 1.0)
- walk_volume(vc, v_pos, v_dir, self.sampler,
- (<void *> idata), NULL, max_t)
- if (j % (10*chunksize)) == 0:
- with gil:
- PyErr_CheckSignals()
- for i in range(Nch):
- im.image[vi, vj, i] = idata.rgba[i]
- free(idata)
- free(v_pos)
- free(v_dir)
- return hit
-
- cdef void setup(self, PartitionedGrid pg):
- return
-
- def __dealloc__(self):
- self.image.image = None
- self.image.vp_pos = None
- self.image.vp_dir = None
- self.image.zbuffer = None
- self.image.camera_data = None
- self.image.image_used = None
- free(self.image)
-
-
-cdef void projection_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef int i
- cdef np.float64_t dl = (exit_t - enter_t)
- cdef int di = (index[0]*vc.dims[1]+index[1])*vc.dims[2]+index[2]
- for i in range(imin(4, vc.n_fields)):
- im.rgba[i] += vc.data[i][di] * dl
-
-cdef struct VolumeRenderAccumulator:
- int n_fits
- int n_samples
- FieldInterpolationTable *fits
- int field_table_ids[6]
- np.float64_t star_coeff
- np.float64_t star_er
- np.float64_t star_sigma_num
- kdtree_utils.kdtree *star_list
- np.float64_t *light_dir
- np.float64_t *light_rgba
- int grey_opacity
-
-
-cdef class ProjectionSampler(ImageSampler):
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = projection_sampler
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void interpolated_projection_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- for j in range(imin(3, vc.n_fields)):
- im.rgba[j] += dvs[j] * dt
- for j in range(3):
- dp[j] += ds[j]
-
-cdef class InterpolatedProjectionSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- n_samples = 10, **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.n_samples = n_samples
- self.supp_data = <void *> self.vra
-
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = interpolated_projection_sampler
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef int cell_offset = index[0] * (vc.dims[1]) * (vc.dims[2]) \
- + index[1] * (vc.dims[2]) + index[2]
- if vc.mask[cell_offset] != 1:
- return
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits,
- vri.fits, vri.field_table_ids, vri.grey_opacity)
- for j in range(3):
- dp[j] += ds[j]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_gradient_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- cdef np.float64_t *grad
- grad = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- eval_gradient(vc.dims, dp, vc.data[0] + offset, grad)
- FIT_eval_transfer_with_light(dt, dvs, grad,
- vri.light_dir, vri.light_rgba,
- im.rgba, vri.n_fits,
- vri.fits, vri.field_table_ids, vri.grey_opacity)
- for j in range(3):
- dp[j] += ds[j]
- free(grad)
-
-cdef class star_kdtree_container:
- cdef kdtree_utils.kdtree *tree
- cdef public np.float64_t sigma
- cdef public np.float64_t coeff
-
- def __init__(self):
- self.tree = kdtree_utils.kd_create(3)
-
- def add_points(self,
- np.ndarray[np.float64_t, ndim=1] pos_x,
- np.ndarray[np.float64_t, ndim=1] pos_y,
- np.ndarray[np.float64_t, ndim=1] pos_z,
- np.ndarray[np.float64_t, ndim=2] star_colors):
- cdef int i
- cdef np.float64_t *pointer = <np.float64_t *> star_colors.data
- for i in range(pos_x.shape[0]):
- kdtree_utils.kd_insert3(self.tree,
- pos_x[i], pos_y[i], pos_z[i], <void *> (pointer + i*3))
-
- def __dealloc__(self):
- kdtree_utils.kd_free(self.tree)
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_stars_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- cdef kdtree_utils.kdres *ballq = NULL
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t slopes[6]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- cdef np.float64_t cell_left[3]
- cdef np.float64_t local_dds[3]
- cdef np.float64_t pos[3]
- cdef int nstars, i, j
- cdef np.float64_t *colors = NULL
- cdef np.float64_t gexp, gaussian, px, py, pz
- px = py = pz = -1
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vc.n_fields):
- slopes[i] = offset_interpolate(vc.dims, dp,
- vc.data[i] + offset)
- cdef np.float64_t temp
- # Now we get the ball-tree result for the stars near our cell center.
- for i in range(3):
- cell_left[i] = index[i] * vc.dds[i] + vc.left_edge[i]
- pos[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- local_dds[i] = v_dir[i] * dt
- ballq = kdtree_utils.kd_nearest_range3(
- vri.star_list, cell_left[0] + vc.dds[0]*0.5,
- cell_left[1] + vc.dds[1]*0.5,
- cell_left[2] + vc.dds[2]*0.5,
- vri.star_er + 0.9*vc.dds[0])
- # ~0.866 + a bit
-
- nstars = kdtree_utils.kd_res_size(ballq)
- for i in range(vc.n_fields):
- temp = slopes[i]
- slopes[i] -= offset_interpolate(vc.dims, dp,
- vc.data[i] + offset)
- slopes[i] *= -1.0/vri.n_samples
- dvs[i] = temp
- for _ in range(vri.n_samples):
- # Now we add the contribution from stars
- kdtree_utils.kd_res_rewind(ballq)
- for i in range(nstars):
- kdtree_utils.kd_res_item3(ballq, &px, &py, &pz)
- colors = <np.float64_t *> kdtree_utils.kd_res_item_data(ballq)
- kdtree_utils.kd_res_next(ballq)
- gexp = (px - pos[0])*(px - pos[0]) \
- + (py - pos[1])*(py - pos[1]) \
- + (pz - pos[2])*(pz - pos[2])
- gaussian = vri.star_coeff * exp(-gexp/vri.star_sigma_num)
- for j in range(3): im.rgba[j] += gaussian*dt*colors[j]
- for i in range(3):
- pos[i] += local_dds[i]
- FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits, vri.fits,
- vri.field_table_ids, vri.grey_opacity)
- for i in range(vc.n_fields):
- dvs[i] += slopes[i]
- kdtree_utils.kd_res_free(ballq)
-
-cdef class VolumeRenderSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- cdef kdtree_utils.kdtree **trees
- cdef object tree_containers
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- tf_obj, n_samples = 10,
- star_list = None, **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- cdef int i
- cdef np.ndarray[np.float64_t, ndim=1] temp
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.fits = <FieldInterpolationTable *> \
- malloc(sizeof(FieldInterpolationTable) * 6)
- self.vra.n_fits = tf_obj.n_field_tables
- assert(self.vra.n_fits <= 6)
- self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
- self.vra.n_samples = n_samples
- self.my_field_tables = []
- for i in range(self.vra.n_fits):
- temp = tf_obj.tables[i].y
- FIT_initialize_table(&self.vra.fits[i],
- temp.shape[0],
- <np.float64_t *> temp.data,
- tf_obj.tables[i].x_bounds[0],
- tf_obj.tables[i].x_bounds[1],
- tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
- tf_obj.weight_table_ids[i])
- self.my_field_tables.append((tf_obj.tables[i],
- tf_obj.tables[i].y))
- for i in range(6):
- self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
- self.supp_data = <void *> self.vra
- cdef star_kdtree_container skdc
- self.tree_containers = star_list
- if star_list is None:
- self.trees = NULL
- else:
- self.trees = <kdtree_utils.kdtree **> malloc(
- sizeof(kdtree_utils.kdtree*) * len(star_list))
- for i in range(len(star_list)):
- skdc = star_list[i]
- self.trees[i] = skdc.tree
-
- cdef void setup(self, PartitionedGrid pg):
- cdef star_kdtree_container star_tree
- if self.trees == NULL:
- self.sampler = volume_render_sampler
- else:
- star_tree = self.tree_containers[pg.parent_grid_id]
- self.vra.star_list = self.trees[pg.parent_grid_id]
- self.vra.star_sigma_num = 2.0*star_tree.sigma**2.0
- self.vra.star_er = 2.326 * star_tree.sigma
- self.vra.star_coeff = star_tree.coeff
- self.sampler = volume_render_stars_sampler
-
- def __dealloc__(self):
- for i in range(self.vra.n_fits):
- free(self.vra.fits[i].d0)
- free(self.vra.fits[i].dy)
- free(self.vra.fits)
- free(self.vra)
-
-cdef class LightSourceRenderSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- tf_obj, n_samples = 10,
- light_dir=[1.,1.,1.],
- light_rgba=[1.,1.,1.,1.],
- **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- cdef int i
- cdef np.ndarray[np.float64_t, ndim=1] temp
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.fits = <FieldInterpolationTable *> \
- malloc(sizeof(FieldInterpolationTable) * 6)
- self.vra.n_fits = tf_obj.n_field_tables
- assert(self.vra.n_fits <= 6)
- self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
- self.vra.n_samples = n_samples
- self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
- self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
- light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
- for i in range(3):
- self.vra.light_dir[i] = light_dir[i]
- for i in range(4):
- self.vra.light_rgba[i] = light_rgba[i]
- self.my_field_tables = []
- for i in range(self.vra.n_fits):
- temp = tf_obj.tables[i].y
- FIT_initialize_table(&self.vra.fits[i],
- temp.shape[0],
- <np.float64_t *> temp.data,
- tf_obj.tables[i].x_bounds[0],
- tf_obj.tables[i].x_bounds[1],
- tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
- tf_obj.weight_table_ids[i])
- self.my_field_tables.append((tf_obj.tables[i],
- tf_obj.tables[i].y))
- for i in range(6):
- self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
- self.supp_data = <void *> self.vra
-
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = volume_render_gradient_sampler
-
- def __dealloc__(self):
- for i in range(self.vra.n_fits):
- free(self.vra.fits[i].d0)
- free(self.vra.fits[i].dy)
- free(self.vra.light_dir)
- free(self.vra.light_rgba)
- free(self.vra.fits)
- free(self.vra)
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef int walk_volume(VolumeContainer *vc,
np.float64_t v_pos[3],
np.float64_t v_dir[3],
- sampler_function *sampler,
+ sampler_function *sample,
void *data,
np.float64_t *return_t = NULL,
np.float64_t max_t = 1.0) nogil:
@@ -1031,7 +139,7 @@
else:
i = 2
exit_t = fmin(tmax[i], max_t)
- sampler(vc, v_pos, v_dir, enter_t, exit_t, cur_ind, data)
+ sample(vc, v_pos, v_dir, enter_t, exit_t, cur_ind, data)
cur_ind[i] += step[i]
enter_t = tmax[i]
tmax[i] += tdelta[i]
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/image_samplers.pxd
--- /dev/null
+++ b/yt/utilities/lib/image_samplers.pxd
@@ -0,0 +1,103 @@
+"""
+Definitions for image samplers
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+cimport kdtree_utils
+from .volume_container cimport VolumeContainer
+from .lenses cimport \
+ calculate_extent_function, \
+ generate_vector_info_function
+from .partitioned_grid cimport PartitionedGrid
+
+from field_interpolation_tables cimport \
+ FieldInterpolationTable
+
+DEF Nch = 4
+
+cdef struct ImageContainer:
+ np.float64_t[:,:,:] vp_pos
+ np.float64_t[:,:,:] vp_dir
+ np.float64_t *center
+ np.float64_t[:,:,:] image
+ np.float64_t[:,:] zbuffer
+ np.int64_t[:,:] image_used
+ np.int64_t[:,:] mesh_lines
+ 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
+
+cdef struct VolumeRenderAccumulator:
+ int n_fits
+ int n_samples
+ FieldInterpolationTable *fits
+ int field_table_ids[6]
+ np.float64_t star_coeff
+ np.float64_t star_er
+ np.float64_t star_sigma_num
+ kdtree_utils.kdtree *star_list
+ np.float64_t *light_dir
+ np.float64_t *light_rgba
+ int grey_opacity
+
+cdef struct ImageAccumulator:
+ np.float64_t rgba[Nch]
+ void *supp_data
+
+cdef class ImageSampler:
+ cdef ImageContainer *image
+ cdef public object acenter, aimage, ax_vec, ay_vec
+ cdef public object azbuffer
+ cdef public object aimage_used
+ cdef public object amesh_lines
+ cdef void *supp_data
+ cdef np.float64_t width[3]
+ cdef public object lens_type
+ cdef calculate_extent_function *extent_function
+ cdef generate_vector_info_function *vector_function
+ cdef void setup(self, PartitionedGrid pg)
+ @staticmethod
+ cdef void sample(VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil
+
+cdef class ProjectionSampler(ImageSampler):
+ pass
+
+cdef class InterpolatedProjectionSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
+
+cdef class VolumeRenderSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
+ cdef kdtree_utils.kdtree **trees
+ cdef object tree_containers
+
+cdef class LightSourceRenderSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/image_samplers.pyx
--- /dev/null
+++ b/yt/utilities/lib/image_samplers.pyx
@@ -0,0 +1,471 @@
+"""
+Image sampler definitions
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from libc.stdlib cimport malloc, calloc, free, abs
+from libc.math cimport exp, floor, log2, \
+ 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,\
+ FIT_eval_transfer_with_light
+cimport lenses
+from .grid_traversal cimport walk_volume
+from .fixed_interpolator cimport \
+ offset_interpolate, \
+ fast_interpolate, \
+ trilinear_interpolate, \
+ eval_gradient, \
+ offset_fill, \
+ vertex_interp
+
+cdef extern from "platform_dep.h":
+ long int lrint(double x) nogil
+
+DEF Nch = 4
+
+from cython.parallel import prange, parallel, threadid
+from vec3_ops cimport dot, subtract, L2_norm, fma
+
+from cpython.exc cimport PyErr_CheckSignals
+
+cdef class ImageSampler:
+ def __init__(self,
+ np.float64_t[:,:,:] vp_pos,
+ np.float64_t[:,:,:] vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ *args, **kwargs):
+ self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
+ cdef np.float64_t[:,:] zbuffer
+ cdef np.int64_t[:,:] image_used
+ cdef np.int64_t[:,:] mesh_lines
+ 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")
+
+ image_used = np.zeros((image.shape[0], image.shape[1]), "int64")
+ mesh_lines = np.zeros((image.shape[0], image.shape[1]), "int64")
+
+ self.lens_type = kwargs.pop("lens_type", None)
+ if self.lens_type == "plane-parallel":
+ self.extent_function = lenses.calculate_extent_plane_parallel
+ self.vector_function = lenses.generate_vector_info_plane_parallel
+ 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 += "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)
+
+ if camera_data is not None and self.lens_type == 'perspective':
+ self.extent_function = lenses.calculate_extent_perspective
+ else:
+ self.extent_function = lenses.calculate_extent_null
+ self.vector_function = lenses.generate_vector_info_null
+
+ # 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
+ # arguments.
+ self.image.vp_pos = vp_pos
+ self.image.vp_dir = vp_dir
+ self.image.image = self.aimage = image
+ self.acenter = center
+ self.image.center = <np.float64_t *> center.data
+ self.ax_vec = x_vec
+ self.image.x_vec = <np.float64_t *> x_vec.data
+ self.ay_vec = y_vec
+ self.image.y_vec = <np.float64_t *> y_vec.data
+ self.image.zbuffer = self.azbuffer = zbuffer
+ self.image.image_used = self.aimage_used = image_used
+ self.image.mesh_lines = self.amesh_lines = mesh_lines
+ self.image.nv[0] = image.shape[0]
+ self.image.nv[1] = image.shape[1]
+ for i in range(4): self.image.bounds[i] = bounds[i]
+ self.image.pdx = (bounds[1] - bounds[0])/self.image.nv[0]
+ self.image.pdy = (bounds[3] - bounds[2])/self.image.nv[1]
+ for i in range(3):
+ self.width[i] = width[i]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __call__(self, PartitionedGrid pg, int num_threads = 0):
+ # This routine will iterate over all of the vectors and cast each in
+ # turn. Might benefit from a more sophisticated intersection check,
+ # like http://courses.csusm.edu/cs697exz/ray_box.htm
+ cdef int vi, vj, hit, i, j
+ cdef np.int64_t iter[4]
+ cdef VolumeContainer *vc = pg.container
+ cdef ImageContainer *im = self.image
+ self.setup(pg)
+ cdef np.float64_t *v_pos
+ cdef np.float64_t *v_dir
+ cdef np.float64_t max_t
+ hit = 0
+ cdef np.int64_t nx, ny, size
+ self.extent_function(self.image, vc, iter)
+ iter[0] = i64clip(iter[0]-1, 0, im.nv[0])
+ iter[1] = i64clip(iter[1]+1, 0, im.nv[0])
+ iter[2] = i64clip(iter[2]-1, 0, im.nv[1])
+ iter[3] = i64clip(iter[3]+1, 0, im.nv[1])
+ nx = (iter[1] - iter[0])
+ ny = (iter[3] - iter[2])
+ size = nx * ny
+ cdef ImageAccumulator *idata
+ cdef np.float64_t width[3]
+ cdef int chunksize = 100
+ for i in range(3):
+ width[i] = self.width[i]
+ with nogil, parallel(num_threads = num_threads):
+ idata = <ImageAccumulator *> malloc(sizeof(ImageAccumulator))
+ idata.supp_data = self.supp_data
+ v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ for j in prange(size, schedule="static", chunksize=chunksize):
+ vj = j % ny
+ vi = (j - vj) / ny + iter[0]
+ vj = vj + iter[2]
+ # Dynamically calculate the position
+ self.vector_function(im, vi, vj, width, v_dir, v_pos)
+ for i in range(Nch):
+ idata.rgba[i] = im.image[vi, vj, i]
+ max_t = fclip(im.zbuffer[vi, vj], 0.0, 1.0)
+ walk_volume(vc, v_pos, v_dir, self.sample,
+ (<void *> idata), NULL, max_t)
+ if (j % (10*chunksize)) == 0:
+ with gil:
+ PyErr_CheckSignals()
+ for i in range(Nch):
+ im.image[vi, vj, i] = idata.rgba[i]
+ free(idata)
+ free(v_pos)
+ free(v_dir)
+ return hit
+
+ cdef void setup(self, PartitionedGrid pg):
+ return
+
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ return
+
+ def ensure_code_unit_params(self, params):
+ for param_name in ['center', 'vp_pos', 'vp_dir', 'width']:
+ param = params[param_name]
+ if hasattr(param, 'in_units'):
+ params[param_name] = param.in_units('code_length')
+ bounds = params['bounds']
+ if hasattr(bounds[0], 'units'):
+ params['bounds'] = tuple(b.in_units('code_length').d for b in bounds)
+
+ return params
+
+ def __dealloc__(self):
+ self.image.image = None
+ self.image.vp_pos = None
+ self.image.vp_dir = None
+ self.image.zbuffer = None
+ self.image.camera_data = None
+ self.image.image_used = None
+ free(self.image)
+
+cdef class ProjectionSampler(ImageSampler):
+
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef int i
+ cdef np.float64_t dl = (exit_t - enter_t)
+ cdef int di = (index[0]*vc.dims[1]+index[1])*vc.dims[2]+index[2]
+ for i in range(imin(4, vc.n_fields)):
+ im.rgba[i] += vc.data[i][di] * dl
+
+
+cdef class InterpolatedProjectionSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ n_samples = 10, **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.n_samples = n_samples
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ for j in range(imin(3, vc.n_fields)):
+ im.rgba[j] += dvs[j] * dt
+ for j in range(3):
+ dp[j] += ds[j]
+
+
+cdef class VolumeRenderSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ tf_obj, n_samples = 10,
+ star_list = None, **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ cdef int i
+ cdef np.ndarray[np.float64_t, ndim=1] temp
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.fits = <FieldInterpolationTable *> \
+ malloc(sizeof(FieldInterpolationTable) * 6)
+ self.vra.n_fits = tf_obj.n_field_tables
+ assert(self.vra.n_fits <= 6)
+ self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
+ self.vra.n_samples = n_samples
+ self.my_field_tables = []
+ for i in range(self.vra.n_fits):
+ temp = tf_obj.tables[i].y
+ FIT_initialize_table(&self.vra.fits[i],
+ temp.shape[0],
+ <np.float64_t *> temp.data,
+ tf_obj.tables[i].x_bounds[0],
+ tf_obj.tables[i].x_bounds[1],
+ tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
+ tf_obj.weight_table_ids[i])
+ self.my_field_tables.append((tf_obj.tables[i],
+ tf_obj.tables[i].y))
+ for i in range(6):
+ self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef int cell_offset = index[0] * (vc.dims[1]) * (vc.dims[2]) \
+ + index[1] * (vc.dims[2]) + index[2]
+ if vc.mask[cell_offset] != 1:
+ return
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits,
+ vri.fits, vri.field_table_ids, vri.grey_opacity)
+ for j in range(3):
+ dp[j] += ds[j]
+
+ def __dealloc__(self):
+ for i in range(self.vra.n_fits):
+ free(self.vra.fits[i].d0)
+ free(self.vra.fits[i].dy)
+ free(self.vra.fits)
+ free(self.vra)
+
+cdef class LightSourceRenderSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ tf_obj, n_samples = 10,
+ light_dir=[1.,1.,1.],
+ light_rgba=[1.,1.,1.,1.],
+ **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ cdef int i
+ cdef np.ndarray[np.float64_t, ndim=1] temp
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.fits = <FieldInterpolationTable *> \
+ malloc(sizeof(FieldInterpolationTable) * 6)
+ self.vra.n_fits = tf_obj.n_field_tables
+ assert(self.vra.n_fits <= 6)
+ self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
+ self.vra.n_samples = n_samples
+ self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
+ self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
+ light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
+ for i in range(3):
+ self.vra.light_dir[i] = light_dir[i]
+ for i in range(4):
+ self.vra.light_rgba[i] = light_rgba[i]
+ self.my_field_tables = []
+ for i in range(self.vra.n_fits):
+ temp = tf_obj.tables[i].y
+ FIT_initialize_table(&self.vra.fits[i],
+ temp.shape[0],
+ <np.float64_t *> temp.data,
+ tf_obj.tables[i].x_bounds[0],
+ tf_obj.tables[i].x_bounds[1],
+ tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
+ tf_obj.weight_table_ids[i])
+ self.my_field_tables.append((tf_obj.tables[i],
+ tf_obj.tables[i].y))
+ for i in range(6):
+ self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ cdef np.float64_t *grad
+ grad = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ eval_gradient(vc.dims, dp, vc.data[0] + offset, grad)
+ FIT_eval_transfer_with_light(dt, dvs, grad,
+ vri.light_dir, vri.light_rgba,
+ im.rgba, vri.n_fits,
+ vri.fits, vri.field_table_ids, vri.grey_opacity)
+ for j in range(3):
+ dp[j] += ds[j]
+ free(grad)
+
+
+ def __dealloc__(self):
+ for i in range(self.vra.n_fits):
+ free(self.vra.fits[i].d0)
+ free(self.vra.fits[i].dy)
+ free(self.vra.light_dir)
+ free(self.vra.light_rgba)
+ free(self.vra.fits)
+ free(self.vra)
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/lenses.pxd
--- /dev/null
+++ b/yt/utilities/lib/lenses.pxd
@@ -0,0 +1,46 @@
+"""
+Definitions for the lens code
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from .volume_container cimport VolumeContainer
+from vec3_ops cimport dot, subtract, L2_norm, fma
+from libc.math cimport exp, floor, log2, \
+ fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
+from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, i64clip
+
+cdef extern from "platform_dep.h":
+ long int lrint(double x) nogil
+
+cdef extern from "limits.h":
+ cdef int SHRT_MAX
+
+cdef struct ImageContainer
+
+ctypedef void calculate_extent_function(ImageContainer *image,
+ VolumeContainer *vc, np.int64_t rv[4]) nogil
+
+ctypedef void generate_vector_info_function(ImageContainer *im,
+ np.int64_t vi, np.int64_t vj,
+ np.float64_t width[2],
+ np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil
+
+cdef generate_vector_info_function generate_vector_info_plane_parallel
+cdef generate_vector_info_function generate_vector_info_null
+cdef calculate_extent_function calculate_extent_plane_parallel
+cdef calculate_extent_function calculate_extent_perspective
+cdef calculate_extent_function calculate_extent_null
diff -r b264176ae1527a32714dbd782c7f5a3e79dc993a -r caf2aa0abab31990427014ab0503c1ba5d87c58b yt/utilities/lib/lenses.pyx
--- /dev/null
+++ b/yt/utilities/lib/lenses.pyx
@@ -0,0 +1,204 @@
+"""
+Functions for computing the extent of lenses and whatnot
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from .image_samplers cimport ImageContainer
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef void calculate_extent_plane_parallel(ImageContainer *image,
+ VolumeContainer *vc, np.int64_t rv[4]) nogil:
+ # We do this for all eight corners
+ cdef np.float64_t temp
+ cdef np.float64_t *edges[2]
+ cdef np.float64_t cx, cy
+ cdef np.float64_t extrema[4]
+ cdef int i, j, k
+ edges[0] = vc.left_edge
+ edges[1] = vc.right_edge
+ extrema[0] = extrema[2] = 1e300; extrema[1] = extrema[3] = -1e300
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ # This should rotate it into the vector plane
+ temp = edges[i][0] * image.x_vec[0]
+ temp += edges[j][1] * image.x_vec[1]
+ temp += edges[k][2] * image.x_vec[2]
+ if temp < extrema[0]: extrema[0] = temp
+ if temp > extrema[1]: extrema[1] = temp
+ temp = edges[i][0] * image.y_vec[0]
+ temp += edges[j][1] * image.y_vec[1]
+ temp += edges[k][2] * image.y_vec[2]
+ if temp < extrema[2]: extrema[2] = temp
+ if temp > extrema[3]: extrema[3] = temp
+ cx = cy = 0.0
+ for i in range(3):
+ cx += image.center[i] * image.x_vec[i]
+ cy += image.center[i] * image.y_vec[i]
+ rv[0] = lrint((extrema[0] - cx - image.bounds[0])/image.pdx)
+ rv[1] = rv[0] + lrint((extrema[1] - extrema[0])/image.pdx)
+ 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]
+
+ cam_width[1] = cam_width[0] * image.nv[1] / image.nv[0]
+
+ 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.
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef void calculate_extent_null(ImageContainer *image,
+ VolumeContainer *vc, np.int64_t rv[4]) nogil:
+ rv[0] = 0
+ rv[1] = image.nv[0]
+ rv[2] = 0
+ rv[3] = image.nv[1]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef void generate_vector_info_plane_parallel(ImageContainer *im,
+ np.int64_t vi, np.int64_t vj,
+ np.float64_t width[2],
+ # Now outbound
+ np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
+ cdef int i
+ cdef np.float64_t px, py
+ px = width[0] * (<np.float64_t>vi)/(<np.float64_t>im.nv[0]-1) - width[0]/2.0
+ py = width[1] * (<np.float64_t>vj)/(<np.float64_t>im.nv[1]-1) - width[1]/2.0
+ # atleast_3d will add to beginning and end
+ v_pos[0] = im.vp_pos[0,0,0]*px + im.vp_pos[0,3,0]*py + im.vp_pos[0,9,0]
+ v_pos[1] = im.vp_pos[0,1,0]*px + im.vp_pos[0,4,0]*py + im.vp_pos[0,10,0]
+ v_pos[2] = im.vp_pos[0,2,0]*px + im.vp_pos[0,5,0]*py + im.vp_pos[0,11,0]
+ for i in range(3): v_dir[i] = im.vp_dir[0,i,0]
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+cdef void generate_vector_info_null(ImageContainer *im,
+ np.int64_t vi, np.int64_t vj,
+ np.float64_t width[2],
+ # Now outbound
+ np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
+ cdef int i
+ for i in range(3):
+ # Here's a funny thing: we use vi here because our *image* will be
+ # flattened. That means that im.nv will be a better one-d offset,
+ # since vp_pos has funny strides.
+ v_pos[i] = im.vp_pos[vi, vj, i]
+ v_dir[i] = im.vp_dir[vi, vj, i]
+
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/yt_analysis/yt/commits/43cf58059a74/
Changeset: 43cf58059a74
Branch: yt
User: MatthewTurk
Date: 2016-07-07 20:28:37+00:00
Summary: Rearranging.
Affected #: 5 files
diff -r caf2aa0abab31990427014ab0503c1ba5d87c58b -r 43cf58059a74820ab85b6c213d6d22237e096ecb yt/utilities/lib/field_interpolation_tables.pxd
--- a/yt/utilities/lib/field_interpolation_tables.pxd
+++ b/yt/utilities/lib/field_interpolation_tables.pxd
@@ -17,6 +17,7 @@
cimport numpy as np
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, fabs
from libc.stdlib cimport malloc
+from libc.math cimport isnormal
cdef struct FieldInterpolationTable:
# Note that we make an assumption about retaining a reference to values
@@ -32,10 +33,6 @@
int weight_table_id
int nbins
-cdef extern from "math.h":
- double expf(double x) nogil
- int isnormal(double x) nogil
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r caf2aa0abab31990427014ab0503c1ba5d87c58b -r 43cf58059a74820ab85b6c213d6d22237e096ecb yt/utilities/lib/image_samplers.pxd
--- a/yt/utilities/lib/image_samplers.pxd
+++ b/yt/utilities/lib/image_samplers.pxd
@@ -21,41 +21,16 @@
from .volume_container cimport VolumeContainer
from .lenses cimport \
calculate_extent_function, \
- generate_vector_info_function
+ generate_vector_info_function, \
+ ImageContainer
from .partitioned_grid cimport PartitionedGrid
-from field_interpolation_tables cimport \
- FieldInterpolationTable
-
DEF Nch = 4
-cdef struct ImageContainer:
- np.float64_t[:,:,:] vp_pos
- np.float64_t[:,:,:] vp_dir
- np.float64_t *center
- np.float64_t[:,:,:] image
- np.float64_t[:,:] zbuffer
- np.int64_t[:,:] image_used
- np.int64_t[:,:] mesh_lines
- 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
-
-cdef struct VolumeRenderAccumulator:
- int n_fits
- int n_samples
- FieldInterpolationTable *fits
- int field_table_ids[6]
- np.float64_t star_coeff
- np.float64_t star_er
- np.float64_t star_sigma_num
- kdtree_utils.kdtree *star_list
- np.float64_t *light_dir
- np.float64_t *light_rgba
- int grey_opacity
+# NOTE: We don't want to import the field_interpolator_tables here, as it
+# breaks a bunch of C++ interop. Maybe some day it won't. So, we just forward
+# declare.
+cdef struct VolumeRenderAccumulator
cdef struct ImageAccumulator:
np.float64_t rgba[Nch]
diff -r caf2aa0abab31990427014ab0503c1ba5d87c58b -r 43cf58059a74820ab85b6c213d6d22237e096ecb yt/utilities/lib/image_samplers.pyx
--- a/yt/utilities/lib/image_samplers.pyx
+++ b/yt/utilities/lib/image_samplers.pyx
@@ -43,6 +43,20 @@
from cpython.exc cimport PyErr_CheckSignals
+cdef struct VolumeRenderAccumulator:
+ int n_fits
+ int n_samples
+ FieldInterpolationTable *fits
+ int field_table_ids[6]
+ np.float64_t star_coeff
+ np.float64_t star_er
+ np.float64_t star_sigma_num
+ kdtree_utils.kdtree *star_list
+ np.float64_t *light_dir
+ np.float64_t *light_rgba
+ int grey_opacity
+
+
cdef class ImageSampler:
def __init__(self,
np.float64_t[:,:,:] vp_pos,
@@ -89,7 +103,8 @@
self.extent_function = lenses.calculate_extent_perspective
else:
self.extent_function = lenses.calculate_extent_null
- self.vector_function = lenses.generate_vector_info_null
+ self.\
+ vector_function = lenses.generate_vector_info_null
# These assignments are so we can track the objects and prevent their
# de-allocation from reference counts. Note that we do this to the
diff -r caf2aa0abab31990427014ab0503c1ba5d87c58b -r 43cf58059a74820ab85b6c213d6d22237e096ecb yt/utilities/lib/lenses.pxd
--- a/yt/utilities/lib/lenses.pxd
+++ b/yt/utilities/lib/lenses.pxd
@@ -29,7 +29,21 @@
cdef extern from "limits.h":
cdef int SHRT_MAX
-cdef struct ImageContainer
+cdef struct ImageContainer:
+ np.float64_t[:,:,:] vp_pos
+ np.float64_t[:,:,:] vp_dir
+ np.float64_t *center
+ np.float64_t[:,:,:] image
+ np.float64_t[:,:] zbuffer
+ np.int64_t[:,:] image_used
+ np.int64_t[:,:] mesh_lines
+ 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
+
ctypedef void calculate_extent_function(ImageContainer *image,
VolumeContainer *vc, np.int64_t rv[4]) nogil
diff -r caf2aa0abab31990427014ab0503c1ba5d87c58b -r 43cf58059a74820ab85b6c213d6d22237e096ecb yt/utilities/lib/mesh_traversal.pyx
--- a/yt/utilities/lib/mesh_traversal.pyx
+++ b/yt/utilities/lib/mesh_traversal.pyx
@@ -21,7 +21,9 @@
cimport pyembree.rtcore_ray as rtcr
cimport pyembree.rtcore_geometry as rtcg
cimport pyembree.rtcore_scene as rtcs
-from grid_traversal cimport ImageSampler, \
+from .image_samplers cimport \
+ ImageSampler
+from .lenses cimport \
ImageContainer
from cython.parallel import prange, parallel, threadid
from yt.visualization.image_writer import apply_colormap
https://bitbucket.org/yt_analysis/yt/commits/8fcfc06227be/
Changeset: 8fcfc06227be
Branch: yt
User: MatthewTurk
Date: 2016-07-07 20:30:43+00:00
Summary: Fixing imports
Affected #: 5 files
diff -r 43cf58059a74820ab85b6c213d6d22237e096ecb -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 yt/analysis_modules/level_sets/contour_finder.py
--- a/yt/analysis_modules/level_sets/contour_finder.py
+++ b/yt/analysis_modules/level_sets/contour_finder.py
@@ -21,7 +21,7 @@
from yt.utilities.lib.contour_finding import \
ContourTree, TileContourTree, link_node_contours, \
update_joins
-from yt.utilities.lib.grid_traversal import \
+from yt.utilities.lib.partitioned_grid import \
PartitionedGrid
def identify_contours(data_source, field, min_val, max_val,
diff -r 43cf58059a74820ab85b6c213d6d22237e096ecb -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 yt/utilities/amr_kdtree/amr_kdtree.py
--- a/yt/utilities/amr_kdtree/amr_kdtree.py
+++ b/yt/utilities/amr_kdtree/amr_kdtree.py
@@ -26,7 +26,7 @@
from yt.utilities.lib.amr_kdtools import Node
from yt.utilities.parallel_tools.parallel_analysis_interface import \
ParallelAnalysisInterface
-from yt.utilities.lib.grid_traversal import PartitionedGrid
+from yt.utilities.lib.partitioned_grid import PartitionedGrid
from yt.utilities.math_utils import periodic_position
steps = np.array([[-1, -1, -1], [-1, -1, 0], [-1, -1, 1],
diff -r 43cf58059a74820ab85b6c213d6d22237e096ecb -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 yt/visualization/volume_rendering/off_axis_projection.py
--- a/yt/visualization/volume_rendering/off_axis_projection.py
+++ b/yt/visualization/volume_rendering/off_axis_projection.py
@@ -17,7 +17,7 @@
from .transfer_functions import ProjectionTransferFunction
from .utils import data_source_or_all
from yt.funcs import mylog, iterable
-from yt.utilities.lib.grid_traversal import \
+from yt.utilities.lib.partitioned_grid import \
PartitionedGrid
from yt.data_objects.api import ImageArray
import numpy as np
diff -r 43cf58059a74820ab85b6c213d6d22237e096ecb -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 yt/visualization/volume_rendering/old_camera.py
--- a/yt/visualization/volume_rendering/old_camera.py
+++ b/yt/visualization/volume_rendering/old_camera.py
@@ -28,9 +28,11 @@
from .transfer_functions import ProjectionTransferFunction
from yt.utilities.lib.grid_traversal import \
- pixelize_healpix, \
- arr_fisheye_vectors, arr_pix2vec_nest, \
- PartitionedGrid, ProjectionSampler, VolumeRenderSampler, \
+ pixelize_healpix, arr_fisheye_vectors, arr_pix2vec_nest
+from yt.utilities.lib.partitioned_grid import \
+ PartitionedGrid
+from yt.utilities.image_samplers import \
+ ProjectionSampler, VolumeRenderSampler, \
LightSourceRenderSampler, InterpolatedProjectionSampler
from yt.utilities.lib.misc_utilities import \
lines
diff -r 43cf58059a74820ab85b6c213d6d22237e096ecb -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 yt/visualization/volume_rendering/utils.py
--- a/yt/visualization/volume_rendering/utils.py
+++ b/yt/visualization/volume_rendering/utils.py
@@ -1,7 +1,7 @@
import numpy as np
from yt.data_objects.static_output import Dataset
from yt.utilities.lib import bounding_volume_hierarchy
-from yt.utilities.lib.grid_traversal import \
+from yt.utilities.lib.image_samplers import \
VolumeRenderSampler, InterpolatedProjectionSampler, ProjectionSampler
from yt.utilities.on_demand_imports import NotAModule
https://bitbucket.org/yt_analysis/yt/commits/6c2f519faf2f/
Changeset: 6c2f519faf2f
Branch: yt
User: MatthewTurk
Date: 2016-07-07 20:39:02+00:00
Summary: Fixing setup to get the right build arguments for fixed_interpolator
Affected #: 3 files
diff -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 -r 6c2f519faf2f8ea61530341d1dd77819418a7c5b setup.py
--- a/setup.py
+++ b/setup.py
@@ -174,6 +174,20 @@
extra_link_args=omp_args,
depends=["yt/utilities/lib/kdtree.h",
"yt/utilities/lib/fixed_interpolator.h"]),
+ Extension("yt.utilities.lib.image_samplers",
+ ["yt/utilities/lib/image_samplers.pyx",
+ "yt/utilities/lib/fixed_interpolator.c"],
+ include_dirs=["yt/utilities/lib/"],
+ libraries=std_libs,
+ extra_compile_args=omp_args,
+ extra_link_args=omp_args,
+ depends=["yt/utilities/lib/fixed_interpolator.h"]),
+ Extension("yt.utilities.lib.partitioned_grid",
+ ["yt/utilities/lib/partitioned_grid.pyx",
+ "yt/utilities/lib/fixed_interpolator.c"],
+ include_dirs=["yt/utilities/lib/"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/fixed_interpolator.h"]),
Extension("yt.utilities.lib.element_mappings",
["yt/utilities/lib/element_mappings.pyx"],
libraries=std_libs),
@@ -186,7 +200,7 @@
"particle_mesh_operations", "depth_first_octree", "fortran_reader",
"interpolators", "misc_utilities", "basic_octree", "image_utilities",
"points_in_volume", "quad_tree", "ray_integrators", "mesh_utilities",
- "amr_kdtools", "image_samplers", "lenses"
+ "amr_kdtools", "lenses",
]
for ext_name in lib_exts:
cython_extensions.append(
diff -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 -r 6c2f519faf2f8ea61530341d1dd77819418a7c5b yt/utilities/lib/partitioned_grid.pxd
--- a/yt/utilities/lib/partitioned_grid.pxd
+++ b/yt/utilities/lib/partitioned_grid.pxd
@@ -17,7 +17,6 @@
import numpy as np
cimport numpy as np
cimport cython
-cimport kdtree_utils
from .volume_container cimport VolumeContainer
cdef class PartitionedGrid:
@@ -27,7 +26,6 @@
cdef public object RightEdge
cdef public int parent_grid_id
cdef VolumeContainer *container
- cdef kdtree_utils.kdtree *star_list
cdef np.float64_t star_er
cdef np.float64_t star_sigma_num
cdef np.float64_t star_coeff
diff -r 8fcfc06227be8ac44e13743e2ad5db8509df80a4 -r 6c2f519faf2f8ea61530341d1dd77819418a7c5b yt/utilities/lib/partitioned_grid.pyx
--- a/yt/utilities/lib/partitioned_grid.pyx
+++ b/yt/utilities/lib/partitioned_grid.pyx
@@ -17,6 +17,7 @@
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, calloc, free, abs
+from .fixed_interpolator cimport offset_interpolate
cdef class PartitionedGrid:
@@ -28,8 +29,7 @@
mask,
np.ndarray[np.float64_t, ndim=1] left_edge,
np.ndarray[np.float64_t, ndim=1] right_edge,
- np.ndarray[np.int64_t, ndim=1] dims,
- star_kdtree_container star_tree = None):
+ np.ndarray[np.int64_t, ndim=1] dims):
# The data is likely brought in via a slice, so we copy it
cdef np.ndarray[np.float64_t, ndim=3] tdata
cdef np.ndarray[np.uint8_t, ndim=3] mask_data
@@ -56,16 +56,6 @@
tdata = data[i]
c.data[i] = <np.float64_t *> tdata.data
c.mask = <np.uint8_t *> mask_data.data
- if star_tree is None:
- self.star_list = NULL
- else:
- self.set_star_tree(star_tree)
-
- def set_star_tree(self, star_kdtree_container star_tree):
- self.star_list = star_tree.tree
- self.star_sigma_num = 2.0*star_tree.sigma**2.0
- self.star_er = 2.326 * star_tree.sigma
- self.star_coeff = star_tree.coeff
def __dealloc__(self):
# The data fields are not owned by the container, they are owned by us!
https://bitbucket.org/yt_analysis/yt/commits/a6f4cd1d07f7/
Changeset: a6f4cd1d07f7
Branch: yt
User: MatthewTurk
Date: 2016-07-07 20:41:18+00:00
Summary: Replacing ensure_code_unit_params
Affected #: 2 files
diff -r 6c2f519faf2f8ea61530341d1dd77819418a7c5b -r a6f4cd1d07f763a1177f9cc296d57d0a3ea15ed3 yt/visualization/volume_rendering/old_camera.py
--- a/yt/visualization/volume_rendering/old_camera.py
+++ b/yt/visualization/volume_rendering/old_camera.py
@@ -31,7 +31,7 @@
pixelize_healpix, arr_fisheye_vectors, arr_pix2vec_nest
from yt.utilities.lib.partitioned_grid import \
PartitionedGrid
-from yt.utilities.image_samplers import \
+from yt.utilities.lib.image_samplers import \
ProjectionSampler, VolumeRenderSampler, \
LightSourceRenderSampler, InterpolatedProjectionSampler
from yt.utilities.lib.misc_utilities import \
diff -r 6c2f519faf2f8ea61530341d1dd77819418a7c5b -r a6f4cd1d07f763a1177f9cc296d57d0a3ea15ed3 yt/visualization/volume_rendering/utils.py
--- a/yt/visualization/volume_rendering/utils.py
+++ b/yt/visualization/volume_rendering/utils.py
@@ -125,3 +125,14 @@
[re[0], re[1], re[2]],
[le[0], re[1], re[2]],
], dtype='float64')
+
+def ensure_code_unit_params(params):
+ for param_name in ['center', 'vp_pos', 'vp_dir', 'width']:
+ param = params[param_name]
+ if hasattr(param, 'in_units'):
+ params[param_name] = param.in_units('code_length')
+ bounds = params['bounds']
+ if hasattr(bounds[0], 'units'):
+ params['bounds'] = tuple(b.in_units('code_length').d for b in bounds)
+ return params
+
https://bitbucket.org/yt_analysis/yt/commits/6a4eea9e6d94/
Changeset: 6a4eea9e6d94
Branch: yt
User: ngoldbaum
Date: 2016-07-16 17:41:31+00:00
Summary: Merged in MatthewTurk/yt (pull request #2264)
Refactor volume rendering Cython code
Affected #: 20 files
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -62,6 +62,11 @@
yt/utilities/lib/marching_cubes.c
yt/utilities/lib/png_writer.h
yt/utilities/lib/write_array.c
+yt/utilities/lib/perftools_wrap.c
+yt/utilities/lib/partitioned_grid.c
+yt/utilities/lib/volume_container.c
+yt/utilities/lib/lenses.c
+yt/utilities/lib/image_samplers.c
syntax: glob
*.pyc
*.pyd
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 setup.py
--- a/setup.py
+++ b/setup.py
@@ -175,6 +175,20 @@
extra_link_args=omp_args,
depends=["yt/utilities/lib/kdtree.h",
"yt/utilities/lib/fixed_interpolator.h"]),
+ Extension("yt.utilities.lib.image_samplers",
+ ["yt/utilities/lib/image_samplers.pyx",
+ "yt/utilities/lib/fixed_interpolator.c"],
+ include_dirs=["yt/utilities/lib/"],
+ libraries=std_libs,
+ extra_compile_args=omp_args,
+ extra_link_args=omp_args,
+ depends=["yt/utilities/lib/fixed_interpolator.h"]),
+ Extension("yt.utilities.lib.partitioned_grid",
+ ["yt/utilities/lib/partitioned_grid.pyx",
+ "yt/utilities/lib/fixed_interpolator.c"],
+ include_dirs=["yt/utilities/lib/"],
+ libraries=std_libs,
+ depends=["yt/utilities/lib/fixed_interpolator.h"]),
Extension("yt.utilities.lib.element_mappings",
["yt/utilities/lib/element_mappings.pyx"],
libraries=std_libs),
@@ -187,7 +201,7 @@
"particle_mesh_operations", "depth_first_octree", "fortran_reader",
"interpolators", "misc_utilities", "basic_octree", "image_utilities",
"points_in_volume", "quad_tree", "ray_integrators", "mesh_utilities",
- "amr_kdtools"
+ "amr_kdtools", "lenses",
]
for ext_name in lib_exts:
cython_extensions.append(
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/analysis_modules/level_sets/contour_finder.py
--- a/yt/analysis_modules/level_sets/contour_finder.py
+++ b/yt/analysis_modules/level_sets/contour_finder.py
@@ -21,7 +21,7 @@
from yt.utilities.lib.contour_finding import \
ContourTree, TileContourTree, link_node_contours, \
update_joins
-from yt.utilities.lib.grid_traversal import \
+from yt.utilities.lib.partitioned_grid import \
PartitionedGrid
def identify_contours(data_source, field, min_val, max_val,
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -22,8 +22,10 @@
from .oct_container cimport OctreeContainer, OctAllocationContainer, Oct
cimport oct_visitors
from .oct_visitors cimport cind
+from yt.utilities.lib.volume_container cimport \
+ VolumeContainer
from yt.utilities.lib.grid_traversal cimport \
- VolumeContainer, sample_function, walk_volume
+ sampler_function, walk_volume
from yt.utilities.lib.bitarray cimport ba_get_value, ba_set_value
cdef extern from "math.h":
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/amr_kdtree/amr_kdtree.py
--- a/yt/utilities/amr_kdtree/amr_kdtree.py
+++ b/yt/utilities/amr_kdtree/amr_kdtree.py
@@ -28,7 +28,7 @@
from yt.utilities.lib.amr_kdtools import Node
from yt.utilities.parallel_tools.parallel_analysis_interface import \
ParallelAnalysisInterface
-from yt.utilities.lib.grid_traversal import PartitionedGrid
+from yt.utilities.lib.partitioned_grid import PartitionedGrid
from yt.utilities.math_utils import periodic_position
steps = np.array([[-1, -1, -1], [-1, -1, 0], [-1, -1, 1],
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/contour_finding.pyx
--- a/yt/utilities/lib/contour_finding.pyx
+++ b/yt/utilities/lib/contour_finding.pyx
@@ -26,8 +26,10 @@
from yt.geometry.oct_visitors cimport \
Oct
from .amr_kdtools cimport Node
-from .grid_traversal cimport VolumeContainer, PartitionedGrid, \
- vc_index, vc_pos_index
+from .partitioned_grid cimport \
+ PartitionedGrid
+from .volume_container cimport \
+ VolumeContainer, vc_index, vc_pos_index
import sys
cdef inline ContourID *contour_create(np.int64_t contour_id,
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/field_interpolation_tables.pxd
--- a/yt/utilities/lib/field_interpolation_tables.pxd
+++ b/yt/utilities/lib/field_interpolation_tables.pxd
@@ -17,8 +17,7 @@
cimport numpy as np
from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, fabs
from libc.stdlib cimport malloc
-
-DEF Nch = 4
+from libc.math cimport isnormal
cdef struct FieldInterpolationTable:
# Note that we make an assumption about retaining a reference to values
@@ -34,10 +33,6 @@
int weight_table_id
int nbins
-cdef extern from "math.h":
- double expf(double x) nogil
- int isnormal(double x) nogil
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/grid_traversal.pxd
--- a/yt/utilities/lib/grid_traversal.pxd
+++ b/yt/utilities/lib/grid_traversal.pxd
@@ -17,22 +17,8 @@
import numpy as np
cimport numpy as np
cimport cython
-cimport kdtree_utils
-
-cdef struct ImageContainer:
- np.float64_t[:,:,:] vp_pos
- np.float64_t[:,:,:] vp_dir
- np.float64_t *center
- np.float64_t[:,:,:] image
- np.float64_t[:,:] zbuffer
- np.int64_t[:,:] image_used
- np.int64_t[:,:] mesh_lines
- 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
+from .image_samplers cimport ImageContainer, ImageSampler
+from .volume_container cimport VolumeContainer, vc_index, vc_pos_index
ctypedef void sampler_function(
VolumeContainer *vc,
@@ -43,82 +29,11 @@
int index[3],
void *data) nogil
-ctypedef void calculate_extent_function(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil
-
-cdef calculate_extent_function calculate_extent_plane_parallel
-
-ctypedef void generate_vector_info_function(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil
-
-cdef generate_vector_info_function generate_vector_info_plane_parallel
-cdef generate_vector_info_function generate_vector_info_null
-
-cdef class ImageSampler:
- cdef ImageContainer *image
- cdef sampler_function *sampler
- cdef public object acenter, aimage, ax_vec, ay_vec
- cdef public object azbuffer
- cdef public object aimage_used
- cdef public object amesh_lines
- cdef void *supp_data
- cdef np.float64_t width[3]
- cdef public object lens_type
- cdef calculate_extent_function *extent_function
- cdef generate_vector_info_function *vector_function
-
- cdef void setup(self, PartitionedGrid pg)
-
-cdef struct VolumeContainer:
- int n_fields
- np.float64_t **data
- # The mask has dimensions one fewer in each direction than data
- np.uint8_t *mask
- np.float64_t left_edge[3]
- np.float64_t right_edge[3]
- np.float64_t dds[3]
- np.float64_t idds[3]
- int dims[3]
-
-cdef class PartitionedGrid:
- cdef public object my_data
- cdef public object source_mask
- cdef public object LeftEdge
- cdef public object RightEdge
- cdef public int parent_grid_id
- cdef VolumeContainer *container
- cdef kdtree_utils.kdtree *star_list
- cdef np.float64_t star_er
- cdef np.float64_t star_sigma_num
- cdef np.float64_t star_coeff
- cdef void get_vector_field(self, np.float64_t pos[3],
- np.float64_t *vel, np.float64_t *vel_mag)
-
-ctypedef void sample_function(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil
-
cdef int walk_volume(VolumeContainer *vc,
np.float64_t v_pos[3],
np.float64_t v_dir[3],
- sample_function *sampler,
+ sampler_function *sampler,
void *data,
np.float64_t *return_t = *,
np.float64_t max_t = *) nogil
-cdef inline int vc_index(VolumeContainer *vc, int i, int j, int k):
- return (i*vc.dims[1]+j)*vc.dims[2]+k
-
-cdef inline int vc_pos_index(VolumeContainer *vc, np.float64_t *spos):
- cdef int index[3]
- cdef int i
- for i in range(3):
- index[i] = <int> ((spos[i] - vc.left_edge[i]) * vc.idds[i])
- return vc_index(vc, index[0], index[1], index[2])
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/grid_traversal.pyx
--- a/yt/utilities/lib/grid_traversal.pyx
+++ b/yt/utilities/lib/grid_traversal.pyx
@@ -17,8 +17,6 @@
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, \
fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
@@ -28,913 +26,23 @@
FIT_eval_transfer_with_light
from fixed_interpolator cimport *
-cdef extern from "platform_dep.h":
- long int lrint(double x) nogil
-
from cython.parallel import prange, parallel, threadid
-from vec3_ops cimport dot, subtract, L2_norm, fma
-
from cpython.exc cimport PyErr_CheckSignals
+from .image_samplers cimport \
+ ImageSampler, \
+ ImageContainer, \
+ VolumeRenderAccumulator
+
DEF Nch = 4
-cdef class PartitionedGrid:
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def __cinit__(self,
- int parent_grid_id, data,
- mask,
- np.ndarray[np.float64_t, ndim=1] left_edge,
- np.ndarray[np.float64_t, ndim=1] right_edge,
- np.ndarray[np.int64_t, ndim=1] dims,
- star_kdtree_container star_tree = None):
- # The data is likely brought in via a slice, so we copy it
- cdef np.ndarray[np.float64_t, ndim=3] tdata
- cdef np.ndarray[np.uint8_t, ndim=3] mask_data
- self.container = NULL
- self.parent_grid_id = parent_grid_id
- self.LeftEdge = left_edge
- self.RightEdge = right_edge
- self.container = <VolumeContainer *> \
- malloc(sizeof(VolumeContainer))
- cdef VolumeContainer *c = self.container # convenience
- cdef int n_fields = len(data)
- c.n_fields = n_fields
- for i in range(3):
- c.left_edge[i] = left_edge[i]
- c.right_edge[i] = right_edge[i]
- c.dims[i] = dims[i]
- c.dds[i] = (c.right_edge[i] - c.left_edge[i])/dims[i]
- c.idds[i] = 1.0/c.dds[i]
- self.my_data = data
- self.source_mask = mask
- mask_data = mask
- c.data = <np.float64_t **> malloc(sizeof(np.float64_t*) * n_fields)
- for i in range(n_fields):
- tdata = data[i]
- c.data[i] = <np.float64_t *> tdata.data
- c.mask = <np.uint8_t *> mask_data.data
- if star_tree is None:
- self.star_list = NULL
- else:
- self.set_star_tree(star_tree)
-
- def set_star_tree(self, star_kdtree_container star_tree):
- self.star_list = star_tree.tree
- self.star_sigma_num = 2.0*star_tree.sigma**2.0
- self.star_er = 2.326 * star_tree.sigma
- self.star_coeff = star_tree.coeff
-
- def __dealloc__(self):
- # The data fields are not owned by the container, they are owned by us!
- # So we don't need to deallocate them.
- if self.container == NULL: return
- if self.container.data != NULL: free(self.container.data)
- free(self.container)
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def integrate_streamline(self, pos, np.float64_t h, mag):
- cdef np.float64_t cmag[1]
- cdef np.float64_t k1[3]
- cdef np.float64_t k2[3]
- cdef np.float64_t k3[3]
- cdef np.float64_t k4[3]
- cdef np.float64_t newpos[3]
- cdef np.float64_t oldpos[3]
- for i in range(3):
- newpos[i] = oldpos[i] = pos[i]
- self.get_vector_field(newpos, k1, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + 0.5*k1[i]*h
-
- if not (self.LeftEdge[0] < newpos[0] and newpos[0] < self.RightEdge[0] and \
- self.LeftEdge[1] < newpos[1] and newpos[1] < self.RightEdge[1] and \
- self.LeftEdge[2] < newpos[2] and newpos[2] < self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k2, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + 0.5*k2[i]*h
-
- if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
- self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
- self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k3, cmag)
- for i in range(3):
- newpos[i] = oldpos[i] + k3[i]*h
-
- if not (self.LeftEdge[0] <= newpos[0] and newpos[0] <= self.RightEdge[0] and \
- self.LeftEdge[1] <= newpos[1] and newpos[1] <= self.RightEdge[1] and \
- self.LeftEdge[2] <= newpos[2] and newpos[2] <= self.RightEdge[2]):
- if mag is not None:
- mag[0] = cmag[0]
- for i in range(3):
- pos[i] = newpos[i]
- return
-
- self.get_vector_field(newpos, k4, cmag)
-
- for i in range(3):
- pos[i] = oldpos[i] + h*(k1[i]/6.0 + k2[i]/3.0 + k3[i]/3.0 + k4[i]/6.0)
-
- if mag is not None:
- for i in range(3):
- newpos[i] = pos[i]
- self.get_vector_field(newpos, k4, cmag)
- mag[0] = cmag[0]
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- cdef void get_vector_field(self, np.float64_t pos[3],
- np.float64_t *vel, np.float64_t *vel_mag):
- cdef np.float64_t dp[3]
- cdef int ci[3]
- cdef VolumeContainer *c = self.container # convenience
-
- for i in range(3):
- ci[i] = (int)((pos[i]-self.LeftEdge[i])/c.dds[i])
- dp[i] = (pos[i] - ci[i]*c.dds[i] - self.LeftEdge[i])/c.dds[i]
-
- cdef int offset = ci[0] * (c.dims[1] + 1) * (c.dims[2] + 1) \
- + ci[1] * (c.dims[2] + 1) + ci[2]
-
- vel_mag[0] = 0.0
- for i in range(3):
- vel[i] = offset_interpolate(c.dims, dp, c.data[i] + offset)
- vel_mag[0] += vel[i]*vel[i]
- vel_mag[0] = np.sqrt(vel_mag[0])
- if vel_mag[0] != 0.0:
- for i in range(3):
- vel[i] /= vel_mag[0]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void calculate_extent_plane_parallel(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil:
- # We do this for all eight corners
- cdef np.float64_t temp
- cdef np.float64_t *edges[2]
- cdef np.float64_t cx, cy
- cdef np.float64_t extrema[4]
- cdef int i, j, k
- edges[0] = vc.left_edge
- edges[1] = vc.right_edge
- extrema[0] = extrema[2] = 1e300; extrema[1] = extrema[3] = -1e300
- for i in range(2):
- for j in range(2):
- for k in range(2):
- # This should rotate it into the vector plane
- temp = edges[i][0] * image.x_vec[0]
- temp += edges[j][1] * image.x_vec[1]
- temp += edges[k][2] * image.x_vec[2]
- if temp < extrema[0]: extrema[0] = temp
- if temp > extrema[1]: extrema[1] = temp
- temp = edges[i][0] * image.y_vec[0]
- temp += edges[j][1] * image.y_vec[1]
- temp += edges[k][2] * image.y_vec[2]
- if temp < extrema[2]: extrema[2] = temp
- if temp > extrema[3]: extrema[3] = temp
- cx = cy = 0.0
- for i in range(3):
- cx += image.center[i] * image.x_vec[i]
- cy += image.center[i] * image.y_vec[i]
- rv[0] = lrint((extrema[0] - cx - image.bounds[0])/image.pdx)
- rv[1] = rv[0] + lrint((extrema[1] - extrema[0])/image.pdx)
- 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]
-
- cam_width[1] = cam_width[0] * image.nv[1] / image.nv[0]
-
- 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.
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void calculate_extent_null(ImageContainer *image,
- VolumeContainer *vc, np.int64_t rv[4]) nogil:
- rv[0] = 0
- rv[1] = image.nv[0]
- rv[2] = 0
- rv[3] = image.nv[1]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void generate_vector_info_plane_parallel(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- # Now outbound
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
- cdef int i
- cdef np.float64_t px, py
- px = width[0] * (<np.float64_t>vi)/(<np.float64_t>im.nv[0]-1) - width[0]/2.0
- py = width[1] * (<np.float64_t>vj)/(<np.float64_t>im.nv[1]-1) - width[1]/2.0
- # atleast_3d will add to beginning and end
- v_pos[0] = im.vp_pos[0,0,0]*px + im.vp_pos[0,3,0]*py + im.vp_pos[0,9,0]
- v_pos[1] = im.vp_pos[0,1,0]*px + im.vp_pos[0,4,0]*py + im.vp_pos[0,10,0]
- v_pos[2] = im.vp_pos[0,2,0]*px + im.vp_pos[0,5,0]*py + im.vp_pos[0,11,0]
- for i in range(3): v_dir[i] = im.vp_dir[0,i,0]
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
-cdef void generate_vector_info_null(ImageContainer *im,
- np.int64_t vi, np.int64_t vj,
- np.float64_t width[2],
- # Now outbound
- np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil:
- cdef int i
- for i in range(3):
- # Here's a funny thing: we use vi here because our *image* will be
- # flattened. That means that im.nv will be a better one-d offset,
- # since vp_pos has funny strides.
- v_pos[i] = im.vp_pos[vi, vj, i]
- v_dir[i] = im.vp_dir[vi, vj, i]
-
-cdef struct ImageAccumulator:
- np.float64_t rgba[Nch]
- void *supp_data
-
-cdef class ImageSampler:
- def __init__(self,
- np.float64_t[:,:,:] vp_pos,
- np.float64_t[:,:,:] vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- *args, **kwargs):
- self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
- cdef np.float64_t[:,:] zbuffer
- cdef np.int64_t[:,:] image_used
- cdef np.int64_t[:,:] mesh_lines
- 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")
-
- image_used = np.zeros((image.shape[0], image.shape[1]), "int64")
- mesh_lines = np.zeros((image.shape[0], image.shape[1]), "int64")
-
- self.lens_type = kwargs.pop("lens_type", None)
- if self.lens_type == "plane-parallel":
- self.extent_function = calculate_extent_plane_parallel
- self.vector_function = generate_vector_info_plane_parallel
- 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 += "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)
-
- if camera_data is not None and self.lens_type == 'perspective':
- self.extent_function = calculate_extent_perspective
- else:
- self.extent_function = calculate_extent_null
- self.vector_function = generate_vector_info_null
-
- self.sampler = NULL
- # 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
- # arguments.
- self.image.vp_pos = vp_pos
- self.image.vp_dir = vp_dir
- self.image.image = self.aimage = image
- self.acenter = center
- self.image.center = <np.float64_t *> center.data
- self.ax_vec = x_vec
- self.image.x_vec = <np.float64_t *> x_vec.data
- self.ay_vec = y_vec
- self.image.y_vec = <np.float64_t *> y_vec.data
- self.image.zbuffer = self.azbuffer = zbuffer
- self.image.image_used = self.aimage_used = image_used
- self.image.mesh_lines = self.amesh_lines = mesh_lines
- self.image.nv[0] = image.shape[0]
- self.image.nv[1] = image.shape[1]
- for i in range(4): self.image.bounds[i] = bounds[i]
- self.image.pdx = (bounds[1] - bounds[0])/self.image.nv[0]
- self.image.pdy = (bounds[3] - bounds[2])/self.image.nv[1]
- for i in range(3):
- self.width[i] = width[i]
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- @cython.cdivision(True)
- def __call__(self, PartitionedGrid pg, int num_threads = 0):
- # This routine will iterate over all of the vectors and cast each in
- # turn. Might benefit from a more sophisticated intersection check,
- # like http://courses.csusm.edu/cs697exz/ray_box.htm
- cdef int vi, vj, hit, i, j
- cdef np.int64_t iter[4]
- cdef VolumeContainer *vc = pg.container
- cdef ImageContainer *im = self.image
- self.setup(pg)
- if self.sampler == NULL: raise RuntimeError
- cdef np.float64_t *v_pos
- cdef np.float64_t *v_dir
- cdef np.float64_t max_t
- hit = 0
- cdef np.int64_t nx, ny, size
- self.extent_function(self.image, vc, iter)
- iter[0] = i64clip(iter[0]-1, 0, im.nv[0])
- iter[1] = i64clip(iter[1]+1, 0, im.nv[0])
- iter[2] = i64clip(iter[2]-1, 0, im.nv[1])
- iter[3] = i64clip(iter[3]+1, 0, im.nv[1])
- nx = (iter[1] - iter[0])
- ny = (iter[3] - iter[2])
- size = nx * ny
- cdef ImageAccumulator *idata
- cdef np.float64_t width[3]
- cdef int chunksize = 100
- for i in range(3):
- width[i] = self.width[i]
- with nogil, parallel(num_threads = num_threads):
- idata = <ImageAccumulator *> malloc(sizeof(ImageAccumulator))
- idata.supp_data = self.supp_data
- v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- for j in prange(size, schedule="static", chunksize=chunksize):
- vj = j % ny
- vi = (j - vj) / ny + iter[0]
- vj = vj + iter[2]
- # Dynamically calculate the position
- self.vector_function(im, vi, vj, width, v_dir, v_pos)
- for i in range(Nch):
- idata.rgba[i] = im.image[vi, vj, i]
- max_t = fclip(im.zbuffer[vi, vj], 0.0, 1.0)
- walk_volume(vc, v_pos, v_dir, self.sampler,
- (<void *> idata), NULL, max_t)
- if (j % (10*chunksize)) == 0:
- with gil:
- PyErr_CheckSignals()
- for i in range(Nch):
- im.image[vi, vj, i] = idata.rgba[i]
- free(idata)
- free(v_pos)
- free(v_dir)
- return hit
-
- cdef void setup(self, PartitionedGrid pg):
- return
-
- def __dealloc__(self):
- self.image.image = None
- self.image.vp_pos = None
- self.image.vp_dir = None
- self.image.zbuffer = None
- self.image.camera_data = None
- self.image.image_used = None
- free(self.image)
-
-
-cdef void projection_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef int i
- cdef np.float64_t dl = (exit_t - enter_t)
- cdef int di = (index[0]*vc.dims[1]+index[1])*vc.dims[2]+index[2]
- for i in range(imin(4, vc.n_fields)):
- im.rgba[i] += vc.data[i][di] * dl
-
-cdef struct VolumeRenderAccumulator:
- int n_fits
- int n_samples
- FieldInterpolationTable *fits
- int field_table_ids[6]
- np.float64_t star_coeff
- np.float64_t star_er
- np.float64_t star_sigma_num
- kdtree_utils.kdtree *star_list
- np.float64_t *light_dir
- np.float64_t *light_rgba
- int grey_opacity
-
-
-cdef class ProjectionSampler(ImageSampler):
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = projection_sampler
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void interpolated_projection_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- for j in range(imin(3, vc.n_fields)):
- im.rgba[j] += dvs[j] * dt
- for j in range(3):
- dp[j] += ds[j]
-
-cdef class InterpolatedProjectionSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- n_samples = 10, **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.n_samples = n_samples
- self.supp_data = <void *> self.vra
-
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = interpolated_projection_sampler
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef int cell_offset = index[0] * (vc.dims[1]) * (vc.dims[2]) \
- + index[1] * (vc.dims[2]) + index[2]
- if vc.mask[cell_offset] != 1:
- return
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits,
- vri.fits, vri.field_table_ids, vri.grey_opacity)
- for j in range(3):
- dp[j] += ds[j]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_gradient_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- cdef np.float64_t *grad
- grad = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vri.n_samples):
- for j in range(vc.n_fields):
- dvs[j] = offset_interpolate(vc.dims, dp,
- vc.data[j] + offset)
- eval_gradient(vc.dims, dp, vc.data[0] + offset, grad)
- FIT_eval_transfer_with_light(dt, dvs, grad,
- vri.light_dir, vri.light_rgba,
- im.rgba, vri.n_fits,
- vri.fits, vri.field_table_ids, vri.grey_opacity)
- for j in range(3):
- dp[j] += ds[j]
- free(grad)
-
-cdef class star_kdtree_container:
- cdef kdtree_utils.kdtree *tree
- cdef public np.float64_t sigma
- cdef public np.float64_t coeff
-
- def __init__(self):
- self.tree = kdtree_utils.kd_create(3)
-
- def add_points(self,
- np.ndarray[np.float64_t, ndim=1] pos_x,
- np.ndarray[np.float64_t, ndim=1] pos_y,
- np.ndarray[np.float64_t, ndim=1] pos_z,
- np.ndarray[np.float64_t, ndim=2] star_colors):
- cdef int i
- cdef np.float64_t *pointer = <np.float64_t *> star_colors.data
- for i in range(pos_x.shape[0]):
- kdtree_utils.kd_insert3(self.tree,
- pos_x[i], pos_y[i], pos_z[i], <void *> (pointer + i*3))
-
- def __dealloc__(self):
- kdtree_utils.kd_free(self.tree)
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void volume_render_stars_sampler(
- VolumeContainer *vc,
- np.float64_t v_pos[3],
- np.float64_t v_dir[3],
- np.float64_t enter_t,
- np.float64_t exit_t,
- int index[3],
- void *data) nogil:
- cdef ImageAccumulator *im = <ImageAccumulator *> data
- cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
- im.supp_data
- cdef kdtree_utils.kdres *ballq = NULL
- # we assume this has vertex-centered data.
- cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
- + index[1] * (vc.dims[2] + 1) + index[2]
- cdef np.float64_t slopes[6]
- cdef np.float64_t dp[3]
- cdef np.float64_t ds[3]
- cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
- cdef np.float64_t dvs[6]
- cdef np.float64_t cell_left[3]
- cdef np.float64_t local_dds[3]
- cdef np.float64_t pos[3]
- cdef int nstars, i, j
- cdef np.float64_t *colors = NULL
- cdef np.float64_t gexp, gaussian, px, py, pz
- px = py = pz = -1
- for i in range(3):
- dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
- dp[i] *= vc.idds[i]
- ds[i] = v_dir[i] * vc.idds[i] * dt
- for i in range(vc.n_fields):
- slopes[i] = offset_interpolate(vc.dims, dp,
- vc.data[i] + offset)
- cdef np.float64_t temp
- # Now we get the ball-tree result for the stars near our cell center.
- for i in range(3):
- cell_left[i] = index[i] * vc.dds[i] + vc.left_edge[i]
- pos[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
- local_dds[i] = v_dir[i] * dt
- ballq = kdtree_utils.kd_nearest_range3(
- vri.star_list, cell_left[0] + vc.dds[0]*0.5,
- cell_left[1] + vc.dds[1]*0.5,
- cell_left[2] + vc.dds[2]*0.5,
- vri.star_er + 0.9*vc.dds[0])
- # ~0.866 + a bit
-
- nstars = kdtree_utils.kd_res_size(ballq)
- for i in range(vc.n_fields):
- temp = slopes[i]
- slopes[i] -= offset_interpolate(vc.dims, dp,
- vc.data[i] + offset)
- slopes[i] *= -1.0/vri.n_samples
- dvs[i] = temp
- for _ in range(vri.n_samples):
- # Now we add the contribution from stars
- kdtree_utils.kd_res_rewind(ballq)
- for i in range(nstars):
- kdtree_utils.kd_res_item3(ballq, &px, &py, &pz)
- colors = <np.float64_t *> kdtree_utils.kd_res_item_data(ballq)
- kdtree_utils.kd_res_next(ballq)
- gexp = (px - pos[0])*(px - pos[0]) \
- + (py - pos[1])*(py - pos[1]) \
- + (pz - pos[2])*(pz - pos[2])
- gaussian = vri.star_coeff * exp(-gexp/vri.star_sigma_num)
- for j in range(3): im.rgba[j] += gaussian*dt*colors[j]
- for i in range(3):
- pos[i] += local_dds[i]
- FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits, vri.fits,
- vri.field_table_ids, vri.grey_opacity)
- for i in range(vc.n_fields):
- dvs[i] += slopes[i]
- kdtree_utils.kd_res_free(ballq)
-
-cdef class VolumeRenderSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- cdef kdtree_utils.kdtree **trees
- cdef object tree_containers
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- tf_obj, n_samples = 10,
- star_list = None, **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- cdef int i
- cdef np.ndarray[np.float64_t, ndim=1] temp
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.fits = <FieldInterpolationTable *> \
- malloc(sizeof(FieldInterpolationTable) * 6)
- self.vra.n_fits = tf_obj.n_field_tables
- assert(self.vra.n_fits <= 6)
- self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
- self.vra.n_samples = n_samples
- self.my_field_tables = []
- for i in range(self.vra.n_fits):
- temp = tf_obj.tables[i].y
- FIT_initialize_table(&self.vra.fits[i],
- temp.shape[0],
- <np.float64_t *> temp.data,
- tf_obj.tables[i].x_bounds[0],
- tf_obj.tables[i].x_bounds[1],
- tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
- tf_obj.weight_table_ids[i])
- self.my_field_tables.append((tf_obj.tables[i],
- tf_obj.tables[i].y))
- for i in range(6):
- self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
- self.supp_data = <void *> self.vra
- cdef star_kdtree_container skdc
- self.tree_containers = star_list
- if star_list is None:
- self.trees = NULL
- else:
- self.trees = <kdtree_utils.kdtree **> malloc(
- sizeof(kdtree_utils.kdtree*) * len(star_list))
- for i in range(len(star_list)):
- skdc = star_list[i]
- self.trees[i] = skdc.tree
-
- cdef void setup(self, PartitionedGrid pg):
- cdef star_kdtree_container star_tree
- if self.trees == NULL:
- self.sampler = volume_render_sampler
- else:
- star_tree = self.tree_containers[pg.parent_grid_id]
- self.vra.star_list = self.trees[pg.parent_grid_id]
- self.vra.star_sigma_num = 2.0*star_tree.sigma**2.0
- self.vra.star_er = 2.326 * star_tree.sigma
- self.vra.star_coeff = star_tree.coeff
- self.sampler = volume_render_stars_sampler
-
- def __dealloc__(self):
- for i in range(self.vra.n_fits):
- free(self.vra.fits[i].d0)
- free(self.vra.fits[i].dy)
- free(self.vra.fits)
- free(self.vra)
-
-cdef class LightSourceRenderSampler(ImageSampler):
- cdef VolumeRenderAccumulator *vra
- cdef public object tf_obj
- cdef public object my_field_tables
- def __cinit__(self,
- np.ndarray vp_pos,
- np.ndarray vp_dir,
- np.ndarray[np.float64_t, ndim=1] center,
- bounds,
- np.ndarray[np.float64_t, ndim=3] image,
- np.ndarray[np.float64_t, ndim=1] x_vec,
- np.ndarray[np.float64_t, ndim=1] y_vec,
- np.ndarray[np.float64_t, ndim=1] width,
- tf_obj, n_samples = 10,
- light_dir=[1.,1.,1.],
- light_rgba=[1.,1.,1.,1.],
- **kwargs):
- ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
- x_vec, y_vec, width, **kwargs)
- cdef int i
- cdef np.ndarray[np.float64_t, ndim=1] temp
- # Now we handle tf_obj
- self.vra = <VolumeRenderAccumulator *> \
- malloc(sizeof(VolumeRenderAccumulator))
- self.vra.fits = <FieldInterpolationTable *> \
- malloc(sizeof(FieldInterpolationTable) * 6)
- self.vra.n_fits = tf_obj.n_field_tables
- assert(self.vra.n_fits <= 6)
- self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
- self.vra.n_samples = n_samples
- self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
- self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
- light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
- for i in range(3):
- self.vra.light_dir[i] = light_dir[i]
- for i in range(4):
- self.vra.light_rgba[i] = light_rgba[i]
- self.my_field_tables = []
- for i in range(self.vra.n_fits):
- temp = tf_obj.tables[i].y
- FIT_initialize_table(&self.vra.fits[i],
- temp.shape[0],
- <np.float64_t *> temp.data,
- tf_obj.tables[i].x_bounds[0],
- tf_obj.tables[i].x_bounds[1],
- tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
- tf_obj.weight_table_ids[i])
- self.my_field_tables.append((tf_obj.tables[i],
- tf_obj.tables[i].y))
- for i in range(6):
- self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
- self.supp_data = <void *> self.vra
-
- cdef void setup(self, PartitionedGrid pg):
- self.sampler = volume_render_gradient_sampler
-
- def __dealloc__(self):
- for i in range(self.vra.n_fits):
- free(self.vra.fits[i].d0)
- free(self.vra.fits[i].dy)
- free(self.vra.light_dir)
- free(self.vra.light_rgba)
- free(self.vra.fits)
- free(self.vra)
-
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef int walk_volume(VolumeContainer *vc,
np.float64_t v_pos[3],
np.float64_t v_dir[3],
- sampler_function *sampler,
+ sampler_function *sample,
void *data,
np.float64_t *return_t = NULL,
np.float64_t max_t = 1.0) nogil:
@@ -1031,7 +139,7 @@
else:
i = 2
exit_t = fmin(tmax[i], max_t)
- sampler(vc, v_pos, v_dir, enter_t, exit_t, cur_ind, data)
+ sample(vc, v_pos, v_dir, enter_t, exit_t, cur_ind, data)
cur_ind[i] += step[i]
enter_t = tmax[i]
tmax[i] += tdelta[i]
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/image_samplers.pxd
--- /dev/null
+++ b/yt/utilities/lib/image_samplers.pxd
@@ -0,0 +1,78 @@
+"""
+Definitions for image samplers
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+cimport kdtree_utils
+from .volume_container cimport VolumeContainer
+from .lenses cimport \
+ calculate_extent_function, \
+ generate_vector_info_function, \
+ ImageContainer
+from .partitioned_grid cimport PartitionedGrid
+
+DEF Nch = 4
+
+# NOTE: We don't want to import the field_interpolator_tables here, as it
+# breaks a bunch of C++ interop. Maybe some day it won't. So, we just forward
+# declare.
+cdef struct VolumeRenderAccumulator
+
+cdef struct ImageAccumulator:
+ np.float64_t rgba[Nch]
+ void *supp_data
+
+cdef class ImageSampler:
+ cdef ImageContainer *image
+ cdef public object acenter, aimage, ax_vec, ay_vec
+ cdef public object azbuffer
+ cdef public object aimage_used
+ cdef public object amesh_lines
+ cdef void *supp_data
+ cdef np.float64_t width[3]
+ cdef public object lens_type
+ cdef calculate_extent_function *extent_function
+ cdef generate_vector_info_function *vector_function
+ cdef void setup(self, PartitionedGrid pg)
+ @staticmethod
+ cdef void sample(VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil
+
+cdef class ProjectionSampler(ImageSampler):
+ pass
+
+cdef class InterpolatedProjectionSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
+
+cdef class VolumeRenderSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
+ cdef kdtree_utils.kdtree **trees
+ cdef object tree_containers
+
+cdef class LightSourceRenderSampler(ImageSampler):
+ cdef VolumeRenderAccumulator *vra
+ cdef public object tf_obj
+ cdef public object my_field_tables
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/image_samplers.pyx
--- /dev/null
+++ b/yt/utilities/lib/image_samplers.pyx
@@ -0,0 +1,486 @@
+"""
+Image sampler definitions
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from libc.stdlib cimport malloc, calloc, free, abs
+from libc.math cimport exp, floor, log2, \
+ 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,\
+ FIT_eval_transfer_with_light
+cimport lenses
+from .grid_traversal cimport walk_volume
+from .fixed_interpolator cimport \
+ offset_interpolate, \
+ fast_interpolate, \
+ trilinear_interpolate, \
+ eval_gradient, \
+ offset_fill, \
+ vertex_interp
+
+cdef extern from "platform_dep.h":
+ long int lrint(double x) nogil
+
+DEF Nch = 4
+
+from cython.parallel import prange, parallel, threadid
+from vec3_ops cimport dot, subtract, L2_norm, fma
+
+from cpython.exc cimport PyErr_CheckSignals
+
+cdef struct VolumeRenderAccumulator:
+ int n_fits
+ int n_samples
+ FieldInterpolationTable *fits
+ int field_table_ids[6]
+ np.float64_t star_coeff
+ np.float64_t star_er
+ np.float64_t star_sigma_num
+ kdtree_utils.kdtree *star_list
+ np.float64_t *light_dir
+ np.float64_t *light_rgba
+ int grey_opacity
+
+
+cdef class ImageSampler:
+ def __init__(self,
+ np.float64_t[:,:,:] vp_pos,
+ np.float64_t[:,:,:] vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ *args, **kwargs):
+ self.image = <ImageContainer *> calloc(sizeof(ImageContainer), 1)
+ cdef np.float64_t[:,:] zbuffer
+ cdef np.int64_t[:,:] image_used
+ cdef np.int64_t[:,:] mesh_lines
+ 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")
+
+ image_used = np.zeros((image.shape[0], image.shape[1]), "int64")
+ mesh_lines = np.zeros((image.shape[0], image.shape[1]), "int64")
+
+ self.lens_type = kwargs.pop("lens_type", None)
+ if self.lens_type == "plane-parallel":
+ self.extent_function = lenses.calculate_extent_plane_parallel
+ self.vector_function = lenses.generate_vector_info_plane_parallel
+ 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 += "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)
+
+ if camera_data is not None and self.lens_type == 'perspective':
+ self.extent_function = lenses.calculate_extent_perspective
+ else:
+ self.extent_function = lenses.calculate_extent_null
+ self.\
+ vector_function = lenses.generate_vector_info_null
+
+ # 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
+ # arguments.
+ self.image.vp_pos = vp_pos
+ self.image.vp_dir = vp_dir
+ self.image.image = self.aimage = image
+ self.acenter = center
+ self.image.center = <np.float64_t *> center.data
+ self.ax_vec = x_vec
+ self.image.x_vec = <np.float64_t *> x_vec.data
+ self.ay_vec = y_vec
+ self.image.y_vec = <np.float64_t *> y_vec.data
+ self.image.zbuffer = self.azbuffer = zbuffer
+ self.image.image_used = self.aimage_used = image_used
+ self.image.mesh_lines = self.amesh_lines = mesh_lines
+ self.image.nv[0] = image.shape[0]
+ self.image.nv[1] = image.shape[1]
+ for i in range(4): self.image.bounds[i] = bounds[i]
+ self.image.pdx = (bounds[1] - bounds[0])/self.image.nv[0]
+ self.image.pdy = (bounds[3] - bounds[2])/self.image.nv[1]
+ for i in range(3):
+ self.width[i] = width[i]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ def __call__(self, PartitionedGrid pg, int num_threads = 0):
+ # This routine will iterate over all of the vectors and cast each in
+ # turn. Might benefit from a more sophisticated intersection check,
+ # like http://courses.csusm.edu/cs697exz/ray_box.htm
+ cdef int vi, vj, hit, i, j
+ cdef np.int64_t iter[4]
+ cdef VolumeContainer *vc = pg.container
+ cdef ImageContainer *im = self.image
+ self.setup(pg)
+ cdef np.float64_t *v_pos
+ cdef np.float64_t *v_dir
+ cdef np.float64_t max_t
+ hit = 0
+ cdef np.int64_t nx, ny, size
+ self.extent_function(self.image, vc, iter)
+ iter[0] = i64clip(iter[0]-1, 0, im.nv[0])
+ iter[1] = i64clip(iter[1]+1, 0, im.nv[0])
+ iter[2] = i64clip(iter[2]-1, 0, im.nv[1])
+ iter[3] = i64clip(iter[3]+1, 0, im.nv[1])
+ nx = (iter[1] - iter[0])
+ ny = (iter[3] - iter[2])
+ size = nx * ny
+ cdef ImageAccumulator *idata
+ cdef np.float64_t width[3]
+ cdef int chunksize = 100
+ for i in range(3):
+ width[i] = self.width[i]
+ with nogil, parallel(num_threads = num_threads):
+ idata = <ImageAccumulator *> malloc(sizeof(ImageAccumulator))
+ idata.supp_data = self.supp_data
+ v_pos = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ v_dir = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ for j in prange(size, schedule="static", chunksize=chunksize):
+ vj = j % ny
+ vi = (j - vj) / ny + iter[0]
+ vj = vj + iter[2]
+ # Dynamically calculate the position
+ self.vector_function(im, vi, vj, width, v_dir, v_pos)
+ for i in range(Nch):
+ idata.rgba[i] = im.image[vi, vj, i]
+ max_t = fclip(im.zbuffer[vi, vj], 0.0, 1.0)
+ walk_volume(vc, v_pos, v_dir, self.sample,
+ (<void *> idata), NULL, max_t)
+ if (j % (10*chunksize)) == 0:
+ with gil:
+ PyErr_CheckSignals()
+ for i in range(Nch):
+ im.image[vi, vj, i] = idata.rgba[i]
+ free(idata)
+ free(v_pos)
+ free(v_dir)
+ return hit
+
+ cdef void setup(self, PartitionedGrid pg):
+ return
+
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ return
+
+ def ensure_code_unit_params(self, params):
+ for param_name in ['center', 'vp_pos', 'vp_dir', 'width']:
+ param = params[param_name]
+ if hasattr(param, 'in_units'):
+ params[param_name] = param.in_units('code_length')
+ bounds = params['bounds']
+ if hasattr(bounds[0], 'units'):
+ params['bounds'] = tuple(b.in_units('code_length').d for b in bounds)
+
+ return params
+
+ def __dealloc__(self):
+ self.image.image = None
+ self.image.vp_pos = None
+ self.image.vp_dir = None
+ self.image.zbuffer = None
+ self.image.camera_data = None
+ self.image.image_used = None
+ free(self.image)
+
+cdef class ProjectionSampler(ImageSampler):
+
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef int i
+ cdef np.float64_t dl = (exit_t - enter_t)
+ cdef int di = (index[0]*vc.dims[1]+index[1])*vc.dims[2]+index[2]
+ for i in range(imin(4, vc.n_fields)):
+ im.rgba[i] += vc.data[i][di] * dl
+
+
+cdef class InterpolatedProjectionSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ n_samples = 10, **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.n_samples = n_samples
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ for j in range(imin(3, vc.n_fields)):
+ im.rgba[j] += dvs[j] * dt
+ for j in range(3):
+ dp[j] += ds[j]
+
+
+cdef class VolumeRenderSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ tf_obj, n_samples = 10,
+ star_list = None, **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ cdef int i
+ cdef np.ndarray[np.float64_t, ndim=1] temp
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.fits = <FieldInterpolationTable *> \
+ malloc(sizeof(FieldInterpolationTable) * 6)
+ self.vra.n_fits = tf_obj.n_field_tables
+ assert(self.vra.n_fits <= 6)
+ self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
+ self.vra.n_samples = n_samples
+ self.my_field_tables = []
+ for i in range(self.vra.n_fits):
+ temp = tf_obj.tables[i].y
+ FIT_initialize_table(&self.vra.fits[i],
+ temp.shape[0],
+ <np.float64_t *> temp.data,
+ tf_obj.tables[i].x_bounds[0],
+ tf_obj.tables[i].x_bounds[1],
+ tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
+ tf_obj.weight_table_ids[i])
+ self.my_field_tables.append((tf_obj.tables[i],
+ tf_obj.tables[i].y))
+ for i in range(6):
+ self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef int cell_offset = index[0] * (vc.dims[1]) * (vc.dims[2]) \
+ + index[1] * (vc.dims[2]) + index[2]
+ if vc.mask[cell_offset] != 1:
+ return
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ FIT_eval_transfer(dt, dvs, im.rgba, vri.n_fits,
+ vri.fits, vri.field_table_ids, vri.grey_opacity)
+ for j in range(3):
+ dp[j] += ds[j]
+
+ def __dealloc__(self):
+ for i in range(self.vra.n_fits):
+ free(self.vra.fits[i].d0)
+ free(self.vra.fits[i].dy)
+ free(self.vra.fits)
+ free(self.vra)
+
+cdef class LightSourceRenderSampler(ImageSampler):
+ def __cinit__(self,
+ np.ndarray vp_pos,
+ np.ndarray vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec,
+ np.ndarray[np.float64_t, ndim=1] width,
+ tf_obj, n_samples = 10,
+ light_dir=[1.,1.,1.],
+ light_rgba=[1.,1.,1.,1.],
+ **kwargs):
+ ImageSampler.__init__(self, vp_pos, vp_dir, center, bounds, image,
+ x_vec, y_vec, width, **kwargs)
+ cdef int i
+ cdef np.ndarray[np.float64_t, ndim=1] temp
+ # Now we handle tf_obj
+ self.vra = <VolumeRenderAccumulator *> \
+ malloc(sizeof(VolumeRenderAccumulator))
+ self.vra.fits = <FieldInterpolationTable *> \
+ malloc(sizeof(FieldInterpolationTable) * 6)
+ self.vra.n_fits = tf_obj.n_field_tables
+ assert(self.vra.n_fits <= 6)
+ self.vra.grey_opacity = getattr(tf_obj, "grey_opacity", 0)
+ self.vra.n_samples = n_samples
+ self.vra.light_dir = <np.float64_t *> malloc(sizeof(np.float64_t) * 3)
+ self.vra.light_rgba = <np.float64_t *> malloc(sizeof(np.float64_t) * 4)
+ light_dir /= np.sqrt(light_dir[0]**2 + light_dir[1]**2 + light_dir[2]**2)
+ for i in range(3):
+ self.vra.light_dir[i] = light_dir[i]
+ for i in range(4):
+ self.vra.light_rgba[i] = light_rgba[i]
+ self.my_field_tables = []
+ for i in range(self.vra.n_fits):
+ temp = tf_obj.tables[i].y
+ FIT_initialize_table(&self.vra.fits[i],
+ temp.shape[0],
+ <np.float64_t *> temp.data,
+ tf_obj.tables[i].x_bounds[0],
+ tf_obj.tables[i].x_bounds[1],
+ tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
+ tf_obj.weight_table_ids[i])
+ self.my_field_tables.append((tf_obj.tables[i],
+ tf_obj.tables[i].y))
+ for i in range(6):
+ self.vra.field_table_ids[i] = tf_obj.field_table_ids[i]
+ self.supp_data = <void *> self.vra
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ @staticmethod
+ cdef void sample(
+ VolumeContainer *vc,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int index[3],
+ void *data) nogil:
+ cdef ImageAccumulator *im = <ImageAccumulator *> data
+ cdef VolumeRenderAccumulator *vri = <VolumeRenderAccumulator *> \
+ im.supp_data
+ # we assume this has vertex-centered data.
+ cdef int offset = index[0] * (vc.dims[1] + 1) * (vc.dims[2] + 1) \
+ + index[1] * (vc.dims[2] + 1) + index[2]
+ cdef np.float64_t dp[3]
+ cdef np.float64_t ds[3]
+ cdef np.float64_t dt = (exit_t - enter_t) / vri.n_samples
+ cdef np.float64_t dvs[6]
+ cdef np.float64_t *grad
+ grad = <np.float64_t *> malloc(3 * sizeof(np.float64_t))
+ for i in range(3):
+ dp[i] = (enter_t + 0.5 * dt) * v_dir[i] + v_pos[i]
+ dp[i] -= index[i] * vc.dds[i] + vc.left_edge[i]
+ dp[i] *= vc.idds[i]
+ ds[i] = v_dir[i] * vc.idds[i] * dt
+ for i in range(vri.n_samples):
+ for j in range(vc.n_fields):
+ dvs[j] = offset_interpolate(vc.dims, dp,
+ vc.data[j] + offset)
+ eval_gradient(vc.dims, dp, vc.data[0] + offset, grad)
+ FIT_eval_transfer_with_light(dt, dvs, grad,
+ vri.light_dir, vri.light_rgba,
+ im.rgba, vri.n_fits,
+ vri.fits, vri.field_table_ids, vri.grey_opacity)
+ for j in range(3):
+ dp[j] += ds[j]
+ free(grad)
+
+
+ def __dealloc__(self):
+ for i in range(self.vra.n_fits):
+ free(self.vra.fits[i].d0)
+ free(self.vra.fits[i].dy)
+ free(self.vra.light_dir)
+ free(self.vra.light_rgba)
+ free(self.vra.fits)
+ free(self.vra)
diff -r e7aabcd611651728be2dca5e817883f99c75a4a5 -r 6a4eea9e6d94b5a2020190a89cf503a6aafec306 yt/utilities/lib/lenses.pxd
--- /dev/null
+++ b/yt/utilities/lib/lenses.pxd
@@ -0,0 +1,60 @@
+"""
+Definitions for the lens code
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+cimport numpy as np
+cimport cython
+from .volume_container cimport VolumeContainer
+from vec3_ops cimport dot, subtract, L2_norm, fma
+from libc.math cimport exp, floor, log2, \
+ fabs, atan, atan2, asin, cos, sin, sqrt, acos, M_PI
+from yt.utilities.lib.fp_utils cimport imax, fmax, imin, fmin, iclip, fclip, i64clip
+
+cdef extern from "platform_dep.h":
+ long int lrint(double x) nogil
+
+cdef extern from "limits.h":
+ cdef int SHRT_MAX
+
+cdef struct ImageContainer:
+ np.float64_t[:,:,:] vp_pos
+ np.float64_t[:,:,:] vp_dir
+ np.float64_t *center
+ np.float64_t[:,:,:] image
+ np.float64_t[:,:] zbuffer
+ np.int64_t[:,:] image_used
+ np.int64_t[:,:] mesh_lines
+ 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
+
+
+ctypedef void calculate_extent_function(ImageContainer *image,
+ VolumeContainer *vc, np.int64_t rv[4]) nogil
+
+ctypedef void generate_vector_info_function(ImageContainer *im,
+ np.int64_t vi, np.int64_t vj,
+ np.float64_t width[2],
+ np.float64_t v_dir[3], np.float64_t v_pos[3]) nogil
+
+cdef generate_vector_info_function generate_vector_info_plane_parallel
+cdef generate_vector_info_function generate_vector_info_null
+cdef calculate_extent_function calculate_extent_plane_parallel
+cdef calculate_extent_function calculate_extent_perspective
+cdef calculate_extent_function calculate_extent_null
This diff is so big that we needed to truncate the remainder.
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