Commits

Matthew Turk committed 2d81eb7

Split off the marching cubes and isocontour code; removed some imports to make
sure we're always using the new refactored versions of things.

  • Participants
  • Parent commits 060e702
  • Branches geometry_handling

Comments (0)

Files changed (7)

yt/utilities/_amr_utils/fixed_interpolator.pxd

+"""
+Fixed interpolator includes
+
+Author: Matthew Turk <matthewturk@gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+  Copyright (C) 2011 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/>.
+"""
+
+cimport numpy as np
+
+
+cdef extern from "FixedInterpolator.h":
+    np.float64_t fast_interpolate(int ds[3], int ci[3], np.float64_t dp[3],
+                                  np.float64_t *data) nogil
+    np.float64_t offset_interpolate(int ds[3], np.float64_t dp[3],
+                                    np.float64_t *data) nogil
+    np.float64_t trilinear_interpolate(int ds[3], int ci[3], np.float64_t dp[3],
+                                       np.float64_t *data) nogil
+    void eval_gradient(int ds[3], np.float64_t dp[3], np.float64_t *data,
+                       np.float64_t grad[3]) nogil
+    void offset_fill(int *ds, np.float64_t *data, np.float64_t *gridval) nogil
+    void vertex_interp(np.float64_t v1, np.float64_t v2, np.float64_t isovalue,
+                       np.float64_t vl[3], np.float64_t dds[3],
+                       np.float64_t x, np.float64_t y, np.float64_t z,
+                       int vind1, int vind2) nogil
+

yt/utilities/_amr_utils/grid_traversal.pyx

 from fp_utils cimport imax, fmax, imin, fmin, iclip, fclip
 from field_interpolation_tables cimport \
     FieldInterpolationTable, FIT_initialize_table, FIT_eval_transfer
+from fixed_interpolator cimport *
 
 from cython.parallel import prange, parallel, threadid
 
                 int index[3],
                 void *data) nogil
 
-cdef extern from "FixedInterpolator.h":
-    np.float64_t fast_interpolate(int ds[3], int ci[3], np.float64_t dp[3],
-                                  np.float64_t *data) nogil
-    np.float64_t offset_interpolate(int ds[3], np.float64_t dp[3],
-                                    np.float64_t *data) nogil
-    np.float64_t trilinear_interpolate(int ds[3], int ci[3], np.float64_t dp[3],
-                                       np.float64_t *data) nogil
-    void eval_gradient(int ds[3], np.float64_t dp[3], np.float64_t *data,
-                       np.float64_t grad[3]) nogil
-    void offset_fill(int *ds, np.float64_t *data, np.float64_t *gridval) nogil
-    void vertex_interp(np.float64_t v1, np.float64_t v2, np.float64_t isovalue,
-                       np.float64_t vl[3], np.float64_t dds[3],
-                       np.float64_t x, np.float64_t y, np.float64_t z,
-                       int vind1, int vind2) nogil
-
 cdef struct VolumeContainer:
     int n_fields
     np.float64_t **data

yt/utilities/_amr_utils/marching_cubes.pyx

