[yt-svn] commit/yt: MatthewTurk: Merged in fbogert/yt-theia/yt-3.0 (pull request #862)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sun Jun 15 05:45:59 PDT 2014
1 new commit in yt:
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