[Yt-svn] yt: Adding camera module and a stripped down, very basic homogen...
hg at spacepope.org
hg at spacepope.org
Fri May 14 16:51:00 PDT 2010
hg Repository: yt
details: yt/rev/1442a00b97e8
changeset: 1670:1442a00b97e8
user: Matthew Turk <matthewturk at gmail.com>
date:
Fri May 14 16:50:53 2010 -0700
description:
Adding camera module and a stripped down, very basic homogenized volume object.
diffstat:
yt/extensions/volume_rendering/__init__.py | 4 +-
yt/extensions/volume_rendering/camera.py | 119 +++++++++++++++++++++++
yt/extensions/volume_rendering/grid_partitioner.py | 86 +++++++++++++++++
3 files changed, 208 insertions(+), 1 deletions(-)
diffs (235 lines):
diff -r fe498f728dae -r 1442a00b97e8 yt/extensions/volume_rendering/__init__.py
--- a/yt/extensions/volume_rendering/__init__.py Wed May 12 23:55:20 2010 -0700
+++ b/yt/extensions/volume_rendering/__init__.py Fri May 14 16:50:53 2010 -0700
@@ -31,9 +31,11 @@
ProjectionTransferFunction
from yt.amr_utils import PartitionedGrid, VectorPlane, \
TransferFunctionProxy
-from grid_partitioner import HomogenizedBrickCollection, \
+from grid_partitioner import HomogenizedVolume, \
+ HomogenizedBrickCollection, \
export_partitioned_grids, \
import_partitioned_grids
from image_handling import export_rgba, import_rgba, \
plot_channel, plot_rgb
from software_sampler import VolumeRendering
+from camera import Camera
diff -r fe498f728dae -r 1442a00b97e8 yt/extensions/volume_rendering/camera.py
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/yt/extensions/volume_rendering/camera.py Fri May 14 16:50:53 2010 -0700
@@ -0,0 +1,119 @@
+"""
+Import the components of the volume rendering extension
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2009 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as na
+from grid_partitioner import HomogenizedVolume
+from TransferFunction import ProjectionTransferFunction
+from yt.funcs import *
+import yt.amr_utils as au
+
+class Camera(object):
+ def __init__(self, center, normal_vector, width,
+ resolution, transfer_function,
+ north_vector = None,
+ volume = None, fields = None,
+ log_fields = None,
+ sub_samples = 5, pf = None):
+ if pf is not None: self.pf = pf
+ if not iterable(resolution):
+ resolution = (resolution, resolution)
+ self.resolution = resolution
+ self.sub_samples = sub_samples
+ if not iterable(width):
+ width = (width, width, width) # front/back, left/right, top/bottom
+ self.width = width
+ self.center = center
+ if fields is None: fields = ["Density"]
+ self.fields = fields
+ if transfer_function is None:
+ transfer_function = ProjectionTransferFunction()
+ self.transfer_function = transfer_function
+ self._setup_normalized_vectors(normal_vector, north_vector)
+ if volume is None:
+ volume = HomogenizedVolume(fields, pf = self.pf,
+ log_fields = log_fields)
+ self.volume = volume
+
+ def _setup_normalized_vectors(self, normal_vector, north_vector):
+ # Now we set up our various vectors
+ normal_vector /= na.sqrt( na.dot(normal_vector, normal_vector))
+ if north_vector is None:
+ vecs = na.identity(3)
+ t = na.cross(normal_vector, vecs).sum(axis=1)
+ ax = t.argmax()
+ north_vector = na.cross(vecs[ax,:], normal_vector).ravel()
+ north_vector /= na.sqrt(na.dot(north_vector, north_vector))
+ east_vector = -na.cross(north_vector, normal_vector).ravel()
+ east_vector /= na.sqrt(na.dot(east_vector, east_vector))
+ self.unit_vectors = [north_vector, east_vector, normal_vector]
+ self.box_vectors = na.array([self.unit_vectors[0]*self.width[0],
+ self.unit_vectors[1]*self.width[1],
+ self.unit_vectors[2]*self.width[2]])
+
+ self.origin = self.center - 0.5*self.width[0]*self.unit_vectors[0] \
+ - 0.5*self.width[1]*self.unit_vectors[1] \
+ - 0.5*self.width[2]*self.unit_vectors[2]
+ self.back_center = self.center - 0.5*self.width[0]*self.unit_vectors[2]
+ self.front_center = self.center + 0.5*self.width[0]*self.unit_vectors[2]
+ self.inv_mat = na.linalg.pinv(self.unit_vectors)
+
+ def get_vector_plane(self, image):
+ # We should move away from pre-generation of vectors like this and into
+ # the usage of on-the-fly generation in the VolumeIntegrator module
+ # We might have a different width and back_center
+ px = na.linspace(-self.width[0]/2.0, self.width[0]/2.0,
+ self.resolution[0])[:,None]
+ py = na.linspace(-self.width[1]/2.0, self.width[1]/2.0,
+ self.resolution[1])[None,:]
+ inv_mat = self.inv_mat
+ bc = self.back_center
+ vectors = na.zeros((self.resolution[0], self.resolution[1], 3),
+ dtype='float64', order='C')
+ vectors[:,:,0] = inv_mat[0,0]*px+inv_mat[0,1]*py+self.back_center[0]
+ vectors[:,:,1] = inv_mat[1,0]*px+inv_mat[1,1]*py+self.back_center[1]
+ vectors[:,:,2] = inv_mat[2,0]*px+inv_mat[2,1]*py+self.back_center[2]
+ bounds = (px.min(), px.max(), py.min(), py.max())
+ vector_plane = au.VectorPlane(vectors, self.box_vectors[2],
+ self.back_center, bounds, image,
+ self.unit_vectors[1],
+ self.unit_vectors[0])
+ return vector_plane
+
+ def snapshot(self):
+ image = na.zeros((self.resolution[0], self.resolution[1], 3),
+ dtype='float64', order='C')
+ vector_plane = self.get_vector_plane(image)
+ tfp = au.TransferFunctionProxy(self.transfer_function) # Reset it every time
+ tfp.ns = self.sub_samples
+ self.volume.initialize_source()
+ pbar = get_pbar("Ray casting",
+ (self.volume.brick_dimensions + 1).prod(axis=-1).sum())
+ total_cells = 0
+ for brick in self.volume.traverse(self.back_center, self.front_center):
+ brick.cast_plane(tfp, vector_plane)
+ total_cells += na.prod(brick.my_data[0].shape)
+ pbar.update(total_cells)
+ pbar.finish()
+ return image
diff -r fe498f728dae -r 1442a00b97e8 yt/extensions/volume_rendering/grid_partitioner.py
--- a/yt/extensions/volume_rendering/grid_partitioner.py Wed May 12 23:55:20 2010 -0700
+++ b/yt/extensions/volume_rendering/grid_partitioner.py Fri May 14 16:50:53 2010 -0700
@@ -32,6 +32,92 @@
from yt.lagos import ParallelAnalysisInterface, only_on_root, parallel_root_only
from yt.parallel_tools import DistributedObjectCollection
+# HISTORICAL NOTE OF SOME IMPORT:
+# The homogenized brick class should and will be removed. The more general,
+# and less fancy-parallel, HomogenizedVolume is the way to go. The HBC is a
+# complicated mechanism that was designed for back when we were going to be
+# passing bricks all around between processors.
+
+class HomogenizedVolume(ParallelAnalysisInterface):
+ bricks = None
+ def __init__(self, fields = "Density", source = None, pf = None,
+ log_fields = None):
+ # Typically, initialized as hanging off a hierarchy. But, not always.
+ if pf is not None: self.pf = pf
+ if source is None: source = self.pf.h.all_data()
+ self.source = source
+ self.fields = ensure_list(fields)
+ if log_fields is not None:
+ log_fields = ensure_list(log_fields)
+ else:
+ log_fields = [self.pf.field_info[field].take_log
+ for field in self.fields]
+ self.log_fields = log_fields
+
+ def traverse(self, back_point, front_point):
+ if self.bricks is None: self.initialize_source()
+ vec = front_point - back_point
+ dist = na.minimum(
+ na.sum((self.brick_left_edges - back_point) * vec, axis=1),
+ na.sum((self.brick_right_edges - back_point) * vec, axis=1))
+ ind = na.argsort(dist)
+ for b in self.bricks[ind]: yield b
+
+ def _partition_grid(self, grid):
+
+ # This is not super efficient, as it re-fills the regions once for each
+ # field.
+ vcds = []
+ for field, log_field in zip(self.fields, self.log_fields):
+ vcd = grid.get_vertex_centered_data(field).astype('float64')
+ if log_field: vcd = na.log10(vcd)
+ vcds.append(vcd)
+
+ GF = GridFaces(grid.Children + [grid])
+ PP = ProtoPrism(grid.id, grid.LeftEdge, grid.RightEdge, GF)
+
+ pgs = []
+ for P in PP.sweep(0):
+ sl = P.get_brick(grid.LeftEdge, grid.dds, grid.child_mask)
+ if len(sl) == 0: continue
+ dd = [d[sl[0][0]:sl[0][1]+1,
+ sl[1][0]:sl[1][1]+1,
+ sl[2][0]:sl[2][1]+1].copy() for d in vcds]
+ pgs.append(PartitionedGrid(grid.id, len(self.fields), dd,
+ P.LeftEdge, P.RightEdge, sl[-1]))
+ return pgs
+
+ def initialize_source(self, source = None):
+ if self.bricks is not None and source is not None:
+ raise NotImplementedError("Sorry, dynamic shuffling of bricks is" +
+ " not yet supported")
+ if self.bricks is not None and source is None: return
+ bricks = []
+ self._preload(self.source._grids, self.fields, self.pf.h.io)
+ pbar = get_pbar("Partitioning ", len(self.source._grids))
+ for i, g in enumerate(self.source._grids):
+ pbar.update(i)
+ bricks += self._partition_grid(g)
+ pbar.finish()
+ bricks = na.array(bricks, dtype='object')
+ NB = len(bricks)
+ # Now we set up our (local for now) hierarchy. Note that to calculate
+ # intersection, we only need to do the left edge & right edge.
+ #
+ # We're going to double up a little bit here in memory.
+ self.brick_left_edges = na.zeros( (NB, 3), dtype='float64')
+ self.brick_right_edges = na.zeros( (NB, 3), dtype='float64')
+ self.brick_parents = na.zeros( NB, dtype='int64')
+ self.brick_dimensions = na.zeros( (NB, 3), dtype='int64')
+ for i,b in enumerate(bricks):
+ self.brick_left_edges[i,:] = b.LeftEdge
+ self.brick_right_edges[i,:] = b.RightEdge
+ self.brick_parents[i] = b.parent_grid_id
+ self.brick_dimensions[i,:] = b.my_data[0].shape
+ # Vertex-centered means we subtract one from the shape
+ self.brick_dimensions -= 1
+ self.bricks = na.array(bricks, dtype='object')
+
class HomogenizedBrickCollection(DistributedObjectCollection):
def __init__(self, source):
# The idea here is that we have two sources -- the global_domain
More information about the yt-svn
mailing list