summaryrefslogtreecommitdiff
path: root/src/silx/math/marchingcubes.pyx
diff options
context:
space:
mode:
Diffstat (limited to 'src/silx/math/marchingcubes.pyx')
-rw-r--r--src/silx/math/marchingcubes.pyx246
1 files changed, 246 insertions, 0 deletions
diff --git a/src/silx/math/marchingcubes.pyx b/src/silx/math/marchingcubes.pyx
new file mode 100644
index 0000000..0409691
--- /dev/null
+++ b/src/silx/math/marchingcubes.pyx
@@ -0,0 +1,246 @@
+# coding: utf-8
+# /*##########################################################################
+#
+# Copyright (c) 2015-2018 European Synchrotron Radiation Facility
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#
+# ###########################################################################*/
+"""This module provides marching cubes implementation.
+
+It provides a :class:`MarchingCubes` class allowing to build an isosurface
+from data provided as a 3D data set or slice by slice.
+"""
+
+__authors__ = ["T. Vincent"]
+__license__ = "MIT"
+__date__ = "16/08/2017"
+
+
+import numpy
+cimport numpy as cnumpy
+cimport cython
+
+cimport silx.math.mc as mc
+
+
+# From numpy_common.pxi to avoid warnings while compiling C code
+# See this thread:
+# https://mail.python.org/pipermail//cython-devel/2012-March/002137.html
+cdef extern from *:
+ bint FALSE "0"
+ void import_array()
+ void import_umath()
+
+if FALSE:
+ import_array()
+ import_umath()
+
+
+cdef class MarchingCubes:
+ """Compute isosurface using marching cubes algorithm.
+
+ It builds a surface from a 3D scalar dataset as a 3D contour at a
+ given value.
+ The resulting surface is not topologically correct.
+
+ See: http://paulbourke.net/geometry/polygonise/
+
+ Lorensen, W. E. and Cline, H. E. Marching cubes: A high resolution 3D
+ surface construction algorithm. Computer Graphics, 21, 4 (July 1987).
+ ACM, 163-169.
+
+ Generated vertex and normal coordinates are in the same order
+ as input array, i.e., (dim 0, dim 1, dim 2).
+
+ Expected indices in memory of a (2, 2, 2) dataset:
+
+ dim 0 (depth)
+ |
+ |
+ 4 +------+ 5
+ /| /|
+ / | / |
+ 6 +------+ 7|
+ | | | |
+ |0 +---|--+ 1 --- dim 2 (width)
+ | / | /
+ |/ |/
+ 2 +------+ 3
+ /
+ /
+ dim 1 (height)
+
+ Example with a 3D data set:
+
+ >>> vertices, normals, indices = MarchingCubes(data, isolevel=1.)
+
+ Example of code for processing a list of images:
+
+ >>> mc = MarchingCubes(isolevel=1.) # Create object with iso-level=1
+ >>> previous_image = images[0]
+ >>> for image in images[1:]:
+ ... mc.process_image(previous_image, image) # Process one slice
+ ... previous_image = image
+
+ >>> vertices = mc.get_vertices() # Array of vertex positions
+ >>> normals = mc.get_normals() # Array of normals
+ >>> triangle_indices = mc.get_indices() # Array of indices of vertices
+
+ :param data: 3D dataset of float32 or None
+ :type data: numpy.ndarray of float32 of dimension 3
+ :param float isolevel: The value for which to generate the isosurface
+ :param bool invert_normals:
+ True (default) for normals oriented in direction of gradient descent
+ :param sampling: Sampling along each dimension (depth, height, width)
+ """
+ cdef mc.MarchingCubes[float, float] * c_mc # Pointer to the C++ instance
+
+ def __cinit__(self, data=None, isolevel=None,
+ invert_normals=True, sampling=(1, 1, 1)):
+ self.c_mc = new mc.MarchingCubes[float, float](isolevel)
+ self.c_mc.invert_normals = bool(invert_normals)
+ self.c_mc.sampling[0] = sampling[0]
+ self.c_mc.sampling[1] = sampling[1]
+ self.c_mc.sampling[2] = sampling[2]
+
+ if data is not None:
+ self.process(data)
+
+ def __dealloc__(self):
+ del self.c_mc
+
+ def __getitem__(self, key):
+ """Allows one to unpack object as a single liner:
+
+ vertices, normals, indices = MarchingCubes(...)
+ """
+ if key == 0:
+ return self.get_vertices()
+ elif key == 1:
+ return self.get_normals()
+ elif key == 2:
+ return self.get_indices()
+ else:
+ raise IndexError("Index out of range")
+
+ def process(self, data):
+ """Compute an isosurface from a 3D scalar field.
+
+ This builds vertices, normals and indices arrays.
+ Vertices and normals coordinates are in the same order as input array,
+ i.e., (dim 0, dim 1, dim 2).
+
+ :param numpy.ndarray data: 3D scalar field
+ """
+ # Make sure data is a 3D contiguous array of native endian float32
+ data = numpy.ascontiguousarray(data, dtype='=f4')
+ assert data.ndim == 3
+ cdef float[:] c_data = numpy.ravel(data)
+ cdef unsigned int depth, height, width
+
+ depth = data.shape[0]
+ height = data.shape[1]
+ width = data.shape[2]
+
+ self.c_mc.process(&c_data[0], depth, height, width)
+
+ def process_slice(self, slice0, slice1):
+ """Process a new slice to build the isosurface.
+
+ :param numpy.ndarray slice0: Slice previously provided as slice1.
+ :param numpy.ndarray slice1: Slice to process.
+ """
+ # Make sure slices are 2D contiguous arrays of native endian float32
+ slice0 = numpy.ascontiguousarray(slice0, dtype='=f4')
+ assert slice0.ndim == 2
+ slice1 = numpy.ascontiguousarray(slice1, dtype='=f4')
+ assert slice1.ndim == 2
+
+ assert slice0.shape[0] == slice1.shape[0]
+ assert slice0.shape[1] == slice1.shape[1]
+
+ cdef float[:] c_slice0 = numpy.ravel(slice0)
+ cdef float[:] c_slice1 = numpy.ravel(slice1)
+
+ if self.c_mc.depth == 0:
+ # Starts a new isosurface, bootstrap with slice size
+ self.c_mc.set_slice_size(slice1.shape[0], slice1.shape[1])
+
+ assert slice1.shape[0] == self.c_mc.height
+ assert slice1.shape[1] == self.c_mc.width
+
+ self.c_mc.process_slice(&c_slice0[0], &c_slice1[0])
+
+ def finish_process(self):
+ """Clear internal cache after processing slice by slice."""
+ self.c_mc.finish_process()
+
+ def reset(self):
+ """Reset internal resources including computed isosurface info."""
+ self.c_mc.reset()
+
+ @cython.embedsignature(False)
+ @property
+ def shape(self):
+ """The shape of the processed scalar field (depth, height, width)."""
+ return self.c_mc.depth, self.c_mc.height, self.c_mc.width
+
+ @cython.embedsignature(False)
+ @property
+ def sampling(self):
+ """The sampling over each dimension (depth, height, width).
+
+ Default: 1, 1, 1
+ """
+ return (self.c_mc.sampling[0],
+ self.c_mc.sampling[1],
+ self.c_mc.sampling[2])
+
+ @cython.embedsignature(False)
+ @property
+ def isolevel(self):
+ """The iso-level at which to generate the isosurface"""
+ return self.c_mc.isolevel
+
+ @cython.embedsignature(False)
+ @property
+ def invert_normals(self):
+ """True to use gradient descent as normals."""
+ return self.c_mc.invert_normals
+
+ def get_vertices(self):
+ """Vertices currently computed (ndarray of dim NbVertices x 3)
+
+ Order is dim0, dim1, dim2 (i.e., z, y, x if dim0 is depth).
+ """
+ return numpy.array(self.c_mc.vertices).reshape(-1, 3)
+
+ def get_normals(self):
+ """Normals currently computed (ndarray of dim NbVertices x 3)
+
+ Order is dim0, dim1, dim2 (i.e., z, y, x if dim0 is depth).
+ """
+ return numpy.array(self.c_mc.normals).reshape(-1, 3)
+
+ def get_indices(self):
+ """Triangle indices currently computed (ndarray of dim NbTriangles x 3)
+ """
+ return numpy.array(self.c_mc.indices,
+ dtype=numpy.uint32).reshape(-1, 3)