wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-05 13:21:12 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
patch-1
scikit-image/skimage/measure/_marching_cubes_lewiner_cy.pyx
T
Almar Klein f96e458624 MC: typed input arg
2016-07-25 12:26:28 +02:00

1435 lines
54 KiB
Cython
Raw Permalink Blame History

# -*- coding: utf-8 -*-
#
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
"""
This is an implementation of the marching cubes algorithm proposed in:
Efficient implementation of Marching Cubes' cases with topological guarantees.
Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan Tavares.
Journal of Graphics Tools 8(2): pp. 1-15 (december 2003)
This algorithm has the advantage that it provides topologically correct
results, and the algorithms implementation is relatively simple. Most
of the magic is in the lookup tables, which are provided as open source.
Originally implemented in C++ by Thomas Lewiner in 2002, ported to Cython
by Almar Klein in 2012. Adapted for scikit-image in 2016.
"""
# Cython specific imports
import numpy as np
cimport numpy as np
import cython
# Enable low level memory management
# from libc.stdlib cimport malloc, free
cdef extern from "stdlib.h": # The cimport does not work on my Linux Laptop
void free(void* ptr)
void* malloc(size_t size)
# Define tiny winy number
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 float [:, :] vertices_ = vertices
cdef float [:] v1, v2, v3
cdef int [:, :] faces_ = faces
cdef int [:] vertices_map1_ = vertices_map1
cdef int [:] 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]):
i1, i2, i3 = faces_[j][0], faces_[j][1], faces_[j][2]
v1, v2, v3 = vertices_[i1], vertices_[i2], vertices_[i3]
if v1[0] == v2[0] and v1[1] == v2[1] and v1[2] == v2[2]:
vertices_map1_[i1] = vertices_map1_[i2] = 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
procedure.
This "struct" keeps track of the current cell location, and the values
of corners of the cube. Gradients for the cube corners are calculated
when needed.
Additionally, it keeps track of the array of vertices, faces and normals.
Notes on vertices
-----------------
The vertices are stored in a C-array that is increased in size with
factors of two if needed. The same applies to the faces and normals.
Notes on faces
--------------
To keep track of the vertices already defined, this class maintains
two faceLayer arrays. faceLayer1 is of the current layer (z-value)
and faceLayer2 is of the next. Both face layers have 4 elements per
cell in that layer, 1 for each unique edge per cell (see
get_index_in_facelayer). These are initialized as -1, and set to the
index in the vertex array when a new vertex is created.
In summary, this allows us to keep track of the already created
vertices without keeping a very big array.
Notes on normals
----------------
The normal is simply defined as the gradient. Each time that a face is
created, we also add the gradient of that vertex position to the
normals array. The gradients are all calculated from the differences between
the 8 corners of the current cube, but because the final value of a normal
was contributed from multiple cells, the normals are quite accurate.
"""
# Reference to LUTS object
cdef LutProvider luts
# Location of cube
cdef int x
cdef int y
cdef int z
# Stepsize
cdef int step
# Values of cube corners (isovalue subtracted)
cdef double v0
cdef double v1
cdef double v2
cdef double v3
cdef double v4
cdef double v5
cdef double v6
cdef double v7
# Small arrays to store the above values in (allowing indexing)
# and also the gradient at these points
cdef double *vv
cdef double *vg
# Max value of the eight corners
cdef double vmax
# Vertex position of center of cube (only calculated if needed)
cdef double v12_x
cdef double v12_y
cdef double v12_z
# And corresponding gradient
cdef double v12_xg
cdef double v12_yg
cdef double v12_zg
cdef int v12_calculated # a boolean
# The index value, our magic 256 bit word
cdef int index
# Dimensions of the total volume
cdef int nx
cdef int ny
cdef int nz
# Arrays with face information
cdef int *faceLayer # The current facelayer (reference-copy of one of the below)
cdef int *faceLayer1 # The actual first face layer
cdef int *faceLayer2 # The actual second face layer
# Stuff to store the output vertices
cdef float *_vertices
cdef float *_normals
cdef float *_values
cdef int _vertexCount
cdef int _vertexMaxCount
# Stuff to store the output faces
cdef int *_faces
cdef int _faceCount
cdef int _faceMaxCount
def __init__(self, LutProvider luts, int nx, int ny, int nz):
self.luts = luts
self.nx, self.ny, self.nz = nx, ny, nz
# Allocate face layers
self.faceLayer1 = <int *>malloc(self.nx*self.ny*4 * sizeof(int))
self.faceLayer2 = <int *>malloc(self.nx*self.ny*4 * sizeof(int))
cdef int i
for i in range(self.nx*self.ny*4):
self.faceLayer1[i] = -1
self.faceLayer2[i] = -1
self.faceLayer = self.faceLayer1
def __cinit__(self):
# Init tiny arrays for vertices and gradients at the vertices
self.vv = <double *>malloc(8 * sizeof(double))
self.vg = <double *>malloc(8*3 * sizeof(double))
# Init face layers
self.faceLayer1 = NULL
self.faceLayer2 = NULL
# Init vertices
self._vertexCount = 0
self._vertexMaxCount = 8
self._vertices = <float *>malloc(self._vertexMaxCount*3 * sizeof(float))
self._normals = <float *>malloc(self._vertexMaxCount*3 * sizeof(float))
self._values = <float *>malloc(self._vertexMaxCount * sizeof(float))
# Clear normals and values
cdef int i, j
for i in range(self._vertexMaxCount):
self._values[i] = 0.0
for j in range(3):
self._normals[i*3+j] = 0.0
# Init faces
self._faceCount = 0
self._faceMaxCount = 8
self._faces = <int *>malloc(self._faceMaxCount * sizeof(int))
def __dealloc__(self):
if self.vv is not NULL:
free(self.vv)
if self.vg is not NULL:
free(self.vg)
if self.faceLayer1 is not NULL:
free(self.faceLayer1)
if self.faceLayer2 is not NULL:
free(self.faceLayer2)
if self._vertices is not NULL:
free(self._vertices)
if self._normals is not NULL:
free(self._normals)
if self._values is not NULL:
free(self._values)
if self._faces is not NULL:
free(self._faces)
cdef void _increase_size_vertices(self):
""" Increase the size of the vertices array by a factor two.
"""
# Allocate new array
cdef int newMaxCount = self._vertexMaxCount * 2
cdef float *newVertices = <float *>malloc(newMaxCount*3 * sizeof(float))
cdef float *newNormals = <float *>malloc(newMaxCount*3 * sizeof(float))
cdef float *newValues = <float *>malloc(newMaxCount * sizeof(float))
# Clear
cdef int i, j
for i in range(self._vertexCount, newMaxCount):
newValues[i] = 0.0
for j in range(3):
newNormals[i*3+j] = 0.0
# Copy
for i in range(self._vertexCount):
newValues[i] = self._values[i]
for j in range(3):
newVertices[i*3+j] = self._vertices[i*3+j]
newNormals[i*3+j] = self._normals[i*3+j]
# Apply
free(self._vertices); self._vertices = newVertices
free(self._normals); self._normals = newNormals
free(self._values); self._values = newValues
self._vertexMaxCount = newMaxCount
cdef void _increase_size_faces(self):
""" Increase the size of the faces array by a factor two.
"""
# Allocate new array
cdef int newMaxCount = self._faceMaxCount * 2
cdef int *newFaces = <int *>malloc(newMaxCount * sizeof(int))
# Copy
cdef int i
for i in range(self._faceCount):
newFaces[i] = self._faces[i]
# Apply
free(self._faces)
self._faces = newFaces
self._faceMaxCount = newMaxCount
## Adding results
cdef int add_vertex(self, float x, float y, float z):
""" Add a vertex to the result. Return index in vertex array.
"""
# Check if array is large enough
if self._vertexCount >= self._vertexMaxCount:
self._increase_size_vertices()
# Add vertex
self._vertices[self._vertexCount*3+0] = x
self._vertices[self._vertexCount*3+1] = y
self._vertices[self._vertexCount*3+2] = z
self._vertexCount += 1
return self._vertexCount -1
cdef void add_gradient(self, int vertexIndex, float gx, float gy, float gz):
""" Add a gradient value to the vertex corresponding to the given index.
"""
self._normals[vertexIndex*3+0] += gx
self._normals[vertexIndex*3+1] += gy
self._normals[vertexIndex*3+2] += gz
cdef void add_gradient_from_index(self, int vertexIndex, int i, float strength):
""" Add a gradient value to the vertex corresponding to the given index.
vertexIndex is the index in the large array of vertices that is returned.
i is the index of the array of vertices 0-7 for the current cell.
"""
self.add_gradient(vertexIndex, self.vg[i*3+0] * strength, self.vg[i*3+1] * strength, self.vg[i*3+2] * strength)
cdef add_face(self, int index):
""" Add a face to the result. Also updates the value.
"""
# Check if array is large enough
if self._faceCount >= self._faceMaxCount:
self._increase_size_faces()
# Add face
self._faces[self._faceCount] = index
self._faceCount += 1
# Also update value
if self.vmax > self._values[index]:
self._values[index] = self.vmax
## Getting results
def get_vertices(self):
""" Get the final vertex array.
"""
vertices = np.empty((self._vertexCount,3), np.float32)
cdef float [:, :] vertices_ = vertices
cdef int i, j
for i in range(self._vertexCount):
for j in range(3):
vertices_[i, j] = self._vertices[i*3+j]
return vertices
def get_normals(self):
""" Get the final normals array.
The normals are normalized to unit length.
"""
normals = np.empty((self._vertexCount,3), np.float32)
cdef float [:, :] normals_ = normals
cdef int i, j
cdef double length, dtmp
for i in range(self._vertexCount):
length = 0.0
for j in range(3):
dtmp = self._normals[i*3+j] # Make it double before taking **2!
length += dtmp*dtmp
if length > 0.0:
length = 1.0 / length**0.5
for j in range(3):
normals_[i,j] = self._normals[i*3+j] * length
return normals
def get_faces(self):
faces = np.empty((self._faceCount,), np.int32)
cdef int [:] faces_ = faces
cdef int i, j
for i in range(self._faceCount):
faces_[i] = self._faces[i]
return faces
def get_values(self):
values = np.empty((self._vertexCount,), np.float32)
cdef float [:] values_ = values
cdef int i, j
for i in range(self._vertexCount):
values_[i] = self._values[i]
return values
## Called from marching cube function
cdef void new_z_value(self):
""" This method should be called each time a new z layer is entered.
We will swap the layers with face information and empty the second.
"""
# Swap layers
self.faceLayer1, self.faceLayer2 = self.faceLayer2, self.faceLayer1
# Empty last
cdef int i
for i in range(self.nx*self.ny*4):
self.faceLayer2[i] = -1
cdef void set_cube(self, double isovalue, int x, int y, int z, int step,
double v0, double v1, double v2, double v3,
double v4, double v5, double v6, double v7):
""" Set the cube to the new location.
Set the values of the cube corners. The isovalue is subtracted
from them, such that in further calculations the isovalue can be
taken as zero.
This method also calculated the magic 256 word to identify the
cases (i.e. cell.index).
"""
# Set location and step
self.x = x
self.y = y
self.z = z
self.step = step
# Set values
self.v0 = v0 - isovalue
self.v1 = v1 - isovalue
self.v2 = v2 - isovalue
self.v3 = v3 - isovalue
self.v4 = v4 - isovalue
self.v5 = v5 - isovalue
self.v6 = v6 - isovalue
self.v7 = v7 - isovalue
# Calculate index
cdef int index = 0
if self.v0 > 0.0: index += 1
if self.v1 > 0.0: index += 2
if self.v2 > 0.0: index += 4
if self.v3 > 0.0: index += 8
if self.v4 > 0.0: index += 16
if self.v5 > 0.0: index += 32
if self.v6 > 0.0: index += 64
if self.v7 > 0.0: index += 128
self.index = index
# Reset c12
self.v12_calculated = 0
cdef void add_triangles(self, Lut lut, int lutIndex, int nt):
""" Add triangles.
The vertices for the triangles are specified in the given
Lut at the specified index. There are nt triangles.
The reason that nt should be given is because it is often known
beforehand.
"""
cdef int i, j
cdef int vi
self.prepare_for_adding_triangles()
for i in range(nt):
for j in range(3):
# Get two sides for each element in this vertex
vi = lut.get2(lutIndex, i*3+j)
self._add_face_from_edge_index(vi)
cdef void add_triangles2(self, Lut lut, int lutIndex, int lutIndex2, int nt):
""" Same as add_triangles, except that now the geometry is in a LUT
with 3 dimensions, and an extra index is provided.
"""
cdef int i, j
cdef int vi
self.prepare_for_adding_triangles()
for i in range(nt):
for j in range(3):
# Get two sides for each element in this vertex
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):
""" Add one face from an edge index. Only adds a face if the
vertex already exists. Otherwise also adds a vertex and applies
interpolation.
"""
# typedefs
cdef int indexInVertexArray, indexInFaceLayer
cdef int dx1, dy1, dz1
cdef int dx2, dy2, dz2
cdef int index1, index2
cdef double tmpf1, tmpf2
cdef double fx, fy, fz, ff
cdef double stp = <double>self.step
# Get index in the face layer and corresponding vertex number
indexInFaceLayer = self.get_index_in_facelayer(vi)
indexInVertexArray = self.faceLayer[indexInFaceLayer]
# If we have the center vertex, we have things pre-calculated,
# otherwise we need to interpolate.
# In both cases we distinguish between having this vertex already
# or not.
if vi == 12: # center vertex
if self.v12_calculated == 0:
self.calculate_center_vertex()
if indexInVertexArray >= 0:
# Vertex already calculated, only need to add face and gradient
self.add_face(indexInVertexArray)
self.add_gradient(indexInVertexArray, self.v12_xg, self.v12_yg, self.v12_zg)
else:
# Add precalculated center vertex position (is interpolated)
indexInVertexArray = self.add_vertex( self.v12_x, self.v12_y, self.v12_z)
# Update face layer
self.faceLayer[indexInFaceLayer] = indexInVertexArray
# Add face and gradient
self.add_face(indexInVertexArray)
self.add_gradient(indexInVertexArray, self.v12_xg, self.v12_yg, self.v12_zg)
else:
# Get relative edge indices for x, y and z
dx1, dx2 = self.luts.EDGESRELX.get2(vi,0), self.luts.EDGESRELX.get2(vi,1)
dy1, dy2 = self.luts.EDGESRELY.get2(vi,0), self.luts.EDGESRELY.get2(vi,1)
dz1, dz2 = self.luts.EDGESRELZ.get2(vi,0), self.luts.EDGESRELZ.get2(vi,1)
# Make two vertex indices
index1 = dz1*4 + dy1*2 + dx1
index2 = dz2*4 + dy2*2 + dx2
# Define strength of both corners
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)
self.add_gradient_from_index(indexInVertexArray, index1, tmpf1)
self.add_gradient_from_index(indexInVertexArray, index2, tmpf2)
else:
# Interpolate by applying a kind of center-of-mass method
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,
<double>self.y + stp*fy/ff,
<double>self.z + stp*fz/ff )
# Update face layer
self.faceLayer[indexInFaceLayer] = indexInVertexArray
# Add face and gradient
self.add_face(indexInVertexArray)
self.add_gradient_from_index(indexInVertexArray, index1, tmpf1)
self.add_gradient_from_index(indexInVertexArray, index2, tmpf2)
# # Create vertex non-interpolated
# self.add_vertex( self.x + 0.5* dx1 + 0.5 * dx2,
# self.y + 0.5* dy1 + 0.5 * dy2,
# self.z + 0.5* dz1 + 0.5 * dz2 )
cdef int get_index_in_facelayer(self, int vi):
"""
Get the index of a vertex position, given the edge on which it lies.
We keep a list of faces so we can reuse vertices. This improves
speed because we need less interpolation, and the result is more
compact and can be visualized better because normals can be
interpolated.
For each cell, we store 4 vertex indices; all other edges can be
represented as the edge of another cell. The fourth is the center vertex.
This method returns -1 if no vertex has been defined yet.
vertices edges edge-indices per cell
* 7 ________ 6 _____6__ ________
* /| /| 7/| /| /| /|
* / | / | / | /5 | / | / |
* 4 /_______ / | /__4____ / 10 /_______ / |
* | | |5 | | 11 | | | | | |
* | 3|__|_____|2 | |__|__2__| | |__|_____|
* | / | / 8 3/ 9 / 2 / | /
* | / | / | / | /1 | 1 | /
* |/_______|/ |/___0___|/ |/___0___|/
* 0 1
*/
"""
# Init indices, both are corrected below
cdef int i = self.nx * self.y + self.x # Index of cube to get vertex at
cdef int j = 0 # Vertex number for that cell
cdef int vi_ = vi
cdef int *faceLayer
# Select either upper or lower half
if vi < 8:
# 8 horizontal edges
if vi < 4:
faceLayer = self.faceLayer1
else:
vi -= 4
faceLayer = self.faceLayer2
# Calculate actual index based on edge
#if vi == 0: pass # no step
if vi == 1: # step in x
i += self.step
j = 1
elif vi == 2: # step in y
i += self.nx * self.step
elif vi == 3: # no step
j = 1
elif vi < 12:
# 4 vertical edges
faceLayer = self.faceLayer1
j = 2
#if vi == 8: pass # no step
if vi == 9: # step in x
i += self.step
elif vi == 10: # step in x and y
i += self.nx * self.step + self.step
elif vi == 11: # step in y
i += self.nx * self.step
else:
# center vertex
faceLayer = self.faceLayer1
j = 3
# Store facelayer and return index
self.faceLayer = faceLayer # Dirty way of returning a value
return 4*i + j
cdef void prepare_for_adding_triangles(self):
""" Calculates some things to help adding the triangles:
array with corner values, max corner value, gradient at each corner.
"""
cdef int i
# Copy values in array so we can index them. Note the misalignment
# because the numbering does not correspond with bitwise OR of xyz.
self.vv[0] = self.v0
self.vv[1] = self.v1
self.vv[2] = self.v3#
self.vv[3] = self.v2#
self.vv[4] = self.v4
self.vv[5] = self.v5
self.vv[6] = self.v7#
self.vv[7] = self.v6#
# Calculate max
cdef double vmin, vmax
vmin, vmax = 0.0, 0.0
for i in range(8):
if self.vv[i] > vmax:
vmax = self.vv[i]
if self.vv[i] < vmin:
vmin = self.vv[i]
self.vmax = vmax-vmin
# Calculate gradients
# Derivatives, selected to always point in same direction.
# Note that many corners have the same components as other points,
# by interpolating and averaging the normals this is solved.
# todo: we can potentially reuse these similar to how we store vertex indices in face layers
self.vg[0*3+0], self.vg[0*3+1], self.vg[0*3+2] = self.v0-self.v1, self.v0-self.v3, self.v0-self.v4
self.vg[1*3+0], self.vg[1*3+1], self.vg[1*3+2] = self.v0-self.v1, self.v1-self.v2, self.v1-self.v5
self.vg[2*3+0], self.vg[2*3+1], self.vg[2*3+2] = self.v3-self.v2, self.v1-self.v2, self.v2-self.v6
self.vg[3*3+0], self.vg[3*3+1], self.vg[3*3+2] = self.v3-self.v2, self.v0-self.v3, self.v3-self.v7
self.vg[4*3+0], self.vg[4*3+1], self.vg[4*3+2] = self.v4-self.v5, self.v4-self.v7, self.v0-self.v4
self.vg[5*3+0], self.vg[5*3+1], self.vg[5*3+2] = self.v4-self.v5, self.v5-self.v6, self.v1-self.v5
self.vg[6*3+0], self.vg[6*3+1], self.vg[6*3+2] = self.v7-self.v6, self.v5-self.v6, self.v2-self.v6
self.vg[7*3+0], self.vg[7*3+1], self.vg[7*3+2] = self.v7-self.v6, self.v4-self.v7, self.v3-self.v7
cdef void calculate_center_vertex(self):
""" Calculate interpolated center vertex and its gradient.
"""
cdef double v0, v1, v2, v3, v4, v5, v6, v7
cdef double fx, fy, fz, ff
fx, fy, fz, ff = 0.0, 0.0, 0.0, 0.0
# Define "strength" of each corner of the cube that we need
v0 = 1.0 / (FLT_EPSILON + dabs(self.v0))
v1 = 1.0 / (FLT_EPSILON + dabs(self.v1))
v2 = 1.0 / (FLT_EPSILON + dabs(self.v2))
v3 = 1.0 / (FLT_EPSILON + dabs(self.v3))
v4 = 1.0 / (FLT_EPSILON + dabs(self.v4))
v5 = 1.0 / (FLT_EPSILON + dabs(self.v5))
v6 = 1.0 / (FLT_EPSILON + dabs(self.v6))
v7 = 1.0 / (FLT_EPSILON + dabs(self.v7))
# Apply a kind of center-of-mass method
fx += 0.0*v0; fy += 0.0*v0; fz += 0.0*v0; ff += v0
fx += 1.0*v1; fy += 0.0*v1; fz += 0.0*v1; ff += v1
fx += 1.0*v2; fy += 1.0*v2; fz += 0.0*v2; ff += v2
fx += 0.0*v3; fy += 1.0*v3; fz += 0.0*v3; ff += v3
fx += 0.0*v4; fy += 0.0*v4; fz += 1.0*v4; ff += v4
fx += 1.0*v5; fy += 0.0*v5; fz += 1.0*v5; ff += v5
fx += 1.0*v6; fy += 1.0*v6; fz += 1.0*v6; ff += v6
fx += 0.0*v7; fy += 1.0*v7; fz += 1.0*v7; ff += v7
# Store
cdef double stp = <double>self.step
self.v12_x = self.x + stp * fx / ff
self.v12_y = self.y + stp * fy / ff
self.v12_z = self.z + stp * fz / ff
# Also pre-calculate gradient of center
# note that prepare_for_adding_triangles() must have been called for
# the gradient data to exist.
self.v12_xg = ( v0*self.vg[0*3+0] + v1*self.vg[1*3+0] + v2*self.vg[2*3+0] + v3*self.vg[3*3+0] +
v4*self.vg[4*3+0] + v5*self.vg[5*3+0] + v6*self.vg[6*3+0] + v7*self.vg[7*3+0] )
self.v12_yg = ( v0*self.vg[0*3+1] + v1*self.vg[1*3+1] + v2*self.vg[2*3+1] + v3*self.vg[3*3+1] +
v4*self.vg[4*3+1] + v5*self.vg[5*3+1] + v6*self.vg[6*3+1] + v7*self.vg[7*3+1] )
self.v12_xg = ( v0*self.vg[0*3+2] + v1*self.vg[1*3+2] + v2*self.vg[2*3+2] + v3*self.vg[3*3+2] +
v4*self.vg[4*3+2] + v5*self.vg[5*3+2] + v6*self.vg[6*3+2] + v7*self.vg[7*3+2] )
# Set flag that this stuff is calculated
self.v12_calculated = 1
cdef class Lut:
""" Representation of a lookup table.
The tables are initially defined as numpy arrays. On initialization,
this class converts them to a C array for fast access.
This class defines functions to look up values using 1, 2 or 3 indices.
"""
cdef char* VALUES
cdef int L0 # Length
cdef int L1 # size of tuple
cdef int L2 # size of tuple in tuple (if any)
def __init__(self, array):
# Get the shape of the LUT
self.L1 = 1
self.L2 = 1
#
self.L0 = array.shape[0]
if array.ndim > 1:
self.L1 = array.shape[1]
if array.ndim > 2:
self.L2 = array.shape[2]
# Copy the contents
array = array.ravel()
cdef int n, N
N = self.L0 * self.L1 * self.L2
self.VALUES = <char *> malloc(N * sizeof(char))
for n in range(N):
self.VALUES[n] = array[n]
def __cinit__(self):
self.VALUES = NULL
def __dealloc__(self):
if self.VALUES is not NULL:
free(self.VALUES)
cdef int get1(self, int i0):
return self.VALUES[i0]
cdef int get2(self, int i0, int i1):
return self.VALUES[i0*self.L1 + i1]
cdef int get3(self, int i0, int i1, int i2):
return self.VALUES[i0*self.L1*self.L2 + i1*self.L2 + i2]
cdef class LutProvider:
""" Class that provides a common interface to the many lookup tables
used by the algorithm.
All the lists of lut names are autogenerated to prevent human error.
"""
cdef Lut EDGESRELX # Edges relative X
cdef Lut EDGESRELY
cdef Lut EDGESRELZ
cdef Lut CASESCLASSIC
cdef Lut CASES
cdef Lut TILING1
cdef Lut TILING2
cdef Lut TILING3_1
cdef Lut TILING3_2
cdef Lut TILING4_1
cdef Lut TILING4_2
cdef Lut TILING5
cdef Lut TILING6_1_1
cdef Lut TILING6_1_2
cdef Lut TILING6_2
cdef Lut TILING7_1
cdef Lut TILING7_2
cdef Lut TILING7_3
cdef Lut TILING7_4_1
cdef Lut TILING7_4_2
cdef Lut TILING8
cdef Lut TILING9
cdef Lut TILING10_1_1
cdef Lut TILING10_1_1_
cdef Lut TILING10_1_2
cdef Lut TILING10_2
cdef Lut TILING10_2_
cdef Lut TILING11
cdef Lut TILING12_1_1
cdef Lut TILING12_1_1_
cdef Lut TILING12_1_2
cdef Lut TILING12_2
cdef Lut TILING12_2_
cdef Lut TILING13_1
cdef Lut TILING13_1_
cdef Lut TILING13_2
cdef Lut TILING13_2_
cdef Lut TILING13_3
cdef Lut TILING13_3_
cdef Lut TILING13_4
cdef Lut TILING13_5_1
cdef Lut TILING13_5_2
cdef Lut TILING14
cdef Lut TEST3
cdef Lut TEST4
cdef Lut TEST6
cdef Lut TEST7
cdef Lut TEST10
cdef Lut TEST12
cdef Lut TEST13
cdef Lut SUBCONFIG13
def __init__(self, EDGESRELX, EDGESRELY, EDGESRELZ, CASESCLASSIC, CASES,
TILING1, TILING2, TILING3_1, TILING3_2, TILING4_1, TILING4_2,
TILING5, TILING6_1_1, TILING6_1_2, TILING6_2, TILING7_1, TILING7_2,
TILING7_3, TILING7_4_1, TILING7_4_2, TILING8, TILING9,
TILING10_1_1, TILING10_1_1_, TILING10_1_2, TILING10_2, TILING10_2_,
TILING11, TILING12_1_1, TILING12_1_1_, TILING12_1_2, TILING12_2,
TILING12_2_, TILING13_1, TILING13_1_, TILING13_2, TILING13_2_,
TILING13_3, TILING13_3_, TILING13_4, TILING13_5_1, TILING13_5_2,
TILING14,
TEST3, TEST4, TEST6, TEST7, TEST10, TEST12, TEST13,
SUBCONFIG13,
):
self.EDGESRELX = Lut(EDGESRELX)
self.EDGESRELY = Lut(EDGESRELY)
self.EDGESRELZ = Lut(EDGESRELZ)
self.CASESCLASSIC = Lut(CASESCLASSIC)
self.CASES = Lut(CASES)
self.TILING1 = Lut(TILING1)
self.TILING2 = Lut(TILING2)
self.TILING3_1 = Lut(TILING3_1)
self.TILING3_2 = Lut(TILING3_2)
self.TILING4_1 = Lut(TILING4_1)
self.TILING4_2 = Lut(TILING4_2)
self.TILING5 = Lut(TILING5)
self.TILING6_1_1 = Lut(TILING6_1_1)
self.TILING6_1_2 = Lut(TILING6_1_2)
self.TILING6_2 = Lut(TILING6_2)
self.TILING7_1 = Lut(TILING7_1)
self.TILING7_2 = Lut(TILING7_2)
self.TILING7_3 = Lut(TILING7_3)
self.TILING7_4_1 = Lut(TILING7_4_1)
self.TILING7_4_2 = Lut(TILING7_4_2)
self.TILING8 = Lut(TILING8)
self.TILING9 = Lut(TILING9)
self.TILING10_1_1 = Lut(TILING10_1_1)
self.TILING10_1_1_ = Lut(TILING10_1_1_)
self.TILING10_1_2 = Lut(TILING10_1_2)
self.TILING10_2 = Lut(TILING10_2)
self.TILING10_2_ = Lut(TILING10_2_)
self.TILING11 = Lut(TILING11)
self.TILING12_1_1 = Lut(TILING12_1_1)
self.TILING12_1_1_ = Lut(TILING12_1_1_)
self.TILING12_1_2 = Lut(TILING12_1_2)
self.TILING12_2 = Lut(TILING12_2)
self.TILING12_2_ = Lut(TILING12_2_)
self.TILING13_1 = Lut(TILING13_1)
self.TILING13_1_ = Lut(TILING13_1_)
self.TILING13_2 = Lut(TILING13_2)
self.TILING13_2_ = Lut(TILING13_2_)
self.TILING13_3 = Lut(TILING13_3)
self.TILING13_3_ = Lut(TILING13_3_)
self.TILING13_4 = Lut(TILING13_4)
self.TILING13_5_1 = Lut(TILING13_5_1)
self.TILING13_5_2 = Lut(TILING13_5_2)
self.TILING14 = Lut(TILING14)
self.TEST3 = Lut(TEST3)
self.TEST4 = Lut(TEST4)
self.TEST6 = Lut(TEST6)
self.TEST7 = Lut(TEST7)
self.TEST10 = Lut(TEST10)
self.TEST12 = Lut(TEST12)
self.TEST13 = Lut(TEST13)
self.SUBCONFIG13 = Lut(SUBCONFIG13)
def marching_cubes(float [:, :, :] im not None, 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.
Returns (vertices, faces, normals, values)
"""
# Get dimemsnions
cdef int Nx, Ny, Nz
Nx, Ny, Nz = im.shape[2], im.shape[1], im.shape[0]
# Create cell to use throughout
cdef Cell cell = Cell(luts, Nx, Ny, Nz)
# Typedef variables
cdef int x, y, z, x_st, y_st, z_st
cdef int nt
cdef int case, config, subconfig
# Unfortunately specifying a step in range() siginificantly degrades
# performance. Therefore we use a while loop.
# we have: max_x = Nx_bound + st + st - 1
# -> Nx_bound = max_allowable_x + 1 - 2 * st
# -> Nx_bound = Nx - 2 * st
assert st > 0
cdef int Nx_bound, Ny_bound, Nz_bound
Nx_bound, Ny_bound, Nz_bound = Nx - 2 * st, Ny - 2 * st, Nz - 2 * st # precalculated index range
z = -st
while z < Nz_bound:
z += st
z_st = z + st
cell.new_z_value() # Indicate that we enter a new layer
y = -st
while y < Ny_bound:
y += st
y_st = y + st
x = -st
while x < Nx_bound:
x += st
x_st = x + st
# Initialize cell
cell.set_cube(isovalue, x, y, z, st,
im[z ,y, x], im[z ,y, x_st], im[z ,y_st, x_st], im[z ,y_st, x],
im[z_st,y, x], im[z_st,y, x_st], im[z_st,y_st, x_st], im[z_st,y_st, x] )
# Do classic!
if classic:
# Determine number of vertices
nt = 0
while luts.CASESCLASSIC.get2(cell.index, 3*nt) != -1:
nt += 1
# Add triangles
if nt > 0:
cell.add_triangles(luts.CASESCLASSIC, cell.index, nt)
else:
# Get case, if non-nul, enter the big switch
case = luts.CASES.get2(cell.index, 0)
if case > 0:
config = luts.CASES.get2(cell.index, 1)
the_big_switch(luts, cell, case, config)
# Done
return cell.get_vertices(), cell.get_faces(), cell.get_normals(), cell.get_values()
cdef void the_big_switch(LutProvider luts, Cell cell, int case, int config):
""" The big switch (i.e. if-statement) that I meticulously ported from
the source code provided by Lewiner et. al.
Together with all the look-up tables, this is where the magic is ...
"""
cdef int subconfig = 0
# Sinatures for tests
#test_face(cell, luts.TESTX.get1(config)):
#test_internal(cell, luts, case, config, subconfig, luts.TESTX.get1(config)):
#cell.add_triangles(luts.TILINGX, config, N)
if case == 1:
cell.add_triangles(luts.TILING1, config, 1)
elif case == 2:
cell.add_triangles(luts.TILING2, config, 2)
elif case == 3:
if test_face(cell, luts.TEST3.get1(config)):
cell.add_triangles(luts.TILING3_2, config, 4)
else:
cell.add_triangles(luts.TILING3_1, config, 2)
elif case == 4 :
if test_internal(cell, luts, case, config, subconfig, luts.TEST4.get1(config)):
cell.add_triangles(luts.TILING4_1, config, 2)
else:
cell.add_triangles(luts.TILING4_2, config, 6)
elif case == 5 :
cell.add_triangles(luts.TILING5, config, 3)
elif case == 6 :
if test_face(cell, luts.TEST6.get2(config,0)):
cell.add_triangles(luts.TILING6_2, config, 5)
else:
if test_internal(cell, luts, case, config, subconfig, luts.TEST6.get2(config,1)):
cell.add_triangles(luts.TILING6_1_1, config, 3)
else:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles(luts.TILING6_1_2, config, 9)
elif case == 7 :
# Get subconfig
if test_face(cell, luts.TEST7.get2(config,0)): subconfig += 1
if test_face(cell, luts.TEST7.get2(config,1)): subconfig += 2
if test_face(cell, luts.TEST7.get2(config,2)): subconfig += 4
# Behavior depends on subconfig
if subconfig == 0: cell.add_triangles(luts.TILING7_1, config, 3)
elif subconfig == 1: cell.add_triangles2(luts.TILING7_2, config, 0, 5)
elif subconfig == 2: cell.add_triangles2(luts.TILING7_2, config, 1, 5)
elif subconfig == 3:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING7_3, config, 0, 9)
elif subconfig == 4: cell.add_triangles2(luts.TILING7_2, config, 2, 5)
elif subconfig == 5:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING7_3, config, 1, 9)
elif subconfig == 6:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING7_3, config, 2, 9)
elif subconfig == 7:
if test_internal(cell, luts, case, config, subconfig, luts.TEST7.get2(config,3)):
cell.add_triangles(luts.TILING7_4_2, config, 9)
else:
cell.add_triangles(luts.TILING7_4_1, config, 5)
elif case == 8 :
cell.add_triangles(luts.TILING8, config, 2)
elif case == 9 :
cell.add_triangles(luts.TILING9, config, 4)
elif case == 10 :
if test_face(cell, luts.TEST10.get2(config,0)):
if test_face(cell, luts.TEST10.get2(config,1)):
cell.add_triangles(luts.TILING10_1_1_, config, 4)
else:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles(luts.TILING10_2, config, 8)
else:
if test_face(cell, luts.TEST10.get2(config,1)):
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles(luts.TILING10_2_, config, 8)
else:
if test_internal(cell, luts, case, config, subconfig, luts.TEST10.get2(config,2)):
cell.add_triangles(luts.TILING10_1_1, config, 4)
else:
cell.add_triangles(luts.TILING10_1_2, config, 8)
elif case == 11 :
cell.add_triangles(luts.TILING11, config, 4)
elif case == 12 :
if test_face(cell, luts.TEST12.get2(config,0)):
if test_face(cell, luts.TEST12.get2(config,1)):
cell.add_triangles(luts.TILING12_1_1_, config, 4)
else:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles(luts.TILING12_2, config, 8)
else:
if test_face(cell, luts.TEST12.get2(config,1)):
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles(luts.TILING12_2_, config, 8)
else:
if test_internal(cell, luts, case, config, subconfig, luts.TEST12.get2(config,2)):
cell.add_triangles(luts.TILING12_1_1, config, 4)
else:
cell.add_triangles(luts.TILING12_1_2, config, 8)
elif case == 13 :
# Calculate subconfig
if test_face(cell, luts.TEST13.get2(config,0)): subconfig += 1
if test_face(cell, luts.TEST13.get2(config,1)): subconfig += 2
if test_face(cell, luts.TEST13.get2(config,2)): subconfig += 4
if test_face(cell, luts.TEST13.get2(config,3)): subconfig += 8
if test_face(cell, luts.TEST13.get2(config,4)): subconfig += 16
if test_face(cell, luts.TEST13.get2(config,5)): subconfig += 32
# Map via LUT
subconfig = luts.SUBCONFIG13.get1(subconfig)
# Behavior depends on subconfig
if subconfig==0: cell.add_triangles(luts.TILING13_1, config, 4)
elif subconfig==1: cell.add_triangles2(luts.TILING13_2, config, 0, 6)
elif subconfig==2: cell.add_triangles2(luts.TILING13_2, config, 1, 6)
elif subconfig==3: cell.add_triangles2(luts.TILING13_2, config, 2, 6)
elif subconfig==4: cell.add_triangles2(luts.TILING13_2, config, 3, 6)
elif subconfig==5: cell.add_triangles2(luts.TILING13_2, config, 4, 6)
elif subconfig==6: cell.add_triangles2(luts.TILING13_2, config, 5, 6)
#
elif subconfig==7:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 0, 10)
elif subconfig==8:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 1, 10)
elif subconfig==9:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 2, 10)
elif subconfig==10:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 3, 10)
elif subconfig==11:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 4, 10)
elif subconfig==12:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 5, 10)
elif subconfig==13:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 6, 10)
elif subconfig==14:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 7, 10)
elif subconfig==15:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 8, 10)
elif subconfig==16:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 9, 10)
elif subconfig==17:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 10, 10)
elif subconfig==18:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3, config, 11, 10)
#
elif subconfig==19:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_4, config, 0, 12)
elif subconfig==20:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_4, config, 1, 12)
elif subconfig==21:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_4, config, 2, 12)
elif subconfig==22:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_4, config, 3, 12)
#
elif subconfig==23:
subconfig = 0 # Note: the original source code sets the subconfig, without apparent reason
if test_internal(cell, luts, case, config, subconfig, luts.TEST13.get2(config,6)):
cell.add_triangles2(luts.TILING13_5_1, config, 0, 6)
else:
cell.add_triangles2(luts.TILING13_5_2, config, 0, 10)
elif subconfig==24:
subconfig = 1
if test_internal(cell, luts, case, config, subconfig, luts.TEST13.get2(config,6)):
cell.add_triangles2(luts.TILING13_5_1, config, 1, 6)
else:
cell.add_triangles2(luts.TILING13_5_2, config, 1, 10)
elif subconfig==25:
subconfig = 2 ;
if test_internal(cell, luts, case, config, subconfig, luts.TEST13.get2(config,6)):
cell.add_triangles2(luts.TILING13_5_1, config, 2, 6)
else:
cell.add_triangles2(luts.TILING13_5_2, config, 2, 10)
elif subconfig==26:
subconfig = 3 ;
if test_internal(cell, luts, case, config, subconfig, luts.TEST13.get2(config,6)):
cell.add_triangles2(luts.TILING13_5_1, config, 3, 6)
else:
cell.add_triangles2(luts.TILING13_5_2, config, 3, 10)
#
elif subconfig==27:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 0, 10)
elif subconfig==28:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 1, 10)
elif subconfig==29:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 2, 10)
elif subconfig==30:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 3, 10)
elif subconfig==31:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 4, 10)
elif subconfig==32:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 5, 10)
elif subconfig==33:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config,6, 10)
elif subconfig==34:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 7, 10)
elif subconfig==35:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 8, 10)
elif subconfig==36:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 9, 10)
elif subconfig==37:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 10, 10)
elif subconfig==38:
#cell.calculate_center_vertex() # v12 needed
cell.add_triangles2(luts.TILING13_3_, config, 11, 10)
#
elif subconfig==39:
cell.add_triangles2(luts.TILING13_2_, config, 0, 6)
elif subconfig==40:
cell.add_triangles2(luts.TILING13_2_, config, 1, 6)
elif subconfig==41:
cell.add_triangles2(luts.TILING13_2_, config, 2, 6)
elif subconfig==42:
cell.add_triangles2(luts.TILING13_2_, config, 3, 6)
elif subconfig==43:
cell.add_triangles2(luts.TILING13_2_, config, 4, 6)
elif subconfig==44:
cell.add_triangles2(luts.TILING13_2_, config, 5, 6)
#
elif subconfig==45:
cell.add_triangles(luts.TILING13_1_, config, 4)
#
else:
print("Marching Cubes: Impossible case 13?" )
elif case == 14 :
cell.add_triangles(luts.TILING14, config, 4)
cdef int test_face(Cell cell, int face):
""" Return True of the face contains part of the surface.
"""
# Get face absolute value
cdef int absFace = face
if face < 0:
absFace *= -1
# Get values of corners A B C D
cdef double A, B, C, D
if absFace == 1:
A, B, C, D = cell.v0, cell.v4, cell.v5, cell.v1
elif absFace == 2:
A, B, C, D = cell.v1, cell.v5, cell.v6, cell.v2
elif absFace == 3:
A, B, C, D = cell.v2, cell.v6, cell.v7, cell.v3
elif absFace == 4:
A, B, C, D = cell.v3, cell.v7, cell.v4, cell.v0
elif absFace == 5:
A, B, C, D = cell.v0, cell.v3, cell.v2, cell.v1
elif absFace == 6:
A, B, C, D = cell.v4, cell.v7, cell.v6, cell.v5
# Return sign
cdef double AC_BD = A*C - B*D
if AC_BD > - FLT_EPSILON and AC_BD < FLT_EPSILON:
return face >= 0
else:
return face * A * AC_BD >= 0; # face and A invert signs
cdef int test_internal(Cell cell, LutProvider luts, int case, int config, int subconfig, int s):
""" Return True of the face contains part of the surface.
"""
# Typedefs
cdef double t, At, Bt, Ct, Dt, a, b
cdef int test = 0
cdef int edge = -1 # reference edge of the triangulation
# Calculate At Bt Ct Dt a b
# Select case 4, 10, 7, 12, 13
At, Bt, Ct, Dt = 0.0, 0.0, 0.0, 0.0
if case==4 or case==10:
a = ( cell.v4 - cell.v0 ) * ( cell.v6 - cell.v2 ) - ( cell.v7 - cell.v3 ) * ( cell.v5 - cell.v1 )
b = cell.v2 * ( cell.v4 - cell.v0 ) + cell.v0 * ( cell.v6 - cell.v2 ) - cell.v1 * ( cell.v7 - cell.v3 ) - cell.v3 * ( cell.v5 - cell.v1 )
t = - b / (2*a + FLT_EPSILON)
if t<0 or t>1: return s>0 ;
At = cell.v0 + ( cell.v4 - cell.v0 ) * t
Bt = cell.v3 + ( cell.v7 - cell.v3 ) * t
Ct = cell.v2 + ( cell.v6 - cell.v2 ) * t
Dt = cell.v1 + ( cell.v5 - cell.v1 ) * t
elif case==6 or case==7 or case==12 or case==13:
# Define edge
if case == 6: edge = luts.TEST6.get2(config, 2)
elif case == 7: edge = luts.TEST7.get2(config, 4)
elif case == 12: edge = luts.TEST12.get2(config, 3)
elif case == 13: edge = luts.TILING13_5_1.get3(config, subconfig, 0)
if edge==0:
t = cell.v0 / ( cell.v0 - cell.v1 + FLT_EPSILON )
At = 0
Bt = cell.v3 + ( cell.v2 - cell.v3 ) * t
Ct = cell.v7 + ( cell.v6 - cell.v7 ) * t
Dt = cell.v4 + ( cell.v5 - cell.v4 ) * t
elif edge==1:
t = cell.v1 / ( cell.v1 - cell.v2 + FLT_EPSILON )
At = 0
Bt = cell.v0 + ( cell.v3 - cell.v0 ) * t
Ct = cell.v4 + ( cell.v7 - cell.v4 ) * t
Dt = cell.v5 + ( cell.v6 - cell.v5 ) * t
elif edge==2:
t = cell.v2 / ( cell.v2 - cell.v3 + FLT_EPSILON )
At = 0
Bt = cell.v1 + ( cell.v0 - cell.v1 ) * t
Ct = cell.v5 + ( cell.v4 - cell.v5 ) * t
Dt = cell.v6 + ( cell.v7 - cell.v6 ) * t
elif edge==3:
t = cell.v3 / ( cell.v3 - cell.v0 + FLT_EPSILON )
At = 0
Bt = cell.v2 + ( cell.v1 - cell.v2 ) * t
Ct = cell.v6 + ( cell.v5 - cell.v6 ) * t
Dt = cell.v7 + ( cell.v4 - cell.v7 ) * t
elif edge==4:
t = cell.v4 / ( cell.v4 - cell.v5 + FLT_EPSILON )
At = 0
Bt = cell.v7 + ( cell.v6 - cell.v7 ) * t
Ct = cell.v3 + ( cell.v2 - cell.v3 ) * t
Dt = cell.v0 + ( cell.v1 - cell.v0 ) * t
elif edge==5:
t = cell.v5 / ( cell.v5 - cell.v6 + FLT_EPSILON )
At = 0
Bt = cell.v4 + ( cell.v7 - cell.v4 ) * t
Ct = cell.v0 + ( cell.v3 - cell.v0 ) * t
Dt = cell.v1 + ( cell.v2 - cell.v1 ) * t
elif edge==6:
t = cell.v6 / ( cell.v6 - cell.v7 + FLT_EPSILON )
At = 0
Bt = cell.v5 + ( cell.v4 - cell.v5 ) * t
Ct = cell.v1 + ( cell.v0 - cell.v1 ) * t
Dt = cell.v2 + ( cell.v3 - cell.v2 ) * t
elif edge==7:
t = cell.v7 / ( cell.v7 - cell.v4 + FLT_EPSILON )
At = 0
Bt = cell.v6 + ( cell.v5 - cell.v6 ) * t
Ct = cell.v2 + ( cell.v1 - cell.v2 ) * t
Dt = cell.v3 + ( cell.v0 - cell.v3 ) * t
elif edge==8:
t = cell.v0 / ( cell.v0 - cell.v4 + FLT_EPSILON )
At = 0
Bt = cell.v3 + ( cell.v7 - cell.v3 ) * t
Ct = cell.v2 + ( cell.v6 - cell.v2 ) * t
Dt = cell.v1 + ( cell.v5 - cell.v1 ) * t
elif edge==9:
t = cell.v1 / ( cell.v1 - cell.v5 + FLT_EPSILON )
At = 0
Bt = cell.v0 + ( cell.v4 - cell.v0 ) * t
Ct = cell.v3 + ( cell.v7 - cell.v3 ) * t
Dt = cell.v2 + ( cell.v6 - cell.v2 ) * t
elif edge==10:
t = cell.v2 / ( cell.v2 - cell.v6 + FLT_EPSILON )
At = 0
Bt = cell.v1 + ( cell.v5 - cell.v1 ) * t
Ct = cell.v0 + ( cell.v4 - cell.v0 ) * t
Dt = cell.v3 + ( cell.v7 - cell.v3 ) * t
elif edge==11:
t = cell.v3 / ( cell.v3 - cell.v7 + FLT_EPSILON )
At = 0
Bt = cell.v2 + ( cell.v6 - cell.v2 ) * t
Ct = cell.v1 + ( cell.v5 - cell.v1 ) * t
Dt = cell.v0 + ( cell.v4 - cell.v0 ) * t
else:
print( "Invalid edge %i." % edge )
else:
print( "Invalid ambiguous case %i." % case )
# Process results
if At >= 0: test += 1
if Bt >= 0: test += 2
if Ct >= 0: test += 4
if Dt >= 0: test += 8
# Determine what to return
if test==0: return s>0
elif test==1: return s>0
elif test==2: return s>0
elif test==3: return s>0
elif test==4: return s>0
elif test==5:
if At * Ct - Bt * Dt < FLT_EPSILON: return s>0
elif test==6: return s>0
elif test==7: return s<0
elif test==8: return s>0
elif test==9: return s>0
elif test==10:
if At * Ct - Bt * Dt >= FLT_EPSILON: return s>0
elif test==11: return s<0
elif test==12: return s>0
elif test==13: return s<0
elif test==14: return s<0
elif test==15: return s<0
else: return s<0
Reference in New Issue View Git Blame Copy Permalink