[Yt-svn] yt: 3 new changesets

hg at spacepope.org hg at spacepope.org
Sat Apr 10 17:52:55 PDT 2010


hg Repository: yt
details:   yt/rev/c900b2dfdc66
changeset: 1545:c900b2dfdc66
user:      Matthew Turk <matthewturk at gmail.com>
date:
Fri Apr 09 11:05:00 2010 -0700
description:
Changing some things with the vertex-centering and +1 nature of the texture shapes for visvis

hg Repository: yt
details:   yt/rev/78ccd7ce364b
changeset: 1546:78ccd7ce364b
user:      Matthew Turk <matthewturk at gmail.com>
date:
Sat Apr 10 17:52:17 2010 -0700
description:
This changes the setup.py in lagos to do a cascading search for HDF5: first it
tries HDF5_DIR, then hdf5.cfg, then if all else fails it tries
ctypes.util.find_library.  ctypes isn't available everywhere so this is wrapped
in an ImportError catch.

hg Repository: yt
details:   yt/rev/9cd931a41384
changeset: 1547:9cd931a41384
user:      Matthew Turk <matthewturk at gmail.com>
date:
Sat Apr 10 17:52:49 2010 -0700
description:
Merge

diffstat:

 scripts/yt_lodgeit.py                              |      1 +
 yt/_amr_utils/VolumeIntegrator.pyx                 |    251 +-
 yt/amr_utils.c                                     |  18107 +++++++++++-------------
 yt/extensions/kdtree/fKD.f90                       |     27 +
 yt/extensions/kdtree/fKD.v                         |      2 +
 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/multi_texture.py    |    306 -
 yt/extensions/volume_rendering/software_sampler.py |     12 +-
 yt/lagos/BaseDataTypes.py                          |      2 +-
 yt/lagos/DerivedQuantities.py                      |     18 +-
 yt/lagos/EnzoFields.py                             |      2 +
 yt/lagos/HaloFinding.py                            |      3 +-
 yt/lagos/ParallelTools.py                          |      4 +-
 yt/lagos/StructureFunctionGenerator.py             |    390 +-
 yt/lagos/setup.py                                  |     49 +-
 yt/physical_constants.py                           |      3 +-
 yt/raven/PlotCollection.py                         |     32 +-
 yt/raven/PlotTypes.py                              |     27 +-
 22 files changed, 8927 insertions(+), 10551 deletions(-)

diffs (truncated from 30702 to 300 lines):

diff -r 7886523c2bd3 -r 9cd931a41384 scripts/yt_lodgeit.py
--- a/scripts/yt_lodgeit.py	Fri Apr 09 07:42:18 2010 -0700
+++ b/scripts/yt_lodgeit.py	Sat Apr 10 17:52:49 2010 -0700
@@ -222,6 +222,7 @@
         data = read_file(sys.stdin)
         if not langopt:
             mime = get_mimetype(data, '') or ''
+        fname = ""
     elif len(filenames) == 1:
         fname = filenames[0]
         data = read_file(open(filenames[0], 'rb'))
diff -r 7886523c2bd3 -r 9cd931a41384 yt/_amr_utils/VolumeIntegrator.pyx
--- a/yt/_amr_utils/VolumeIntegrator.pyx	Fri Apr 09 07:42:18 2010 -0700
+++ b/yt/_amr_utils/VolumeIntegrator.pyx	Sat Apr 10 17:52:49 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):



More information about the yt-svn mailing list