[Yt-svn] yt-commit r1531 - in trunk: tests yt yt/_amr_utils yt/extensions/volume_rendering yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Nov 20 19:54:18 PST 2009
Author: mturk
Date: Fri Nov 20 19:54:16 2009
New Revision: 1531
URL: http://yt.enzotools.org/changeset/1531
Log:
Moving all of the utility functions to amr_utils. Eases compilation on Kraken
for CNL. Thanks to Britton for a lot of this...
Added:
trunk/yt/_amr_utils/
trunk/yt/_amr_utils/CICDeposit.pyx
trunk/yt/_amr_utils/DepthFirstOctree.pyx
trunk/yt/_amr_utils/FixedInterpolator.c
trunk/yt/_amr_utils/FixedInterpolator.h
trunk/yt/_amr_utils/Interpolators.pyx
trunk/yt/_amr_utils/PointsInVolume.pyx
trunk/yt/_amr_utils/RayIntegrators.pyx
trunk/yt/_amr_utils/VolumeIntegrator.pyx
trunk/yt/amr_utils.pyx
Removed:
trunk/yt/extensions/volume_rendering/FixedInterpolator.c
trunk/yt/extensions/volume_rendering/FixedInterpolator.h
trunk/yt/extensions/volume_rendering/VolumeIntegrator.c
trunk/yt/extensions/volume_rendering/VolumeIntegrator.pyx
trunk/yt/lagos/DepthFirstOctree.c
trunk/yt/lagos/DepthFirstOctree.pyx
trunk/yt/lagos/Interpolators.c
trunk/yt/lagos/Interpolators.pyx
trunk/yt/lagos/PointsInVolume.c
trunk/yt/lagos/PointsInVolume.pyx
Modified:
trunk/tests/test_lagos.py
trunk/yt/extensions/volume_rendering/__init__.py
trunk/yt/extensions/volume_rendering/setup.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/HelperFunctions.py
trunk/yt/lagos/__init__.py
trunk/yt/lagos/setup.py
trunk/yt/setup.py
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Fri Nov 20 19:54:16 2009
@@ -454,7 +454,7 @@
class TestSliceDataType(DataTypeTestingBase, LagosTestingBase, unittest.TestCase):
def setUp(self):
DataTypeTestingBase.setUp(self)
- self.data = self.hierarchy.slice(0,0.5)
+ self.data = self.hierarchy.slice(0,0.5, center=[0.5, 0.5, 0.5])
class TestCuttingPlane(DataTypeTestingBase, LagosTestingBase, unittest.TestCase):
def setUp(self):
Added: trunk/yt/_amr_utils/CICDeposit.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/CICDeposit.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,80 @@
+"""
+Simle integrators for the radiative transfer equation
+
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: CASA/University of Colorado
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 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/>.
+"""
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+def CICDeposit_3(np.ndarray[np.float64_t, ndim=1] posx,
+ np.ndarray[np.float64_t, ndim=1] posy,
+ np.ndarray[np.float64_t, ndim=1] posz,
+ np.ndarray[np.float32_t, ndim=1] mass,
+ np.int64_t npositions,
+ np.ndarray[np.float32_t, ndim=3] field,
+ np.ndarray[np.float64_t, ndim=1] leftEdge,
+ np.ndarray[np.int32_t, ndim=1] gridDimension,
+ np.float64_t cellSize):
+
+ cdef int i1, j1, k1, n
+ cdef double xpos, ypos, zpos
+ cdef double fact, edge0, edge1, edge2
+ cdef double le0, le1, le2
+ cdef float dx, dy, dz, dx2, dy2, dz2
+
+ edge0 = (<float> gridDimension[0]) - 0.5001
+ edge1 = (<float> gridDimension[1]) - 0.5001
+ edge2 = (<float> gridDimension[2]) - 0.5001
+ fact = 1.0 / cellSize
+
+ le0 = leftEdge[0]
+ le1 = leftEdge[1]
+ le2 = leftEdge[2]
+
+ for n in range(npositions):
+
+ # Compute the position of the central cell
+ xpos = fmin(fmax((posx[n] - le0)*fact, 0.5001), edge0)
+ ypos = fmin(fmax((posy[n] - le1)*fact, 0.5001), edge1)
+ zpos = fmin(fmax((posz[n] - le2)*fact, 0.5001), edge2)
+
+ i1 = <int> (xpos + 0.5)
+ j1 = <int> (ypos + 0.5)
+ k1 = <int> (zpos + 0.5)
+
+ # Compute the weights
+ dx = (<float> i1) + 0.5 - xpos
+ dy = (<float> j1) + 0.5 - ypos
+ dz = (<float> k1) + 0.5 - zpos
+ dx2 = 1.0 - dx
+ dy2 = 1.0 - dy
+ dz2 = 1.0 - dz
+
+ # Interpolate from field into sumfield
+ field[i1-1,j1-1,k1-1] += mass[n] * dx * dy * dz
+ field[i1 ,j1-1,k1-1] += mass[n] * dx2 * dy * dz
+ field[i1-1,j1 ,k1-1] += mass[n] * dx * dy2 * dz
+ field[i1 ,j1 ,k1-1] += mass[n] * dx2 * dy2 * dz
+ field[i1-1,j1-1,k1 ] += mass[n] * dx * dy * dz2
+ field[i1 ,j1-1,k1 ] += mass[n] * dx2 * dy * dz2
+ field[i1-1,j1 ,k1 ] += mass[n] * dx * dy2 * dz2
+ field[i1 ,j1 ,k1 ] += mass[n] * dx2 * dy2 * dz2
Added: trunk/yt/_amr_utils/DepthFirstOctree.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/DepthFirstOctree.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,167 @@
+"""
+This is a recursive function to return a depth-first octree
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 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 np
+cimport numpy as np
+cimport cython
+
+cdef class position:
+ cdef public int output_pos, refined_pos
+ def __cinit__(self):
+ self.output_pos = 0
+ self.refined_pos = 0
+
+cdef class OctreeGrid:
+ cdef public object child_indices, fields, left_edges, dimensions, dx
+ cdef public int level
+ def __cinit__(self,
+ np.ndarray[np.int32_t, ndim=3] child_indices,
+ np.ndarray[np.float64_t, ndim=4] fields,
+ np.ndarray[np.float64_t, ndim=1] left_edges,
+ np.ndarray[np.int32_t, ndim=1] dimensions,
+ np.ndarray[np.float64_t, ndim=1] dx,
+ int level):
+ self.child_indices = child_indices
+ self.fields = fields
+ self.left_edges = left_edges
+ self.dimensions = dimensions
+ self.dx = dx
+ self.level = level
+
+cdef class OctreeGridList:
+ cdef public object grids
+ def __cinit__(self, grids):
+ self.grids = grids
+
+ def __getitem__(self, int item):
+ return self.grids[item]
+
+ at cython.boundscheck(False)
+def RecurseOctreeDepthFirst(int i_i, int j_i, int k_i,
+ int i_f, int j_f, int k_f,
+ position curpos, int gi,
+ np.ndarray[np.float64_t, ndim=2] output,
+ np.ndarray[np.int32_t, ndim=1] refined,
+ OctreeGridList grids):
+ cdef int i, i_off, j, j_off, k, k_off, ci, fi
+ cdef int child_i, child_j, child_k
+ cdef OctreeGrid child_grid
+ cdef OctreeGrid grid = grids[gi-1]
+ cdef np.ndarray[np.int32_t, ndim=3] child_indices = grid.child_indices
+ cdef np.ndarray[np.int32_t, ndim=1] dimensions = grid.dimensions
+ cdef np.ndarray[np.float64_t, ndim=4] fields = grid.fields
+ cdef np.ndarray[np.float64_t, ndim=1] leftedges = grid.left_edges
+ cdef np.float64_t dx = grid.dx[0]
+ cdef np.float64_t child_dx
+ cdef np.ndarray[np.float64_t, ndim=1] child_leftedges
+ cdef np.float64_t cx, cy, cz
+ for k_off in range(k_f):
+ k = k_off + k_i
+ cz = (leftedges[2] + k*dx)
+ for j_off in range(j_f):
+ j = j_off + j_i
+ cy = (leftedges[1] + j*dx)
+ for i_off in range(i_f):
+ i = i_off + i_i
+ cx = (leftedges[0] + i*dx)
+ ci = grid.child_indices[i,j,k]
+ if ci == -1:
+ for fi in range(fields.shape[0]):
+ output[curpos.output_pos,fi] = fields[fi,i,j,k]
+ refined[curpos.refined_pos] = 0
+ curpos.output_pos += 1
+ curpos.refined_pos += 1
+ else:
+ refined[curpos.refined_pos] = 1
+ curpos.refined_pos += 1
+ child_grid = grids[ci-1]
+ child_dx = child_grid.dx[0]
+ child_leftedges = child_grid.left_edges
+ child_i = int((cx - child_leftedges[0])/child_dx)
+ child_j = int((cy - child_leftedges[1])/child_dx)
+ child_k = int((cz - child_leftedges[2])/child_dx)
+ s = RecurseOctreeDepthFirst(child_i, child_j, child_k, 2, 2, 2,
+ curpos, ci, output, refined, grids)
+ return s
+
+ at cython.boundscheck(False)
+def RecurseOctreeByLevels(int i_i, int j_i, int k_i,
+ int i_f, int j_f, int k_f,
+ np.ndarray[np.int32_t, ndim=1] curpos,
+ int gi,
+ np.ndarray[np.float64_t, ndim=2] output,
+ np.ndarray[np.int32_t, ndim=2] genealogy,
+ np.ndarray[np.float64_t, ndim=2] corners,
+ OctreeGridList grids):
+ cdef np.int32_t i, i_off, j, j_off, k, k_off, ci, fi
+ cdef int child_i, child_j, child_k
+ cdef OctreeGrid child_grid
+ cdef OctreeGrid grid = grids[gi-1]
+ cdef int level = grid.level
+ cdef np.ndarray[np.int32_t, ndim=3] child_indices = grid.child_indices
+ cdef np.ndarray[np.int32_t, ndim=1] dimensions = grid.dimensions
+ cdef np.ndarray[np.float64_t, ndim=4] fields = grid.fields
+ cdef np.ndarray[np.float64_t, ndim=1] leftedges = grid.left_edges
+ cdef np.float64_t dx = grid.dx[0]
+ cdef np.float64_t child_dx
+ cdef np.ndarray[np.float64_t, ndim=1] child_leftedges
+ cdef np.float64_t cx, cy, cz
+ cdef int cp
+ for i_off in range(i_f):
+ i = i_off + i_i
+ cx = (leftedges[0] + i*dx)
+ if i_f > 2: print k, cz
+ for j_off in range(j_f):
+ j = j_off + j_i
+ cy = (leftedges[1] + j*dx)
+ for k_off in range(k_f):
+ k = k_off + k_i
+ cz = (leftedges[2] + k*dx)
+ cp = curpos[level]
+ corners[cp, 0] = cx
+ corners[cp, 1] = cy
+ corners[cp, 2] = cz
+ genealogy[curpos[level], 2] = level
+ # always output data
+ for fi in range(fields.shape[0]):
+ output[cp,fi] = fields[fi,i,j,k]
+ ci = child_indices[i,j,k]
+ if ci > -1:
+ child_grid = grids[ci-1]
+ child_dx = child_grid.dx[0]
+ child_leftedges = child_grid.left_edges
+ child_i = int((cx-child_leftedges[0])/child_dx)
+ child_j = int((cy-child_leftedges[1])/child_dx)
+ child_k = int((cz-child_leftedges[2])/child_dx)
+ # set current child id to id of next cell to examine
+ genealogy[cp, 0] = curpos[level+1]
+ # set next parent id to id of current cell
+ genealogy[curpos[level+1]:curpos[level+1]+8, 1] = cp
+ s = RecurseOctreeByLevels(child_i, child_j, child_k, 2, 2, 2,
+ curpos, ci, output, genealogy,
+ corners, grids)
+ curpos[level] += 1
+ return s
+
Added: trunk/yt/_amr_utils/FixedInterpolator.c
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/FixedInterpolator.c Fri Nov 20 19:54:16 2009
@@ -0,0 +1,74 @@
+/************************************************************************
+* 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/>.
+*
+************************************************************************/
+
+
+//
+// A small, tiny, itty bitty module for computation-intensive interpolation
+// that I can't seem to make fast in Cython
+//
+
+#include "FixedInterpolator.h"
+
+#define VINDEX(A,B,C) data[((((A)+ci[0])*(ds[1]+1)+((B)+ci[1]))*(ds[2]+1)+ci[2]+(C))]
+// (((C*ds[1])+B)*ds[0]+A)
+
+npy_float64 fast_interpolate(int *ds, int *ci, npy_float64 *dp,
+ npy_float64 *data)
+{
+ int i;
+ npy_float64 dv, dm[3];
+ for(i=0;i<3;i++)dm[i] = (1.0 - dp[i]);
+ dv = 0.0;
+ dv += VINDEX(0,0,0) * (dm[0]*dm[1]*dm[2]);
+ dv += VINDEX(0,0,1) * (dm[0]*dm[1]*dp[2]);
+ dv += VINDEX(0,1,0) * (dm[0]*dp[1]*dm[2]);
+ dv += VINDEX(0,1,1) * (dm[0]*dp[1]*dp[2]);
+ dv += VINDEX(1,0,0) * (dp[0]*dm[1]*dm[2]);
+ dv += VINDEX(1,0,1) * (dp[0]*dm[1]*dp[2]);
+ dv += VINDEX(1,1,0) * (dp[0]*dp[1]*dm[2]);
+ dv += VINDEX(1,1,1) * (dp[0]*dp[1]*dp[2]);
+ /*assert(dv < -20);*/
+ return dv;
+}
+
+npy_float64 trilinear_interpolate(int *ds, int *ci, npy_float64 *dp,
+ npy_float64 *data)
+{
+ /* dims is one less than the dimensions of the array */
+ int i;
+ npy_float64 dm[3], vz[4];
+ //dp is the distance to the plane. dm is val, dp = 1-val
+ for(i=0;i<3;i++)dm[i] = (1.0 - dp[i]);
+
+ //First interpolate in z
+ vz[0] = dm[2]*VINDEX(0,0,0) + dp[2]*VINDEX(0,0,1);
+ vz[1] = dm[2]*VINDEX(0,1,0) + dp[2]*VINDEX(0,1,1);
+ vz[2] = dm[2]*VINDEX(1,0,0) + dp[2]*VINDEX(1,0,1);
+ vz[3] = dm[2]*VINDEX(1,1,0) + dp[2]*VINDEX(1,1,1);
+
+ //Then in y
+ vz[0] = dm[1]*vz[0] + dp[1]*vz[1];
+ vz[1] = dm[1]*vz[2] + dp[1]*vz[3];
+
+ //Then in x
+ vz[0] = dm[0]*vz[0] + dp[0]*vz[1];
+ /*assert(dv < -20);*/
+ return vz[0];
+}
Added: trunk/yt/_amr_utils/FixedInterpolator.h
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/FixedInterpolator.h Fri Nov 20 19:54:16 2009
@@ -0,0 +1,40 @@
+/************************************************************************
+* 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/>.
+*
+************************************************************************/
+
+
+//
+// A small, tiny, itty bitty module for computation-intensive interpolation
+// that I can't seem to make fast in Cython
+//
+
+#include "Python.h"
+
+#include <stdio.h>
+#include <math.h>
+#include <signal.h>
+#include <ctype.h>
+
+#include "numpy/ndarrayobject.h"
+
+npy_float64 fast_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
+ npy_float64 *data);
+
+npy_float64 trilinear_interpolate(int ds[3], int ci[3], npy_float64 dp[3],
+ npy_float64 *data);
Added: trunk/yt/_amr_utils/Interpolators.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/Interpolators.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,116 @@
+"""
+Simple interpolators
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 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 np
+cimport numpy as np
+cimport cython
+
+ at cython.boundscheck(False)
+def UnilinearlyInterpolate(np.ndarray[np.float64_t, ndim=1] table,
+ np.ndarray[np.float64_t, ndim=1] x_vals,
+ np.ndarray[np.float64_t, ndim=1] x_bins,
+ np.ndarray[np.int32_t, ndim=1] x_is,
+ np.ndarray[np.float64_t, ndim=1] output):
+ cdef double x, xp, xm
+ cdef int i, x_i, y_i
+ for i in range(x_vals.shape[0]):
+ x_i = x_is[i]
+ x = x_vals[i]
+ dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+ xp = (x - x_bins[x_i]) * dx_inv
+ xm = (x_bins[x_i+1] - x) * dx_inv
+ output[i] = table[x_i ] * (xm) \
+ + table[x_i+1] * (xp)
+
+ at cython.boundscheck(False)
+def BilinearlyInterpolate(np.ndarray[np.float64_t, ndim=2] table,
+ np.ndarray[np.float64_t, ndim=1] x_vals,
+ np.ndarray[np.float64_t, ndim=1] y_vals,
+ np.ndarray[np.float64_t, ndim=1] x_bins,
+ np.ndarray[np.float64_t, ndim=1] y_bins,
+ np.ndarray[np.int32_t, ndim=1] x_is,
+ np.ndarray[np.int32_t, ndim=1] y_is,
+ np.ndarray[np.float64_t, ndim=1] output):
+ cdef double x, xp, xm
+ cdef double y, yp, ym
+ cdef double dx_inv, dy_inv
+ cdef int i, x_i, y_i
+ for i in range(x_vals.shape[0]):
+ x_i = x_is[i]
+ y_i = y_is[i]
+ x = x_vals[i]
+ y = y_vals[i]
+ dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+ dy_inv = 1.0 / (y_bins[y_i+1] - y_bins[y_i])
+ xp = (x - x_bins[x_i]) * dx_inv
+ yp = (y - y_bins[y_i]) * dy_inv
+ xm = (x_bins[x_i+1] - x) * dx_inv
+ ym = (y_bins[y_i+1] - y) * dy_inv
+ output[i] = table[x_i , y_i ] * (xm*ym) \
+ + table[x_i+1, y_i ] * (xp*ym) \
+ + table[x_i , y_i+1] * (xm*yp) \
+ + table[x_i+1, y_i+1] * (xp*yp)
+
+ at cython.boundscheck(False)
+def TrilinearlyInterpolate(np.ndarray[np.float64_t, ndim=3] table,
+ np.ndarray[np.float64_t, ndim=1] x_vals,
+ np.ndarray[np.float64_t, ndim=1] y_vals,
+ np.ndarray[np.float64_t, ndim=1] z_vals,
+ np.ndarray[np.float64_t, ndim=1] x_bins,
+ np.ndarray[np.float64_t, ndim=1] y_bins,
+ np.ndarray[np.float64_t, ndim=1] z_bins,
+ np.ndarray[np.int_t, ndim=1] x_is,
+ np.ndarray[np.int_t, ndim=1] y_is,
+ np.ndarray[np.int_t, ndim=1] z_is,
+ np.ndarray[np.float64_t, ndim=1] output):
+ cdef double x, xp, xm
+ cdef double y, yp, ym
+ cdef double z, zp, zm
+ cdef double dx_inv, dy_inv, dz_inv
+ cdef int i, x_i, y_i, z_i
+ for i in range(x_vals.shape[0]):
+ x_i = x_is[i]
+ y_i = y_is[i]
+ z_i = z_is[i]
+ x = x_vals[i]
+ y = y_vals[i]
+ z = z_vals[i]
+ dx_inv = 1.0 / (x_bins[x_i+1] - x_bins[x_i])
+ dy_inv = 1.0 / (y_bins[y_i+1] - y_bins[y_i])
+ dz_inv = 1.0 / (z_bins[z_i+1] - z_bins[z_i])
+ xp = (x - x_bins[x_i]) * dx_inv
+ yp = (y - y_bins[y_i]) * dy_inv
+ zp = (z - z_bins[z_i]) * dz_inv
+ xm = (x_bins[x_i+1] - x) * dx_inv
+ ym = (y_bins[y_i+1] - y) * dy_inv
+ zm = (z_bins[z_i+1] - z) * dz_inv
+ output[i] = table[x_i ,y_i ,z_i ] * (xm*ym*zm) \
+ + table[x_i+1,y_i ,z_i ] * (xp*ym*zm) \
+ + table[x_i ,y_i+1,z_i ] * (xm*yp*zm) \
+ + table[x_i ,y_i ,z_i+1] * (xm*ym*zp) \
+ + table[x_i+1,y_i ,z_i+1] * (xp*ym*zp) \
+ + table[x_i ,y_i+1,z_i+1] * (xm*yp*zp) \
+ + table[x_i+1,y_i+1,z_i ] * (xp*yp*zm) \
+ + table[x_i+1,y_i+1,z_i+1] * (xp*yp*zp)
Added: trunk/yt/_amr_utils/PointsInVolume.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/PointsInVolume.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,73 @@
+"""
+Checks for points contained in a volume
+
+Author: John Wise <jwise77 at gmail.com>
+Affiliation: Princeton
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2009 John. 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 np
+cimport numpy as np
+cimport cython
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def PointsInVolume(np.ndarray[np.float64_t, ndim=2] points,
+ np.ndarray[np.int8_t, ndim=1] pmask, # pixel mask
+ np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] right_edge,
+ np.ndarray[np.int32_t, ndim=3] mask,
+ float dx):
+ cdef np.ndarray[np.int8_t, ndim=1] \
+ valid = np.zeros(points.shape[0], dtype='int8')
+ cdef int i, dim, count
+ cdef int ex
+ cdef double dx_inv
+ cdef unsigned int idx[3]
+ count = 0
+ dx_inv = 1.0 / dx
+ for i in xrange(points.shape[0]):
+ if pmask[i] == 0:
+ continue
+ ex = 1
+ for dim in xrange(3):
+ if points[i,dim] < left_edge[dim] or points[i,dim] > right_edge[dim]:
+ valid[i] = ex = 0
+ break
+ if ex == 1:
+ for dim in xrange(3):
+ idx[dim] = <unsigned int> \
+ ((points[i,dim] - left_edge[dim]) * dx_inv)
+ if mask[idx[0], idx[1], idx[2]] == 1:
+ valid[i] = 1
+ count += 1
+
+ cdef np.ndarray[np.int32_t, ndim=1] result = np.empty(count, dtype='int32')
+ count = 0
+ for i in xrange(points.shape[0]):
+ if valid[i] == 1 and pmask[i] == 1:
+ result[count] = i
+ count += 1
+
+ return result
+
+
+
Added: trunk/yt/_amr_utils/RayIntegrators.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/RayIntegrators.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,374 @@
+"""
+Simle integrators for the radiative transfer equation
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 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 np
+cimport numpy as np
+cimport cython
+from stdlib cimport malloc, free, abs
+
+cdef extern from "math.h":
+ double exp(double x)
+ float expf(float x)
+
+ at cython.boundscheck(False)
+def Transfer3D(np.ndarray[np.float64_t, ndim=3] i_s,
+ np.ndarray[np.float64_t, ndim=4] o_s,
+ np.ndarray[np.float64_t, ndim=4] e,
+ np.ndarray[np.float64_t, ndim=4] a,
+ int imin, int imax, int jmin, int jmax,
+ int kmin, int kmax, int istride, int jstride,
+ np.float64_t dx):
+ """
+ This function accepts an incoming slab (*i_s*), a buffer
+ for an outgoing set of values at every point in the grid (*o_s*),
+ an emission array (*e*), an absorption array (*a*), and dimensions of
+ the grid (*imin*, *imax*, *jmin*, *jmax*, *kmin*, *kmax*) as well
+ as strides in the *i* and *j* directions, and a *dx* of the grid being
+ integrated.
+ """
+ cdef int i, ii
+ cdef int j, jj
+ cdef int k, kk
+ cdef int n, nn
+ nn = o_s.shape[3] # This might be slow
+ cdef np.float64_t *temp = <np.float64_t *>malloc(sizeof(np.float64_t) * nn)
+ for i in range((imax-imin)*istride):
+ ii = i + imin*istride
+ for j in range((jmax-jmin)*jstride):
+ jj = j + jmin*jstride
+ # Not sure about the ordering of the loops here
+ for n in range(nn):
+ temp[n] = i_s[ii,jj,n]
+ for k in range(kmax-kmin):
+ kk = k + kmin#*kstride, which doesn't make any sense
+ for n in range(nn):
+ o_s[i,j,k,n] = temp[n] + dx*(e[i,j,k,n] - temp[n]*a[i,j,k,n])
+ temp[n] = o_s[i,j,k,n]
+ for n in range(nn):
+ i_s[ii,jj,n] = temp[n]
+ free(temp)
+
+ at cython.boundscheck(False)
+def TransferShells(np.ndarray[np.float64_t, ndim=3] i_s,
+ np.ndarray[np.float64_t, ndim=3] data,
+ np.ndarray[np.float64_t, ndim=2] shells):
+ """
+ This function accepts an incoming slab (*i_s*), a buffer of *data*,
+ and a list of shells specified as [ (value, tolerance, r, g, b), ... ].
+ """
+ cdef int i, ii
+ cdef int j, jj
+ cdef int k, kk
+ cdef int n, nn
+ cdef np.float64_t dist
+ ii = data.shape[0]
+ jj = data.shape[1]
+ kk = data.shape[2]
+ nn = shells.shape[0]
+ cdef float rgba[4]
+ cdef float alpha
+ for i in range(ii):
+ for j in range(jj):
+ # Not sure about the ordering of the loops here
+ for k in range(kk):
+ for n in range(nn):
+ dist = shells[n, 0] - data[i,j,k]
+ if dist < 0: dist *= -1.0
+ if dist < shells[n,1]:
+ dist = exp(-dist/8.0)
+ rgba[0] = shells[n,2]
+ rgba[1] = shells[n,3]
+ rgba[2] = shells[n,4]
+ rgba[3] = shells[n,5]
+ alpha = i_s[i,j,3]
+ dist *= dist # This might improve appearance
+ i_s[i,j,0] += (1.0 - alpha)*rgba[0]*dist*rgba[3]
+ i_s[i,j,1] += (1.0 - alpha)*rgba[1]*dist*rgba[3]
+ i_s[i,j,2] += (1.0 - alpha)*rgba[2]*dist*rgba[3]
+ i_s[i,j,3] += (1.0 - alpha)*rgba[3]*dist*rgba[3]
+ break
+
+ at cython.boundscheck(False)
+def Transfer1D(float i_s,
+ np.ndarray[np.float_t, ndim=1] o_s,
+ np.ndarray[np.float_t, ndim=1] e,
+ np.ndarray[np.float_t, ndim=1] a,
+ np.ndarray[np.float_t, ndim=1] dx,
+ int imin, int imax):
+ cdef int i
+ for i in range(imin, imax):
+ o_s[i] = i_s + dx[i]*(e[i] - i_s*a[i])
+ i_s = o_s[i]
+ return i_s
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def VoxelTraversal(np.ndarray[np.int_t, ndim=3] grid_mask,
+ np.ndarray[np.float64_t, ndim=3] grid_t,
+ np.ndarray[np.float64_t, ndim=3] grid_dt,
+ np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] right_edge,
+ np.ndarray[np.float64_t, ndim=1] dx,
+ np.ndarray[np.float64_t, ndim=1] u,
+ np.ndarray[np.float64_t, ndim=1] v):
+ # We're roughly following Amanatides & Woo
+ # Find the first place the ray hits the grid on its path
+ # Do left edge then right edge in each dim
+ cdef int i, x, y
+ cdef np.float64_t tl, tr, intersect_t, enter_t, exit_t, dt_tolerance
+ cdef np.ndarray[np.int64_t, ndim=1] step = np.empty((3,), dtype=np.int64)
+ cdef np.ndarray[np.int64_t, ndim=1] cur_ind = np.empty((3,), dtype=np.int64)
+ cdef np.ndarray[np.float64_t, ndim=1] tdelta = np.empty((3,), dtype=np.float64)
+ cdef np.ndarray[np.float64_t, ndim=1] tmax = np.empty((3,), dtype=np.float64)
+ cdef np.ndarray[np.float64_t, ndim=1] intersect = np.empty((3,), dtype=np.float64)
+ intersect_t = 1
+ dt_tolerance = 1e-6
+ # recall p = v * t + u
+ # where p is position, v is our vector, u is the start point
+ for i in range(3):
+ # As long as we're iterating, set some other stuff, too
+ if(v[i] < 0): step[i] = -1
+ else: step[i] = 1
+ x = (i+1)%3
+ y = (i+2)%3
+ tl = (left_edge[i] - u[i])/v[i]
+ tr = (right_edge[i] - u[i])/v[i]
+ if (left_edge[x] <= (u[x] + tl*v[x]) <= right_edge[x]) and \
+ (left_edge[y] <= (u[y] + tl*v[y]) <= right_edge[y]) and \
+ (0.0 <= tl < intersect_t):
+ intersect_t = tl
+ if (left_edge[x] <= (u[x] + tr*v[x]) <= right_edge[x]) and \
+ (left_edge[y] <= (u[y] + tr*v[y]) <= right_edge[y]) and \
+ (0.0 <= tr < intersect_t):
+ intersect_t = tr
+ # if fully enclosed
+ if (left_edge[0] <= u[0] <= right_edge[0]) and \
+ (left_edge[1] <= u[1] <= right_edge[1]) and \
+ (left_edge[2] <= u[2] <= right_edge[2]):
+ intersect_t = 0.0
+ if not (0 <= intersect_t <= 1): return
+ # Now get the indices of the intersection
+ intersect = u + intersect_t * v
+ for i in range(3):
+ cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i])
+ tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i]
+ if cur_ind[i] == grid_mask.shape[i] and step[i] < 0:
+ cur_ind[i] = grid_mask.shape[i] - 1
+ if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i]
+ if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i]
+ tdelta[i] = (dx[i]/v[i])
+ if tdelta[i] < 0: tdelta[i] *= -1
+ # The variable intersect contains the point we first pierce the grid
+ enter_t = intersect_t
+ while 1:
+ if (not (0 <= cur_ind[0] < grid_mask.shape[0])) or \
+ (not (0 <= cur_ind[1] < grid_mask.shape[1])) or \
+ (not (0 <= cur_ind[2] < grid_mask.shape[2])):
+ break
+ # Note that we are calculating t on the fly, but we get *negative* t
+ # values from what they should be.
+ # If we've reached t = 1, we are done.
+ grid_mask[cur_ind[0], cur_ind[1], cur_ind[2]] = 1
+ if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0):
+ grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = 1.0 - enter_t
+ break
+ if tmax[0] < tmax[1]:
+ if tmax[0] < tmax[2]:
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+ grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[0] - enter_t
+ enter_t = tmax[0]
+ tmax[0] += tdelta[0]
+ cur_ind[0] += step[0]
+ else:
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+ grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
+ else:
+ if tmax[1] < tmax[2]:
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+ grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[1] - enter_t
+ enter_t = tmax[1]
+ tmax[1] += tdelta[1]
+ cur_ind[1] += step[1]
+ else:
+ grid_t[cur_ind[0], cur_ind[1], cur_ind[2]] = enter_t
+ grid_dt[cur_ind[0], cur_ind[1], cur_ind[2]] = tmax[2] - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
+ return
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def PlaneVoxelIntegration(np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] right_edge,
+ np.ndarray[np.float64_t, ndim=1] dx,
+ np.ndarray[np.float64_t, ndim=2] ug,
+ np.ndarray[np.float64_t, ndim=1] v,
+ np.ndarray[np.float64_t, ndim=2] image,
+ np.ndarray[np.float64_t, ndim=3] data,
+ np.ndarray[np.float64_t, ndim=2] shells):
+ # We're roughly following Amanatides & Woo on a ray-by-ray basis
+ # Note that for now it's just shells, but this can and should be
+ # generalized to transfer functions
+ cdef int i, x, y, vi
+ intersect_t = 1
+ dt_tolerance = 1e-6
+ cdef int nv = ug.shape[0]
+ cdef int nshells = shells.shape[0]
+ cdef np.ndarray[np.float64_t, ndim=1] u = np.empty((3,), dtype=np.float64)
+ # Copy things into temporary location for passing between functions
+ for vi in range(nv):
+ for i in range(3): u[i] = ug[vi, i]
+ integrate_ray(u, v, left_edge, right_edge, dx,
+ nshells, vi, data, shells, image)
+
+ at cython.wraparound(False)
+ at cython.boundscheck(False)
+def integrate_ray(np.ndarray[np.float64_t, ndim=1] u,
+ np.ndarray[np.float64_t, ndim=1] v,
+ np.ndarray[np.float64_t, ndim=1] left_edge,
+ np.ndarray[np.float64_t, ndim=1] right_edge,
+ np.ndarray[np.float64_t, ndim=1] dx,
+ int nshells, int ind,
+ np.ndarray[np.float64_t, ndim=3] data,
+ np.ndarray[np.float64_t, ndim=2] shells,
+ np.ndarray[np.float64_t, ndim=2] image):
+ cdef int step[3], x, y, i, n
+ cdef np.float64_t intersect_t = 1
+ cdef np.float64_t dt_tolerance = 1e-6
+ cdef np.float64_t tl, tr, enter_t, exit_t
+ cdef np.int64_t cur_ind[3]
+ cdef np.float64_t tdelta[3]
+ cdef np.float64_t tmax[3]
+ cdef np.float64_t intersect[3]
+ cdef np.float64_t dt, dv
+ cdef np.float64_t dist, alpha
+ cdef np.float64_t one = 1.0
+ cdef int dims[3]
+ cdef np.float64_t rgba[4], temp_x, temp_y
+ for i in range(3):
+ # As long as we're iterating, set some other stuff, too
+ dims[i] = data.shape[i]
+ if(v[i] < 0): step[i] = -1
+ else: step[i] = 1
+ x = (i+1)%3
+ y = (i+2)%3
+ tl = (left_edge[i] - u[i])/v[i]
+ tr = (right_edge[i] - u[i])/v[i]
+ temp_x = (u[x] + tl*v[x])
+ temp_y = (u[y] + tl*v[y])
+ if (left_edge[x] <= temp_x) and (temp_x <= right_edge[x]) and \
+ (left_edge[y] <= temp_y) and (temp_y <= right_edge[y]) and \
+ (0.0 <= tl) and (tl < intersect_t):
+ intersect_t = tl
+ temp_x = (u[x] + tr*v[x])
+ temp_y = (u[y] + tr*v[y])
+ if (left_edge[x] <= temp_x) and (temp_x <= right_edge[x]) and \
+ (left_edge[y] <= temp_y) and (temp_y <= right_edge[y]) and \
+ (0.0 <= tr) and (tr < intersect_t):
+ intersect_t = tr
+ # if fully enclosed
+ if (left_edge[0] <= u[0] <= right_edge[0]) and \
+ (left_edge[1] <= u[1] <= right_edge[1]) and \
+ (left_edge[2] <= u[2] <= right_edge[2]):
+ intersect_t = 0.0
+ if not (0 <= intersect_t <= 1):
+ #print "Returning: intersect_t ==", intersect_t
+ return
+ # Now get the indices of the intersection
+ for i in range(3): intersect[i] = u[i] + intersect_t * v[i]
+ cdef int ncells = 0
+ for i in range(3):
+ cur_ind[i] = np.floor((intersect[i] + 1e-8*dx[i] - left_edge[i])/dx[i])
+ tmax[i] = (((cur_ind[i]+step[i])*dx[i])+left_edge[i]-u[i])/v[i]
+ if cur_ind[i] == dims[i] and step[i] < 0:
+ cur_ind[i] = dims[i] - 1
+ if step[i] > 0: tmax[i] = (((cur_ind[i]+1)*dx[i])+left_edge[i]-u[i])/v[i]
+ if step[i] < 0: tmax[i] = (((cur_ind[i]+0)*dx[i])+left_edge[i]-u[i])/v[i]
+ tdelta[i] = (dx[i]/v[i])
+ if tdelta[i] < 0: tdelta[i] *= -1
+ # The variable intersect contains the point we first pierce the grid
+ enter_t = intersect_t
+ if (not (0 <= cur_ind[0] < dims[0])) or \
+ (not (0 <= cur_ind[1] < dims[1])) or \
+ (not (0 <= cur_ind[2] < dims[2])):
+ #print "Returning: cur_ind", cur_ind[0], cur_ind[1], cur_ind[2]
+ #print " dims: ", dims[0], dims[1], dims[2]
+ #print " intersect:", intersect[0], intersect[1], intersect[2]
+ #print " intersect:", intersect_t
+ #print " u :", u[0], u[1], u[2]
+ #
+ return
+ #print cur_ind[0], dims[0], cur_ind[1], dims[1], cur_ind[2], dims[2]
+ dv = data[cur_ind[0], cur_ind[1], cur_ind[2]]
+ #dt = 1e300
+ while 1:
+ if image[ind,3] >= 1.0: break
+ if (not (0 <= cur_ind[0] < dims[0])) or \
+ (not (0 <= cur_ind[1] < dims[1])) or \
+ (not (0 <= cur_ind[2] < dims[2])):
+ break
+ # Do our transfer here
+ for n in range(nshells):
+ dist = shells[n, 0] - dv
+ if dist < shells[n,1]:
+ dist = exp(-dist/8.0)
+ alpha = (1.0 - shells[n,5])*shells[n,5]#*dt
+ image[ind,0] += alpha*shells[n,2]*dist
+ image[ind,1] += alpha*shells[n,3]*dist
+ image[ind,2] += alpha*shells[n,4]*dist
+ image[ind,3] += alpha*shells[n,5]*dist
+ #image[ind,i] += rgba[i]*dist*rgba[3]/dt
+ #print rgba[i], image[ind,i], dist, dt
+ break
+ if (tmax[0] > 1.0) and (tmax[1] > 1.0) and (tmax[2] > 1.0):
+ dt = 1.0 - enter_t
+ break
+ if tmax[0] < tmax[1]:
+ if tmax[0] < tmax[2]:
+ dt = tmax[0] - enter_t
+ enter_t = tmax[0]
+ tmax[0] += tdelta[0]
+ cur_ind[0] += step[0]
+ else:
+ dt = tmax[2] - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
+ else:
+ if tmax[1] < tmax[2]:
+ dt = tmax[1] - enter_t
+ enter_t = tmax[1]
+ tmax[1] += tdelta[1]
+ cur_ind[1] += step[1]
+ else:
+ dt = tmax[2] - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ cur_ind[2] += step[2]
+ dv = data[cur_ind[0], cur_ind[1], cur_ind[2]]
Added: trunk/yt/_amr_utils/VolumeIntegrator.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/_amr_utils/VolumeIntegrator.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,384 @@
+"""
+Simple integrators for the radiative transfer equation
+
+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 np
+cimport numpy as np
+cimport cython
+from stdlib cimport malloc, free, abs
+
+cdef inline int imax(int i0, int i1):
+ if i0 > i1: return i0
+ return i1
+
+cdef inline np.float64_t fmax(np.float64_t f0, np.float64_t f1):
+ if f0 > f1: return f0
+ return f1
+
+cdef inline int imin(int i0, int i1):
+ if i0 < i1: return i0
+ return i1
+
+cdef inline np.float64_t fmin(np.float64_t f0, np.float64_t f1):
+ if f0 < f1: return f0
+ return f1
+
+cdef inline int iclip(int i, int a, int b):
+ return imin(imax(i, a), b)
+
+cdef inline np.float64_t fclip(np.float64_t f,
+ np.float64_t a, np.float64_t b):
+ return fmin(fmax(f, a), b)
+
+cdef extern from "math.h":
+ double exp(double x)
+ float expf(float x)
+ double floor(double x)
+ double ceil(double x)
+ double fmod(double x, double y)
+ double log2(double x)
+
+cdef extern from "FixedInterpolator.h":
+ np.float64_t fast_interpolate(int *ds, int *ci, np.float64_t *dp,
+ np.float64_t *data)
+cdef extern from "FixedInterpolator.h":
+ np.float64_t trilinear_interpolate(int *ds, int *ci, np.float64_t *dp,
+ np.float64_t *data)
+
+cdef class VectorPlane
+
+cdef class TransferFunctionProxy:
+ cdef np.float64_t x_bounds[2]
+ cdef np.float64_t *vs[4]
+ cdef int nbins
+ cdef np.float64_t dbin
+ 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
+
+ cdef void eval_transfer(self, np.float64_t dt, np.float64_t dv,
+ np.float64_t *rgba):
+ cdef int i
+ cdef int bin_id
+ cdef np.float64_t tf, trgba[4], bv, dx, dy, dd,ta
+ dx = self.dbin
+
+ # get source alpha first
+ # First locate our points
+ bin_id = iclip(<int> floor((dv - self.x_bounds[0]) / dx),
+ 0, self.nbins-2)
+ # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+ bv = self.vs[3][bin_id] # This is x0
+ dy = self.vs[3][bin_id+1]-bv # dy
+ dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
+ # This is our final value for transfer function on the entering face
+ tf = bv+dd*(dy/dx)
+ ta = tf # Store the source alpha
+ for i in range(3):
+ # Recall that linear interpolation is y0 + (x-x0) * dx/dy
+ bv = self.vs[i][bin_id] # This is x0
+ dy = self.vs[i][bin_id+1]-bv # dy
+ dd = dv-(self.x_bounds[0] + bin_id * dx) # x - x0
+ # This is our final value for transfer function on the entering face
+ tf = bv+dd*(dy/dx)
+ # alpha blending
+ rgba[i] += (1. - rgba[3])*ta*tf*dt
+ #update alpha
+ rgba[3] += (1. - rgba[3])*ta*dt
+
+ # We should really do some alpha blending.
+ # Front to back blending is defined as:
+ # dst.rgb = dst.rgb + (1 - dst.a) * src.a * src.rgb
+ # dst.a = dst.a + (1 - dst.a) * src.a
+
+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
+ cdef public object ax_vec, ay_vec
+ cdef np.float64_t *x_vec, *y_vec
+
+ def __cinit__(self,
+ np.ndarray[np.float64_t, ndim=3] vp_pos,
+ np.ndarray[np.float64_t, ndim=1] vp_dir,
+ np.ndarray[np.float64_t, ndim=1] center,
+ bounds,
+ np.ndarray[np.float64_t, ndim=3] image,
+ np.ndarray[np.float64_t, ndim=1] x_vec,
+ np.ndarray[np.float64_t, ndim=1] y_vec):
+ cdef int i, j
+ self.avp_pos = vp_pos
+ self.avp_dir = vp_dir
+ self.acenter = center
+ self.aimage = image
+ self.ax_vec = x_vec
+ self.ay_vec = y_vec
+ self.vp_pos = <np.float64_t *> vp_pos.data
+ self.vp_dir = <np.float64_t *> vp_dir.data
+ self.center = <np.float64_t *> center.data
+ self.image = <np.float64_t *> image.data
+ self.x_vec = <np.float64_t *> x_vec.data
+ self.y_vec = <np.float64_t *> y_vec.data
+ self.nv = vp_pos.shape[0]
+ for i in range(4): self.bounds[i] = bounds[i]
+ self.pdx = (self.bounds[1] - self.bounds[0])/self.nv
+ self.pdy = (self.bounds[3] - self.bounds[2])/self.nv
+
+ cdef void get_start_stop(self, np.float64_t *ex, int *rv):
+ # Extrema need to be re-centered
+ cdef np.float64_t cx, cy
+ cx = cy = 0.0
+ for i in range(3):
+ cx += self.center[i] * self.x_vec[i]
+ cy += self.center[i] * self.y_vec[i]
+ rv[0] = <int> floor((ex[0] - cx - self.bounds[0])/self.pdx)
+ rv[1] = rv[0] + <int> ceil((ex[1] - ex[0])/self.pdx)
+ rv[2] = <int> floor((ex[2] - cy - self.bounds[2])/self.pdy)
+ rv[3] = rv[2] + <int> ceil((ex[3] - ex[2])/self.pdy)
+
+ cdef inline void copy_into(self, np.float64_t *fv, np.float64_t *tv,
+ int i, int j, int nk):
+ # We know the first two dimensions of our from-vector, and our
+ # to-vector is flat and 'ni' long
+ cdef int k
+ for k in range(nk):
+ tv[k] = fv[(((k*self.nv)+j)*self.nv+i)]
+
+ cdef inline void copy_back(self, np.float64_t *fv, np.float64_t *tv,
+ int i, int j, int nk):
+ cdef int k
+ for k in range(nk):
+ tv[(((k*self.nv)+j)*self.nv+i)] = fv[k]
+
+cdef class PartitionedGrid:
+ cdef public object my_data
+ cdef public object LeftEdge
+ cdef public object RightEdge
+ 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 public np.float64_t min_dds
+ cdef int ns
+ cdef int dims[3]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def __cinit__(self,
+ 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
+ self.LeftEdge = left_edge
+ self.RightEdge = right_edge
+ for i in range(3):
+ self.left_edge[i] = left_edge[i]
+ self.right_edge[i] = right_edge[i]
+ self.dims[i] = dims[i]
+ self.dds[i] = (self.right_edge[i] - self.left_edge[i])/dims[i]
+ self.my_data = data
+ self.data = <np.float64_t*> data.data
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def cast_plane(self, TransferFunctionProxy tf, VectorPlane vp):
+ # This routine will iterate over all of the vectors and cast each in
+ # turn. Might benefit from a more sophisticated intersection check,
+ # like http://courses.csusm.edu/cs697exz/ray_box.htm
+ cdef int vi, vj, hit, i, ni, nj, nn
+ self.ns = 5 #* (1 + <int> log2(self.dds[0] / self.min_dds))
+ cdef int iter[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)
+ for i in range(4): iter[i] = iclip(iter[i], 0, vp.nv)
+ hit = 0
+ for vj in range(iter[0], iter[1]):
+ for vi 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, 4)
+ self.integrate_ray(v_pos, vp.vp_dir, rgba, tf)
+ vp.copy_back(rgba, vp.image, vi, vj, 4)
+ return hit
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef void calculate_extent(self, VectorPlane vp,
+ np.float64_t extrema[4]):
+ # We do this for all eight corners
+ cdef np.float64_t *edges[2], temp
+ edges[0] = self.left_edge
+ edges[1] = self.right_edge
+ extrema[0] = extrema[2] = 1e300; extrema[1] = extrema[3] = -1e300
+ cdef int i, j, k
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ # This should rotate it into the vector plane
+ temp = edges[i][0] * vp.x_vec[0]
+ temp += edges[j][1] * vp.x_vec[1]
+ temp += edges[k][2] * vp.x_vec[2]
+ if temp < extrema[0]: extrema[0] = temp
+ if temp > extrema[1]: extrema[1] = temp
+ temp = edges[i][0] * vp.y_vec[0]
+ temp += edges[j][1] * vp.y_vec[1]
+ temp += edges[k][2] * vp.y_vec[2]
+ if temp < extrema[2]: extrema[2] = temp
+ if temp > extrema[3]: extrema[3] = temp
+ #print extrema[0], extrema[1], extrema[2], extrema[3]
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ cdef int integrate_ray(self, np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t rgba[4],
+ TransferFunctionProxy tf):
+ cdef int cur_ind[3], step[3], x, y, i, n, flat_ind, hit, direction
+ cdef np.float64_t intersect_t = 1.0
+ cdef np.float64_t intersect[3], tmax[3], tdelta[3]
+ cdef np.float64_t enter_t, dist, alpha, dt
+ cdef np.float64_t tr, tl, temp_x, temp_y, dv
+ for i in range(3):
+ if (v_dir[i] < 0):
+ step[i] = -1
+ else:
+ step[i] = 1
+ x = (i+1) % 3
+ y = (i+2) % 3
+ tl = (self.left_edge[i] - v_pos[i])/v_dir[i]
+ tr = (self.right_edge[i] - v_pos[i])/v_dir[i]
+ temp_x = (v_pos[x] + tl*v_dir[x])
+ temp_y = (v_pos[y] + tl*v_dir[y])
+ if self.left_edge[x] <= temp_x and temp_x <= self.right_edge[x] and \
+ self.left_edge[y] <= temp_y and temp_y <= self.right_edge[y] and \
+ 0.0 <= tl and tl < intersect_t:
+ direction = i
+ intersect_t = tl
+ temp_x = (v_pos[x] + tr*v_dir[x])
+ temp_y = (v_pos[y] + tr*v_dir[y])
+ if self.left_edge[x] <= temp_x and temp_x <= self.right_edge[x] and \
+ self.left_edge[y] <= temp_y and temp_y <= self.right_edge[y] and \
+ 0.0 <= tr and tr < intersect_t:
+ direction = i
+ intersect_t = tr
+ if self.left_edge[0] <= v_pos[0] and v_pos[0] <= self.right_edge[0] and \
+ self.left_edge[1] <= v_pos[1] and v_pos[1] <= self.right_edge[1] and \
+ self.left_edge[2] <= v_pos[2] and v_pos[2] <= self.right_edge[2]:
+ intersect_t = 0.0
+ if not ((0.0 <= intersect_t) and (intersect_t < 1.0)): return 0
+ for i in range(3):
+ intersect[i] = v_pos[i] + intersect_t * v_dir[i]
+ cur_ind[i] = <int> floor((intersect[i] +
+ step[i]*1e-8*self.dds[i] -
+ self.left_edge[i])/self.dds[i])
+ tmax[i] = (((cur_ind[i]+step[i])*self.dds[i])+
+ self.left_edge[i]-v_pos[i])/v_dir[i]
+ if cur_ind[i] == self.dims[i] and step[i] < 0:
+ cur_ind[i] = self.dims[i] - 1
+ if cur_ind[i] < 0 or cur_ind[i] >= self.dims[i]: return 0
+ if step[i] > 0:
+ tmax[i] = (((cur_ind[i]+1)*self.dds[i])
+ +self.left_edge[i]-v_pos[i])/v_dir[i]
+ if step[i] < 0:
+ tmax[i] = (((cur_ind[i]+0)*self.dds[i])
+ +self.left_edge[i]-v_pos[i])/v_dir[i]
+ tdelta[i] = (self.dds[i]/v_dir[i])
+ if tdelta[i] < 0: tdelta[i] *= -1
+ # We have to jumpstart our calculation
+ enter_t = intersect_t
+ while 1:
+ # dims here is one less than the dimensions of the data,
+ # but we are tracing on the grid, not on the data...
+ if (not (0 <= cur_ind[0] < self.dims[0])) or \
+ (not (0 <= cur_ind[1] < self.dims[1])) or \
+ (not (0 <= cur_ind[2] < self.dims[2])):
+ break
+ hit += 1
+ if tmax[0] < tmax[1]:
+ if tmax[0] < tmax[2]:
+ self.sample_values(v_pos, v_dir, enter_t, tmax[0], cur_ind,
+ rgba, tf)
+ cur_ind[0] += step[0]
+ dt = fmin(tmax[0], 1.0) - enter_t
+ enter_t = tmax[0]
+ tmax[0] += tdelta[0]
+ else:
+ self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+ rgba, tf)
+ cur_ind[2] += step[2]
+ dt = fmin(tmax[2], 1.0) - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ else:
+ if tmax[1] < tmax[2]:
+ self.sample_values(v_pos, v_dir, enter_t, tmax[1], cur_ind,
+ rgba, tf)
+ cur_ind[1] += step[1]
+ dt = fmin(tmax[1], 1.0) - enter_t
+ enter_t = tmax[1]
+ tmax[1] += tdelta[1]
+ else:
+ self.sample_values(v_pos, v_dir, enter_t, tmax[2], cur_ind,
+ rgba, tf)
+ cur_ind[2] += step[2]
+ dt = fmin(tmax[2], 1.0) - enter_t
+ enter_t = tmax[2]
+ tmax[2] += tdelta[2]
+ if enter_t > 1.0: break
+ return hit
+
+ cdef void sample_values(self,
+ np.float64_t v_pos[3],
+ np.float64_t v_dir[3],
+ np.float64_t enter_t,
+ np.float64_t exit_t,
+ int ci[3],
+ np.float64_t *rgba,
+ TransferFunctionProxy tf):
+ cdef np.float64_t cp[3], dp[3], temp, dt, t, dv
+ cdef int dti, i
+ dt = (exit_t - enter_t) / (self.ns-1) # five samples, so divide by four
+ for dti in range(self.ns - 1):
+ t = enter_t + dt * dti
+ for i in range(3):
+ cp[i] = v_pos[i] + t * v_dir[i]
+ dp[i] = fclip(fmod(cp[i], self.dds[i])/self.dds[i], 0, 1.0)
+ dv = trilinear_interpolate(self.dims, ci, dp, self.data)
+ tf.eval_transfer(dt, dv, rgba)
Added: trunk/yt/amr_utils.pyx
==============================================================================
--- (empty file)
+++ trunk/yt/amr_utils.pyx Fri Nov 20 19:54:16 2009
@@ -0,0 +1,43 @@
+"""
+Container file to hold all our Cython routines. This is to avoid problems with
+static linking.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008 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/>.
+"""
+
+#cython embedsignature=True
+#cython cdivision=True
+#cython always_allow_keywords=True
+
+# Set up some imports
+import numpy as np
+cimport numpy as np
+cimport cython
+
+# We include all of our files
+
+include "_amr_utils/DepthFirstOctree.pyx"
+include "_amr_utils/Interpolators.pyx"
+include "_amr_utils/PointsInVolume.pyx"
+include "_amr_utils/RayIntegrators.pyx"
+include "_amr_utils/VolumeIntegrator.pyx"
+include "_amr_utils/CICDeposit.pyx"
\ No newline at end of file
Modified: trunk/yt/extensions/volume_rendering/__init__.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/__init__.py (original)
+++ trunk/yt/extensions/volume_rendering/__init__.py Fri Nov 20 19:54:16 2009
@@ -26,7 +26,7 @@
import numpy as na
from TransferFunction import TransferFunction, ColorTransferFunction
-from VolumeIntegrator import PartitionedGrid, VectorPlane, \
+from yt.amr_utils import PartitionedGrid, VectorPlane, \
TransferFunctionProxy
from grid_partitioner import partition_all_grids, partition_grid, \
export_partitioned_grids, \
Modified: trunk/yt/extensions/volume_rendering/setup.py
==============================================================================
--- trunk/yt/extensions/volume_rendering/setup.py (original)
+++ trunk/yt/extensions/volume_rendering/setup.py Fri Nov 20 19:54:16 2009
@@ -4,12 +4,11 @@
import os.path
+os.system("cython -a yt/extensions/volume_rendering/VolumeIntegrator.pyx")
+
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
config = Configuration('volume_rendering',parent_package,top_path)
config.make_config_py() # installs __config__.py
config.make_svn_version_py()
- config.add_extension("VolumeIntegrator",
- ["yt/extensions/volume_rendering/VolumeIntegrator.c",
- "yt/extensions/volume_rendering/FixedInterpolator.c",])
return config
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Fri Nov 20 19:54:16 2009
@@ -28,7 +28,6 @@
data_object_registry = {}
from yt.lagos import *
-import PointsInVolume as PV
def restore_grid_state(func):
"""
@@ -1168,7 +1167,7 @@
def _get_point_indices(self, grid):
if self._pixelmask.max() == 0: return []
- k = PV.PointsInVolume(self._coord, self._pixelmask,
+ k = amr_utils.PointsInVolume(self._coord, self._pixelmask,
grid.LeftEdge, grid.RightEdge,
grid.child_mask, just_one(grid['dx']))
return k
Modified: trunk/yt/lagos/HelperFunctions.py
==============================================================================
--- trunk/yt/lagos/HelperFunctions.py (original)
+++ trunk/yt/lagos/HelperFunctions.py Fri Nov 20 19:54:16 2009
@@ -25,7 +25,6 @@
"""
from yt.lagos import *
-import Interpolators as IT
class UnilinearFieldInterpolator:
def __init__(self, table, boundaries, field_names):
@@ -46,7 +45,7 @@
raise ValueError
my_vals = na.zeros(x_vals.shape, dtype='float64')
- IT.UnilinearlyInterpolate(self.table, x_vals, self.x_bins, x_i, my_vals)
+ amr_utils.UnilinearlyInterpolate(self.table, x_vals, self.x_bins, x_i, my_vals)
return my_vals
class BilinearFieldInterpolator:
@@ -77,7 +76,7 @@
y_i = na.minimum(na.maximum(y_i,0), len(self.y_bins)-2)
my_vals = na.zeros(x_vals.shape, dtype='float64')
- IT.BilinearlyInterpolate(self.table,
+ amr_utils.BilinearlyInterpolate(self.table,
x_vals, y_vals, self.x_bins, self.y_bins,
x_i, y_i, my_vals)
return my_vals.reshape(orig_shape)
@@ -115,7 +114,7 @@
z_i = na.minimum(na.maximum(z_i,0), len(self.z_bins)-2)
my_vals = na.zeros(x_vals.shape, dtype='float64')
- IT.TrilinearlyInterpolate(self.table,
+ amr_utils.TrilinearlyInterpolate(self.table,
x_vals, y_vals, z_vals,
self.x_bins, self.y_bins, self.z_bins,
x_i, y_i, z_i, my_vals)
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Fri Nov 20 19:54:16 2009
@@ -55,6 +55,8 @@
import time
+import yt.amr_utils as amr_utils
+
from Cosmology import *
from EnzoCosmology import *
Modified: trunk/yt/lagos/setup.py
==============================================================================
--- trunk/yt/lagos/setup.py (original)
+++ trunk/yt/lagos/setup.py Fri Nov 20 19:54:16 2009
@@ -20,12 +20,9 @@
config.make_config_py() # installs __config__.py
config.make_svn_version_py()
config.add_extension("PointCombine", "yt/lagos/PointCombine.c", libraries=["m"])
- config.add_extension("RTIntegrator", "yt/lagos/RTIntegrator.c", libraries=["m"])
- config.add_extension("Interpolators", "yt/lagos/Interpolators.c")
- config.add_extension("PointsInVolume", "yt/lagos/PointsInVolume.c")
- #config.add_extension("DepthFirstOctree", "yt/lagos/DepthFirstOctree.c")
config.add_subpackage("hop")
config.add_subpackage("fof")
+ config.add_subpackage("parallelHOP")
H5dir = check_for_hdf5()
if H5dir is not None:
include_dirs=[os.path.join(H5dir,"include")]
@@ -35,6 +32,6 @@
libraries=["m","hdf5"],
library_dirs=library_dirs, include_dirs=include_dirs)
# Uncomment the next two lines if you want particle_density support
- #config.add_extension("cic_deposit", ["yt/lagos/enzo_routines/cic_deposit.pyf",
- # "yt/lagos/enzo_routines/cic_deposit.f"])
+ config.add_extension("cic_deposit", ["yt/lagos/enzo_routines/cic_deposit.pyf",
+ "yt/lagos/enzo_routines/cic_deposit.f"])
return config
Modified: trunk/yt/setup.py
==============================================================================
--- trunk/yt/setup.py (original)
+++ trunk/yt/setup.py Fri Nov 20 19:54:16 2009
@@ -9,6 +9,10 @@
config.add_subpackage('fido')
config.add_subpackage('reason')
config.add_subpackage('extensions')
+ config.add_extension("amr_utils",
+ ["yt/amr_utils.c", "yt/_amr_utils/FixedInterpolator.c"],
+ include_dirs=["yt/_amr_utils/"],
+ libraries=["m"])
config.make_config_py()
config.make_svn_version_py()
return config
More information about the yt-svn
mailing list