#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 np.import_array() # Enable low level memory management from libc.stdlib cimport malloc, free # 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 a0]]] 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 = malloc(self.nx*self.ny*4 * sizeof(int)) self.faceLayer2 = malloc(self.nx*self.ny*4 * sizeof(int)) if (self.faceLayer1 is NULL or self.faceLayer2 is NULL or self.vv is NULL or self.vg is NULL or self._vertices is NULL or self._normals is NULL or self._values is NULL or self._faces is NULL): raise MemoryError() 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 = malloc(8 * sizeof(double)) self.vg = malloc(8*3 * sizeof(double)) # Init face layers self.faceLayer1 = NULL self.faceLayer2 = NULL # Init vertices self._vertexCount = 0 self._vertexMaxCount = 8 self._vertices = malloc(self._vertexMaxCount*3 * sizeof(float)) self._normals = malloc(self._vertexMaxCount*3 * sizeof(float)) self._values = malloc(self._vertexMaxCount * sizeof(float)) # Clear normals and values cdef int i, j if self._values is not NULL and self._normals is not NULL: 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 = malloc(self._faceMaxCount * sizeof(int)) def __dealloc__(self): free(self.vv) free(self.vg) free(self.faceLayer1) free(self.faceLayer2) free(self._vertices) free(self._normals) free(self._values) 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 = malloc(newMaxCount*3 * sizeof(float)) cdef float *newNormals = malloc(newMaxCount*3 * sizeof(float)) cdef float *newValues = malloc(newMaxCount * sizeof(float)) if newVertices is NULL or newNormals is NULL or newValues is NULL: free(newVertices) free(newNormals) free(newValues) raise MemoryError() # 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 = malloc(newMaxCount * sizeof(int)) if newFaces is NULL: raise MemoryError() # 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 = 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 += dx1 * tmpf1; fy += dy1 * tmpf1; fz += dz1 * tmpf1; ff += tmpf1 fx += dx2 * tmpf2; fy += dy2 * tmpf2; fz += dz2 * tmpf2; ff += tmpf2 # Add vertex indexInVertexArray = self.add_vertex( self.x + stp*fx/ff, self.y + stp*fy/ff, 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 = 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_zg = ( 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 signed 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 = malloc(N * sizeof(signed char)) if self.VALUES is NULL: raise MemoryError() 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, np.ndarray[np.npy_bool, ndim=3, cast=True] mask=None): """ marching_cubes(im, double isovalue, LutProvider luts, int st=1, int classic=0) Main entry to apply marching cubes. Masked version of marching cubes. This function will check a masking array (same size as im) to decide if the algorithm must be computed for a given voxel. This adds a small overhead that rapidly gets compensated by the fewer computed 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 cdef bint no_mask = mask is None # Unfortunately specifying a step in range() significantly 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 if no_mask or mask[z_st, y_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