+"""
+Marching cubes implementation
+
+Author: Matthew Turk <matthewturk@gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+  Copyright (C) 2011 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/>.
+"""
+
+cimport numpy as np
+cimport cython
+import numpy as np
+from fp_utils cimport imax, fmax, imin, fmin, iclip, fclip
+from stdlib cimport malloc, free, abs
+from fixed_interpolator cimport *
+
+cdef struct Triangle:
+    Triangle *next
+    np.float64_t p[3][3]
+    np.float64_t val
+
+cdef struct TriangleCollection:
+    int count
+    Triangle *first
+    Triangle *current
+
+cdef Triangle *AddTriangle(Triangle *self,
+                    np.float64_t p0[3], np.float64_t p1[3], np.float64_t p2[3]):
+    cdef Triangle *nn = <Triangle *> malloc(sizeof(Triangle))
+    if self != NULL:
+        self.next = nn
+    cdef int i
+    for i in range(3):
+        nn.p[0][i] = p0[i]
+    for i in range(3):
+        nn.p[1][i] = p1[i]
+    for i in range(3):
+        nn.p[2][i] = p2[i]
+    nn.next = NULL
+    return nn
+
+cdef int CountTriangles(Triangle *first):
+    cdef int count = 0
+    cdef Triangle *this = first
+    while this != NULL:
+        count += 1
+        this = this.next
+    return count
+
+cdef void FillTriangleValues(np.ndarray[np.float64_t, ndim=1] values,
+                             Triangle *first):
+    cdef Triangle *this = first
+    cdef Triangle *last
+    cdef int i = 0
+    while this != NULL:
+        values[i] = this.val
+        i += 1
+        last = this
+        this = this.next
+
+cdef void WipeTriangles(Triangle *first):
+    cdef Triangle *this = first
+    cdef Triangle *last
+    while this != NULL:
+        last = this
+        this = this.next
+        free(last)
+
+cdef void FillAndWipeTriangles(np.ndarray[np.float64_t, ndim=2] vertices,
+                               Triangle *first):
+    cdef int count = 0
+    cdef Triangle *this = first
+    cdef Triangle *last
+    cdef int i, j
+    while this != NULL:
+        for i in range(3):
+            for j in range(3):
+                vertices[count, j] = this.p[i][j]
+            count += 1 # Do it at the end because it's an index
+        last = this
+        this = this.next
+        free(last)
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+cdef int march_cubes(
+                 np.float64_t gv[8], np.float64_t isovalue,
+                 np.float64_t dds[3],
+                 np.float64_t x, np.float64_t y, np.float64_t z,
+                 TriangleCollection *triangles):
+    cdef int *edge_table=[
+    0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
+    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+    0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
+    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+    0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
+    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+    0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
+    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+    0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
+    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
+    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
+    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
+    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
+    0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
+    0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
+    0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
+    0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
+    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
+    0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
+    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
+    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
+    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
+    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
+    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
+    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0   ]
+
+    cdef int **tri_table = \
+    [[-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1],
+    [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1],
+    [3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1],
+    [3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1],
+    [9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1],
+    [9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
+    [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1],
+    [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1],
+    [9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
+    [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1],
+    [3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1],
+    [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1],
+    [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1],
+    [4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1],
+    [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
+    [5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1],
+    [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1],
+    [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1],
+    [0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1],
+    [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1],
+    [10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1],
+    [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1],
+    [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1],
+    [5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1],
+    [9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1],
+    [0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1],
+    [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1],
+    [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1],
+    [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1],
+    [2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1],
+    [7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1],
+    [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1],
+    [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1],
+    [11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1],
+    [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1],
+    [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1],
+    [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1],
+    [11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
+    [1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1],
+    [9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1],
+    [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1],
+    [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
+    [0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1],
+    [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1],
+    [6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1],
+    [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1],
+    [6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1],
+    [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1],
+    [1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1],
+    [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1],
+    [6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1],
+    [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1],
+    [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1],
+    [3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1],
+    [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1],
+    [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1],
+    [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1],
+    [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1],
+    [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1],
+    [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1],
+    [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1],
+    [10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1],
+    [10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1],
+    [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1],
+    [1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1],
+    [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1],
+    [0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1],
+    [10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1],
+    [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1],
+    [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1],
+    [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1],
+    [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1],
+    [3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1],
+    [6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1],
+    [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1],
+    [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1],
+    [10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1],
+    [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1],
+    [7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1],
+    [7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1],
+    [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1],
+    [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1],
+    [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1],
+    [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1],
+    [0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1],
+    [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
+    [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
+    [2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1],
+    [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1],
+    [7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1],
+    [2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1],
+    [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1],
+    [10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1],
+    [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1],
+    [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1],
+    [7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1],
+    [6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1],
+    [8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1],
+    [6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1],
+    [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1],
+    [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1],
+    [8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1],
+    [0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1],
+    [1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1],
+    [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1],
+    [10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1],
+    [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1],
+    [10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1],
+    [5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
+    [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1],
+    [9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1],
+    [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1],
+    [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1],
+    [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1],
+    [7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1],
+    [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1],
+    [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1],
+    [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1],
+    [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1],
+    [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1],
+    [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1],
+    [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1],
+    [6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1],
+    [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1],
+    [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1],
+    [6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1],
+    [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1],
+    [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1],
+    [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1],
+    [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1],
+    [9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1],
+    [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1],
+    [1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1],
+    [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1],
+    [0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1],
+    [5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1],
+    [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1],
+    [11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1],
+    [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1],
+    [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1],
+    [2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1],
+    [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1],
+    [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1],
+    [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1],
+    [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1],
+    [9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1],
+    [9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1],
+    [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1],
+    [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1],
+    [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1],
+    [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1],
+    [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1],
+    [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1],
+    [9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1],
+    [5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1],
+    [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1],
+    [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1],
+    [8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1],
+    [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1],
+    [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1],
+    [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1],
+    [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1],
+    [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1],
+    [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1],
+    [11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1],
+    [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1],
+    [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1],
+    [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1],
+    [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1],
+    [1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1],
+    [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1],
+    [4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1],
+    [0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1],
+    [3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1],
+    [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1],
+    [0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1],
+    [9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1],
+    [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],
+    [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]]
+    cdef np.float64_t vertlist[12][3]
+    cdef int cubeindex = 0
+    cdef int n
+    cdef int nt = 0
+    for n in range(8):
+        if gv[n] < isovalue:
+            cubeindex |= (1 << n)
+    if edge_table[cubeindex] == 0:
+        return 0
+    if (edge_table[cubeindex] & 1): # 0,0,0 with 1,0,0
+        vertex_interp(gv[0], gv[1], isovalue, vertlist[0],
+                      dds, x, y, z, 0, 1)
+    if (edge_table[cubeindex] & 2): # 1,0,0 with 1,1,0
+        vertex_interp(gv[1], gv[2], isovalue, vertlist[1],
+                      dds, x, y, z, 1, 2)
+    if (edge_table[cubeindex] & 4): # 1,1,0 with 0,1,0
+        vertex_interp(gv[2], gv[3], isovalue, vertlist[2],
+                      dds, x, y, z, 2, 3)
+    if (edge_table[cubeindex] & 8): # 0,1,0 with 0,0,0
+        vertex_interp(gv[3], gv[0], isovalue, vertlist[3],
+                      dds, x, y, z, 3, 0)
+    if (edge_table[cubeindex] & 16): # 0,0,1 with 1,0,1
+        vertex_interp(gv[4], gv[5], isovalue, vertlist[4],
+                      dds, x, y, z, 4, 5)
+    if (edge_table[cubeindex] & 32): # 1,0,1 with 1,1,1
+        vertex_interp(gv[5], gv[6], isovalue, vertlist[5],
+                      dds, x, y, z, 5, 6)
+    if (edge_table[cubeindex] & 64): # 1,1,1 with 0,1,1
+        vertex_interp(gv[6], gv[7], isovalue, vertlist[6],
+                      dds, x, y, z, 6, 7)
+    if (edge_table[cubeindex] & 128): # 0,1,1 with 0,0,1
+        vertex_interp(gv[7], gv[4], isovalue, vertlist[7],
+                      dds, x, y, z, 7, 4)
+    if (edge_table[cubeindex] & 256): # 0,0,0 with 0,0,1
+        vertex_interp(gv[0], gv[4], isovalue, vertlist[8],
+                      dds, x, y, z, 0, 4)
+    if (edge_table[cubeindex] & 512): # 1,0,0 with 1,0,1
+        vertex_interp(gv[1], gv[5], isovalue, vertlist[9],
+                      dds, x, y, z, 1, 5)
+    if (edge_table[cubeindex] & 1024): # 1,1,0 with 1,1,1
+        vertex_interp(gv[2], gv[6], isovalue, vertlist[10],
+                      dds, x, y, z, 2, 6)
+    if (edge_table[cubeindex] & 2048): # 0,1,0 with 0,1,1
+        vertex_interp(gv[3], gv[7], isovalue, vertlist[11],
+                      dds, x, y, z, 3, 7)
+    n = 0
+    while 1:
+        triangles.current = AddTriangle(triangles.current,
+                    vertlist[tri_table[cubeindex][n  ]],
+                    vertlist[tri_table[cubeindex][n+1]],
+                    vertlist[tri_table[cubeindex][n+2]])
+        triangles.count += 1
+        nt += 1
+        if triangles.first == NULL:
+            triangles.first = triangles.current
+        n += 3
+        if tri_table[cubeindex][n] == -1: break
+    return nt
+    
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+def march_cubes_grid(np.float64_t isovalue,
+                     np.ndarray[np.float64_t, ndim=3] values,
+                     np.ndarray[np.int32_t, ndim=3] mask,
+                     np.ndarray[np.float64_t, ndim=1] left_edge,
+                     np.ndarray[np.float64_t, ndim=1] dxs,
+                     obj_sample = None):
+    cdef int dims[3]
+    cdef int i, j, k, n, m, nt
+    cdef int offset
+    cdef np.float64_t gv[8], pos[3], point[3], idds[3]
+    cdef np.float64_t *intdata = NULL
+    cdef np.float64_t *sdata = NULL
+    cdef np.float64_t x, y, z, do_sample
+    cdef np.ndarray[np.float64_t, ndim=3] sample
+    cdef np.ndarray[np.float64_t, ndim=1] sampled
+    cdef TriangleCollection triangles
+    cdef Triangle *last, *current
+    if obj_sample is not None:
+        sample = obj_sample
+        sdata = <np.float64_t *> sample.data
+        do_sample = 1
+    else:
+        do_sample = 0
+    for i in range(3):
+        dims[i] = values.shape[i] - 1
+        idds[i] = 1.0 / dxs[i]
+    triangles.first = triangles.current = NULL
+    last = current = NULL
+    triangles.count = 0
+    cdef np.float64_t *data = <np.float64_t *> values.data
+    cdef np.float64_t *dds = <np.float64_t *> dxs.data
+    pos[0] = left_edge[0]
+    for i in range(dims[0]):
+        pos[1] = left_edge[1]
+        for j in range(dims[1]):
+            pos[2] = left_edge[2]
+            for k in range(dims[2]):
+                if mask[i,j,k] == 1:
+                    offset = i * (dims[1] + 1) * (dims[2] + 1) \
+                           + j * (dims[2] + 1) + k
+                    intdata = data + offset
+                    offset_fill(dims, intdata, gv)
+                    nt = march_cubes(gv, isovalue, dds, pos[0], pos[1], pos[2],
+                                &triangles)
+                    if do_sample == 1 and nt > 0:
+                        # At each triangle's center, sample our secondary field
+                        if last == NULL and triangles.first != NULL:
+                            current = triangles.first
+                            last = NULL
+                        elif last != NULL:
+                            current = last.next
+                        while current != NULL:
+                            for n in range(3):
+                                point[n] = 0.0
+                            for n in range(3):
+                                for m in range(3):
+                                    point[m] += (current.p[n][m]-pos[m])*idds[m]
+                            for n in range(3):
+                                point[n] /= 3.0
+                            current.val = offset_interpolate(dims, point,
+                                                             sdata + offset)
+                            last = current
+                            if current.next == NULL: break
+                            current = current.next
+                pos[2] += dds[2]
+            pos[1] += dds[1]
+        pos[0] += dds[0]
+    # Hallo, we are all done.
+    cdef np.ndarray[np.float64_t, ndim=2] vertices 
+    vertices = np.zeros((triangles.count*3,3), dtype='float64')
+    if do_sample == 1:
+        sampled = np.zeros(triangles.count, dtype='float64')
+        FillTriangleValues(sampled, triangles.first)
+        FillAndWipeTriangles(vertices, triangles.first)
+        return vertices, sampled
+    FillAndWipeTriangles(vertices, triangles.first)
+    return vertices
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+def march_cubes_grid_flux(
+                     np.float64_t isovalue,
+                     np.ndarray[np.float64_t, ndim=3] values,
+                     np.ndarray[np.float64_t, ndim=3] v1,
+                     np.ndarray[np.float64_t, ndim=3] v2,
+                     np.ndarray[np.float64_t, ndim=3] v3,
+                     np.ndarray[np.float64_t, ndim=3] flux_field,
+                     np.ndarray[np.int32_t, ndim=3] mask,
+                     np.ndarray[np.float64_t, ndim=1] left_edge,
+                     np.ndarray[np.float64_t, ndim=1] dxs):
+    cdef int dims[3]
+    cdef int i, j, k, n, m
+    cdef int offset
+    cdef np.float64_t gv[8]
+    cdef np.float64_t *intdata = NULL
+    cdef TriangleCollection triangles
+    cdef Triangle *current = NULL
+    cdef Triangle *last = NULL
+    cdef np.float64_t *data = <np.float64_t *> values.data
+    cdef np.float64_t *v1data = <np.float64_t *> v1.data
+    cdef np.float64_t *v2data = <np.float64_t *> v2.data
+    cdef np.float64_t *v3data = <np.float64_t *> v3.data
+    cdef np.float64_t *fdata = <np.float64_t *> flux_field.data
+    cdef np.float64_t *dds = <np.float64_t *> dxs.data
+    cdef np.float64_t flux = 0.0
+    cdef np.float64_t center[3], point[3], wval, temp, area, s
+    cdef np.float64_t cell_pos[3], fv[3], idds[3], normal[3]
+    for i in range(3):
+        dims[i] = values.shape[i] - 1
+        idds[i] = 1.0 / dds[i]
+    triangles.first = triangles.current = NULL
+    triangles.count = 0
+    cell_pos[0] = left_edge[0]
+    for i in range(dims[0]):
+        cell_pos[1] = left_edge[1]
+        for j in range(dims[1]):
+            cell_pos[2] = left_edge[2]
+            for k in range(dims[2]):
+                if mask[i,j,k] == 1:
+                    offset = i * (dims[1] + 1) * (dims[2] + 1) \
+                           + j * (dims[2] + 1) + k
+                    intdata = data + offset
+                    offset_fill(dims, intdata, gv)
+                    march_cubes(gv, isovalue, dds,
+                                cell_pos[0], cell_pos[1], cell_pos[2],
+                                &triangles)
+                    # Now our triangles collection has a bunch.  We now
+                    # calculate fluxes for each.
+                    if last == NULL and triangles.first != NULL:
+                        current = triangles.first
+                        last = NULL
+                    elif last != NULL:
+                        current = last.next
+                    while current != NULL:
+                        # Calculate the center of the triangle
+                        wval = 0.0
+                        for n in range(3):
+                            center[n] = 0.0
+                        for n in range(3):
+                            for m in range(3):
+                                point[m] = (current.p[n][m]-cell_pos[m])*idds[m]
+                            # Now we calculate the value at this point
+                            temp = offset_interpolate(dims, point, intdata)
+                            #print "something", temp, point[0], point[1], point[2]
+                            wval += temp
+                            for m in range(3):
+                                center[m] += temp * point[m]
+                        # Now we divide by our normalizing factor
+                        for n in range(3):
+                            center[n] /= wval
+                        # We have our center point of the triangle, in 0..1
+                        # coordinates.  So now we interpolate our three
+                        # fields.
+                        fv[0] = offset_interpolate(dims, center, v1data + offset)
+                        fv[1] = offset_interpolate(dims, center, v2data + offset)
+                        fv[2] = offset_interpolate(dims, center, v3data + offset)
+                        # We interpolate again the actual value data
+                        wval = offset_interpolate(dims, center, fdata + offset)
+                        # Now we have our flux vector and our field value!
+                        # We just need a normal vector with which we can
+                        # dot it.  The normal should be equal to the gradient
+                        # in the center of the triangle, or thereabouts.
+                        eval_gradient(dims, center, intdata, normal)
+                        temp = 0.0
+                        for n in range(3):
+                            temp += normal[n]*normal[n]
+                        # Take the negative, to ensure it points inwardly
+                        temp = -(temp**0.5)
+                        # Dump this somewhere for now
+                        temp = wval * (fv[0] * normal[0] +
+                                       fv[1] * normal[1] +
+                                       fv[2] * normal[2])/temp
+                        # Now we need the area of the triangle.  This will take
+                        # a lot of time to calculate compared to the rest.
+                        # We use Heron's formula.
+                        for n in range(3):
+                            fv[n] = 0.0
+                        for n in range(3):
+                            fv[0] += (current.p[0][n] - current.p[2][n])**2.0
+                            fv[1] += (current.p[1][n] - current.p[0][n])**2.0
+                            fv[2] += (current.p[2][n] - current.p[1][n])**2.0
+                        s = 0.0
+                        for n in range(3):
+                            fv[n] = fv[n]**0.5
+                            s += 0.5 * fv[n]
+                        area = (s*(s-fv[0])*(s-fv[1])*(s-fv[2]))
+                        area = area**0.5
+                        flux += temp*area
+                        last = current
+                        if current.next == NULL: break
+                        current = current.next
+                cell_pos[2] += dds[2]
+            cell_pos[1] += dds[1]
+        cell_pos[0] += dds[0]
+    # Hallo, we are all done.
+    WipeTriangles(triangles.first)
+    return flux
+

yt/utilities/_amr_utils/setup.py

     config.add_extension("Interpolators", 
                 ["yt/utilities/_amr_utils/Interpolators.pyx"],
                 libraries=["m"], depends=["yt/utilities/_amr_utils/fp_utils.pxd"])
+    config.add_extension("marching_cubes", 
+                ["yt/utilities/_amr_utils/marching_cubes.pyx",
+                 "yt/utilities/_amr_utils/FixedInterpolator.c"],
+                libraries=["m"],
+                depends=["yt/utilities/_amr_utils/fp_utils.pxd",
+                         "yt/utilities/_amr_utils/fixed_interpolator.pxd",
+                         "yt/utilities/_amr_utils/FixedInterpolator.h",
+                ])
     config.add_extension("misc_utilities", 
                 ["yt/utilities/_amr_utils/misc_utilities.pyx"],
                 libraries=["m"], depends=["yt/utilities/_amr_utils/fp_utils.pxd"])
           )
     config.add_extension("grid_traversal", 
                ["yt/utilities/_amr_utils/grid_traversal.pyx",
-                "yt/utilities/_amr_utils/FixedInterpolator.c"],
+                "yt/utilities/_amr_utils/FixedInterpolator.c",
+                "yt/utilities/_amr_utils/kdtree.c"],
                include_dirs=["yt/utilities/_amr_utils/"],
                libraries=["m"], 
                extra_compile_args=['-fopenmp'],
                extra_link_args=['-fopenmp'],
                depends = ["yt/utilities/_amr_utils/VolumeIntegrator.pyx",
                           "yt/utilities/_amr_utils/fp_utils.pxd",
+                          "yt/utilities/_amr_utils/kdtree.h",
                           "yt/utilities/_amr_utils/FixedInterpolator.h",
+                          "yt/utilities/_amr_utils/fixed_interpolator.pxd",
                           ]
           )
     if os.environ.get("GPERFTOOLS", "no").upper() != "NO":

yt/utilities/amr_utils.py

 from ._amr_utils.PointsInVolume import *
 from ._amr_utils.QuadTree import *
 from ._amr_utils.RayIntegrators import *
-from ._amr_utils.VolumeIntegrator import *
+from ._amr_utils.grid_traversal import *
+from ._amr_utils.marching_cubes import *
+#from ._amr_utils.VolumeIntegrator import *

yt/visualization/volume_rendering/api.py

                              PlanckTransferFunction, \
                              MultiVariateTransferFunction, \
                              ProjectionTransferFunction
-from yt.utilities.amr_utils import PartitionedGrid, VectorPlane, \
-    TransferFunctionProxy
 from grid_partitioner import HomogenizedVolume, \
                              export_partitioned_grids, \
                              import_partitioned_grids

yt/visualization/volume_rendering/camera.py

 from .grid_partitioner import HomogenizedVolume
 from .transfer_functions import ProjectionTransferFunction
 
-from yt.utilities.amr_utils import TransferFunctionProxy, VectorPlane, \
-    arr_vec2pix_nest, arr_pix2vec_nest, AdaptiveRaySource, \
-    arr_ang2pix_nest
+#from yt.utilities.amr_utils import \
+#    arr_vec2pix_nest, arr_pix2vec_nest, AdaptiveRaySource, \
+#    arr_ang2pix_nest
 from yt.visualization.image_writer import write_bitmap
 from yt.data_objects.data_containers import data_object_registry
 from yt.utilities.parallel_tools.parallel_analysis_interface import \