wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-11 11:25:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

mc deal with degenerate triangles

Browse Source
This commit is contained in:
almar
2016-07-25 12:26:28 +02:00
committed by Almar Klein
parent e2ebc0d68f
commit 8d9ad26488
2 changed files with 68 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/measure/_marching_cubes_lewiner.py
+11 -3
View File
@@ -14,7 +14,7 @@ from . import _marching_cubes_lewiner_cy
def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
step_size=1, use_classic=False):
step_size=1, allow_degenerate=True, use_classic=False):
"""
Lewiner marching cubes algorithm to find surfaces in 3d volumetric data
@@ -35,6 +35,10 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
Step size in voxels. Default 1. Larger steps yield faster but
coarser results. The result will always be topologically correct
though.
allow_degenerate : bool
Whether to allow degenerate triangles in the end-result. Default True.
If False, degenerate triangles are removed, making the algorithm
about twice as slow.
use_classic : bool
If given and True, the classic marching cubes by Lorensen (1987)
is used. This option is included for reference purposes. Note
@@ -80,7 +84,7 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
raise ValueError('Input volume should be a 3D numpy array.')
if volume.shape[0] < 2 or volume.shape[1] < 2 or volume.shape[2] < 2:
raise ValueError("Input array must be at least 2x2x2.")
volume = np.array(volume, dtype=np.float32, order="C", copy=False)
volume = np.ascontiguousarray(volume, np.float32) # no copy if not necessary
# Check/convert other inputs:
# level
@@ -119,7 +123,11 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
if spacing != (1, 1, 1):
vertices = vertices * np.r_[spacing]
return vertices, faces, normals, values
if allow_degenerate:
return vertices, faces, normals, values
else:
fun = _marching_cubes_lewiner_cy.remove_degenerate_faces
return fun(vertices, faces, normals, values)
def _toArray(args):
skimage/measure/_marching_cubes_lewiner_cy.pyx
+57 -9
View File
@@ -1,6 +1,11 @@
# -*- coding: utf-8 -*-
# Copyright (C) 2012, Almar Klein
# Copyright (C) 2002, Thomas Lewiner
#
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
"""
This is an implementation of the marching cubes algorithm proposed in:
@@ -43,11 +48,55 @@ cdef double FLT_EPSILON = np.spacing(1.0) #0.0000001
# Define abs function for doubles
cdef inline double dabs(double a): return a if a>=0 else -a
cdef inline int imin(int a, int b): return a if a<b else b
# todo: allow dynamic isovalue?
# todo: can we disable Cython from checking for zero division? Sometimes we know that it never happens!
def remove_degenerate_faces(vertices, faces, *arrays):
vertices_map0 = np.arange(len(vertices), dtype=np.int32)
vertices_map1 = vertices_map0.copy()
faces_ok = np.ones(len(faces), dtype=np.int32)
cdef np.ndarray[FLOAT32_T, ndim=2] vertices_ = vertices
cdef np.ndarray[FLOAT32_T, ndim=1] v1, v2, v3
cdef np.ndarray[INT32_T, ndim=2] faces_ = faces
cdef np.ndarray[INT32_T, ndim=1] face_
cdef np.ndarray[INT32_T, ndim=1] vertices_map1_ = vertices_map1
cdef np.ndarray[INT32_T, ndim=1] faces_ok_ = faces_ok
cdef int j, i1, i2, i3
# Iterate over all faces. When we encounter a degenerate triangle,
# we update the vertex map, i.e. we merge the corresponding vertices.
for j in range(faces_.shape[0]):
face_ = faces_[j]
i1, i2, i3 = face_[0], face_[1], face_[2]
v1, v2, v3 = vertices_[face_[0]], vertices_[face_[1]], vertices_[face_[2]]
if v1[0] == v2[0] and v1[1] == v2[1] and v1[2] == v2[2]:
vertices_map1_[i1] = vertices_map1_[i1] = imin(vertices_map1_[i1], vertices_map1_[i2])
faces_ok_[j] = 0
if v1[0] == v3[0] and v1[1] == v3[1] and v1[2] == v3[2]:
vertices_map1_[i1] = vertices_map1_[i3] = imin(vertices_map1_[i1], vertices_map1_[i3])
faces_ok_[j] = 0
if v2[0] == v3[0] and v2[1] == v3[1] and v2[2] == v3[2]:
vertices_map1_[i2] = vertices_map1_[i3] = imin(vertices_map1_[i2], vertices_map1_[i3])
faces_ok_[j] = 0
# Create mask and mapping to new vertex indices
vertices_ok = vertices_map1 == vertices_map0
vertices_map2 = np.cumsum(vertices_ok) - 1
# Apply selection and mapping
faces2 = vertices_map2[vertices_map1[faces[faces_ok>0]]]
vertices2 = vertices[vertices_ok]
arrays2 = [arr[vertices_ok] for arr in arrays]
return (vertices2, faces2) + tuple(arrays2)
cdef class Cell:
""" Class to keep track of some stuff during the whole cube marching
@@ -449,7 +498,6 @@ cdef class Cell:
vi = lut.get3(lutIndex, lutIndex2, i*3+j)
self._add_face_from_edge_index(vi)
## Used internally
cdef void _add_face_from_edge_index(self, int vi):
@@ -505,6 +553,8 @@ cdef class Cell:
tmpf1 = 1.0 / (FLT_EPSILON + dabs(self.vv[index1]))
tmpf2 = 1.0 / (FLT_EPSILON + dabs(self.vv[index2]))
# print('indexInVertexArray', self.x, self.y, self.z, '-', vi, indexInVertexArray, indexInFaceLayer)
if indexInVertexArray >= 0:
# Vertex already calculated, only need to add face and gradient
self.add_face(indexInVertexArray)
@@ -516,6 +566,7 @@ cdef class Cell:
fx, fy, fz, ff = 0.0, 0.0, 0.0, 0.0
fx += <double>dx1 * tmpf1; fy += <double>dy1 * tmpf1; fz += <double>dz1 * tmpf1; ff += tmpf1
fx += <double>dx2 * tmpf2; fy += <double>dy2 * tmpf2; fz += <double>dz2 * tmpf2; ff += tmpf2
# Add vertex
indexInVertexArray = self.add_vertex(
<double>self.x + stp*fx/ff,
@@ -550,7 +601,7 @@ cdef class Cell:
This method returns -1 if no vertex has been defined yet.
vertices edes edge-indices per cell
vertices edges edge-indices per cell
* 7 ________ 6 _____6__ ________
* /| /| 7/| /| /| /|
* / | / | / | /5 | / | / |
@@ -560,7 +611,7 @@ cdef class Cell:
* | / | / 8 3/ 9 / 2 / | /
* | / | / | / | /1 | 1 | /
* |/_______|/ |/___0___|/ |/___0___|/
* 0 1 0 1
* 0 1
*/
"""
@@ -590,8 +641,8 @@ cdef class Cell:
elif vi == 3: # no step
j = 1
elif vi<12:
# four vertical edges
elif vi < 12:
# 4 vertical edges
faceLayer = self.faceLayer1
j = 2
@@ -889,9 +940,6 @@ cdef class LutProvider:
self.SUBCONFIG13 = Lut(SUBCONFIG13)
@cython.boundscheck(False)
@cython.wraparound(False)
def marching_cubes(im, double isovalue, LutProvider luts, int st=1, int classic=0):
""" marching_cubes(im, double isovalue, LutProvider luts, int st=1, int classic=0)
This is the main entry to apply marching cubes.