[Yt-svn] yt: 3 new changesets

hg at spacepope.org hg at spacepope.org
Sun Apr 25 10:07:56 PDT 2010


hg Repository: yt
details:   yt/rev/0736cbd95c33
changeset: 1601:0736cbd95c33
user:      mturk
date:
Sun Apr 25 10:00:23 2010 -0700
description:
[svn r1696] Backport from hg:

 * New multivariate volume rendering.  This allows for having multiple
   different ways of controlling (independently) the emission, absorption, and
   colors of fields in a volume rendering -- for instance, Planck emission,
   density weighted, with approximate rayleigh scattering.
 * New off-axis projections in the volume rendering module
 * New PNG writer to directly output PNG files, rather than using matplotlib
 * Support for the Tiger code
 * Colored logging output (optional)
 * Some fixes for the fixed res projections
 * John's fix for setting specific marker types in the particle plots

The off-axis projections are currently not "fixed" with units, but that should
be coming soon.  The direct PNG imaging does not (and probably will not, at
least for a while) support colorbars or any other compositing of information.

hg Repository: yt
details:   yt/rev/609855779aa9
changeset: 1602:609855779aa9
user:      mturk
date:
Sun Apr 25 10:01:12 2010 -0700
description:
[svn r1697] Re-cythonizing

hg Repository: yt
details:   yt/rev/f763594fd26e
changeset: 1603:f763594fd26e
user:      Matthew Turk <matthewturk at gmail.com>
date:
Sun Apr 25 10:07:45 2010 -0700
description:
Merging with trunk

diffstat:

 .hgignore                                          |     3 +
 scripts/fbranch                                    |     5 +-
 scripts/fbury                                      |     5 +-
 scripts/fdigup                                     |     5 +-
 scripts/fido                                       |     5 +-
 scripts/fimport                                    |     5 +-
 yt/_amr_utils/VolumeIntegrator.pyx                 |   263 +-
 yt/_amr_utils/fortran_reader.pyx                   |    70 +
 yt/_amr_utils/png_writer.pyx                       |   171 +
 yt/amr_utils.c                                     |  5750 +++++++++++++++----------
 yt/amr_utils.pyx                                   |     2 +
 yt/arraytypes.py                                   |    11 +
 yt/config.py                                       |     1 +
 yt/extensions/EnzoSimulation.py                    |    65 +-
 yt/extensions/MergerTree.py                        |   840 +++
 yt/extensions/StarAnalysis.py                      |    24 +-
 yt/extensions/_colormap_data.py                    |   797 +++
 yt/extensions/image_writer.py                      |    77 +
 yt/extensions/kdtree/fKD.f90                       |    27 +
 yt/extensions/kdtree/fKD.v                         |     2 +
 yt/extensions/volume_rendering/TransferFunction.py |   119 +-
 yt/extensions/volume_rendering/UBVRI.py            |   105 +
 yt/extensions/volume_rendering/__init__.py         |     5 +-
 yt/extensions/volume_rendering/grid_partitioner.py |    28 +-
 yt/extensions/volume_rendering/image_handling.py   |     2 +-
 yt/extensions/volume_rendering/multi_texture.py    |   302 +
 yt/extensions/volume_rendering/software_sampler.py |    14 +-
 yt/fido/share_data.py                              |    83 +
 yt/funcs.py                                        |     9 +-
 yt/lagos/BaseDataTypes.py                          |    21 +-
 yt/lagos/BaseGridType.py                           |    24 +
 yt/lagos/DataReadingFuncs.py                       |    26 +
 yt/lagos/FieldInfoContainer.py                     |     9 +
 yt/lagos/HaloFinding.py                            |    93 +-
 yt/lagos/HierarchyType.py                          |    76 +
 yt/lagos/OutputTypes.py                            |    77 +-
 yt/lagos/ParallelTools.py                          |    49 +-
 yt/lagos/StructureFunctionGenerator.py             |   746 +++
 yt/lagos/__init__.py                               |     2 +
 yt/lagos/kd.py                                     |     2 +-
 yt/lagos/parallelHOP/parallelHOP.py                |    31 +-
 yt/lagos/setup.py                                  |     4 +-
 yt/logger.py                                       |    26 +
 yt/mods.py                                         |     3 +-
 yt/physical_constants.py                           |     3 +-
 yt/raven/Callbacks.py                              |     5 +-
 yt/raven/ColorMaps.py                              |     9 +
 yt/reason/plot_editors.py                          |    10 +-
 yt/reason/reason_v2.py                             |     5 +-
 yt/reason/tvtk_interface.py                        |     5 +-
 yt/setup.py                                        |    62 +-
 51 files changed, 7505 insertions(+), 2578 deletions(-)

