[Yt-svn] yt: 3 new changesets
hg at spacepope.org
hg at spacepope.org
Tue Apr 6 16:44:20 PDT 2010
hg Repository: yt
details: yt/rev/1f8362113835
changeset: 1517:1f8362113835
user: Matthew Turk <matthewturk at gmail.com>
date:
Fri Apr 02 08:24:41 2010 -0700
description:
Typo, missed anorm as a used parameter
hg Repository: yt
details: yt/rev/a28d29a7491a
changeset: 1518:a28d29a7491a
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Apr 06 16:41:12 2010 -0700
description:
Adding "projload" convenience command -- which may go away ! -- to
convenience.py. This lets you load a projection affiliated with a parameter
file off disk without having to instantiate all the affiliated grid objects.
hg Repository: yt
details: yt/rev/515bb4db647d
changeset: 1519:515bb4db647d
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Apr 06 16:44:16 2010 -0700
description:
Merging
diffstat:
yt/_amr_utils/VolumeIntegrator.pyx | 251 +-
yt/amr_utils.c | 18107 +++++++++++-------------
yt/convenience.py | 25 +
yt/extensions/image_panner/pan_and_scan_widget.py | 9 +-
yt/extensions/image_panner/vm_panner.py | 4 +-
yt/extensions/volume_rendering/TransferFunction.py | 103 +-
yt/extensions/volume_rendering/UBVRI.py | 105 -
yt/extensions/volume_rendering/__init__.py | 4 +-
yt/extensions/volume_rendering/grid_partitioner.py | 28 +-
yt/extensions/volume_rendering/image_handling.py | 2 +-
yt/extensions/volume_rendering/software_sampler.py | 12 +-
yt/lagos/HaloFinding.py | 60 +-
yt/lagos/ParallelTools.py | 4 +-
yt/lagos/parallelHOP/parallelHOP.py | 27 +-
yt/mods.py | 2 +-
yt/physical_constants.py | 3 +-
16 files changed, 8666 insertions(+), 10080 deletions(-)
diffs (truncated from 29850 to 300 lines):
diff -r 3da92a5ea545 -r 515bb4db647d yt/_amr_utils/VolumeIntegrator.pyx
--- a/yt/_amr_utils/VolumeIntegrator.pyx Fri Apr 02 08:02:09 2010 -0700
+++ b/yt/_amr_utils/VolumeIntegrator.pyx Tue Apr 06 16:44:16 2010 -0700
@@ -73,147 +73,81 @@
cdef class VectorPlane
-cdef struct FieldInterpolationTable:
- # Note that we make an assumption about retaining a reference to values
- # externally.
- np.float64_t *values
- np.float64_t bounds[2]
- np.float64_t dbin
- np.float64_t idbin
- int field_id
- int weight_field_id
- int weight_table_id
- int nbins
+cdef class TransferFunctionProxy:
+ cdef np.float64_t x_bounds[2]
+ cdef np.float64_t *vs[4]
+ cdef int nbins
+ cdef public int ns
+ cdef np.float64_t dbin, idbin
+ cdef np.float64_t light_color[3]
+ cdef np.float64_t light_dir[3]
+ cdef int use_light
+ cdef public object tf_obj
+ def __cinit__(self, tf_obj):
+ self.tf_obj = tf_obj
+ cdef np.ndarray[np.float64_t, ndim=1] temp
+ temp = tf_obj.red.y
+ self.vs[0] = <np.float64_t *> temp.data
+ temp = tf_obj.green.y
+ self.vs[1] = <np.float64_t *> temp.data
+ temp = tf_obj.blue.y
+ self.vs[2] = <np.float64_t *> temp.data
+ temp = tf_obj.alpha.y
+ self.vs[3] = <np.float64_t *> temp.data
+ self.x_bounds[0] = tf_obj.x_bounds[0]
+ self.x_bounds[1] = tf_obj.x_bounds[1]
+ self.nbins = tf_obj.nbins
+ self.dbin = (self.x_bounds[1] - self.x_bounds[0])/self.nbins
+ self.idbin = 1.0/self.dbin
+ self.light_color[0] = tf_obj.light_color[0]
+ self.light_color[1] = tf_obj.light_color[1]
+ self.light_color[2] = tf_obj.light_color[2]
+ self.light_dir[0] = tf_obj.light_dir[0]
+ self.light_dir[1] = tf_obj.light_dir[1]
+ self.light_dir[2] = tf_obj.light_dir[2]
+ cdef np.float64_t normval = 0.0
+ for i in range(3): normval += self.light_dir[i]**2
+ normval = normval**0.5
+ for i in range(3): self.light_dir[i] /= normval
+ self.use_light = tf_obj.use_light
-cdef void FIT_initialize_table(FieldInterpolationTable *fit, int nbins,
- np.float64_t *values, np.float64_t bounds1, np.float64_t bounds2,
- int field_id, int weight_field_id = -1, int weight_table_id = -1):
- fit.bounds[0] = bounds1; fit.bounds[1] = bounds2
- fit.nbins = nbins
- fit.dbin = (fit.bounds[1] - fit.bounds[0])/fit.nbins
- fit.idbin = 1.0/fit.dbin
- # Better not pull this out from under us, yo
- fit.values = values
- fit.field_id = field_id
- fit.weight_field_id = weight_field_id
- fit.weight_table_id = weight_table_id
-
-cdef np.float64_t FIT_get_value(FieldInterpolationTable *fit,
- np.float64_t *dvs):
- cdef np.float64_t bv, dy, dd, tf
- cdef int bin_id
- bin_id = <int> ((dvs[fit.field_id] - fit.bounds[0]) * fit.idbin)
- dd = dvs[fit.field_id] - (fit.bounds[0] + bin_id * fit.dbin) # x - x0
- if bin_id > fit.nbins - 2 or bin_id < 0: return 0.0
- bv = fit.values[bin_id]
- dy = fit.values[bin_id + 1] - bv
- if fit.weight_field_id != -1:
- return dvs[fit.weight_field_id] * (bv + dd*dy*fit.idbin)
- return (bv + dd*dy*fit.idbin)
-
-cdef class TransferFunctionProxy:
- cdef int n_fields
- cdef int n_field_tables
- cdef public int ns
-
- # These are the field tables and their affiliated storage.
- # We have one field_id for every table. Note that a single field can
- # correspond to multiple tables, and each field table will only have
- # interpolate called once.
- cdef FieldInterpolationTable field_tables[6]
- cdef np.float64_t istorage[6]
-
- # Here are the field tables that correspond to each of the six channels.
- # We have three emission channels, three absorption channels.
- cdef int field_table_ids[6]
-
- # We store a reference to the transfer function object and to the field
- # interpolation tables
- cdef public object tf_obj
- cdef public object my_field_tables
-
- def __cinit__(self, tf_obj):
- # We have N fields. We have 6 channels. We have M field tables.
- # The idea is that we can have multiple channels corresponding to the
- # same field table. So, we create storage for the outputs from all the
- # field tables. We need to know which field value to pass in to the
- # field table, and we need to know which table to use for each of the
- # six channels.
- cdef int i
- cdef np.ndarray[np.float64_t, ndim=1] temp
- cdef FieldInterpolationTable fit
-
- self.tf_obj = tf_obj
-
- self.n_field_tables = tf_obj.n_field_tables
- for i in range(6): self.istorage[i] = 0.0
-
- self.my_field_tables = []
- for i in range(self.n_field_tables):
- temp = tf_obj.tables[i].y
- FIT_initialize_table(&self.field_tables[i],
- temp.shape[0],
- <np.float64_t *> temp.data,
- tf_obj.tables[i].x_bounds[0],
- tf_obj.tables[i].x_bounds[1],
- tf_obj.field_ids[i], tf_obj.weight_field_ids[i],
- tf_obj.weight_table_ids[i])
- self.my_field_tables.append((tf_obj.tables[i],
- tf_obj.tables[i].y))
- self.field_tables[i].field_id = tf_obj.field_ids[i]
- self.field_tables[i].weight_field_id = tf_obj.weight_field_ids[i]
- print "Field table", i, "corresponds to",
- print self.field_tables[i].field_id,
- print "(Weighted with ", self.field_tables[i].weight_field_id,
- print ")"
-
- for i in range(6):
- self.field_table_ids[i] = tf_obj.field_table_ids[i]
- print "Channel", i, "corresponds to", self.field_table_ids[i]
-
@cython.boundscheck(False)
@cython.wraparound(False)
- cdef void eval_transfer(self, np.float64_t dt, np.float64_t *dvs,
- np.float64_t *rgba, np.float64_t *grad):
- cdef int i, fid, use
- cdef np.float64_t ta, tf, trgba[6], dot_prod
- # This very quick pass doesn't hurt us too badly, and it helps for
- # early-cutoff. We check all the field tables, because we want to be
- # able to attenuate even in the presence of no emissivity.
- use = 0
- for i in range(self.n_field_tables):
- fid = self.field_tables[i].field_id
- if (dvs[fid] >= self.field_tables[i].bounds[0]) and \
- (dvs[fid] <= self.field_tables[i].bounds[1]):
- use = 1
- break
- for i in range(self.n_field_tables):
- self.istorage[i] = FIT_get_value(&self.field_tables[i], dvs)
- # We have to do this after the interpolation
- for i in range(self.n_field_tables):
- fid = self.field_tables[i].weight_table_id
- if fid != -1: self.istorage[i] *= self.istorage[fid]
- for i in range(6):
- trgba[i] = self.istorage[self.field_table_ids[i]]
- #print i, trgba[i],
- #print
- # A few words on opacity. We're going to be integrating equation 1.23
- # from Rybicki & Lightman. dI_\nu / ds = -\alpha_\nu I_\nu + j_\nu
- # \alpha_nu = \kappa \rho , but we leave that up to the input
- # transfer function.
- # SOoooooOOOooo, the upshot is that we are doing a rectangular
- # integration here:
- # I_{i+1} = ds * C_i + (1.0 - ds*alpha_i) * I_i
- for i in range(3):
- # This is the new way: alpha corresponds to opacity of a given
- # slice. Previously it was ill-defined, but represented some
- # measure of emissivity.
- ta = fmax((1.0 - dt*trgba[i+3]), 0.0)
- rgba[i ] = dt*trgba[i ] + ta * rgba[i ]
- rgba[i+3] = dt*trgba[i+3] + ta * rgba[i+3]
- # This is the old way:
- #rgba[i ] += trgba[i] * (1.0 - rgba[i+3])*dt*trgba[i+3]
- #rgba[i+3] += trgba[i] * (1.0 - rgba[i+3])*dt*trgba[i+3]
+ cdef void interpolate(self, np.float64_t dv, np.float64_t *trgba):
+ cdef int bin_id, channel
+ cdef np.float64_t bv, dy, dd, tf
+ bin_id = <int> ((dv - self.x_bounds[0]) * self.idbin)
+ # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+ dd = dv-(self.x_bounds[0] + bin_id * self.dbin) # x - x0
+ if bin_id > self.nbins - 2 or bin_id < 0:
+ for channel in range(4): trgba[channel] = 0.0
+ return
+ for channel in range(4):
+ bv = self.vs[channel][bin_id] # This is x0
+ dy = self.vs[channel][bin_id+1]-bv # dy
+ # This is our final value for transfer function on the entering face
+ trgba[channel] = bv+dd*dy*self.idbin
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void eval_transfer(self, np.float64_t dt, np.float64_t dv,
+ np.float64_t *rgba, np.float64_t *grad):
+ cdef int i
+ cdef np.float64_t ta, tf, trgba[4], dot_prod
+ self.interpolate(dv, trgba)
+ # get source alpha first
+ # First locate our points
+ dot_prod = 0.0
+ if self.use_light:
+ for i in range(3):
+ dot_prod += self.light_dir[i] * grad[i]
+ dot_prod = fmax(0.0, dot_prod)
+ for i in range(3):
+ trgba[i] += dot_prod*self.light_color[i]
+ # alpha blending
+ ta = (1.0 - rgba[3])*dt*trgba[3]
+ for i in range(4):
+ rgba[i] += ta*trgba[i]
cdef class VectorPlane:
cdef public object avp_pos, avp_dir, acenter, aimage
@@ -283,26 +217,24 @@
cdef public object my_data
cdef public object LeftEdge
cdef public object RightEdge
- cdef np.float64_t *data[6]
- cdef np.float64_t dvs[6]
+ cdef np.float64_t *data
cdef np.float64_t left_edge[3]
cdef np.float64_t right_edge[3]
cdef np.float64_t dds[3]
cdef np.float64_t idds[3]
cdef int dims[3]
cdef public int parent_grid_id
- cdef public int n_fields
@cython.boundscheck(False)
@cython.wraparound(False)
def __cinit__(self,
- int parent_grid_id, int n_fields, data,
+ int parent_grid_id,
+ np.ndarray[np.float64_t, ndim=3] data,
np.ndarray[np.float64_t, ndim=1] left_edge,
np.ndarray[np.float64_t, ndim=1] right_edge,
np.ndarray[np.int64_t, ndim=1] dims):
# The data is likely brought in via a slice, so we copy it
cdef int i, j, k, size
- cdef np.ndarray[np.float64_t, ndim=3] tdata
self.parent_grid_id = parent_grid_id
self.LeftEdge = left_edge
self.RightEdge = right_edge
@@ -313,10 +245,7 @@
self.dds[i] = (self.right_edge[i] - self.left_edge[i])/dims[i]
self.idds[i] = 1.0/self.dds[i]
self.my_data = data
- self.n_fields = n_fields
- for i in range(n_fields):
- tdata = data[i]
- self.data[i] = <np.float64_t *> tdata.data
+ self.data = <np.float64_t*> data.data
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -326,7 +255,7 @@
# like http://courses.csusm.edu/cs697exz/ray_box.htm
cdef int vi, vj, hit, i, ni, nj, nn
cdef int iter[4]
- cdef np.float64_t v_pos[3], v_dir[3], rgba[6], extrema[4]
+ cdef np.float64_t v_pos[3], v_dir[3], rgba[4], extrema[4]
self.calculate_extent(vp, extrema)
vp.get_start_stop(extrema, iter)
iter[0] = iclip(iter[0], 0, vp.nv[0])
@@ -337,9 +266,9 @@
for vi in range(iter[0], iter[1]):
for vj in range(iter[2], iter[3]):
vp.copy_into(vp.vp_pos, v_pos, vi, vj, 3)
- vp.copy_into(vp.image, rgba, vi, vj, 6)
+ vp.copy_into(vp.image, rgba, vi, vj, 4)
self.integrate_ray(v_pos, vp.vp_dir, rgba, tf)
- vp.copy_back(rgba, vp.image, vi, vj, 6)
+ vp.copy_back(rgba, vp.image, vi, vj, 4)
return hit
@cython.boundscheck(False)
@@ -494,11 +423,12 @@
for dti in range(tf.ns):
for i in range(3):
dp[i] += ds[i]
- for i in range(self.n_fields):
- self.dvs[i] = trilinear_interpolate(self.dims, ci, dp, self.data[i])
- #if (dv < tf.x_bounds[0]) or (dv > tf.x_bounds[1]):
- # continue
- tf.eval_transfer(dt, self.dvs, rgba, grad)
+ dv = trilinear_interpolate(self.dims, ci, dp, self.data)
+ if (dv < tf.x_bounds[0]) or (dv > tf.x_bounds[1]):
+ continue
+ if tf.use_light == 1:
+ eval_gradient(self.dims, ci, dp, self.data, grad)
More information about the yt-svn
mailing list