[yt-svn] commit/yt: 7 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sun Jun 15 05:45:57 PDT 2014
7 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/e43c9e8bb769/
Changeset: e43c9e8bb769
Branch: yt-3.0
User: Alex Bogert
Date: 2014-04-24 22:54:28
Summary: begin merging theia into yt-3.0
Affected #: 21 files
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
@@ -0,0 +1,94 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+#include <helpers/cuda.cu>
+
+texture<float , cudaTextureType3D, cudaReadModeElementType> volume;
+texture<float4, cudaTextureType1D, cudaReadModeElementType> transfer;
+
+
+extern "C"
+__global__ void front_to_back(uint *buffer, float *actual_matrix, int buffer_w,
+ int buffer_h, int max_steps, float density,
+ float offset, float scale,
+ float brightness, float step_size,
+ ) {
+
+ const float opacity_threshold = 0.95f;
+
+ float3x4 matrix = mat_array_to_3x4(actual_matrix);
+
+ int x = BLOCK_X;
+ int y = BLOCK_Y;
+
+ if ((x >= buffer_w) || (y >= buffer_h)) return;
+
+ float u = COORD_STD(x, buffer_w);
+ float v = COORD_STD(y, buffer_h);
+
+ // calculate eye ray in world space
+ Ray eye_ray = make_ray_from_eye(matrix, u, v);
+
+ // find intersection with box
+ float tnear, tfar;
+ int hit = intersect_std_box(eye_ray, &tnear, &tfar);
+
+ if (!hit) {
+ buffer[y*buffer_w + x] = 0;
+ return;
+ }
+
+ if (tnear < 0.0f) tnear = 0.0f; // clamp to near plane
+
+ // march along ray from front to back, accumulating color
+ float4 sum = make_float4(0.0f);
+ float t = tnear;
+ float3 pos = eye_ray.o + eye_ray.d*tnear;
+ float3 step = eye_ray.d*step_size;
+
+ // Cast Rays
+ float4 col = make_float4(0.0,0.0,0.0,0.0);
+ float top = (1.0/scale) + offset;
+ for (int i=0; i<=max_steps; i++) {
+ float sample = tex3D(volume, pos.x*0.5f+0.5f, pos.y*0.5f+0.5f, pos.z*0.5f+0.5f);
+
+
+ if (sample > offset) {
+ col = tex1D(transfer, (sample-offset)*scale);
+ }
+ else {
+ col = make_float4(0.0,0.0,0.0,0.0);
+ }
+
+ col.w *= density;
+
+ // pre-multiply alpha
+ col.x *= col.w;
+ col.y *= col.w;
+ col.z *= col.w;
+
+ // "over" operator for front-to-back blending
+ sum = sum + col*(1.0f - sum.w);
+
+ // exit early if opaque
+ if (sum.w > opacity_threshold)
+ break;
+
+ // Increment step size and position
+ t += step_size;
+ pos += step;
+
+ // If ray has cast too far, end
+ if (t > tfar) break;
+ }
+
+ sum *= brightness;
+
+ buffer[SCREEN_XY(x,y,buffer_w)] = rgbaFloatToInt(sum);
+}
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -0,0 +1,167 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+import pyRGBA.theia.helpers.cuda as cdh
+import numpy as np
+
+from pyRGBA.theia.transfer.linear_transfer import LinearTransferFunction
+
+from pyRGBA.theia.volumes.cube import Unigrid
+from pyRGBA.theia.surfaces.array_surface import ArraySurface
+import os
+
+class FrontToBackRaycaster:
+ r"""
+ This is a basic camera intended to be a base class.
+ The camera provies utilities for controling the view point onto the data.
+
+
+ Parameters
+ ----------
+ pos : numpy array
+
+ Example:
+
+ from pyRGBA.cameras.camera import Camera
+
+ cam = Camera()
+
+ cam.zoom(10.0)
+ cam.rotateX(np.pi)
+ cam.rotateY(np.pi)
+ cam.rotateZ(np.pi)
+ cam.translateX(-0.5)
+ cam.translateY(-0.5)
+ cam.translateZ(-0.5)
+
+ view = cam.get_matrix()
+
+ """
+
+ def __init__(self, matrix = None, surface = None,
+ volume = None, transfer = None,
+ size = (864, 480)) :
+
+ # kernel module
+ self.kernel_module = cdh.source_file(
+ self.base_directory(__file__, "front_to_back.cu"))
+
+ #kernel function definitions
+ self.cast_rays_identifier = self.kernel_module.get_function("front_to_back")
+ self.transfer_identifier = self.kernel_module.get_texref("transfer")
+ self.volume_identifier = self.kernel_module.get_texref("volume")
+
+ self.set_matrix(matrix)
+
+ if (surface == None) :
+ self.set_surface(ArraySurface(size = size))
+ else :
+ self.set_surface(surface)
+
+
+ self.volume = None
+ self.set_transfer(transfer)
+
+ self.set_sample_size()
+ self.set_max_samples()
+ self.set_density_scale()
+ self.set_brightness()
+
+ def __call__(self):
+ self.cast()
+
+ def get_surface(self) :
+ return self.surface.get_array()
+
+ def get_sample_size(self) :
+ return self.sample_size
+
+ def get_max_samples(self) :
+ return self.max_samples
+
+ def get_density_scale(self) :
+ return self.density_scale
+
+ def get_brightness(self):
+ return self.brightness
+
+ def set_sample_size(self, size = 0.01) :
+ self.sample_size = size
+
+ def set_max_samples(self, max = 5000) :
+ self.max_samples = max
+
+ def set_density_scale(self, scale = 0.05) :
+ self.density_scale = scale
+
+ def set_brightness(self, brightness = 1.0):
+ self.brightness = brightness
+
+ def cast(self):
+ w, h = self.surface.bounds
+
+ self.cast_rays_identifier(np.uint64(self.surface.device_ptr),
+ drv.In(self.matrix),
+ np.int32(w), np.int32(h),
+ np.int32(self.max_samples),
+ np.float32(self.density_scale),
+ np.float32(self.transfer.offset),
+ np.float32(self.transfer.scale),
+ np.float32(self.brightness),
+ np.float32(self.sample_size),
+ block=self.block,
+ grid=self.grid
+ )
+
+
+ def set_matrix(self, matrix):
+ self.matrix = matrix
+
+ def set_surface(self, surface = None, block_size = 32):
+ if surface == None:
+ self.surface = None
+ return
+
+ self.surface = surface
+ self.grid = cdh.block_estimate(block_size, self.surface.bounds)
+ self.block = (block_size, block_size, 1)
+
+ def send_volume_to_gpu(self, volume = None) :
+ if (volume != None) :
+ self.set_volume(Unigrid(array = volume, allocate = True))
+
+ def set_volume(self, volume):
+ if volume == None:
+ self.volume = None
+ return
+
+ self.volume = volume
+ self.volume_identifier.set_flags(self.volume.flags)
+ self.volume_identifier.set_filter_mode(self.volume.filter_mode)
+ self.volume_identifier.set_address_mode(0, self.volume.address_mode)
+ self.volume_identifier.set_address_mode(1, self.volume.address_mode)
+ self.volume_identifier.set_array(self.volume.cuda_volume_array)
+
+ def set_transfer(self, transfer):
+ if transfer == None:
+ self.transfer = None
+ return
+
+ self.transfer = transfer
+ self.transfer_identifier.set_flags(self.transfer.flags)
+ self.transfer_identifier.set_filter_mode(self.transfer.filter_mode)
+ self.transfer_identifier.set_address_mode(0, self.transfer.address_mode)
+ self.transfer_identifier.set_array(self.transfer.cuda_transfer_array)
+
+ def base_directory(self, dir, file):
+ base = os.path.dirname(dir)
+ src = os.path.join(base, file)
+ return src
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/camera.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/camera.py
@@ -0,0 +1,139 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+#Matix helpers for camera control
+import pyRGBA.theia.helpers.matrix_transformation as mth
+
+import numpy as np
+
+class Camera:
+ r"""
+ This is a basic camera intended to be a base class.
+ The camera provies utilities for controling the view point onto the data.
+
+
+ Parameters
+ ----------
+ pos : numpy array
+ The 3d position vector of the camera
+ rot : numpy array
+ The 3d rotation vector of the camera
+ scale : float
+ The scalar acting on camera position to zoom
+
+ Example:
+
+ from pyRGBA.cameras.camera import Camera
+
+ cam = Camera()
+
+ cam.zoom(10.0)
+ cam.rotateX(np.pi)
+ cam.rotateY(np.pi)
+ cam.rotateZ(np.pi)
+ cam.translateX(-0.5)
+ cam.translateY(-0.5)
+ cam.translateZ(-0.5)
+
+ view = cam.get_matrix()
+
+ """
+ def __init__(self, rot = np.array([0.0, 0.0, 0.0]), scale = 1.0,
+ pos = np.array([0.0, 0.0, 0.0]), ):
+ #Values to control camera
+ self.rot = rot
+ self.scale = scale
+ self.pos = pos
+
+ #TODO : users should have control over perspective frustrum
+ self.stdmatrix = mth.perspective( [0.0, 0.0, 0.25] )
+
+
+ def zoom(self, percentage = 10.0):
+ r"""This will zoom the camera by percentage specified.
+
+ Parameters
+ ----------
+ percentage : float
+ If percentage is postive the camera will zoom in, if negative
+ then the camera will zoom out. Cannot exceed 100%.
+ """
+ if (percentage > 100.0) :
+ percentage = 100.0
+ self.scale = ((100.0 - percentage)/100.0)*self.scale
+
+ def rotateX(self, degrees = 0.1):
+ r"""This will rotate the camera around its X axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[0] += degrees
+
+ def rotateY(self, degrees = 0.1):
+ r"""This will rotate the camera around its Y axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[1] += degrees
+
+ def rotateZ(self, degrees = 0.1):
+ r"""This will rotate the camera around its Z axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[2] += degrees
+
+ def translateX(self, dist = 0.1):
+ r"""This will move the camera's position along its X axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[0] += dist
+
+ def translateY(self, dist = 0.1):
+ r"""This will move the camera's position along its Y axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[1] += dist
+
+ def translateZ(self, dist = 0.1):
+ r"""This will move the camera's position along its Z axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[2] += dist
+
+ def get_matrix(self):
+ r"""This will calculate the curent modelview matrix
+ from the camera's position, rotation, and zoom scale.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ return mth.rotate_x(self.rot[0]) * mth.rotate_y(self.rot[1]) * mth.rotate_z(self.rot[2]) * mth.scale((self.scale, self.scale, self.scale)) * mth.perspective([0.0, 0.0, 0.25]) *mth.translate((self.pos[0], self.pos[1], self.pos[2]))
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/helpers/cuda.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.cu
@@ -0,0 +1,189 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+
+#include <helper_cuda.h>
+#include <helper_math.h>
+
+#define BLOCK_X (blockIdx.x*blockDim.x + threadIdx.x)
+#define BLOCK_Y (blockIdx.y*blockDim.y + threadIdx.y)
+#define SCREEN_XY(a, b, w) (b*w+a)
+#define COORD_STD(x, w) ((x / (float) w)*2.0f-1.0f)
+
+#define COLOR_BLACK make_float3(0.0f, 0.0f, 0.0f)
+
+typedef struct {
+ float4 m[3];
+} float3x4;
+
+typedef struct {
+ float4 m[4];
+} float4x4;
+
+
+struct Ray {
+ float3 o; // origin
+ float3 d; // direction
+};
+
+
+
+__device__ float4 mul(const float3x4 &M, const float4 &v);
+__device__ float3 mul(const float3x4 &M, const float3 &v);
+__device__ float3x4 mat_array_to_3x4(float *m);
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v);
+__device__ int intersect_box(Ray r, float3 boxmin, float3 boxmax,
+ float *tnear, float *tfar);
+__device__ int intersect_std_box(Ray r, float *tnear, float *tfar);
+
+__device__ float mag(const float3 &v);
+__device__ float avg(const float3 &v);
+__device__ uint rgbaFloatToInt(float4 rgba);
+__device__ uint rgbFloatToInt(float3 rgb);
+__device__ float4 intToRGBAFloat(uint irgba);
+__device__ float3 intToRGBFloat(uint irgba);
+
+
+
+__device__ uint rgbFloatToInt(float3 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ return (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ uint rgbaFloatToInt(float4 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ rgba.w = __saturatef(rgba.w);
+ return (uint(rgba.w*255)<<24) | (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ float4 intToRGBAFloat(uint irgba)
+{
+ float4 rgba;
+ rgba.w = (irgba >> 24 & 0xFF)/255;
+ rgba.z = (irgba >> 16 & 0xFF)/255;
+ rgba.y = (irgba >> 8 & 0xFF)/255;
+ rgba.x = (irgba & 0xFF)/255;
+ return rgba;
+}
+
+__device__ float3 intToRGBFloat(uint irgba)
+{
+ float3 rgb;
+ rgb.z = (irgba >> 16 & 0xFF)/255;
+ rgb.y = (irgba >> 8 & 0xFF)/255;
+ rgb.x = (irgba & 0xFF)/255;
+ return rgb;
+}
+
+
+__device__ float3x4 mat_array_to_3x4(float *m) {
+ float3x4 matrix;
+ matrix.m[0] = make_float4(m[0], m[1], m[2], m[3]);
+ matrix.m[1] = make_float4(m[4], m[5], m[6], m[7]);
+ matrix.m[2] = make_float4(m[8], m[9], m[10], m[11]);
+ return matrix;
+}
+
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v) {
+ Ray eyeRay;
+ eyeRay.o = make_float3(mul(m, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
+ eyeRay.d = normalize(make_float3(u, v, -2.0f));
+ eyeRay.d = mul(m, eyeRay.d);
+ return eyeRay;
+}
+
+__device__
+int intersect_box(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
+{
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+int intersect_std_box(Ray r, float *tnear, float *tfar)
+{
+ const float3 boxmin = make_float3(-1.0f, -1.0f, -1.0f);
+ const float3 boxmax = make_float3( 1.0f, 1.0f, 1.0f);
+
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+float4 mul(const float3x4 &M, const float4 &v)
+{
+ float4 r;
+ r.x = dot(v, M.m[0]);
+ r.y = dot(v, M.m[1]);
+ r.z = dot(v, M.m[2]);
+ r.w = 1.0f;
+ return r;
+}
+
+// transform vector by matrix with translation
+
+__device__
+float3 mul(const float3x4 &M, const float3 &v)
+{
+ float3 r;
+ r.x = dot(v, make_float3(M.m[0]));
+ r.y = dot(v, make_float3(M.m[1]));
+ r.z = dot(v, make_float3(M.m[2]));
+ return r;
+}
+
+
+__device__
+float mag(const float3 &v) {
+ return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+}
+
+
+__device__
+float avg(const float3 &v) {
+ return (v.x + v.y + v.z) / 3.0;
+}
+
+
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/helpers/cuda.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.py
@@ -0,0 +1,58 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+from pycuda.compiler import SourceModule
+import pycuda.gpuarray as gpuarray
+import pycuda.autoinit
+import math
+import os
+
+def sample_path():
+ #return '/pfs/sw/cuda-5.5/samples/common/inc'
+ return os.environ.get('PYRGBA_CUDA_SAMPLES')
+
+def cuda_helpers_path():
+ return os.path.dirname(__file__) + "/../"
+
+def source_file(path):
+ return SourceModule(open(path).read(),
+ include_dirs=[sample_path(),cuda_helpers_path()],
+ no_extern_c=True, keep=True)
+
+def block_estimate(thread_size, shape):
+ (x, y) = shape
+ return (int(math.ceil(float(x) / thread_size)), int(math.ceil(float(y) / thread_size)), 1)
+
+def np3d_to_device_array(np_array, allow_surface_bind=True):
+ d, h, w = np_array.shape
+
+ descr = drv.ArrayDescriptor3D()
+ descr.width = w
+ descr.height = h
+ descr.depth = d
+ descr.format = drv.dtype_to_array_format(np_array.dtype)
+ descr.num_channels = 1
+ descr.flags = 0
+
+ if allow_surface_bind:
+ descr.flags = drv.array3d_flags.SURFACE_LDST
+
+ device_array = drv.Array(descr)
+
+ copy = drv.Memcpy3D()
+ copy.set_src_host(np_array)
+ copy.set_dst_array(device_array)
+ copy.width_in_bytes = copy.src_pitch = np_array.strides[1]
+ copy.src_height = copy.height = h
+ copy.depth = d
+
+ copy()
+
+ return device_array
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/helpers/cuda/matrix_helper.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda/matrix_helper.cu
@@ -0,0 +1,265 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+
+#include <helper_cuda.h>
+#include <helper_math.h>
+
+
+typedef float3 vec3;
+typedef float4 vec4;
+typedef float2 vec2;
+
+
+typedef struct {
+ vec4 m[4];
+} mat4x4;
+
+typedef struct {
+ vec3 m[3];
+} mat3x3;
+
+typedef struct {
+ vec2 m[2];
+} mat2x2;
+
+typedef struct {
+ vec3 m[4];
+} mat3x4;
+
+typedef struct {
+ vec4 m[3];
+} mat4x3;
+
+typedef struct {
+ vec4 m[2];
+} mat4x2;
+
+typedef struct {
+ vec2 m[4];
+} mat2x4;
+
+typedef mat4 mat4x4;
+typedef mat3 mat3x3;
+typedef mat2 mat2x2;
+
+
+
+// Column Major Always
+// All vectors are k x 1, and cannot be at the front of a multiplication
+
+
+// -----------------------------------------------------------------------------
+// Create Matrices, identity by default
+// -----------------------------------------------------------------------------
+__device__ __host__ mat4 make_mat4() {
+ mat4 tmp;
+ tmp.m[0] = make_float4(1.0f, 0.0f, 0.0f, 0.0f);
+ tmp.m[1] = make_float4(0.0f, 1.0f, 0.0f, 0.0f);
+ tmp.m[2] = make_float4(0.0f, 0.0f, 1.0f, 0.0f);
+ tmp.m[3] = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
+ return tmp;
+}
+
+__device__ __host__ mat3 make_mat3() {
+ mat3 tmp;
+ tmp.m[0] = make_float3(1.0f, 0.0f, 0.0f);
+ tmp.m[1] = make_float3(0.0f, 1.0f, 0.0f);
+ tmp.m[2] = make_float3(0.0f, 0.0f, 1.0f);
+ return tmp;
+}
+
+__device__ __host__ mat2 make_mat2() {
+ mat3 tmp;
+ tmp.m[0] = make_float2(1.0f, 0.0f);
+ tmp.m[1] = make_float2(0.0f, 1.0f);
+ return tmp;
+}
+
+
+// -----------------------------------------------------------------------------
+// Get Transpose vectors, get column vector
+// -----------------------------------------------------------------------------
+// 4s
+__device__ __host__ vec4 mvec(mat4x4 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ else if(col == 2) tmp = make_float4(a.m[0].z, a.m[1].z, a.m[2].z, a.m[3].z);
+ else if(col == 3) tmp = make_float4(a.m[0].w, a.m[1].w, a.m[2].w, a.m[3].w);
+ return tmp;
+}
+
+__device__ __host__ vec3 mvec(mat4x3 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ else if(col == 2) tmp = make_float3(a.m[0].z, a.m[1].z, a.m[2].z);
+ else if(col == 3) tmp = make_float3(a.m[0].w, a.m[1].w, a.m[2].w);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat4x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ else if(col == 2) tmp = make_float2(a.m[0].z, a.m[1].z);
+ else if(col == 3) tmp = make_float2(a.m[0].w, a.m[1].w);
+ return tmp;
+}
+// 3s
+
+__device__ __host__ vec3 mvec(mat3x4 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ else if(col == 2) tmp = make_float4(a.m[0].z, a.m[1].z, a.m[2].z, a.m[3].z);
+ return tmp;
+}
+
+__device__ __host__ vec4 mvec(mat3x3 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ else if(col == 2) tmp = make_float3(a.m[0].z, a.m[1].z, a.m[2].z);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat3x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ else if(col == 2) tmp = make_float2(a.m[0].z, a.m[1].z);
+ return tmp;
+}
+// 2s
+
+__device__ __host__ vec2 mvec(mat2x4 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ return tmp;
+}
+
+__device__ __host__ vec3 mvec(mat2x3 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat2x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ return tmp;
+}
+
+
+// -----------------------------------------------------------------------------
+// Transpose Matrix
+// Cannot transpose anything other than square matrices
+// -----------------------------------------------------------------------------
+__device__ __host__ mat4 transpose(mat4 a) {
+ mat4 tmp;
+ for(int i = 0; i < 4; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+__device__ __host__ mat3 transpose(mat3 a) {
+ mat3 tmp;
+ for(int i = 0; i < 3; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+__device__ __host__ mat2 transpose(mat2 a) {
+ mat2 tmp;
+ for(int i = 0; i < 2; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+// -----------------------------------------------------------------------------
+// Dot Products
+// -----------------------------------------------------------------------------
+inline __device__ __host__ float dot(vec4 a, vec4 b) {
+ return (a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w);
+}
+
+inline __device__ __host__ float dot(vec3 a, vec3 b) {
+ return (a.x*b.x + a.y*b.y + a.z*b.z);
+}
+
+inline __device__ __host__ float dot(vec2 a, vec2 b) {
+ return (a.x*b.x + a.y*b.y);
+}
+
+// -----------------------------------------------------------------------------
+// Scale Matrix
+// -----------------------------------------------------------------------------
+
+
+// -----------------------------------------------------------------------------
+// Multiply Matrix
+// -----------------------------------------------------------------------------
+//4s
+__device__ __host__ mat4x4 operator*(mat4x4 a, mat4x4 b) {
+ mat4x4 tmp;
+ for(int i = 0; i < 4; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ mat4x3 operator*(mat4x4 a, mat4x3 b) {
+ mat4x3 tmp;
+ for(int i = 0; i < 3; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ mat4x2 operator*(mat4x4 a, mat4x2 b) {
+ mat4x2 tmp;
+ for(int i = 0; i < 2; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ vec4 operator*(mat4x4 a, vec4 b) {
+ vec4 tmp;
+ tmp.x = dot(a.m[0], b);
+ tmp.y = dot(a.m[1], b);
+ tmp.z = dot(a.m[2], b);
+ tmp.w = dot(a.m[3], b);
+ return tmp;
+}
+
+__device__ __host__ mat3x4 operator*(mat3x3 a, mat3x4 b);
+__device__ __host__ mat3x3 operator*(mat3x3 a, mat3x3 b);
+__device__ __host__ mat3x2 operator*(mat3x3 a, mat3x2 b);
+__device__ __host__ vec3 operator*(mat3x3 a, vec3 b);
+
+__device__ __host__ mat2x4 operator*(mat2x2 a, mat2x4 b);
+__device__ __host__ mat2x3 operator*(mat2x2 a, mat2x3 b);
+__device__ __host__ mat2x2 operator*(mat2x2 a, mat2x2 b);
+__device__ __host__ vec2 operator*(mat2x2 a, vec2 b);
+
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/helpers/matrix_transformation.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/matrix_transformation.py
@@ -0,0 +1,80 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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
+from math import sin, cos
+
+def identity():
+ return np.matrix([[1.0, 0.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0, 0.0],
+ [0.0, 0.0, 1.0, 0.0],
+ [0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+def standard():
+ return np.matrix([[1.0, 0.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0, 0.0],
+ [0.0, 0.0, 1.0, 4.0],
+ [0.0, 0.0, -1.0, 0.0]]).astype(np.float32)
+
+def perspective(eye):
+ [ex, ey, ez] = eye
+ a = -(ex/ez)
+ b = -(ey/ez)
+ c = 1.0/ez
+ return np.matrix([[1.0, 0.0, 0.0, 0.0],
+ [0.0, 1.0, 0.0, 0.0],
+ [ a, b, 1.0, c],
+ [0.0, 0.0, 0.0, 0.0]]).astype(np.float32)
+
+def rotate(k, vec):
+ [l, m, n] = vec
+ return np.matrix([[l*l*(1.0-cos(k))+cos(k),
+ l*m*(1.0-cos(k))+n*sin(k),
+ l*n*(1.0-cos(k))-m*sin(k),
+ 0.0],
+ [m*l*(1.0-cos(k))-n*sin(k),
+ m*m*(1.0-cos(k))+cos(k),
+ m*n*(1.0-cos(k))+l*sin(k),
+ 0.0],
+ [n*l*(1.0-cos(k))+m*sin(k),
+ n*m*(1.0-cos(k))-l*sin(k),
+ n*n*(1.0-cos(k))+cos(k),
+ 0.0],
+ [0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+
+def rotate_x(k):
+ return np.matrix([[1.0, 0.0, 0.0, 0.0],
+ [0.0, cos(k), -sin(k), 0.0],
+ [0.0, sin(k), cos(k), 0.0],
+ [0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+
+def rotate_y(k):
+ return np.matrix([[ cos(k), 0.0, sin(k), 0.0],
+ [ 0.0, 1.0, 0.0, 0.0],
+ [-sin(k), 0.0, cos(k), 0.0],
+ [ 0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+
+def rotate_z(k):
+ return np.matrix([[cos(k), -sin(k), 0.0, 0.0],
+ [sin(k), cos(k), 0.0, 0.0],
+ [ 0.0, 0.0, 1.0, 0.0],
+ [ 0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+
+def scale(vec):
+ [x, y, z] = vec
+ return np.matrix([[x, 0.0, 0.0, 0.0],
+ [0.0, y, 0.0, 0.0],
+ [0.0, 0.0, z, 0.0],
+ [0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
+
+def translate(vec):
+ [x, y, z] = vec
+ return np.matrix([[1.0, 0.0, 0.0, x],
+ [0.0, 1.0, 0.0, y],
+ [0.0, 0.0, 1.0, z],
+ [0.0, 0.0, 0.0, 1.0]]).astype(np.float32)
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/scene.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/scene.py
@@ -0,0 +1,95 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+from pyRGBA.theia.source import TheiaSource
+from pyRGBA.theia.camera import Camera
+
+class TheiaScene:
+ r"""
+ TheiaScene
+
+ This is the highest level entry point into Theia.
+ This class is a container for a Theia source and
+ a camera for controlling a view onto the source.
+
+ Call :
+ update()
+ To perform the volume rendering.
+ Call :
+ get_results()
+ To get the resulting image array.
+
+
+ Parameters
+ ----------
+ ts : TheiaSource
+ The container for the raycasting algorithm and the volumetric data
+ volume :
+ If a TheiaSource is not specified a volme can be used to initialize
+ the scene. However, rendering will only occur when a raycaster
+ is also provided.
+
+ camera :
+ A camera for controling view point. Must be able to provide a 4x4
+ matrix to the raycaster via get_matrix()
+ If no camera is provided the default Theia Camera is used.
+
+ Example:
+
+ from pyRGBA.theia.scene import TheiaScene
+ from pyRGBA.algorithms.front_to_back import FrontToBackRaycaster
+ import numpy as np
+
+ #load 3D numpy array of float32
+ volume = np.load("/home/bogert/log_densities_1024.npy")
+
+ scene = TheiaScene( volume = volume, raycaster = FrontToBackRaycaster() )
+
+ scene.camera.rotateX(1.0)
+ scene.update()
+
+ surface = scene.get_results()
+ #now do something with surface
+
+ """
+ def __init__(self, ts = None, camera = None, volume = None, raycaster = None) :
+ #Values to control screen resolution
+ if ts != None :
+ self.source = ts
+ else :
+ self.source = TheiaSource(volume = volume, raycaster = raycaster)
+
+ if camera == None :
+ self.camera = Camera()
+ else :
+ self.camera = camera
+
+ def get_raycaster(self):
+ """
+ This is for programmers who would prefer calling :
+ scene.get_raycaster()
+ instead of :
+ scene.source.raycaster
+ """
+ return self.source.raycaster
+
+ def update(self):
+ """
+ This will send the camera's modelview matrix to
+ the raycaster and call the raycaster's cast function
+ to perform the volume rendering.
+ """
+ self.source.update(self.camera.get_matrix())
+
+ def get_results(self):
+ """
+ This will cause the results of the volume rendering to be
+ copied back from the GPU and then returned as a 2d numpy array.
+ """
+ return self.source.get_results()
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/source.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/source.py
@@ -0,0 +1,129 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+class TheiaSource:
+ r"""
+ TheiaSource
+
+ This is a wrapper class to create an interface for consistent access to
+ raycasters and the volume they act on, as well as the returned results.
+
+
+ The raycaster shold be a GPU based raycasting algorithm to act on
+ the volume. The raycaster class must provide the following functions:
+ send_volume_to_gpu(volume)
+ Accepts a volume and moves the data onto the graphics card
+ memory.
+ set_matrix(camera.get_matrix())
+ Accepts a modelview matrix (4x4 numpy matrix) from the
+ Theia Camera class and sets the view on the volume
+ cast()
+ This must cause the ray caster to act on the volume
+ in the graphics cards memory, sampling the volume and
+ filling the surface with results
+
+ get_surface()
+ This must return a 2d numpy array with the results of
+ the last raycast
+
+
+ Any UI's using a TheiaSource can call the functions:
+
+ update()
+ This will call the raycaster's cast function
+
+ get_results()
+ This returns a 2d numpy array with the results of the latest update
+
+ Parameters
+ ----------
+ volume :
+ The data the camera will volume render.
+ raycaster :
+ A ray casting algorithm conforming to the requirements
+ detailed above
+
+ Example:
+
+ from pyRGBA.theia.source import TheiaSource
+ from pyRGBA.algorithms.front_to_back import FrontToBackRaycaster
+ import numpy as np
+
+ #load a 3d numpy array of float32
+ volume = np.load("/home/bogert/log_densities_1024.npy")
+ raycaster = FrontToBackRaycaster()
+
+ ts = TheiaSource(volume = volume, raycaster = raycaster)
+
+ ts.update()
+
+ image_array = ts.get_results()
+
+ """
+ def __init__(self, volume = None, raycaster = None):
+ #Values to control screen resolution
+ self.raycaster = None
+ self.volume = None
+
+ self.set_volume(volume)
+ self.set_raycaster(raycaster)
+
+ def set_volume(self, volume = None):
+ r"""This will set the TheiaSource volume to
+ the volume supplied. If the TheiaSource
+ has a raycaster object it will send the
+ volume to the gpu.
+
+ Parameters
+ ----------
+ volume :
+ The volumetric data to send to the raycaster
+ """
+ self.volume = volume
+ if (self.volume != None and self.raycaster != None) :
+ self.raycaster.send_volume_to_gpu(self.volume)
+
+ def set_raycaster(self, raycaster = None):
+ r"""This will set the TheiaSource raycaster to
+ the raycaster supplied. If the TheiaSource
+ has a volume set previously it will send the
+ volume to the gpu via the raycaster.
+
+ Parameters
+ ----------
+ raycaster :
+ The raycaster used to render the volume.
+ """
+ self.raycaster = raycaster
+ if (self.volume != None and self.raycaster != None) :
+ self.raycaster.send_volume_to_gpu(self.volume)
+
+ def update(self, matrix = None):
+ r"""This will update the raycaster with any
+ matrix information if given and then
+ call the raycaster's cast() function to
+ render the volume.
+
+ Parameters
+ ----------
+ matrix : 4x4 numpy matrix
+ The modelview matrix used to define the
+ view onto the volume
+ """
+ if (matrix != None) :
+ self.raycaster.set_matrix(matrix)
+ self.raycaster.cast()
+
+ def get_results(self):
+ r"""This will return the results of the
+ raycasting after the update function
+ is called as a 2d numpy array.
+
+ """
+ return self.raycaster.get_surface()
+
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/surfaces/array_surface.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/surfaces/array_surface.py
@@ -0,0 +1,112 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+
+import numpy as np
+
+from pyRGBA.theia.surfaces.surface import Surface
+
+
+"""
+An array surface is a type of surface object that is bound to a specified
+n-dimensional numpy array. It contains functions that allow data to be
+readily copied back and forth between host and device.
+"""
+
+class ArraySurface(Surface):
+ def __init__(self, array=None, size=None, copy=True) :
+ """
+ Creates a new ArraySurface. May optionally be given initial values.
+
+ Parameters
+ ----------
+ array: An ndarray to initialize local_array to. Will copy to
+ device by default.
+
+ size: A shape that would serve as a basis for an ndarray.
+ this argument is not necessary if array argument is specified
+ as the shape can be inferred from it.
+
+ copy: When True, will copy local array to device.
+ """
+ Surface.__init__(self)
+
+ self.local_array = None
+ self.bounds = size
+
+ if (array != None):
+ self.set_array(array, copy=copy)
+ else:
+ if size != None:
+ self.set_array(np.zeros(size, dtype=np.uint32), copy=copy)
+
+
+ def device_to_local(self):
+ """
+ Copies device memory to local host memory.
+ """
+
+ self.local_array = drv.from_device_like(self.device_ptr, self.local_array)
+
+
+ def local_to_device(self, early_free=True):
+ """
+ Copies locally set array data to device memory.
+ Parameters
+ ----------
+ early_free: boolean
+ When True, will automatically free all device memory
+ before allocation, preventing double-allocation for a short time.
+ """
+
+ if self.device_ptr != None and early_free == True:
+ self.device_ptr.free()
+
+ self.device_ptr = drv.to_device(self.local_array)
+
+
+ def set_array(self, array, copy=True):
+ """
+ Sets the local array in the Surface to the first argument, array, a
+ a numpy array.
+ Parameters
+ ----------
+ array : 2D numpy array
+ Container defining the surface.
+ copy: boolean
+ When True, will immediately copy local data to device.
+ """
+ self.local_array = drv.pagelocked_empty_like(array)
+ self.bounds = array.shape
+ if copy == True:
+ self.local_to_device()
+
+
+ def get_array(self, copy=True):
+ """
+ Gets a copy of the local array and returns it.
+
+ Parameters
+ ----------
+
+ copy: boolean
+ When True, will immediately copy device data
+ into local memory.
+ """
+
+ if copy == True:
+ self.device_to_local()
+ return self.local_array
+
+
+ def get_bounds(self):
+ """
+ Returns the shape of the local ndarray.
+ """
+ return self.bounds
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/surfaces/surface.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/surfaces/surface.py
@@ -0,0 +1,24 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+
+class Surface:
+ """
+ Containe for holding image array, results of volume rendering.
+
+ Parameters
+ ----------
+ device_ptr: uint64
+ A pycuda allocation object representing a point to device memory.
+ bounds:
+ A tuple representing the shape of the underlying data, often
+ numpy.shape.
+ """
+ def __init__(self):
+ self.device_ptr = None
+ self.bounds = None
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/transfer/helper.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/transfer/helper.py
@@ -0,0 +1,56 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+
+"""Transfer Function Helper Module
+This module supports transfer functions that are defined in this directory. General
+functions should be placed here that might be utilized by multiple transfer functions.
+"""
+
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+
+
+def yt_to_rgba_transfer(yt_tf):
+ """Takes a YT Transfer function (ColorTransferFunction) and forms a NumPy Array.
+ This is commonly used to hack together a contiguous array from a Yt
+ transfer function.
+ """
+ return (yt_tf.red.y, yt_tf.green.y, yt_tf.blue.y, yt_tf.alpha.y)
+
+def bounds_to_scale(bounds):
+ """Defines a transfer function scale factor given a transfer function.
+ """
+ mi, ma = bounds
+
+ return (1.0)/(ma - mi)
+
+def bounds_to_offset(bounds):
+ """Defines a transfer function offset given a transfer function.
+ """
+ mi, ma = bounds
+
+ return mi
+
+def make_yt_transfer(bounds = (0.0, 1.0), colormap = "algae", bins = 1000, scale = 1.0, scale_func = None):
+ """Defines a transfer function offset given a transfer function.
+ """
+ mi, ma = bounds
+ transfer = ColorTransferFunction( (mi, ma), bins)
+ if scale_func == None :
+ transfer.map_to_colormap(mi, ma, colormap=colormap, scale=1.0)
+ else :
+ transfer.map_to_colormap(mi, ma, colormap=colormap, scale_func = scale_func)
+
+
+ return transfer
+
+def base_directory(file):
+
+ dir = os.path.dirname(file)
+ src = os.path.join(dir, src)
+ return src
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
@@ -0,0 +1,88 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+
+"""Linear Transfer Function
+This module defines a linear-type transfer function. This type
+of transfer function is constructed from a contiguous array of
+RGBA values. In this case, the array input should represent a
+numpy array.
+"""
+
+import pycuda.driver as drv
+import pyRGBA.theia.transfer.helper as tfh
+
+import numpy as np
+
+class LinearTransferFunction:
+ def __init__(self, arrays = None,
+ filter_mode = drv.filter_mode.LINEAR,
+ address_mode = drv.address_mode.CLAMP,
+ flags = drv.TRSF_NORMALIZED_COORDINATES,
+ range = (0.0, 1.0)):
+
+ (low, up) = range
+ self.scale = 1.0/(up - low)
+ self.offset = low
+ if (low < 0 and up < 0):
+ self.scale *= -1.0
+
+
+ self.range = range
+ self.flags = flags
+ self.filter_mode = filter_mode
+ self.address_mode = address_mode
+ self.interleaved = None
+
+ if arrays != None:
+ self.arrays_to_transfer_function(arrays)
+
+
+ def set_filter_mode(self, mode):
+ """Sets the Texture Filter Mode
+
+ """
+ self.filter_mode = mode
+
+ def set_address_mode(self, mode):
+ self.address_mode = mode
+
+ def set_flags(self, flags):
+ self.flags = flags
+
+ def set_range(self, range):
+ self.range = range
+
+ self.offset = tfh.bounds_to_offset(range)
+ self.scale = tfh.bounds_to_scale(range)
+
+
+ def arrays_to_transfer_function(self, arrays):
+
+ (r_array, g_array, b_array, a_array) = arrays
+ (size, ) = r_array.shape
+
+
+ self.interleaved = np.asarray(np.dstack((r_array.reshape(1, size),
+ g_array.reshape(1, size),
+ b_array.reshape(1, size),
+ a_array.reshape(1, size))),
+ dtype=np.float32, order='F')
+ self.cuda_transfer_array = drv.make_multichannel_2d_array(np.asarray(self.interleaved.transpose((2,1,0)),
+ dtype=np.float32, order='F'), order='F')
+
+ def yt_to_function(self, tf):
+ size = tf.nbins
+ self.interleaved = np.asarray(np.dstack((tf.red.y.reshape(1, size),
+ tf.green.y.reshape(1, size),
+ tf.blue.y.reshape(1, size),
+ tf.alpha.y.reshape(1, size))),
+ dtype=np.float32, order='F')
+
+ self.cuda_transfer_array = drv.make_multichannel_2d_array(np.asarray(self.interleaved.transpose((2,1,0)),
+ dtype=np.float32, order='F'), order='F')
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/volumes/cube.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/volumes/cube.py
@@ -0,0 +1,60 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+import pyRGBA.theia.helpers.cuda as cdh
+
+from pyRGBA.theia.volumes.volume import Volume
+
+import numpy as np
+
+class Unigrid(Volume):
+ def __init__(self, array = None, allocate = False,
+ flags = drv.TRSF_NORMALIZED_COORDINATES,
+ filter_mode = drv.filter_mode.LINEAR,
+ address_mode = drv.address_mode.CLAMP):
+ Volume.__init__(self)
+
+ self.flags = flags
+ self.filter_mode = filter_mode
+ self.address_mode = address_mode
+
+ self.volume = None
+ self.cuda_volume_array = None
+
+ if array != None:
+ self.set_array(array)
+ if allocate == True:
+ self.allocate_and_copy_to_gpu()
+
+
+ def set_filter_mode(self, mode):
+ self.filter_mode = mode
+
+ def set_address_mode(self, mode):
+ self.address_mode = mode
+
+ def set_flags(self, flags):
+ self.flags = flags
+
+ def set_array(self, np_array):
+ self.volume = np_array
+
+ def allocate_and_copy_to_gpu(self):
+ if self.cuda_volume_array:
+ self.deallocate_from_gpu()
+
+ self.cuda_volume_array = cdh.np3d_to_device_array(self.volume)
+
+ def __del__(self) :
+ self.deallocate_from_gpu()
+
+ def deallocate_from_gpu(self):
+ self.cuda_volume_array.free()
+ self.cuda_volume_array = None
diff -r 7bd133d91d3934b8f4617f73b01cd102fb23af8b -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e yt/visualization/volume_rendering/theia/volumes/volume.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/volumes/volume.py
@@ -0,0 +1,12 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+
+class Volume(object):
+ def __init__(self):
+ pass
https://bitbucket.org/yt_analysis/yt/commits/92d6f4ea86f7/
Changeset: 92d6f4ea86f7
Branch: yt-3.0
User: fbogert
Date: 2014-05-01 19:12:37
Summary: more comments
Affected #: 11 files
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
--- a/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
@@ -8,21 +8,25 @@
#include <helpers/cuda.cu>
+//A 3d texture holding scalar values to be volume rendered
texture<float , cudaTextureType3D, cudaReadModeElementType> volume;
+//A 1d texture holding color transfer function values to act on raycaster results
texture<float4, cudaTextureType1D, cudaReadModeElementType> transfer;
extern "C"
-__global__ void front_to_back(uint *buffer, float *actual_matrix, int buffer_w,
+__global__ void front_to_back(uint *buffer, float *modelviewmatrix, int buffer_w,
int buffer_h, int max_steps, float density,
float offset, float scale,
float brightness, float step_size,
) {
+ //Rays will terminate when opacity_threshold is reached
+ const float opacity_threshold = 0.95f;
- const float opacity_threshold = 0.95f;
+ //We don't use the far and near clipping information from the modelviewmatrix
+ float3x4 matrix = mat_array_to_3x4(modelviewmatrix);
- float3x4 matrix = mat_array_to_3x4(actual_matrix);
-
+ //The X,Y coordinate of the surface (image) array
int x = BLOCK_X;
int y = BLOCK_Y;
@@ -38,6 +42,7 @@
float tnear, tfar;
int hit = intersect_std_box(eye_ray, &tnear, &tfar);
+ // If the ray misses the box set the coloro to black
if (!hit) {
buffer[y*buffer_w + x] = 0;
return;
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- a/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -8,13 +8,13 @@
import pycuda.driver as drv
-import pyRGBA.theia.helpers.cuda as cdh
+import yt.visualization.volume_rendering.theia.helpers.cuda as cdh
import numpy as np
-from pyRGBA.theia.transfer.linear_transfer import LinearTransferFunction
+from yt.visualization.volume_rendering.theia.transfer.linear_transfer import LinearTransferFunction
-from pyRGBA.theia.volumes.cube import Unigrid
-from pyRGBA.theia.surfaces.array_surface import ArraySurface
+from yt.visualization.volume_rendering.theia.volumes.unigrid import Unigrid
+from yt.visualization.volume_rendering.theia.surfaces.array_surface import ArraySurface
import os
class FrontToBackRaycaster:
@@ -29,7 +29,7 @@
Example:
- from pyRGBA.cameras.camera import Camera
+ from yt.visualization.volume_rendering.cameras.camera import Camera
cam = Camera()
@@ -76,34 +76,73 @@
def __call__(self):
self.cast()
-
+ """
+ Parameters
+ ----------
+ """
def get_surface(self) :
return self.surface.get_array()
+ """
+ Parameters
+ ----------
+ """
def get_sample_size(self) :
return self.sample_size
+ """
+ Parameters
+ ----------
+ """
def get_max_samples(self) :
return self.max_samples
+ """
+ Parameters
+ ----------
+ """
def get_density_scale(self) :
return self.density_scale
+ """
+ Parameters
+ ----------
+ """
def get_brightness(self):
return self.brightness
+ """
+ Parameters
+ ----------
+ """
def set_sample_size(self, size = 0.01) :
self.sample_size = size
+ """
+ Parameters
+ ----------
+ """
def set_max_samples(self, max = 5000) :
self.max_samples = max
+ """
+ Parameters
+ ----------
+ """
def set_density_scale(self, scale = 0.05) :
self.density_scale = scale
+ """
+ Parameters
+ ----------
+ """
def set_brightness(self, brightness = 1.0):
self.brightness = brightness
+ """
+ Parameters
+ ----------
+ """
def cast(self):
w, h = self.surface.bounds
@@ -121,9 +160,17 @@
)
+ """
+ Parameters
+ ----------
+ """
def set_matrix(self, matrix):
self.matrix = matrix
+ """
+ Parameters
+ ----------
+ """
def set_surface(self, surface = None, block_size = 32):
if surface == None:
self.surface = None
@@ -133,10 +180,18 @@
self.grid = cdh.block_estimate(block_size, self.surface.bounds)
self.block = (block_size, block_size, 1)
+ """
+ Parameters
+ ----------
+ """
def send_volume_to_gpu(self, volume = None) :
if (volume != None) :
self.set_volume(Unigrid(array = volume, allocate = True))
+ """
+ Parameters
+ ----------
+ """
def set_volume(self, volume):
if volume == None:
self.volume = None
@@ -149,6 +204,10 @@
self.volume_identifier.set_address_mode(1, self.volume.address_mode)
self.volume_identifier.set_array(self.volume.cuda_volume_array)
+ """
+ Parameters
+ ----------
+ """
def set_transfer(self, transfer):
if transfer == None:
self.transfer = None
@@ -160,6 +219,18 @@
self.transfer_identifier.set_address_mode(0, self.transfer.address_mode)
self.transfer_identifier.set_array(self.transfer.cuda_transfer_array)
+ """
+ Attach the base directory path to the desired source file.
+ Parameters
+ ----------
+ dir : string
+ Directory where source file is located
+
+ file : string
+ Source file name
+
+
+ """
def base_directory(self, dir, file):
base = os.path.dirname(dir)
src = os.path.join(base, file)
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/camera.py
--- a/yt/visualization/volume_rendering/theia/camera.py
+++ b/yt/visualization/volume_rendering/theia/camera.py
@@ -7,7 +7,7 @@
#-----------------------------------------------------------------------------
#Matix helpers for camera control
-import pyRGBA.theia.helpers.matrix_transformation as mth
+import yt.visualization.volume_rendering.theia.helpers.matrix_transformation as mth
import numpy as np
@@ -28,7 +28,7 @@
Example:
- from pyRGBA.cameras.camera import Camera
+ from yt.visualization.volume_rendering.cameras.camera import Camera
cam = Camera()
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/helpers/cuda.py
--- a/yt/visualization/volume_rendering/theia/helpers/cuda.py
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.py
@@ -16,7 +16,7 @@
def sample_path():
#return '/pfs/sw/cuda-5.5/samples/common/inc'
- return os.environ.get('PYRGBA_CUDA_SAMPLES')
+ return os.environ.get('CUDA_SAMPLES')
def cuda_helpers_path():
return os.path.dirname(__file__) + "/../"
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/scene.py
--- a/yt/visualization/volume_rendering/theia/scene.py
+++ b/yt/visualization/volume_rendering/theia/scene.py
@@ -6,8 +6,8 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
-from pyRGBA.theia.source import TheiaSource
-from pyRGBA.theia.camera import Camera
+from yt.visualization.volume_rendering.theia.source import TheiaSource
+from yt.visualization.volume_rendering.theia.camera import Camera
class TheiaScene:
r"""
@@ -41,8 +41,8 @@
Example:
- from pyRGBA.theia.scene import TheiaScene
- from pyRGBA.algorithms.front_to_back import FrontToBackRaycaster
+ from yt.visualization.volume_rendering.theia.scene import TheiaScene
+ from yt.visualization.volume_rendering.algorithms.front_to_back import FrontToBackRaycaster
import numpy as np
#load 3D numpy array of float32
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/source.py
--- a/yt/visualization/volume_rendering/theia/source.py
+++ b/yt/visualization/volume_rendering/theia/source.py
@@ -50,8 +50,8 @@
Example:
- from pyRGBA.theia.source import TheiaSource
- from pyRGBA.algorithms.front_to_back import FrontToBackRaycaster
+ from yt.visualization.volume_rendering.theia.source import TheiaSource
+ from yt.visualization.volume_rendering.algorithms.front_to_back import FrontToBackRaycaster
import numpy as np
#load a 3d numpy array of float32
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/surfaces/array_surface.py
--- a/yt/visualization/volume_rendering/theia/surfaces/array_surface.py
+++ b/yt/visualization/volume_rendering/theia/surfaces/array_surface.py
@@ -10,7 +10,7 @@
import numpy as np
-from pyRGBA.theia.surfaces.surface import Surface
+from yt.visualization.volume_rendering.theia.surfaces.surface import Surface
"""
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
--- a/yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
+++ b/yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
@@ -15,7 +15,7 @@
"""
import pycuda.driver as drv
-import pyRGBA.theia.transfer.helper as tfh
+import yt.visualizatoin.volume_rendering.theia.transfer.helper as tfh
import numpy as np
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/volumes/cube.py
--- a/yt/visualization/volume_rendering/theia/volumes/cube.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#-----------------------------------------------------------------------------
-# Copyright (c) 2014, 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 pycuda.driver as drv
-import pyRGBA.theia.helpers.cuda as cdh
-
-from pyRGBA.theia.volumes.volume import Volume
-
-import numpy as np
-
-class Unigrid(Volume):
- def __init__(self, array = None, allocate = False,
- flags = drv.TRSF_NORMALIZED_COORDINATES,
- filter_mode = drv.filter_mode.LINEAR,
- address_mode = drv.address_mode.CLAMP):
- Volume.__init__(self)
-
- self.flags = flags
- self.filter_mode = filter_mode
- self.address_mode = address_mode
-
- self.volume = None
- self.cuda_volume_array = None
-
- if array != None:
- self.set_array(array)
- if allocate == True:
- self.allocate_and_copy_to_gpu()
-
-
- def set_filter_mode(self, mode):
- self.filter_mode = mode
-
- def set_address_mode(self, mode):
- self.address_mode = mode
-
- def set_flags(self, flags):
- self.flags = flags
-
- def set_array(self, np_array):
- self.volume = np_array
-
- def allocate_and_copy_to_gpu(self):
- if self.cuda_volume_array:
- self.deallocate_from_gpu()
-
- self.cuda_volume_array = cdh.np3d_to_device_array(self.volume)
-
- def __del__(self) :
- self.deallocate_from_gpu()
-
- def deallocate_from_gpu(self):
- self.cuda_volume_array.free()
- self.cuda_volume_array = None
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/volumes/unigrid.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/volumes/unigrid.py
@@ -0,0 +1,60 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+import yt.visualization.volume_rendering.theia.helpers.cuda as cdh
+
+from yt.visualization.volume_rendering.theia.volumes.volume import Volume
+
+import numpy as np
+
+class Unigrid(Volume):
+ def __init__(self, array = None, allocate = False,
+ flags = drv.TRSF_NORMALIZED_COORDINATES,
+ filter_mode = drv.filter_mode.LINEAR,
+ address_mode = drv.address_mode.CLAMP):
+ Volume.__init__(self)
+
+ self.flags = flags
+ self.filter_mode = filter_mode
+ self.address_mode = address_mode
+
+ self.volume = None
+ self.cuda_volume_array = None
+
+ if array != None:
+ self.set_array(array)
+ if allocate == True:
+ self.allocate_and_copy_to_gpu()
+
+
+ def set_filter_mode(self, mode):
+ self.filter_mode = mode
+
+ def set_address_mode(self, mode):
+ self.address_mode = mode
+
+ def set_flags(self, flags):
+ self.flags = flags
+
+ def set_array(self, np_array):
+ self.volume = np_array
+
+ def allocate_and_copy_to_gpu(self):
+ if self.cuda_volume_array:
+ self.deallocate_from_gpu()
+
+ self.cuda_volume_array = cdh.np3d_to_device_array(self.volume)
+
+ def __del__(self) :
+ self.deallocate_from_gpu()
+
+ def deallocate_from_gpu(self):
+ self.cuda_volume_array.free()
+ self.cuda_volume_array = None
diff -r e43c9e8bb76981e5cb8a912946370ab45d8a2d2e -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 yt/visualization/volume_rendering/theia/volumes/volume.py
--- a/yt/visualization/volume_rendering/theia/volumes/volume.py
+++ b/yt/visualization/volume_rendering/theia/volumes/volume.py
@@ -6,7 +6,10 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+"""
+A Base Class for all volumes.
+"""
class Volume(object):
def __init__(self):
pass
https://bitbucket.org/yt_analysis/yt/commits/faac00a4ab44/
Changeset: faac00a4ab44
Branch: yt-3.0
User: Alex Bogert
Date: 2014-05-01 20:22:01
Summary: Final push
Affected #: 3 files
diff -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
--- a/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
@@ -18,7 +18,7 @@
__global__ void front_to_back(uint *buffer, float *modelviewmatrix, int buffer_w,
int buffer_h, int max_steps, float density,
float offset, float scale,
- float brightness, float step_size,
+ float brightness, float step_size
) {
//Rays will terminate when opacity_threshold is reached
const float opacity_threshold = 0.95f;
diff -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- a/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -12,41 +12,43 @@
import numpy as np
from yt.visualization.volume_rendering.theia.transfer.linear_transfer import LinearTransferFunction
+from yt.visualization.volume_rendering.theia.transfer.helper import *
+
from yt.visualization.volume_rendering.theia.volumes.unigrid import Unigrid
from yt.visualization.volume_rendering.theia.surfaces.array_surface import ArraySurface
import os
class FrontToBackRaycaster:
- r"""
- This is a basic camera intended to be a base class.
- The camera provies utilities for controling the view point onto the data.
+ r"""
+ This is the python wrapper for the CUDA raycaster. This
+ raycaster can act on unigrid data only. Use yt to map
+ AMR data to unigrid prior to using this algorithm.
- Parameters
- ----------
- pos : numpy array
+ Parameters
+ ----------
+ matrix : 4x4 numpy.matriix
+ ModelViewMatrix defines view onto volume
- Example:
+ surface: yt.visualization.volume_rendering.theia.surface.array_surface.ArraySurface
+ The surface to contain the image information from
+ the results of the raycasting
- from yt.visualization.volume_rendering.cameras.camera import Camera
+ volume: 3D numpy array Float32
+ Scalar values to be volume rendered
- cam = Camera()
-
- cam.zoom(10.0)
- cam.rotateX(np.pi)
- cam.rotateY(np.pi)
- cam.rotateZ(np.pi)
- cam.translateX(-0.5)
- cam.translateY(-0.5)
- cam.translateZ(-0.5)
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Transfer function defining how the raycasting results are
+ to be colored.
- view = cam.get_matrix()
-
- """
+ size : Tuple of 2 Floats
+ If no surface is specified creates an ArraySurface of the specified size.
+
+ """
def __init__(self, matrix = None, surface = None,
- volume = None, transfer = None,
+ volume = None, tf = None,
size = (864, 480)) :
# kernel module
@@ -54,9 +56,12 @@
self.base_directory(__file__, "front_to_back.cu"))
#kernel function definitions
- self.cast_rays_identifier = self.kernel_module.get_function("front_to_back")
- self.transfer_identifier = self.kernel_module.get_texref("transfer")
- self.volume_identifier = self.kernel_module.get_texref("volume")
+ self.cast_rays_identifier = \
+ self.kernel_module.get_function("front_to_back")
+ self.transfer_identifier = \
+ self.kernel_module.get_texref("transfer")
+ self.volume_identifier = \
+ self.kernel_module.get_texref("volume")
self.set_matrix(matrix)
@@ -67,46 +72,58 @@
self.volume = None
- self.set_transfer(transfer)
+ self.set_transfer(tf)
self.set_sample_size()
self.set_max_samples()
self.set_density_scale()
self.set_brightness()
+ """
+ Envoke the ray caster to cast rays
+ """
def __call__(self):
self.cast()
"""
- Parameters
+ Returns
----------
+ A 2d numpy array containing the image of the
+ volumetric rendering
"""
def get_surface(self) :
return self.surface.get_array()
"""
- Parameters
+ Returns
----------
+ The sample size the ray caster is
+ configured to take.
"""
def get_sample_size(self) :
return self.sample_size
"""
- Parameters
+ Returns
----------
+ The Max number of samples per ray the ray caster is
+ configured to take.
"""
def get_max_samples(self) :
return self.max_samples
"""
- Parameters
+ Returns
----------
+ The Global density scalar.
"""
def get_density_scale(self) :
return self.density_scale
"""
- Parameters
+ Returns
----------
+ The Global brightness scalar.
+
"""
def get_brightness(self):
return self.brightness
@@ -114,6 +131,9 @@
"""
Parameters
----------
+ size : scalra Float
+ The distance between each sample, a smaller size will result in
+ more samples and can cause performance loss.
"""
def set_sample_size(self, size = 0.01) :
self.sample_size = size
@@ -121,6 +141,9 @@
"""
Parameters
----------
+ max : scalar Int
+ The limit on the number of samples the ray caster will
+ take per ray.
"""
def set_max_samples(self, max = 5000) :
self.max_samples = max
@@ -128,6 +151,8 @@
"""
Parameters
----------
+ scale : scalar Float
+ Global multiplier on volume data
"""
def set_density_scale(self, scale = 0.05) :
self.density_scale = scale
@@ -135,13 +160,15 @@
"""
Parameters
----------
+ brightness : scalar Float
+ Global multiplier on returned alpha values
"""
def set_brightness(self, brightness = 1.0):
self.brightness = brightness
"""
- Parameters
- ----------
+ Causes the ray caster to act on the volume data
+ and puts the results in the surface array.
"""
def cast(self):
w, h = self.surface.bounds
@@ -161,15 +188,22 @@
"""
+ Set the ModelViewMatrix, this does not send data to the GPU
Parameters
----------
+ matrix : 4x4 numpy.matrix
+ ModelViewMatrix
+
"""
def set_matrix(self, matrix):
self.matrix = matrix
"""
+ Setup the image array for ray casting results
Parameters
----------
+ surface : yt.visualization..volume_rendering.theia.surfaces.ArraySurface
+ Surface to contain results of the ray casting
"""
def set_surface(self, surface = None, block_size = 32):
if surface == None:
@@ -181,8 +215,13 @@
self.block = (block_size, block_size, 1)
"""
+ This function will convert a 3d numpy array into a unigrid volume
+ and move the data to the gpu texture memory
Parameters
----------
+ volume : 3D numpy array float32
+ Contains scalar volumes to be acted on by the ray caster
+
"""
def send_volume_to_gpu(self, volume = None) :
if (volume != None) :
@@ -191,8 +230,11 @@
"""
Parameters
----------
+ volume : yt.visualization..volume_rendering.theia.volumes.Unigrid
+ Contains scalar volumes to be acted on by the ray caster
+
"""
- def set_volume(self, volume):
+ def set_volume(self, volume = None):
if volume == None:
self.volume = None
return
@@ -207,13 +249,16 @@
"""
Parameters
----------
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Used to color the results of the raycasting
+
"""
- def set_transfer(self, transfer):
- if transfer == None:
+ def set_transfer(self, tf = None):
+ if tf == None:
self.transfer = None
return
- self.transfer = transfer
+ self.transfer = LinearTransferFunction(arrays = yt_to_rgba_transfer(tf), range = tf.x_bounds)
self.transfer_identifier.set_flags(self.transfer.flags)
self.transfer_identifier.set_filter_mode(self.transfer.filter_mode)
self.transfer_identifier.set_address_mode(0, self.transfer.address_mode)
diff -r 92d6f4ea86f79e4fd1bdd30f3afe54247c5f4fe6 -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
--- a/yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
+++ b/yt/visualization/volume_rendering/theia/transfer/linear_transfer.py
@@ -15,7 +15,7 @@
"""
import pycuda.driver as drv
-import yt.visualizatoin.volume_rendering.theia.transfer.helper as tfh
+import yt.visualization.volume_rendering.theia.transfer.helper as tfh
import numpy as np
https://bitbucket.org/yt_analysis/yt/commits/944adbdb415a/
Changeset: 944adbdb415a
Branch: yt-3.0
User: Alex Bogert
Date: 2014-05-26 20:00:23
Summary: added doc strings and renamed density_scale as opacity
Affected #: 6 files
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 doc/source/cookbook/amrkdtree_to_uniformgrid.py
--- /dev/null
+++ b/doc/source/cookbook/amrkdtree_to_uniformgrid.py
@@ -0,0 +1,16 @@
+import numpy as np
+import yt
+
+ds = yt.load("~/pfs/galaxy/new_tests/feedback_8bz/DD0021/DD0021")
+
+lmax = ds.parameters['MaximumRefinementLevel']
+domain_center = (ds.domain_right_edge - ds.domain_left_edge)/2
+cell_size = pf.domain_width/(pf.domain_dimensions*2**lmax)
+left_edge = domain_center - 512*cell_size
+
+ncells = 1024
+cgrid = pf.h.covering_grid(lmax, left_edge, np.array([ncells,]*3))
+density_map = cgrid["density"].astype(dtype="float32")
+print density_map.min()
+print density_map.max()
+np.save("/pfs/goldbaum/galaxy/new_tests/feedback_8bz/gas_density_DD0021_log_densities.npy", density_map)
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 doc/source/cookbook/ffmpeg_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/ffmpeg_volume_rendering.py
@@ -0,0 +1,59 @@
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import subprocess as sp
+
+import numpy as np
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, (v-mi)/(ma-mi) + 0.0)
+
+bolshoi = "/home/bogert/log_densities_1024.npy"
+density_grid = np.load(bolshoi)
+
+ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (1920,1080) ))
+
+mi, ma = 0.0, 3.6
+bins = 5000
+tf = ColorTransferFunction( (mi, ma), bins)
+tf.map_to_colormap(0.5, ma, colormap="spring", scale_func = scale_func)
+ts.source.raycaster.set_transfer(tf)
+
+ts.source.raycaster.set_density_scale(0.03)
+ts.source.raycaster.set_brightness(2.3)
+ts.camera.zoom(30.0)
+
+FFMPEG_BIN = "ffmpeg"
+
+pipe = sp.Popen([ FFMPEG_BIN,
+ '-y', # (optional) overwrite the output file if it already exists
+ '-f', 'rawvideo',
+ '-vcodec','rawvideo',
+ '-s', '1920x1080', # size of one frame
+ '-pix_fmt', 'rgba',
+ '-r', '9.97', # frames per second
+ '-i', '-', # The input comes from a pipe
+ '-an', # Tells FFMPEG not to expect any audio
+ #'-vcodec', 'dnxhd', '-b:v', '220M',
+ '-vcodec', 'libx264', '-preset', 'ultrafast', '-qp', '0',# '-b:v', '220M',
+ '-pix_fmt', 'yuv420p',
+ 'bolshoiplanck2.mkv' ],
+ stdin=sp.PIPE,stdout=sp.PIPE)
+
+
+#ts.source.surface.bounds = (1920,1080)
+for k in range (0,50) :
+ ts.update()
+ ts.camera.rotateX(0.01)
+ ts.camera.rotateZ(0.01)
+ ts.camera.rotateY(0.01)
+ ts.camera.zoom(0.01)
+
+ array = ts.source.get_results()
+
+ array.tofile(pipe.stdin)
+
+pipe.terminate()
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 doc/source/cookbook/opengl_stereo_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_stereo_volume_rendering.py
@@ -0,0 +1,370 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1920
+height = 1080
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+(rpbo, rpycuda_pbo) = [None]*2
+
+#create 2 PBO for stereo scopic rendering
+def create_PBO(w, h):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ num_texels = w*h
+ array = np.zeros((num_texels, 3),np.float32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+ rpbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, rpbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpycuda_pbo = cuda_gl.RegisteredBuffer(long(rpbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+ glBindBuffer(GL_ARRAY_BUFFER, long(rpbo))
+ glDeleteBuffers(1, long(rpbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpbo,rpycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ w, h, 0, GL_RGB, GL_FLOAT, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity, eyesep
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+ #change the transfer scale
+ elif args[0] == '-':
+ eyesep -= 0.01
+ elif args[0] == '=':
+ eyesep += 0.01
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ #process right eye
+ process_image(eye = False)
+ display_image(eye = False)
+
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, rpycuda_pbo, eyesep, c_tbrightness, c_tdensity
+ """ Use PyCuda """
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ if (eye) :
+ ts.camera.translateX(-eyesep)
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(eyesep)
+ else :
+ ts.camera.translateX(eyesep)
+ dest_mapping = rpycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(-eyesep)
+
+
+def process_image(eye = True):
+ global output_texture, pbo, rpbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process(eye)
+ # download texture from PBO
+ if (eye) :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+ else :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(rpbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ if (eye) :
+ glDrawBuffer(GL_BACK_LEFT)
+ else :
+ glDrawBuffer(GL_BACK_RIGHT)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1920, 1080)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH | GLUT_STEREO)
+ glutInitWindowSize(*initial_size)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 doc/source/cookbook/opengl_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_volume_rendering.py
@@ -0,0 +1,327 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1920
+height = 1080
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+
+def create_PBO(w, h):
+ global pbo, pycuda_pbo
+ num_texels = w*h
+ array = np.zeros((num_texels, 3),np.float32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ w, h, 0, GL_RGB, GL_FLOAT, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity, eyesep
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+ #change the transfer scale
+ elif args[0] == '-':
+ eyesep -= 0.01
+ elif args[0] == '=':
+ eyesep += 0.01
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, eyesep, c_tbrightness, c_tdensity
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+
+
+def process_image(eye = True):
+ global output_texture, pbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process()
+ # download texture from PBO
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1920, 1080)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH )
+ glutInitWindowSize(*initial_size)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 doc/source/visualizing/volume_rendering.rst
--- a/doc/source/visualizing/volume_rendering.rst
+++ b/doc/source/visualizing/volume_rendering.rst
@@ -466,3 +466,99 @@
your homogenized volume to then be passed in to the camera. A sample usage is shown
in :ref:`cookbook-amrkdtree_downsampling`.
+Hardware Volume Rendering on NVidia Graphics cards
+--------------------------------------------------
+.. versionadded:: 3.0
+As of yt 3.0, Theia has been added to yt.
+
+Theia is a hardware volume renderer that takes advantage of NVidias CUDA language
+to peform ray casting with GPUs instead of the CPU.
+
+Only unigrid rendering is supported, but yt provides a grid mapping function
+ to get unigrid data from amr or sph formats :
+ :ref:`cookbook-amrkdtree_to_uniformgrid`.
+
+System Requirements
+-------------------
+.. versionadded:: 3.0
+
+Nvidia graphics card - The memory limit of the graphics card sets the limit
+ on the size of the data source.
+
+CUDA 5 or later and
+
+The environment variable CUDA_SAMPLES must be set pointing to
+the common/inc samples shipped with CUDA. The following shows an example
+in bash with CUDA 5.5 installed in /usr/local :
+
+export CUDA_SAMPLES=/usr/local/cuda-5.5/samples/common/inc
+
+PyCUDA must also be installed to use Theia.
+
+PyCUDA can be installed following these instructions :
+
+ git clone --recursive http://git.tiker.net/trees/pycuda.git
+
+ python configure.py
+ python setup.py install
+
+
+Tutorial
+--------
+.. versionadded:: 3.0
+
+Currently rendering only works on uniform grids. Here is an example
+on a 1024 cube of float32 scalars.
+
+.. code-block:: python
+ from yt.visualization.volume_rendering.theia.scene import TheiaScene
+ from yt.visualization.volume_rendering.algorithms.front_to_back import FrontToBackRaycaster
+ import numpy as np
+
+ #load 3D numpy array of float32
+ volume = np.load("/home/bogert/log_densities_1024.npy")
+
+ scene = TheiaScene( volume = volume, raycaster = FrontToBackRaycaster() )
+
+ scene.camera.rotateX(1.0)
+ scene.update()
+
+ surface = scene.get_results()
+ #surface now contains an image array 2x2 int32 rbga values
+
+.. _the-theiascene-interface:
+
+The TheiaScene Interface
+--------------------
+.. versionadded:: 3.0
+
+A TheiaScene object has been created to provide a high level entry point for
+controlling the raycaster's view onto the data. The class
+:class:`~yt.visualization.volume_rendering.theia.TheiaScene` encapsulates
+ a Camera object and a TheiaSource that intern encapsulates
+a volume. The :class:`~yt.visualization.volume_rendering.theia.Camera`
+provides controls for rotating, translating, and zooming into the volume.
+Using the :class:`~yt.visualization.volume_rendering.theia.TheiaSource`
+automatically transfers the volume to the graphic's card texture memory.
+
+
+.. warning:: The keyword *no_ghost* has been set to True by default
+ for speed considerations. However, because this turns off ghost
+ zones, there may be artifacts at grid boundaries. If a higher quality
+ rendering is required, use *no_ghost = False*.
+
+.. _camera_movement:
+
+Example Cookbooks
+---------------
+
+OpenGL Example for interactive volume rendering:
+:ref:`cookbook-opengl_volume_rendering`.
+
+OpenGL Stereoscopic Example :
+.. warning:: Frame rate will suffer significantly from stereoscopic rendering.
+ ~2x slower since the volume must be rendered twice.
+:ref:`cookbook-opengl_stereo_volume_rendering`.
+
+Pseudo-Realtime video rendering with ffmpeg :
+:ref:`cookbook-ffmpeg_volume_rendering`.
diff -r faac00a4ab44f455f9099c5eb8d6f08bd15eb570 -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- a/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -76,7 +76,7 @@
self.set_sample_size()
self.set_max_samples()
- self.set_density_scale()
+ self.set_opacity()
self.set_brightness()
"""
@@ -116,8 +116,8 @@
----------
The Global density scalar.
"""
- def get_density_scale(self) :
- return self.density_scale
+ def get_opacity(self) :
+ return self.opacity
"""
Returns
@@ -154,8 +154,8 @@
scale : scalar Float
Global multiplier on volume data
"""
- def set_density_scale(self, scale = 0.05) :
- self.density_scale = scale
+ def set_opacity(self, scale = 0.05) :
+ self.opacity = scale
"""
Parameters
@@ -177,7 +177,7 @@
drv.In(self.matrix),
np.int32(w), np.int32(h),
np.int32(self.max_samples),
- np.float32(self.density_scale),
+ np.float32(self.opacity),
np.float32(self.transfer.offset),
np.float32(self.transfer.scale),
np.float32(self.brightness),
https://bitbucket.org/yt_analysis/yt/commits/61fbae4a3e5c/
Changeset: 61fbae4a3e5c
Branch: yt-3.0
User: Alex Bogert
Date: 2014-05-26 20:07:19
Summary: saw some copy pasta
Affected #: 1 file
diff -r 944adbdb415aae22a8c253cd45df1f3d44cc2a82 -r 61fbae4a3e5cf16d84fec02222ecb8f90e85d8f3 doc/source/visualizing/volume_rendering.rst
--- a/doc/source/visualizing/volume_rendering.rst
+++ b/doc/source/visualizing/volume_rendering.rst
@@ -541,14 +541,6 @@
Using the :class:`~yt.visualization.volume_rendering.theia.TheiaSource`
automatically transfers the volume to the graphic's card texture memory.
-
-.. warning:: The keyword *no_ghost* has been set to True by default
- for speed considerations. However, because this turns off ghost
- zones, there may be artifacts at grid boundaries. If a higher quality
- rendering is required, use *no_ghost = False*.
-
-.. _camera_movement:
-
Example Cookbooks
---------------
https://bitbucket.org/yt_analysis/yt/commits/e68a8cd972f2/
Changeset: e68a8cd972f2
Branch: yt-3.0
User: Alex Bogert
Date: 2014-06-06 20:38:28
Summary: Addressed Cameron's comments. Opengl example is now working, haven't gotten to the 3d example yet.
Affected #: 4 files
diff -r 61fbae4a3e5cf16d84fec02222ecb8f90e85d8f3 -r e68a8cd972f210230cd8d1ca5a316d76bee3d204 doc/source/cookbook/amrkdtree_to_uniformgrid.py
--- a/doc/source/cookbook/amrkdtree_to_uniformgrid.py
+++ b/doc/source/cookbook/amrkdtree_to_uniformgrid.py
@@ -1,16 +1,33 @@
import numpy as np
import yt
+#This is an example of how to map an amr data set
+#to a uniform grid. In this case the highest
+#level of refinement is mapped into a 1024x1024x1024 cube
+
+#first the amr data is loaded
ds = yt.load("~/pfs/galaxy/new_tests/feedback_8bz/DD0021/DD0021")
+#next we get the maxium refinement level
lmax = ds.parameters['MaximumRefinementLevel']
+
+#calculate the center of the domain
domain_center = (ds.domain_right_edge - ds.domain_left_edge)/2
+
+#determine the cellsize in the highest refinement level
cell_size = pf.domain_width/(pf.domain_dimensions*2**lmax)
+
+#calculate the left edge of the new grid
left_edge = domain_center - 512*cell_size
+#the number of cells per side of the new grid
ncells = 1024
+
+#ask yt for the specified covering grid
cgrid = pf.h.covering_grid(lmax, left_edge, np.array([ncells,]*3))
+
+#get a map of the density into the new grid
density_map = cgrid["density"].astype(dtype="float32")
-print density_map.min()
-print density_map.max()
+
+#save the file as a numpy array for convenient future processing
np.save("/pfs/goldbaum/galaxy/new_tests/feedback_8bz/gas_density_DD0021_log_densities.npy", density_map)
diff -r 61fbae4a3e5cf16d84fec02222ecb8f90e85d8f3 -r e68a8cd972f210230cd8d1ca5a316d76bee3d204 doc/source/cookbook/ffmpeg_volume_rendering.py
--- a/doc/source/cookbook/ffmpeg_volume_rendering.py
+++ b/doc/source/cookbook/ffmpeg_volume_rendering.py
@@ -1,59 +1,99 @@
+#This is an example of how to make videos of
+#uniform grid data using Theia and ffmpeg
+
+#The Scene object to hold the ray caster and view camera
from yt.visualization.volume_rendering.theia.scene import TheiaScene
+
+#GPU based raycasting algorithm to use
from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+#These will be used to define how to color the data
from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
from yt.visualization.color_maps import *
+#This will be used to launch ffmpeg
import subprocess as sp
+#Of course we need numpy for math magic
import numpy as np
+#Opacity scaling function
def scale_func(v, mi, ma):
return np.minimum(1.0, (v-mi)/(ma-mi) + 0.0)
+#load the uniform grid from a numpy array file
bolshoi = "/home/bogert/log_densities_1024.npy"
density_grid = np.load(bolshoi)
+#Set the TheiaScene to use the density_grid and
+#setup the raycaster for a resulting 1080p image
ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (1920,1080) ))
+#the min and max values in the data to color
mi, ma = 0.0, 3.6
+
+#setup colortransferfunction
bins = 5000
tf = ColorTransferFunction( (mi, ma), bins)
tf.map_to_colormap(0.5, ma, colormap="spring", scale_func = scale_func)
+
+#pass the transfer function to the ray caster
ts.source.raycaster.set_transfer(tf)
-ts.source.raycaster.set_density_scale(0.03)
+#Initial configuration for start of video
+#set initial opacity and brightness values
+#then zoom into the center of the data 30%
+ts.source.raycaster.set_opacity(0.03)
ts.source.raycaster.set_brightness(2.3)
ts.camera.zoom(30.0)
-FFMPEG_BIN = "ffmpeg"
+#path to ffmpeg executable
+FFMPEG_BIN = "/usr/local/bin/ffmpeg"
pipe = sp.Popen([ FFMPEG_BIN,
'-y', # (optional) overwrite the output file if it already exists
- '-f', 'rawvideo',
+ #This must be set to rawvideo because the image is an array
+ '-f', 'rawvideo',
+ #This must be set to rawvideo because the image is an array
'-vcodec','rawvideo',
- '-s', '1920x1080', # size of one frame
+ #The size of the image array and resulting video
+ '-s', '1920x1080',
+ #This must be rgba to match array format (uint32)
'-pix_fmt', 'rgba',
- '-r', '9.97', # frames per second
- '-i', '-', # The input comes from a pipe
- '-an', # Tells FFMPEG not to expect any audio
- #'-vcodec', 'dnxhd', '-b:v', '220M',
- '-vcodec', 'libx264', '-preset', 'ultrafast', '-qp', '0',# '-b:v', '220M',
+ #frame rate of video
+ '-r', '29.97',
+ #Indicate that the input to ffmpeg comes from a pipe
+ '-i', '-',
+ # Tells FFMPEG not to expect any audio
+ '-an',
+ #Setup video encoder
+ #Use any encoder you life available from ffmpeg
+ '-vcodec', 'libx264', '-preset', 'ultrafast', '-qp', '0',
'-pix_fmt', 'yuv420p',
+ #Name of the output
'bolshoiplanck2.mkv' ],
stdin=sp.PIPE,stdout=sp.PIPE)
-#ts.source.surface.bounds = (1920,1080)
-for k in range (0,50) :
+#Now we loop and produce 500 frames
+for k in range (0,500) :
+ #update the scene resulting in a new image
ts.update()
+
+ #get the image array from the ray caster
+ array = ts.source.get_results()
+
+ #send the image array to ffmpeg
+ array.tofile(pipe.stdin)
+
+ #rotate the scene by 0.01 rads in x,y & z
ts.camera.rotateX(0.01)
ts.camera.rotateZ(0.01)
ts.camera.rotateY(0.01)
+
+ #zoom in 0.01% for a total of a 5% zoom
ts.camera.zoom(0.01)
- array = ts.source.get_results()
- array.tofile(pipe.stdin)
-
+#Close the pipe to ffmpeg
pipe.terminate()
diff -r 61fbae4a3e5cf16d84fec02222ecb8f90e85d8f3 -r e68a8cd972f210230cd8d1ca5a316d76bee3d204 doc/source/cookbook/opengl_volume_rendering.py
--- a/doc/source/cookbook/opengl_volume_rendering.py
+++ b/doc/source/cookbook/opengl_volume_rendering.py
@@ -32,8 +32,8 @@
rightButton = False
#Screen width and height
-width = 1920
-height = 1080
+width = 1024
+height = 1024
eyesep = 0.1
@@ -42,7 +42,7 @@
def create_PBO(w, h):
global pbo, pycuda_pbo
num_texels = w*h
- array = np.zeros((num_texels, 3),np.float32)
+ array = np.zeros((w,h,3),np.uint32)
pbo = glGenBuffers(1)
glBindBuffer(GL_ARRAY_BUFFER, pbo)
@@ -68,8 +68,8 @@
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
# buffer data
- glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
- w, h, 0, GL_RGB, GL_FLOAT, None)
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ w, h, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8, None)
#consistent with C initPixelBuffer()
def destroy_texture():
@@ -156,19 +156,19 @@
factor = 0.001
if leftButton == True:
- ts.camera.rotateX( - deltaY * factor)
- ts.camera.rotateY( - deltaX * factor)
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
if middleButton == True:
- ts.camera.translateX( deltaX* 2.0 * factor)
- ts.camera.translateY( - deltaY* 2.0 * factor)
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
if rightButton == True:
- ts.camera.scale += deltaY * factor
+ ts.camera.scale += deltaY * factor
oldMousePos[0], oldMousePos[1] = x, y
glutPostRedisplay( )
def keyPressed(*args):
- global c_tbrightness, c_tdensity, eyesep
+ global c_tbrightness, c_tdensity
# If escape is pressed, kill everything.
if args[0] == '\033':
print 'Closing..'
@@ -188,12 +188,6 @@
elif args[0] == '\'':
c_tdensity += 0.001
- #change the transfer scale
- elif args[0] == '-':
- eyesep -= 0.01
- elif args[0] == '=':
- eyesep += 0.01
-
def idle():
glutPostRedisplay()
@@ -221,10 +215,11 @@
(dev_ptr, size) = dest_mapping.device_ptr_and_size()
ts.get_raycaster().surface.device_ptr = dev_ptr
ts.update()
+ # ts.get_raycaster().cast()
dest_mapping.unmap()
-def process_image(eye = True):
+def process_image():
global output_texture, pbo, width, height
""" copy image and process using CUDA """
# run the Cuda kernel
@@ -233,9 +228,8 @@
glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
glBindTexture(GL_TEXTURE_2D, output_texture)
- glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
- width, height, 0,
- GL_RGB, GL_FLOAT, None)
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ width, height, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8_REV, None)
def display_image(eye = True):
global width, height
@@ -277,11 +271,11 @@
#note we may need to init cuda_gl here and pass it to camera
def main():
global window, ts, width, height
- (width, height) = (1920, 1080)
+ (width, height) = (1024, 1024)
glutInit(sys.argv)
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH )
- glutInitWindowSize(*initial_size)
+ glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)
window = glutCreateWindow("Stereo Volume Rendering")
@@ -315,6 +309,7 @@
ts.get_raycaster().set_sample_size(0.01)
ts.get_raycaster().set_max_samples(5000)
+ ts.update()
glutMainLoop()
diff -r 61fbae4a3e5cf16d84fec02222ecb8f90e85d8f3 -r e68a8cd972f210230cd8d1ca5a316d76bee3d204 doc/source/visualizing/volume_rendering.rst
--- a/doc/source/visualizing/volume_rendering.rst
+++ b/doc/source/visualizing/volume_rendering.rst
@@ -469,7 +469,6 @@
Hardware Volume Rendering on NVidia Graphics cards
--------------------------------------------------
.. versionadded:: 3.0
-As of yt 3.0, Theia has been added to yt.
Theia is a hardware volume renderer that takes advantage of NVidias CUDA language
to peform ray casting with GPUs instead of the CPU.
https://bitbucket.org/yt_analysis/yt/commits/b6b78bc7de4c/
Changeset: b6b78bc7de4c
Branch: yt-3.0
User: MatthewTurk
Date: 2014-06-15 14:45:51
Summary: Merged in fbogert/yt-theia/yt-3.0 (pull request #862)
CUDA Volume Renderer
Affected #: 26 files
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa doc/source/cookbook/amrkdtree_to_uniformgrid.py
--- /dev/null
+++ b/doc/source/cookbook/amrkdtree_to_uniformgrid.py
@@ -0,0 +1,33 @@
+import numpy as np
+import yt
+
+#This is an example of how to map an amr data set
+#to a uniform grid. In this case the highest
+#level of refinement is mapped into a 1024x1024x1024 cube
+
+#first the amr data is loaded
+ds = yt.load("~/pfs/galaxy/new_tests/feedback_8bz/DD0021/DD0021")
+
+#next we get the maxium refinement level
+lmax = ds.parameters['MaximumRefinementLevel']
+
+#calculate the center of the domain
+domain_center = (ds.domain_right_edge - ds.domain_left_edge)/2
+
+#determine the cellsize in the highest refinement level
+cell_size = pf.domain_width/(pf.domain_dimensions*2**lmax)
+
+#calculate the left edge of the new grid
+left_edge = domain_center - 512*cell_size
+
+#the number of cells per side of the new grid
+ncells = 1024
+
+#ask yt for the specified covering grid
+cgrid = pf.h.covering_grid(lmax, left_edge, np.array([ncells,]*3))
+
+#get a map of the density into the new grid
+density_map = cgrid["density"].astype(dtype="float32")
+
+#save the file as a numpy array for convenient future processing
+np.save("/pfs/goldbaum/galaxy/new_tests/feedback_8bz/gas_density_DD0021_log_densities.npy", density_map)
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa doc/source/cookbook/ffmpeg_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/ffmpeg_volume_rendering.py
@@ -0,0 +1,99 @@
+#This is an example of how to make videos of
+#uniform grid data using Theia and ffmpeg
+
+#The Scene object to hold the ray caster and view camera
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+
+#GPU based raycasting algorithm to use
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+
+#These will be used to define how to color the data
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+#This will be used to launch ffmpeg
+import subprocess as sp
+
+#Of course we need numpy for math magic
+import numpy as np
+
+#Opacity scaling function
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, (v-mi)/(ma-mi) + 0.0)
+
+#load the uniform grid from a numpy array file
+bolshoi = "/home/bogert/log_densities_1024.npy"
+density_grid = np.load(bolshoi)
+
+#Set the TheiaScene to use the density_grid and
+#setup the raycaster for a resulting 1080p image
+ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (1920,1080) ))
+
+#the min and max values in the data to color
+mi, ma = 0.0, 3.6
+
+#setup colortransferfunction
+bins = 5000
+tf = ColorTransferFunction( (mi, ma), bins)
+tf.map_to_colormap(0.5, ma, colormap="spring", scale_func = scale_func)
+
+#pass the transfer function to the ray caster
+ts.source.raycaster.set_transfer(tf)
+
+#Initial configuration for start of video
+#set initial opacity and brightness values
+#then zoom into the center of the data 30%
+ts.source.raycaster.set_opacity(0.03)
+ts.source.raycaster.set_brightness(2.3)
+ts.camera.zoom(30.0)
+
+#path to ffmpeg executable
+FFMPEG_BIN = "/usr/local/bin/ffmpeg"
+
+pipe = sp.Popen([ FFMPEG_BIN,
+ '-y', # (optional) overwrite the output file if it already exists
+ #This must be set to rawvideo because the image is an array
+ '-f', 'rawvideo',
+ #This must be set to rawvideo because the image is an array
+ '-vcodec','rawvideo',
+ #The size of the image array and resulting video
+ '-s', '1920x1080',
+ #This must be rgba to match array format (uint32)
+ '-pix_fmt', 'rgba',
+ #frame rate of video
+ '-r', '29.97',
+ #Indicate that the input to ffmpeg comes from a pipe
+ '-i', '-',
+ # Tells FFMPEG not to expect any audio
+ '-an',
+ #Setup video encoder
+ #Use any encoder you life available from ffmpeg
+ '-vcodec', 'libx264', '-preset', 'ultrafast', '-qp', '0',
+ '-pix_fmt', 'yuv420p',
+ #Name of the output
+ 'bolshoiplanck2.mkv' ],
+ stdin=sp.PIPE,stdout=sp.PIPE)
+
+
+#Now we loop and produce 500 frames
+for k in range (0,500) :
+ #update the scene resulting in a new image
+ ts.update()
+
+ #get the image array from the ray caster
+ array = ts.source.get_results()
+
+ #send the image array to ffmpeg
+ array.tofile(pipe.stdin)
+
+ #rotate the scene by 0.01 rads in x,y & z
+ ts.camera.rotateX(0.01)
+ ts.camera.rotateZ(0.01)
+ ts.camera.rotateY(0.01)
+
+ #zoom in 0.01% for a total of a 5% zoom
+ ts.camera.zoom(0.01)
+
+
+#Close the pipe to ffmpeg
+pipe.terminate()
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa doc/source/cookbook/opengl_stereo_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_stereo_volume_rendering.py
@@ -0,0 +1,370 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1920
+height = 1080
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+(rpbo, rpycuda_pbo) = [None]*2
+
+#create 2 PBO for stereo scopic rendering
+def create_PBO(w, h):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ num_texels = w*h
+ array = np.zeros((num_texels, 3),np.float32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+ rpbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, rpbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpycuda_pbo = cuda_gl.RegisteredBuffer(long(rpbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo, rpbo, rpycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+ glBindBuffer(GL_ARRAY_BUFFER, long(rpbo))
+ glDeleteBuffers(1, long(rpbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ rpbo,rpycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ w, h, 0, GL_RGB, GL_FLOAT, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity, eyesep
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+ #change the transfer scale
+ elif args[0] == '-':
+ eyesep -= 0.01
+ elif args[0] == '=':
+ eyesep += 0.01
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ #process right eye
+ process_image(eye = False)
+ display_image(eye = False)
+
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, rpycuda_pbo, eyesep, c_tbrightness, c_tdensity
+ """ Use PyCuda """
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ if (eye) :
+ ts.camera.translateX(-eyesep)
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(eyesep)
+ else :
+ ts.camera.translateX(eyesep)
+ dest_mapping = rpycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ dest_mapping.unmap()
+ ts.camera.translateX(-eyesep)
+
+
+def process_image(eye = True):
+ global output_texture, pbo, rpbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process(eye)
+ # download texture from PBO
+ if (eye) :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+ else :
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(rpbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB,
+ width, height, 0,
+ GL_RGB, GL_FLOAT, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ if (eye) :
+ glDrawBuffer(GL_BACK_LEFT)
+ else :
+ glDrawBuffer(GL_BACK_RIGHT)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1920, 1080)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH | GLUT_STEREO)
+ glutInitWindowSize(*initial_size)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa doc/source/cookbook/opengl_volume_rendering.py
--- /dev/null
+++ b/doc/source/cookbook/opengl_volume_rendering.py
@@ -0,0 +1,322 @@
+from OpenGL.GL import *
+from OpenGL.GLUT import *
+from OpenGL.GLU import *
+from OpenGL.GL.ARB.vertex_buffer_object import *
+
+import sys, time
+import numpy as np
+import pycuda.driver as cuda_driver
+import pycuda.gl as cuda_gl
+
+from yt.visualization.volume_rendering.theia.scene import TheiaScene
+from yt.visualization.volume_rendering.theia.algorithms.front_to_back import FrontToBackRaycaster
+from yt.visualization.volume_rendering.transfer_functions import ColorTransferFunction
+from yt.visualization.color_maps import *
+
+import numexpr as ne
+
+window = None # Number of the glut window.
+rot_enabled = True
+
+#Theia Scene
+ts = None
+
+#RAY CASTING values
+c_tbrightness = 1.0
+c_tdensity = 0.05
+
+output_texture = None # pointer to offscreen render target
+
+leftButton = False
+middleButton = False
+rightButton = False
+
+#Screen width and height
+width = 1024
+height = 1024
+
+eyesep = 0.1
+
+(pbo, pycuda_pbo) = [None]*2
+
+def create_PBO(w, h):
+ global pbo, pycuda_pbo
+ num_texels = w*h
+ array = np.zeros((w,h,3),np.uint32)
+
+ pbo = glGenBuffers(1)
+ glBindBuffer(GL_ARRAY_BUFFER, pbo)
+ glBufferData(GL_ARRAY_BUFFER, array, GL_DYNAMIC_DRAW)
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pycuda_pbo = cuda_gl.RegisteredBuffer(long(pbo))
+
+def destroy_PBO(self):
+ global pbo, pycuda_pbo
+ glBindBuffer(GL_ARRAY_BUFFER, long(pbo))
+ glDeleteBuffers(1, long(pbo));
+ glBindBuffer(GL_ARRAY_BUFFER, 0)
+ pbo,pycuda_pbo = [None]*2
+
+#consistent with C initPixelBuffer()
+def create_texture(w,h):
+ global output_texture
+ output_texture = glGenTextures(1)
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+ # set basic parameters
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
+ glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST)
+ # buffer data
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ w, h, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8, None)
+
+#consistent with C initPixelBuffer()
+def destroy_texture():
+ global output_texture
+ glDeleteTextures(output_texture);
+ output_texture = None
+
+def init_gl(w = 512 , h = 512):
+ Width, Height = (w, h)
+
+ glClearColor(0.1, 0.1, 0.5, 1.0)
+ glDisable(GL_DEPTH_TEST)
+
+ #matrix functions
+ glViewport(0, 0, Width, Height)
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+
+ #matrix functions
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+ glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
+
+def resize(Width, Height):
+ global width, height
+ (width, height) = Width, Height
+ glViewport(0, 0, Width, Height) # Reset The Current Viewport And Perspective Transformation
+ glMatrixMode(GL_PROJECTION)
+ glLoadIdentity()
+ gluPerspective(60.0, Width/float(Height), 0.1, 10.0)
+
+
+def do_tick():
+ global time_of_last_titleupdate, frame_counter, frames_per_second
+ if ((time.clock () * 1000.0) - time_of_last_titleupdate >= 1000.):
+ frames_per_second = frame_counter # Save The FPS
+ frame_counter = 0 # Reset The FPS Counter
+ szTitle = "%d FPS" % (frames_per_second )
+ glutSetWindowTitle ( szTitle )
+ time_of_last_titleupdate = time.clock () * 1000.0
+ frame_counter += 1
+
+oldMousePos = [ 0, 0 ]
+def mouseButton( button, mode, x, y ):
+ """Callback function (mouse button pressed or released).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively"""
+
+ global leftButton, middleButton, rightButton, oldMousePos
+
+ if button == GLUT_LEFT_BUTTON:
+ if mode == GLUT_DOWN:
+ leftButton = True
+ else:
+ leftButton = False
+
+ if button == GLUT_MIDDLE_BUTTON:
+ if mode == GLUT_DOWN:
+ middleButton = True
+ else:
+ middleButton = False
+
+ if button == GLUT_RIGHT_BUTTON:
+ if mode == GLUT_DOWN:
+ rightButton = True
+ else:
+ rightButton = False
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def mouseMotion( x, y ):
+ """Callback function (mouse moved while button is pressed).
+
+ The current and old mouse positions are stored in
+ a global renderParam and a global list respectively.
+ The global translation vector is updated according to
+ the movement of the mouse pointer."""
+
+ global ts, leftButton, middleButton, rightButton, oldMousePos
+ deltaX = x - oldMousePos[ 0 ]
+ deltaY = y - oldMousePos[ 1 ]
+
+ factor = 0.001
+
+ if leftButton == True:
+ ts.camera.rotateX( - deltaY * factor)
+ ts.camera.rotateY( - deltaX * factor)
+ if middleButton == True:
+ ts.camera.translateX( deltaX* 2.0 * factor)
+ ts.camera.translateY( - deltaY* 2.0 * factor)
+ if rightButton == True:
+ ts.camera.scale += deltaY * factor
+
+ oldMousePos[0], oldMousePos[1] = x, y
+ glutPostRedisplay( )
+
+def keyPressed(*args):
+ global c_tbrightness, c_tdensity
+ # If escape is pressed, kill everything.
+ if args[0] == '\033':
+ print 'Closing..'
+ destroy_PBOs()
+ destroy_texture()
+ exit()
+
+ #change the brightness of the scene
+ elif args[0] == ']':
+ c_tbrightness += 0.025
+ elif args[0] == '[':
+ c_tbrightness -= 0.025
+
+ #change the density scale
+ elif args[0] == ';':
+ c_tdensity -= 0.001
+ elif args[0] == '\'':
+ c_tdensity += 0.001
+
+def idle():
+ glutPostRedisplay()
+
+def display():
+ try:
+ #process left eye
+ process_image()
+ display_image()
+
+ glutSwapBuffers()
+
+ except:
+ from traceback import print_exc
+ print_exc()
+ from os import _exit
+ _exit(0)
+
+def process(eye = True):
+ global ts, pycuda_pbo, eyesep, c_tbrightness, c_tdensity
+
+ ts.get_raycaster().set_opacity(c_tdensity)
+ ts.get_raycaster().set_brightness(c_tbrightness)
+
+ dest_mapping = pycuda_pbo.map()
+ (dev_ptr, size) = dest_mapping.device_ptr_and_size()
+ ts.get_raycaster().surface.device_ptr = dev_ptr
+ ts.update()
+ # ts.get_raycaster().cast()
+ dest_mapping.unmap()
+
+
+def process_image():
+ global output_texture, pbo, width, height
+ """ copy image and process using CUDA """
+ # run the Cuda kernel
+ process()
+ # download texture from PBO
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, np.uint64(pbo))
+ glBindTexture(GL_TEXTURE_2D, output_texture)
+
+ glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA,
+ width, height, 0, GL_RGBA, GL_UNSIGNED_INT_8_8_8_8_REV, None)
+
+def display_image(eye = True):
+ global width, height
+ """ render a screen sized quad """
+ glDisable(GL_DEPTH_TEST)
+ glDisable(GL_LIGHTING)
+ glEnable(GL_TEXTURE_2D)
+ glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE)
+
+ #matix functions should be moved
+ glMatrixMode(GL_PROJECTION)
+ glPushMatrix()
+ glLoadIdentity()
+ glOrtho(-1.0, 1.0, -1.0, 1.0, -1.0, 1.0)
+ glMatrixMode( GL_MODELVIEW)
+ glLoadIdentity()
+ glViewport(0, 0, width, height)
+
+ glBegin(GL_QUADS)
+ glTexCoord2f(0.0, 0.0)
+ glVertex3f(-1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 0.0)
+ glVertex3f(1.0, -1.0, 0.5)
+ glTexCoord2f(1.0, 1.0)
+ glVertex3f(1.0, 1.0, 0.5)
+ glTexCoord2f(0.0, 1.0)
+ glVertex3f(-1.0, 1.0, 0.5)
+ glEnd()
+
+ glMatrixMode(GL_PROJECTION)
+ glPopMatrix()
+
+ glDisable(GL_TEXTURE_2D)
+ glBindTexture(GL_TEXTURE_2D, 0)
+ glBindBuffer(GL_PIXEL_PACK_BUFFER, 0)
+ glBindBuffer(GL_PIXEL_UNPACK_BUFFER, 0)
+
+
+#note we may need to init cuda_gl here and pass it to camera
+def main():
+ global window, ts, width, height
+ (width, height) = (1024, 1024)
+
+ glutInit(sys.argv)
+ glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH )
+ glutInitWindowSize(width, height)
+ glutInitWindowPosition(0, 0)
+ window = glutCreateWindow("Stereo Volume Rendering")
+
+
+ glutDisplayFunc(display)
+ glutIdleFunc(idle)
+ glutReshapeFunc(resize)
+ glutMouseFunc( mouseButton )
+ glutMotionFunc( mouseMotion )
+ glutKeyboardFunc(keyPressed)
+ init_gl(width, height)
+
+ # create texture for blitting to screen
+ create_texture(width, height)
+
+ import pycuda.gl.autoinit
+ import pycuda.gl
+ cuda_gl = pycuda.gl
+
+ create_PBO(width, height)
+ # ----- Load and Set Volume Data -----
+
+ density_grid = np.load("/home/bogert/dd150_log_densities.npy")
+
+ mi, ma= 21.5, 24.5
+ bins = 5000
+ tf = ColorTransferFunction( (mi, ma), bins)
+ tf.map_to_colormap(mi, ma, colormap="algae", scale_func = scale_func)
+
+ ts = TheiaScene(volume = density_grid, raycaster = FrontToBackRaycaster(size = (width, height), tf = tf))
+
+ ts.get_raycaster().set_sample_size(0.01)
+ ts.get_raycaster().set_max_samples(5000)
+ ts.update()
+
+ glutMainLoop()
+
+def scale_func(v, mi, ma):
+ return np.minimum(1.0, np.abs((v)-ma)/np.abs(mi-ma) + 0.0)
+
+# Print message to console, and kick off the main to get it rolling.
+if __name__ == "__main__":
+ print "Hit ESC key to quit, 'a' to toggle animation, and 'e' to toggle cuda"
+ main()
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa doc/source/visualizing/volume_rendering.rst
--- a/doc/source/visualizing/volume_rendering.rst
+++ b/doc/source/visualizing/volume_rendering.rst
@@ -466,3 +466,90 @@
your homogenized volume to then be passed in to the camera. A sample usage is shown
in :ref:`cookbook-amrkdtree_downsampling`.
+Hardware Volume Rendering on NVidia Graphics cards
+--------------------------------------------------
+.. versionadded:: 3.0
+
+Theia is a hardware volume renderer that takes advantage of NVidias CUDA language
+to peform ray casting with GPUs instead of the CPU.
+
+Only unigrid rendering is supported, but yt provides a grid mapping function
+ to get unigrid data from amr or sph formats :
+ :ref:`cookbook-amrkdtree_to_uniformgrid`.
+
+System Requirements
+-------------------
+.. versionadded:: 3.0
+
+Nvidia graphics card - The memory limit of the graphics card sets the limit
+ on the size of the data source.
+
+CUDA 5 or later and
+
+The environment variable CUDA_SAMPLES must be set pointing to
+the common/inc samples shipped with CUDA. The following shows an example
+in bash with CUDA 5.5 installed in /usr/local :
+
+export CUDA_SAMPLES=/usr/local/cuda-5.5/samples/common/inc
+
+PyCUDA must also be installed to use Theia.
+
+PyCUDA can be installed following these instructions :
+
+ git clone --recursive http://git.tiker.net/trees/pycuda.git
+
+ python configure.py
+ python setup.py install
+
+
+Tutorial
+--------
+.. versionadded:: 3.0
+
+Currently rendering only works on uniform grids. Here is an example
+on a 1024 cube of float32 scalars.
+
+.. code-block:: python
+ from yt.visualization.volume_rendering.theia.scene import TheiaScene
+ from yt.visualization.volume_rendering.algorithms.front_to_back import FrontToBackRaycaster
+ import numpy as np
+
+ #load 3D numpy array of float32
+ volume = np.load("/home/bogert/log_densities_1024.npy")
+
+ scene = TheiaScene( volume = volume, raycaster = FrontToBackRaycaster() )
+
+ scene.camera.rotateX(1.0)
+ scene.update()
+
+ surface = scene.get_results()
+ #surface now contains an image array 2x2 int32 rbga values
+
+.. _the-theiascene-interface:
+
+The TheiaScene Interface
+--------------------
+.. versionadded:: 3.0
+
+A TheiaScene object has been created to provide a high level entry point for
+controlling the raycaster's view onto the data. The class
+:class:`~yt.visualization.volume_rendering.theia.TheiaScene` encapsulates
+ a Camera object and a TheiaSource that intern encapsulates
+a volume. The :class:`~yt.visualization.volume_rendering.theia.Camera`
+provides controls for rotating, translating, and zooming into the volume.
+Using the :class:`~yt.visualization.volume_rendering.theia.TheiaSource`
+automatically transfers the volume to the graphic's card texture memory.
+
+Example Cookbooks
+---------------
+
+OpenGL Example for interactive volume rendering:
+:ref:`cookbook-opengl_volume_rendering`.
+
+OpenGL Stereoscopic Example :
+.. warning:: Frame rate will suffer significantly from stereoscopic rendering.
+ ~2x slower since the volume must be rendered twice.
+:ref:`cookbook-opengl_stereo_volume_rendering`.
+
+Pseudo-Realtime video rendering with ffmpeg :
+:ref:`cookbook-ffmpeg_volume_rendering`.
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.cu
@@ -0,0 +1,99 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+#include <helpers/cuda.cu>
+
+//A 3d texture holding scalar values to be volume rendered
+texture<float , cudaTextureType3D, cudaReadModeElementType> volume;
+//A 1d texture holding color transfer function values to act on raycaster results
+texture<float4, cudaTextureType1D, cudaReadModeElementType> transfer;
+
+
+extern "C"
+__global__ void front_to_back(uint *buffer, float *modelviewmatrix, int buffer_w,
+ int buffer_h, int max_steps, float density,
+ float offset, float scale,
+ float brightness, float step_size
+ ) {
+ //Rays will terminate when opacity_threshold is reached
+ const float opacity_threshold = 0.95f;
+
+ //We don't use the far and near clipping information from the modelviewmatrix
+ float3x4 matrix = mat_array_to_3x4(modelviewmatrix);
+
+ //The X,Y coordinate of the surface (image) array
+ int x = BLOCK_X;
+ int y = BLOCK_Y;
+
+ if ((x >= buffer_w) || (y >= buffer_h)) return;
+
+ float u = COORD_STD(x, buffer_w);
+ float v = COORD_STD(y, buffer_h);
+
+ // calculate eye ray in world space
+ Ray eye_ray = make_ray_from_eye(matrix, u, v);
+
+ // find intersection with box
+ float tnear, tfar;
+ int hit = intersect_std_box(eye_ray, &tnear, &tfar);
+
+ // If the ray misses the box set the coloro to black
+ if (!hit) {
+ buffer[y*buffer_w + x] = 0;
+ return;
+ }
+
+ if (tnear < 0.0f) tnear = 0.0f; // clamp to near plane
+
+ // march along ray from front to back, accumulating color
+ float4 sum = make_float4(0.0f);
+ float t = tnear;
+ float3 pos = eye_ray.o + eye_ray.d*tnear;
+ float3 step = eye_ray.d*step_size;
+
+ // Cast Rays
+ float4 col = make_float4(0.0,0.0,0.0,0.0);
+ float top = (1.0/scale) + offset;
+ for (int i=0; i<=max_steps; i++) {
+ float sample = tex3D(volume, pos.x*0.5f+0.5f, pos.y*0.5f+0.5f, pos.z*0.5f+0.5f);
+
+
+ if (sample > offset) {
+ col = tex1D(transfer, (sample-offset)*scale);
+ }
+ else {
+ col = make_float4(0.0,0.0,0.0,0.0);
+ }
+
+ col.w *= density;
+
+ // pre-multiply alpha
+ col.x *= col.w;
+ col.y *= col.w;
+ col.z *= col.w;
+
+ // "over" operator for front-to-back blending
+ sum = sum + col*(1.0f - sum.w);
+
+ // exit early if opaque
+ if (sum.w > opacity_threshold)
+ break;
+
+ // Increment step size and position
+ t += step_size;
+ pos += step;
+
+ // If ray has cast too far, end
+ if (t > tfar) break;
+ }
+
+ sum *= brightness;
+
+ buffer[SCREEN_XY(x,y,buffer_w)] = rgbaFloatToInt(sum);
+}
+
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/algorithms/front_to_back.py
@@ -0,0 +1,283 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+import yt.visualization.volume_rendering.theia.helpers.cuda as cdh
+import numpy as np
+
+from yt.visualization.volume_rendering.theia.transfer.linear_transfer import LinearTransferFunction
+from yt.visualization.volume_rendering.theia.transfer.helper import *
+
+
+from yt.visualization.volume_rendering.theia.volumes.unigrid import Unigrid
+from yt.visualization.volume_rendering.theia.surfaces.array_surface import ArraySurface
+import os
+
+class FrontToBackRaycaster:
+ r"""
+ This is the python wrapper for the CUDA raycaster. This
+ raycaster can act on unigrid data only. Use yt to map
+ AMR data to unigrid prior to using this algorithm.
+
+
+ Parameters
+ ----------
+ matrix : 4x4 numpy.matriix
+ ModelViewMatrix defines view onto volume
+
+ surface: yt.visualization.volume_rendering.theia.surface.array_surface.ArraySurface
+ The surface to contain the image information from
+ the results of the raycasting
+
+ volume: 3D numpy array Float32
+ Scalar values to be volume rendered
+
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Transfer function defining how the raycasting results are
+ to be colored.
+
+ size : Tuple of 2 Floats
+ If no surface is specified creates an ArraySurface of the specified size.
+
+ """
+
+ def __init__(self, matrix = None, surface = None,
+ volume = None, tf = None,
+ size = (864, 480)) :
+
+ # kernel module
+ self.kernel_module = cdh.source_file(
+ self.base_directory(__file__, "front_to_back.cu"))
+
+ #kernel function definitions
+ self.cast_rays_identifier = \
+ self.kernel_module.get_function("front_to_back")
+ self.transfer_identifier = \
+ self.kernel_module.get_texref("transfer")
+ self.volume_identifier = \
+ self.kernel_module.get_texref("volume")
+
+ self.set_matrix(matrix)
+
+ if (surface == None) :
+ self.set_surface(ArraySurface(size = size))
+ else :
+ self.set_surface(surface)
+
+
+ self.volume = None
+ self.set_transfer(tf)
+
+ self.set_sample_size()
+ self.set_max_samples()
+ self.set_opacity()
+ self.set_brightness()
+
+ """
+ Envoke the ray caster to cast rays
+ """
+ def __call__(self):
+ self.cast()
+ """
+ Returns
+ ----------
+ A 2d numpy array containing the image of the
+ volumetric rendering
+ """
+ def get_surface(self) :
+ return self.surface.get_array()
+
+ """
+ Returns
+ ----------
+ The sample size the ray caster is
+ configured to take.
+ """
+ def get_sample_size(self) :
+ return self.sample_size
+
+ """
+ Returns
+ ----------
+ The Max number of samples per ray the ray caster is
+ configured to take.
+ """
+ def get_max_samples(self) :
+ return self.max_samples
+
+ """
+ Returns
+ ----------
+ The Global density scalar.
+ """
+ def get_opacity(self) :
+ return self.opacity
+
+ """
+ Returns
+ ----------
+ The Global brightness scalar.
+
+ """
+ def get_brightness(self):
+ return self.brightness
+
+ """
+ Parameters
+ ----------
+ size : scalra Float
+ The distance between each sample, a smaller size will result in
+ more samples and can cause performance loss.
+ """
+ def set_sample_size(self, size = 0.01) :
+ self.sample_size = size
+
+ """
+ Parameters
+ ----------
+ max : scalar Int
+ The limit on the number of samples the ray caster will
+ take per ray.
+ """
+ def set_max_samples(self, max = 5000) :
+ self.max_samples = max
+
+ """
+ Parameters
+ ----------
+ scale : scalar Float
+ Global multiplier on volume data
+ """
+ def set_opacity(self, scale = 0.05) :
+ self.opacity = scale
+
+ """
+ Parameters
+ ----------
+ brightness : scalar Float
+ Global multiplier on returned alpha values
+ """
+ def set_brightness(self, brightness = 1.0):
+ self.brightness = brightness
+
+ """
+ Causes the ray caster to act on the volume data
+ and puts the results in the surface array.
+ """
+ def cast(self):
+ w, h = self.surface.bounds
+
+ self.cast_rays_identifier(np.uint64(self.surface.device_ptr),
+ drv.In(self.matrix),
+ np.int32(w), np.int32(h),
+ np.int32(self.max_samples),
+ np.float32(self.opacity),
+ np.float32(self.transfer.offset),
+ np.float32(self.transfer.scale),
+ np.float32(self.brightness),
+ np.float32(self.sample_size),
+ block=self.block,
+ grid=self.grid
+ )
+
+
+ """
+ Set the ModelViewMatrix, this does not send data to the GPU
+ Parameters
+ ----------
+ matrix : 4x4 numpy.matrix
+ ModelViewMatrix
+
+ """
+ def set_matrix(self, matrix):
+ self.matrix = matrix
+
+ """
+ Setup the image array for ray casting results
+ Parameters
+ ----------
+ surface : yt.visualization..volume_rendering.theia.surfaces.ArraySurface
+ Surface to contain results of the ray casting
+ """
+ def set_surface(self, surface = None, block_size = 32):
+ if surface == None:
+ self.surface = None
+ return
+
+ self.surface = surface
+ self.grid = cdh.block_estimate(block_size, self.surface.bounds)
+ self.block = (block_size, block_size, 1)
+
+ """
+ This function will convert a 3d numpy array into a unigrid volume
+ and move the data to the gpu texture memory
+ Parameters
+ ----------
+ volume : 3D numpy array float32
+ Contains scalar volumes to be acted on by the ray caster
+
+ """
+ def send_volume_to_gpu(self, volume = None) :
+ if (volume != None) :
+ self.set_volume(Unigrid(array = volume, allocate = True))
+
+ """
+ Parameters
+ ----------
+ volume : yt.visualization..volume_rendering.theia.volumes.Unigrid
+ Contains scalar volumes to be acted on by the ray caster
+
+ """
+ def set_volume(self, volume = None):
+ if volume == None:
+ self.volume = None
+ return
+
+ self.volume = volume
+ self.volume_identifier.set_flags(self.volume.flags)
+ self.volume_identifier.set_filter_mode(self.volume.filter_mode)
+ self.volume_identifier.set_address_mode(0, self.volume.address_mode)
+ self.volume_identifier.set_address_mode(1, self.volume.address_mode)
+ self.volume_identifier.set_array(self.volume.cuda_volume_array)
+
+ """
+ Parameters
+ ----------
+ tf : yt.visualization.volume_rendering.transfer_functions.ColorTransferFunction
+ Used to color the results of the raycasting
+
+ """
+ def set_transfer(self, tf = None):
+ if tf == None:
+ self.transfer = None
+ return
+
+ self.transfer = LinearTransferFunction(arrays = yt_to_rgba_transfer(tf), range = tf.x_bounds)
+ self.transfer_identifier.set_flags(self.transfer.flags)
+ self.transfer_identifier.set_filter_mode(self.transfer.filter_mode)
+ self.transfer_identifier.set_address_mode(0, self.transfer.address_mode)
+ self.transfer_identifier.set_array(self.transfer.cuda_transfer_array)
+
+ """
+ Attach the base directory path to the desired source file.
+ Parameters
+ ----------
+ dir : string
+ Directory where source file is located
+
+ file : string
+ Source file name
+
+
+ """
+ def base_directory(self, dir, file):
+ base = os.path.dirname(dir)
+ src = os.path.join(base, file)
+ return src
+
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/camera.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/camera.py
@@ -0,0 +1,139 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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.
+#-----------------------------------------------------------------------------
+
+#Matix helpers for camera control
+import yt.visualization.volume_rendering.theia.helpers.matrix_transformation as mth
+
+import numpy as np
+
+class Camera:
+ r"""
+ This is a basic camera intended to be a base class.
+ The camera provies utilities for controling the view point onto the data.
+
+
+ Parameters
+ ----------
+ pos : numpy array
+ The 3d position vector of the camera
+ rot : numpy array
+ The 3d rotation vector of the camera
+ scale : float
+ The scalar acting on camera position to zoom
+
+ Example:
+
+ from yt.visualization.volume_rendering.cameras.camera import Camera
+
+ cam = Camera()
+
+ cam.zoom(10.0)
+ cam.rotateX(np.pi)
+ cam.rotateY(np.pi)
+ cam.rotateZ(np.pi)
+ cam.translateX(-0.5)
+ cam.translateY(-0.5)
+ cam.translateZ(-0.5)
+
+ view = cam.get_matrix()
+
+ """
+ def __init__(self, rot = np.array([0.0, 0.0, 0.0]), scale = 1.0,
+ pos = np.array([0.0, 0.0, 0.0]), ):
+ #Values to control camera
+ self.rot = rot
+ self.scale = scale
+ self.pos = pos
+
+ #TODO : users should have control over perspective frustrum
+ self.stdmatrix = mth.perspective( [0.0, 0.0, 0.25] )
+
+
+ def zoom(self, percentage = 10.0):
+ r"""This will zoom the camera by percentage specified.
+
+ Parameters
+ ----------
+ percentage : float
+ If percentage is postive the camera will zoom in, if negative
+ then the camera will zoom out. Cannot exceed 100%.
+ """
+ if (percentage > 100.0) :
+ percentage = 100.0
+ self.scale = ((100.0 - percentage)/100.0)*self.scale
+
+ def rotateX(self, degrees = 0.1):
+ r"""This will rotate the camera around its X axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[0] += degrees
+
+ def rotateY(self, degrees = 0.1):
+ r"""This will rotate the camera around its Y axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[1] += degrees
+
+ def rotateZ(self, degrees = 0.1):
+ r"""This will rotate the camera around its Z axis.
+
+ Parameters
+ ----------
+ degrees : float
+ The amount of rotation in specified in radians.
+ """
+ self.rot[2] += degrees
+
+ def translateX(self, dist = 0.1):
+ r"""This will move the camera's position along its X axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[0] += dist
+
+ def translateY(self, dist = 0.1):
+ r"""This will move the camera's position along its Y axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[1] += dist
+
+ def translateZ(self, dist = 0.1):
+ r"""This will move the camera's position along its Z axis.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ self.pos[2] += dist
+
+ def get_matrix(self):
+ r"""This will calculate the curent modelview matrix
+ from the camera's position, rotation, and zoom scale.
+
+ Parameters
+ ----------
+ dist : float
+ The amount of movement. **This unit is unknown!**
+ """
+ return mth.rotate_x(self.rot[0]) * mth.rotate_y(self.rot[1]) * mth.rotate_z(self.rot[2]) * mth.scale((self.scale, self.scale, self.scale)) * mth.perspective([0.0, 0.0, 0.25]) *mth.translate((self.pos[0], self.pos[1], self.pos[2]))
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/helpers/cuda.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.cu
@@ -0,0 +1,189 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+
+#include <helper_cuda.h>
+#include <helper_math.h>
+
+#define BLOCK_X (blockIdx.x*blockDim.x + threadIdx.x)
+#define BLOCK_Y (blockIdx.y*blockDim.y + threadIdx.y)
+#define SCREEN_XY(a, b, w) (b*w+a)
+#define COORD_STD(x, w) ((x / (float) w)*2.0f-1.0f)
+
+#define COLOR_BLACK make_float3(0.0f, 0.0f, 0.0f)
+
+typedef struct {
+ float4 m[3];
+} float3x4;
+
+typedef struct {
+ float4 m[4];
+} float4x4;
+
+
+struct Ray {
+ float3 o; // origin
+ float3 d; // direction
+};
+
+
+
+__device__ float4 mul(const float3x4 &M, const float4 &v);
+__device__ float3 mul(const float3x4 &M, const float3 &v);
+__device__ float3x4 mat_array_to_3x4(float *m);
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v);
+__device__ int intersect_box(Ray r, float3 boxmin, float3 boxmax,
+ float *tnear, float *tfar);
+__device__ int intersect_std_box(Ray r, float *tnear, float *tfar);
+
+__device__ float mag(const float3 &v);
+__device__ float avg(const float3 &v);
+__device__ uint rgbaFloatToInt(float4 rgba);
+__device__ uint rgbFloatToInt(float3 rgb);
+__device__ float4 intToRGBAFloat(uint irgba);
+__device__ float3 intToRGBFloat(uint irgba);
+
+
+
+__device__ uint rgbFloatToInt(float3 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ return (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ uint rgbaFloatToInt(float4 rgba)
+{
+ rgba.x = __saturatef(rgba.x); // clamp to [0.0, 1.0]
+ rgba.y = __saturatef(rgba.y);
+ rgba.z = __saturatef(rgba.z);
+ rgba.w = __saturatef(rgba.w);
+ return (uint(rgba.w*255)<<24) | (uint(rgba.z*255)<<16) | (uint(rgba.y*255)<<8) | uint(rgba.x*255);
+}
+
+__device__ float4 intToRGBAFloat(uint irgba)
+{
+ float4 rgba;
+ rgba.w = (irgba >> 24 & 0xFF)/255;
+ rgba.z = (irgba >> 16 & 0xFF)/255;
+ rgba.y = (irgba >> 8 & 0xFF)/255;
+ rgba.x = (irgba & 0xFF)/255;
+ return rgba;
+}
+
+__device__ float3 intToRGBFloat(uint irgba)
+{
+ float3 rgb;
+ rgb.z = (irgba >> 16 & 0xFF)/255;
+ rgb.y = (irgba >> 8 & 0xFF)/255;
+ rgb.x = (irgba & 0xFF)/255;
+ return rgb;
+}
+
+
+__device__ float3x4 mat_array_to_3x4(float *m) {
+ float3x4 matrix;
+ matrix.m[0] = make_float4(m[0], m[1], m[2], m[3]);
+ matrix.m[1] = make_float4(m[4], m[5], m[6], m[7]);
+ matrix.m[2] = make_float4(m[8], m[9], m[10], m[11]);
+ return matrix;
+}
+
+__device__ Ray make_ray_from_eye(float3x4 m, float u, float v) {
+ Ray eyeRay;
+ eyeRay.o = make_float3(mul(m, make_float4(0.0f, 0.0f, 0.0f, 1.0f)));
+ eyeRay.d = normalize(make_float3(u, v, -2.0f));
+ eyeRay.d = mul(m, eyeRay.d);
+ return eyeRay;
+}
+
+__device__
+int intersect_box(Ray r, float3 boxmin, float3 boxmax, float *tnear, float *tfar)
+{
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+int intersect_std_box(Ray r, float *tnear, float *tfar)
+{
+ const float3 boxmin = make_float3(-1.0f, -1.0f, -1.0f);
+ const float3 boxmax = make_float3( 1.0f, 1.0f, 1.0f);
+
+ // compute intersection of ray with all six bbox planes
+ float3 invR = make_float3(1.0f) / r.d;
+ float3 tbot = invR * (boxmin - r.o);
+ float3 ttop = invR * (boxmax - r.o);
+
+ // re-order intersections to find smallest and largest on each axis
+ float3 tmin = fminf(ttop, tbot);
+ float3 tmax = fmaxf(ttop, tbot);
+
+ // find the largest tmin and the smallest tmax
+ float largest_tmin = fmaxf(fmaxf(tmin.x, tmin.y), fmaxf(tmin.x, tmin.z));
+ float smallest_tmax = fminf(fminf(tmax.x, tmax.y), fminf(tmax.x, tmax.z));
+
+ *tnear = largest_tmin;
+ *tfar = smallest_tmax;
+
+ return smallest_tmax > largest_tmin;
+}
+
+__device__
+float4 mul(const float3x4 &M, const float4 &v)
+{
+ float4 r;
+ r.x = dot(v, M.m[0]);
+ r.y = dot(v, M.m[1]);
+ r.z = dot(v, M.m[2]);
+ r.w = 1.0f;
+ return r;
+}
+
+// transform vector by matrix with translation
+
+__device__
+float3 mul(const float3x4 &M, const float3 &v)
+{
+ float3 r;
+ r.x = dot(v, make_float3(M.m[0]));
+ r.y = dot(v, make_float3(M.m[1]));
+ r.z = dot(v, make_float3(M.m[2]));
+ return r;
+}
+
+
+__device__
+float mag(const float3 &v) {
+ return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+}
+
+
+__device__
+float avg(const float3 &v) {
+ return (v.x + v.y + v.z) / 3.0;
+}
+
+
+
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/helpers/cuda.py
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda.py
@@ -0,0 +1,58 @@
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, 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 pycuda.driver as drv
+from pycuda.compiler import SourceModule
+import pycuda.gpuarray as gpuarray
+import pycuda.autoinit
+import math
+import os
+
+def sample_path():
+ #return '/pfs/sw/cuda-5.5/samples/common/inc'
+ return os.environ.get('CUDA_SAMPLES')
+
+def cuda_helpers_path():
+ return os.path.dirname(__file__) + "/../"
+
+def source_file(path):
+ return SourceModule(open(path).read(),
+ include_dirs=[sample_path(),cuda_helpers_path()],
+ no_extern_c=True, keep=True)
+
+def block_estimate(thread_size, shape):
+ (x, y) = shape
+ return (int(math.ceil(float(x) / thread_size)), int(math.ceil(float(y) / thread_size)), 1)
+
+def np3d_to_device_array(np_array, allow_surface_bind=True):
+ d, h, w = np_array.shape
+
+ descr = drv.ArrayDescriptor3D()
+ descr.width = w
+ descr.height = h
+ descr.depth = d
+ descr.format = drv.dtype_to_array_format(np_array.dtype)
+ descr.num_channels = 1
+ descr.flags = 0
+
+ if allow_surface_bind:
+ descr.flags = drv.array3d_flags.SURFACE_LDST
+
+ device_array = drv.Array(descr)
+
+ copy = drv.Memcpy3D()
+ copy.set_src_host(np_array)
+ copy.set_dst_array(device_array)
+ copy.width_in_bytes = copy.src_pitch = np_array.strides[1]
+ copy.src_height = copy.height = h
+ copy.depth = d
+
+ copy()
+
+ return device_array
diff -r ea182fc68faeecd6c8b189da3a145811f9cb150a -r b6b78bc7de4c9cce495f58cefaabcee2dfe3bbaa yt/visualization/volume_rendering/theia/helpers/cuda/matrix_helper.cu
--- /dev/null
+++ b/yt/visualization/volume_rendering/theia/helpers/cuda/matrix_helper.cu
@@ -0,0 +1,265 @@
+//-----------------------------------------------------------------------------
+// Copyright (c) 2014, 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.
+//-----------------------------------------------------------------------------
+
+
+#include <helper_cuda.h>
+#include <helper_math.h>
+
+
+typedef float3 vec3;
+typedef float4 vec4;
+typedef float2 vec2;
+
+
+typedef struct {
+ vec4 m[4];
+} mat4x4;
+
+typedef struct {
+ vec3 m[3];
+} mat3x3;
+
+typedef struct {
+ vec2 m[2];
+} mat2x2;
+
+typedef struct {
+ vec3 m[4];
+} mat3x4;
+
+typedef struct {
+ vec4 m[3];
+} mat4x3;
+
+typedef struct {
+ vec4 m[2];
+} mat4x2;
+
+typedef struct {
+ vec2 m[4];
+} mat2x4;
+
+typedef mat4 mat4x4;
+typedef mat3 mat3x3;
+typedef mat2 mat2x2;
+
+
+
+// Column Major Always
+// All vectors are k x 1, and cannot be at the front of a multiplication
+
+
+// -----------------------------------------------------------------------------
+// Create Matrices, identity by default
+// -----------------------------------------------------------------------------
+__device__ __host__ mat4 make_mat4() {
+ mat4 tmp;
+ tmp.m[0] = make_float4(1.0f, 0.0f, 0.0f, 0.0f);
+ tmp.m[1] = make_float4(0.0f, 1.0f, 0.0f, 0.0f);
+ tmp.m[2] = make_float4(0.0f, 0.0f, 1.0f, 0.0f);
+ tmp.m[3] = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
+ return tmp;
+}
+
+__device__ __host__ mat3 make_mat3() {
+ mat3 tmp;
+ tmp.m[0] = make_float3(1.0f, 0.0f, 0.0f);
+ tmp.m[1] = make_float3(0.0f, 1.0f, 0.0f);
+ tmp.m[2] = make_float3(0.0f, 0.0f, 1.0f);
+ return tmp;
+}
+
+__device__ __host__ mat2 make_mat2() {
+ mat3 tmp;
+ tmp.m[0] = make_float2(1.0f, 0.0f);
+ tmp.m[1] = make_float2(0.0f, 1.0f);
+ return tmp;
+}
+
+
+// -----------------------------------------------------------------------------
+// Get Transpose vectors, get column vector
+// -----------------------------------------------------------------------------
+// 4s
+__device__ __host__ vec4 mvec(mat4x4 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ else if(col == 2) tmp = make_float4(a.m[0].z, a.m[1].z, a.m[2].z, a.m[3].z);
+ else if(col == 3) tmp = make_float4(a.m[0].w, a.m[1].w, a.m[2].w, a.m[3].w);
+ return tmp;
+}
+
+__device__ __host__ vec3 mvec(mat4x3 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ else if(col == 2) tmp = make_float3(a.m[0].z, a.m[1].z, a.m[2].z);
+ else if(col == 3) tmp = make_float3(a.m[0].w, a.m[1].w, a.m[2].w);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat4x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ else if(col == 2) tmp = make_float2(a.m[0].z, a.m[1].z);
+ else if(col == 3) tmp = make_float2(a.m[0].w, a.m[1].w);
+ return tmp;
+}
+// 3s
+
+__device__ __host__ vec3 mvec(mat3x4 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ else if(col == 2) tmp = make_float4(a.m[0].z, a.m[1].z, a.m[2].z, a.m[3].z);
+ return tmp;
+}
+
+__device__ __host__ vec4 mvec(mat3x3 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ else if(col == 2) tmp = make_float3(a.m[0].z, a.m[1].z, a.m[2].z);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat3x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ else if(col == 2) tmp = make_float2(a.m[0].z, a.m[1].z);
+ return tmp;
+}
+// 2s
+
+__device__ __host__ vec2 mvec(mat2x4 a, int col) {
+ vec4 tmp;
+ if (col == 0) tmp = make_float4(a.m[0].x, a.m[1].x, a.m[2].x, a.m[3].x);
+ else if(col == 1) tmp = make_float4(a.m[0].y, a.m[1].y, a.m[2].y, a.m[3].y);
+ return tmp;
+}
+
+__device__ __host__ vec3 mvec(mat2x3 a, int col) {
+ vec3 tmp;
+ if (col == 0) tmp = make_float3(a.m[0].x, a.m[1].x, a.m[2].x);
+ else if(col == 1) tmp = make_float3(a.m[0].y, a.m[1].y, a.m[2].y);
+ return tmp;
+}
+
+__device__ __host__ vec2 mvec(mat2x2 a, int col) {
+ vec2 tmp;
+ if (col == 0) tmp = make_float2(a.m[0].x, a.m[1].x);
+ else if(col == 1) tmp = make_float2(a.m[0].y, a.m[1].y);
+ return tmp;
+}
+
+
+// -----------------------------------------------------------------------------
+// Transpose Matrix
+// Cannot transpose anything other than square matrices
+// -----------------------------------------------------------------------------
+__device__ __host__ mat4 transpose(mat4 a) {
+ mat4 tmp;
+ for(int i = 0; i < 4; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+__device__ __host__ mat3 transpose(mat3 a) {
+ mat3 tmp;
+ for(int i = 0; i < 3; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+__device__ __host__ mat2 transpose(mat2 a) {
+ mat2 tmp;
+ for(int i = 0; i < 2; i++)
+ tmp.m[i] = mvec(a, i);
+ return tmp
+}
+
+// -----------------------------------------------------------------------------
+// Dot Products
+// -----------------------------------------------------------------------------
+inline __device__ __host__ float dot(vec4 a, vec4 b) {
+ return (a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w);
+}
+
+inline __device__ __host__ float dot(vec3 a, vec3 b) {
+ return (a.x*b.x + a.y*b.y + a.z*b.z);
+}
+
+inline __device__ __host__ float dot(vec2 a, vec2 b) {
+ return (a.x*b.x + a.y*b.y);
+}
+
+// -----------------------------------------------------------------------------
+// Scale Matrix
+// -----------------------------------------------------------------------------
+
+
+// -----------------------------------------------------------------------------
+// Multiply Matrix
+// -----------------------------------------------------------------------------
+//4s
+__device__ __host__ mat4x4 operator*(mat4x4 a, mat4x4 b) {
+ mat4x4 tmp;
+ for(int i = 0; i < 4; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ mat4x3 operator*(mat4x4 a, mat4x3 b) {
+ mat4x3 tmp;
+ for(int i = 0; i < 3; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ mat4x2 operator*(mat4x4 a, mat4x2 b) {
+ mat4x2 tmp;
+ for(int i = 0; i < 2; i++){
+ tmp.m[i].x = dot(a.m[i], mvec(b, 0));
+ tmp.m[i].y = dot(a.m[i], mvec(b, 1));
+ tmp.m[i].z = dot(a.m[i], mvec(b, 2));
+ tmp.m[i].w = dot(a.m[i], mvec(b, 3));
+ }
+ return tmp;
+}
+
+__device__ __host__ vec4 operator*(mat4x4 a, vec4 b) {
+ vec4 tmp;
+ tmp.x = dot(a.m[0], b);
+ tmp.y = dot(a.m[1], b);
+ tmp.z = dot(a.m[2], b);
+ tmp.w = dot(a.m[3], b);
+ return tmp;
+}
+
+__device__ __host__ mat3x4 operator*(mat3x3 a, mat3x4 b);
+__device__ __host__ mat3x3 operator*(mat3x3 a, mat3x3 b);
+__device__ __host__ mat3x2 operator*(mat3x3 a, mat3x2 b);
+__device__ __host__ vec3 operator*(mat3x3 a, vec3 b);
+
+__device__ __host__ mat2x4 operator*(mat2x2 a, mat2x4 b);
+__device__ __host__ mat2x3 operator*(mat2x2 a, mat2x3 b);
+__device__ __host__ mat2x2 operator*(mat2x2 a, mat2x2 b);
+__device__ __host__ vec2 operator*(mat2x2 a, vec2 b);
+
+
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