diffs (truncated from 18732 to 300 lines):

diff -r 0cfe6fa8de89 -r f763594fd26e .hgignore
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Sun Apr 25 10:07:45 2010 -0700
@@ -0,0 +1,3 @@
+syntax: glob
+*.pyc
+.*.swp
diff -r 0cfe6fa8de89 -r f763594fd26e scripts/fbranch
--- a/scripts/fbranch	Fri Apr 23 17:17:44 2010 -0700
+++ b/scripts/fbranch	Sun Apr 25 10:07:45 2010 -0700
@@ -1,1 +1,4 @@
-basic.py
\ No newline at end of file
+#!python2.5
+import yt.fido
+
+yt.fido.runAction()
diff -r 0cfe6fa8de89 -r f763594fd26e scripts/fbury
--- a/scripts/fbury	Fri Apr 23 17:17:44 2010 -0700
+++ b/scripts/fbury	Sun Apr 25 10:07:45 2010 -0700
@@ -1,1 +1,4 @@
-basic.py
\ No newline at end of file
+#!python2.5
+import yt.fido
+
+yt.fido.runAction()
diff -r 0cfe6fa8de89 -r f763594fd26e scripts/fdigup
--- a/scripts/fdigup	Fri Apr 23 17:17:44 2010 -0700
+++ b/scripts/fdigup	Sun Apr 25 10:07:45 2010 -0700
@@ -1,1 +1,4 @@
-basic.py
\ No newline at end of file
+#!python2.5
+import yt.fido
+
+yt.fido.runAction()
diff -r 0cfe6fa8de89 -r f763594fd26e scripts/fido
--- a/scripts/fido	Fri Apr 23 17:17:44 2010 -0700
+++ b/scripts/fido	Sun Apr 25 10:07:45 2010 -0700
@@ -1,1 +1,4 @@
-basic.py
\ No newline at end of file
+#!python2.5
+from yt.mods import *
+
+fido.runAction()
diff -r 0cfe6fa8de89 -r f763594fd26e scripts/fimport
--- a/scripts/fimport	Fri Apr 23 17:17:44 2010 -0700
+++ b/scripts/fimport	Sun Apr 25 10:07:45 2010 -0700
@@ -1,1 +1,4 @@
-basic.py
\ No newline at end of file
+#!python2.5
+from yt.mods import *
+
+fido.runAction()
diff -r 0cfe6fa8de89 -r f763594fd26e yt/_amr_utils/VolumeIntegrator.pyx
--- a/yt/_amr_utils/VolumeIntegrator.pyx	Fri Apr 23 17:17:44 2010 -0700
+++ b/yt/_amr_utils/VolumeIntegrator.pyx	Sun Apr 25 10:07:45 2010 -0700
@@ -73,87 +73,160 @@
 
 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
+    int pass_through
+
+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,
+              int pass_through = 0):
+    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
+    fit.pass_through = pass_through
+
+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
+    if fit.pass_through == 1: return dvs[fit.field_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 np.float64_t x_bounds[2]
-    cdef np.float64_t *vs[4]
-    cdef int nbins
+    cdef int n_fields
+    cdef int n_field_tables
     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
+
+    # 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
-        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
 
+        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],
+                      tf_obj.tables[i].pass_through)
+            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 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 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 class VectorPlane:
     cdef public object avp_pos, avp_dir, acenter, aimage
     cdef np.float64_t *vp_pos, *vp_dir, *center, *image,
     cdef np.float64_t pdx, pdy, bounds[4]
     cdef int nv[2]
+    cdef int vp_strides[3]
+    cdef int im_strides[3]
     cdef public object ax_vec, ay_vec
     cdef np.float64_t *x_vec, *y_vec
 
@@ -183,6 +256,9 @@
         for i in range(4): self.bounds[i] = bounds[i]
         self.pdx = (self.bounds[1] - self.bounds[0])/self.nv[0]
         self.pdy = (self.bounds[3] - self.bounds[2])/self.nv[1]
+        for i in range(3):
+            self.vp_strides[i] = vp_pos.strides[i] / 8
+            self.im_strides[i] = image.strides[i] / 8
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
@@ -200,41 +276,45 @@
         rv[3] = rv[2] + lrint((ex[3] - ex[2])/self.pdy)



More information about the yt-svn mailing list