diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index 264cf162..c03587d4 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -1,6 +1,7 @@ from ._find_contours import find_contours from ._marching_cubes import (marching_cubes, mesh_surface_area, correct_mesh_orientation) +from ._marching_cubes_lewiner import marching_cubes_lewiner from ._regionprops import regionprops, perimeter from .simple_metrics import compare_mse, compare_nrmse, compare_psnr from ._structural_similarity import compare_ssim, structural_similarity @@ -29,6 +30,7 @@ __all__ = ['find_contours', 'moments_normalized', 'moments_hu', 'marching_cubes', + 'marching_cubes_lewiner', 'mesh_surface_area', 'correct_mesh_orientation', 'profile_line', diff --git a/skimage/measure/_marching_cubes_lewiner.py b/skimage/measure/_marching_cubes_lewiner.py new file mode 100644 index 00000000..12cb5100 --- /dev/null +++ b/skimage/measure/_marching_cubes_lewiner.py @@ -0,0 +1,173 @@ +import sys +import base64 + +import numpy as np + +if sys.version_info[0] == 3: + base64decode = base64.decodebytes +else: + base64decode = base64.decodestring + +from . import _marching_cubes_lewiner_luts as mcluts +from . import _marching_cubes_lewiner_cy + + + +def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.), + step_size=1, use_classic=False): + """ + Lewiner marching cubes algorithm to find surfaces in 3d volumetric data + + In contrast to ``marching_cubes()``, this algorithm resolves + ambiguities and guarantees topologically correct results. + + Parameters + ---------- + volume : (M, N, P) array (the data is internally converted to + float32 if necessary) + level : float + Contour value to search for isosurfaces in `volume`. If not + given or None, the average of the min and max of vol is used. + spacing : length-3 tuple of floats + Voxel spacing in spatial dimensions corresponding to numpy array + indexing dimensions (M, N, P) as in `volume`. + step_size : int + Step size in voxels. Default 1. Larger steps yield faster but + coarser results. The result will always be topologically correct + though. + use_classic : bool + If given and True, the classic marching cubes by Lorensen (1987) + is used. This option is included for reference purposes. Note + that this algorithm has ambiguities and is not guaranteed to + produce a topologically correct result. + + Returns + ------- + verts : (V, 3) array + Spatial coordinates for V unique mesh vertices. Coordinate order + matches input `volume` (M, N, P). + faces : (F, 3) array + Define triangular faces via referencing vertex indices from ``verts``. + This algorithm specifically outputs triangles, so each face has + exactly three indices. + normals : (V, 3) array + The normal direction at each vertex, as calculated from the + data. + values : (V, ) array + Gives a measure for the maximum value of the data in the local region + near each vertex. This can be used by visualization tools to apply + a colormap to the mesh. + + Notes about the algorithm + ------------------------- + + This is an implementation of: + + 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) + + The algorithm is an improved version of Chernyaev's Marching Cubes 33 + algorithm, originally written in C++. It is an efficient algorithm + that relies on heavy use of lookup tables to handle the many different + cases. This keeps the algorithm relatively easy. The current algorithm + is a port of Lewiner's algorithm and written in Cython. + """ + + # Check volume and ensure its in the format that the alg needs + if not isinstance(volume, np.ndarray) or (volume.ndim != 3): + 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) + + # Check/convert other inputs: + # level + if level is None: + level = 0.5 * (volume.min() + volume.max()) + else: + level = float(level) + if level < volume.min() or level > volume.max(): + raise ValueError("Surface level must be within volume data range.") + # spacing + if len(spacing) != 3: + raise ValueError("`spacing` must consist of three floats.") + # step_size + step_size = int(step_size) + if step_size < 1: + raise ValueError('step_size must be at least one.') + # use_classic + use_classic = bool(use_classic) + + # Get LutProvider class (reuse if possible) + L = _getMCLuts() + + # Apply algorithm + func = _marching_cubes_lewiner_cy.marching_cubes + vertices, faces , normals, values = func(volume, level, L, step_size, use_classic) + + if not len(vertices): + raise RuntimeError('No surface found at the given iso value.') + + # Output in z-y-x order, as is common in skimage + vertices = np.fliplr(vertices) + normals = np.fliplr(normals) + + # Finishing touches to output + faces.shape = -1, 3 + if spacing != (1, 1, 1): + vertices = vertices * np.r_[spacing] + + return vertices, faces, normals, values + + +def _toArray(args): + shape, text = args + byts = base64decode(text.encode('utf-8')) + ar = np.frombuffer(byts, dtype='int8') + ar.shape = shape + return ar + + +# Map an edge-index to two relative pixel positions. The ege index +# represents a point that lies somewhere in between these pixels. +# Linear interpolation should be used to determine where it is exactly. +# 0 +# 3 1 -> 0x +# 2 xx +EDGETORELATIVEPOSX = np.array([ [0,1],[1,1],[1,0],[0,0], [0,1],[1,1],[1,0],[0,0], [0,0],[1,1],[1,1],[0,0] ], 'int8') +EDGETORELATIVEPOSY = np.array([ [0,0],[0,1],[1,1],[1,0], [0,0],[0,1],[1,1],[1,0], [0,0],[0,0],[1,1],[1,1] ], 'int8') +EDGETORELATIVEPOSZ = np.array([ [0,0],[0,0],[0,0],[0,0], [1,1],[1,1],[1,1],[1,1], [0,1],[0,1],[0,1],[0,1] ], 'int8') + + +def _getMCLuts(): + """ Kind of lazy obtaining of the luts. + """ + if not hasattr(mcluts, 'THE_LUTS'): + + mcluts.THE_LUTS = _marching_cubes_lewiner_cy.LutProvider( + EDGETORELATIVEPOSX, EDGETORELATIVEPOSY, EDGETORELATIVEPOSZ, + + _toArray(mcluts.CASESCLASSIC), _toArray(mcluts.CASES), + + _toArray(mcluts.TILING1), _toArray(mcluts.TILING2), _toArray(mcluts.TILING3_1), _toArray(mcluts.TILING3_2), + _toArray(mcluts.TILING4_1), _toArray(mcluts.TILING4_2), _toArray(mcluts.TILING5), _toArray(mcluts.TILING6_1_1), + _toArray(mcluts.TILING6_1_2), _toArray(mcluts.TILING6_2), _toArray(mcluts.TILING7_1), + _toArray(mcluts.TILING7_2), _toArray(mcluts.TILING7_3), _toArray(mcluts.TILING7_4_1), + _toArray(mcluts.TILING7_4_2), _toArray(mcluts.TILING8), _toArray(mcluts.TILING9), + _toArray(mcluts.TILING10_1_1), _toArray(mcluts.TILING10_1_1_), _toArray(mcluts.TILING10_1_2), + _toArray(mcluts.TILING10_2), _toArray(mcluts.TILING10_2_), _toArray(mcluts.TILING11), + _toArray(mcluts.TILING12_1_1), _toArray(mcluts.TILING12_1_1_), _toArray(mcluts.TILING12_1_2), + _toArray(mcluts.TILING12_2), _toArray(mcluts.TILING12_2_), _toArray(mcluts.TILING13_1), + _toArray(mcluts.TILING13_1_), _toArray(mcluts.TILING13_2), _toArray(mcluts.TILING13_2_), + _toArray(mcluts.TILING13_3), _toArray(mcluts.TILING13_3_), _toArray(mcluts.TILING13_4), + _toArray(mcluts.TILING13_5_1), _toArray(mcluts.TILING13_5_2), _toArray(mcluts.TILING14), + + _toArray(mcluts.TEST3), _toArray(mcluts.TEST4), _toArray(mcluts.TEST6), + _toArray(mcluts.TEST7), _toArray(mcluts.TEST10), _toArray(mcluts.TEST12), + _toArray(mcluts.TEST13), _toArray(mcluts.SUBCONFIG13), + ) + + return mcluts.THE_LUTS + diff --git a/skimage/measure/_marching_cubes_lewiner_cy.pyx b/skimage/measure/_marching_cubes_lewiner_cy.pyx new file mode 100644 index 00000000..04fe49f2 --- /dev/null +++ b/skimage/measure/_marching_cubes_lewiner_cy.pyx @@ -0,0 +1,1416 @@ +# -*- coding: utf-8 -*- +# Copyright (C) 2012, Almar Klein +# Copyright (C) 2002, Thomas Lewiner + +""" +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) + +I selected this algorithm because it provides topologically correct results, +and because the algorithms implementation is relatively simple. Most of +the magic is in the lookup tables, which are provided as open source. + +This code is distributed under the terms of the (new) BSD License. + +""" + +# 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) + +# Type defs, we support float32 and float64 +ctypedef np.float32_t FLOAT32_T +ctypedef np.float64_t FLOAT64_T +ctypedef np.int32_t INT32_T +FLOAT32 = np.float32 +FLOAT64 = np.float64 +INT32 = np.float32 + + +# 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 + + +# todo: allow dynamic isovalue? +# todo: can we disable Cython from checking for zero division? Sometimes we know that it never happens! + + +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)) + 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 + 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): + 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 = malloc(newMaxCount*3 * sizeof(float)) + cdef float *newNormals = malloc(newMaxCount*3 * sizeof(float)) + cdef float *newValues = 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 = 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),'float32') + cdef np.ndarray[FLOAT32_T, ndim=2] 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),'float32') + cdef np.ndarray[FLOAT32_T, ndim=2] 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,),'int32') + cdef np.ndarray[INT32_T, ndim=1] 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,),'float32') + cdef np.ndarray[FLOAT32_T, ndim=1] 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])) + + 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 edes 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 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: + # four 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_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 = 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) + + + +@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. + Returns (vertices, faces, normals, values) + """ + + # Typdedef image + cdef np.ndarray[FLOAT32_T, ndim=3] im_ = im + + # 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 + cdef int nt + cdef int case, config, subconfig + + # Unfortunately, using a step in the for loops degregades performance. + # Quick and dirty trick below ... + + if st == 1: + + for z in range(0,Nz-1): + cell.new_z_value() # Indicate that we enter a new layer + for y in range(0,Ny-1): + for x in range(0,Nx-1): + + # 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) + + else: + + for z in range(0,Nz-st,st): + cell.new_z_value() # Indicate that we enter a new layer + for y in range(0,Ny-st,st): + for x in range(0,Nx-st,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 diff --git a/skimage/measure/_marching_cubes_lewiner_luts.py b/skimage/measure/_marching_cubes_lewiner_luts.py new file mode 100644 index 00000000..59243662 --- /dev/null +++ b/skimage/measure/_marching_cubes_lewiner_luts.py @@ -0,0 +1,529 @@ +# -*- coding: utf-8 -*- +# Copyright (C) 2012, Almar Klein +# Copyright (C) 2002, Thomas Lewiner + +# This file was auto-generated from LookUpTable.h by createluts.py. + +#static const char casesClassic[256][16] +CASESCLASSIC = (256, 16), """ +/////////////////////wAIA/////////////////8AAQn/////////////////AQgDCQgB//// +/////////wECCv////////////////8ACAMBAgr/////////////CQIKAAIJ/////////////wII +AwIKCAoJCP////////8DCwL/////////////////AAsCCAsA/////////////wEJAAIDC/////// +//////8BCwIBCQsJCAv/////////AwoBCwoD/////////////wAKAQAICggLCv////////8DCQAD +CwkLCgn/////////CQgKCggL/////////////wQHCP////////////////8EAwAHAwT///////// +////AAEJCAQH/////////////wQBCQQHAQcDAf////////8BAgoIBAf/////////////AwQHAwAE +AQIK/////////wkCCgkAAggEB/////////8CCgkCCQcCBwMHCQT/////CAQHAwsC//////////// +/wsEBwsCBAIABP////////8JAAEIBAcCAwv/////////BAcLCQQLCQsCCQIB/////wMKAQMLCgcI +BP////////8BCwoBBAsBAAQHCwT/////BAcICQALCQsKCwAD/////wQHCwQLCQkLCv////////8J +BQT/////////////////CQUEAAgD/////////////wAFBAEFAP////////////8IBQQIAwUDAQX/ +////////AQIKCQUE/////////////wMACAECCgQJBf////////8FAgoFBAIEAAL/////////AgoF +AwIFAwUEAwQI/////wkFBAIDC/////////////8ACwIACAsECQX/////////AAUEAAEFAgML//// +/////wIBBQIFCAIICwQIBf////8KAwsKAQMJBQT/////////BAkFAAgBCAoBCAsK/////wUEAAUA +CwULCgsAA/////8FBAgFCAoKCAv/////////CQcIBQcJ/////////////wkDAAkFAwUHA/////// +//8ABwgAAQcBBQf/////////AQUDAwUH/////////////wkHCAkFBwoBAv////////8KAQIJBQAF +AwAFBwP/////CAACCAIFCAUHCgUC/////wIKBQIFAwMFB/////////8HCQUHCAkDCwL///////// +CQUHCQcCCQIAAgcL/////wIDCwABCAEHCAEFB/////8LAgELAQcHAQX/////////CQUICAUHCgED +CgML/////wUHAAUACQcLAAEACgsKAP8LCgALAAMKBQAIAAcFBwD/CwoFBwsF/////////////woG +Bf////////////////8ACAMFCgb/////////////CQABBQoG/////////////wEIAwEJCAUKBv// +//////8BBgUCBgH/////////////AQYFAQIGAwAI/////////wkGBQkABgACBv////////8FCQgF +CAIFAgYDAgj/////AgMLCgYF/////////////wsACAsCAAoGBf////////8AAQkCAwsFCgb///// +////BQoGAQkCCQsCCQgL/////wYDCwYFAwUBA/////////8ACAsACwUABQEFCwb/////AwsGAAMG +AAYFAAUJ/////wYFCQYJCwsJCP////////8FCgYEBwj/////////////BAMABAcDBgUK//////// +/wEJAAUKBggEB/////////8KBgUBCQcBBwMHCQT/////BgECBgUBBAcI/////////wECBQUCBgMA +BAMEB/////8IBAcJAAUABgUAAgb/////BwMJBwkEAwIJBQkGAgYJ/wMLAgcIBAoGBf////////8F +CgYEBwIEAgACBwv/////AAEJBAcIAgMLBQoG/////wkCAQkLAgkECwcLBAUKBv8IBAcDCwUDBQEF +Cwb/////BQELBQsGAQALBwsEAAQL/wAFCQAGBQADBgsGAwgEB/8GBQkGCQsEBwkHCwn/////CgQJ +BgQK/////////////wQKBgQJCgAIA/////////8KAAEKBgAGBAD/////////CAMBCAEGCAYEBgEK +/////wEECQECBAIGBP////////8DAAgBAgkCBAkCBgT/////AAIEBAIG/////////////wgDAggC +BAQCBv////////8KBAkKBgQLAgP/////////AAgCAggLBAkKBAoG/////wMLAgABBgAGBAYBCv// +//8GBAEGAQoECAECAQsICwH/CQYECQMGCQEDCwYD/////wgLAQgBAAsGAQkBBAYEAf8DCwYDBgAA +BgT/////////BgQICwYI/////////////wcKBgcICggJCv////////8ABwMACgcACQoGBwr///// +CgYHAQoHAQcIAQgA/////woGBwoHAQEHA/////////8BAgYBBggBCAkIBgf/////AgYJAgkBBgcJ +AAkDBwMJ/wcIAAcABgYAAv////////8HAwIGBwL/////////////AgMLCgYICggJCAYH/////wIA +BwIHCwAJBwYHCgkKB/8BCAABBwgBCgcGBwoCAwv/CwIBCwEHCgYBBgcB/////wgJBggGBwkBBgsG +AwEDBv8ACQELBgf/////////////BwgABwAGAwsACwYA/////wcLBv////////////////8HBgv/ +////////////////AwAICwcG/////////////wABCQsHBv////////////8IAQkIAwELBwb///// +////CgECBgsH/////////////wECCgMACAYLB/////////8CCQACCgkGCwf/////////BgsHAgoD +CggDCgkI/////wcCAwYCB/////////////8HAAgHBgAGAgD/////////AgcGAgMHAAEJ//////// +/wEGAgEIBgEJCAgHBv////8KBwYKAQcBAwf/////////CgcGAQcKAQgHAQAI/////wADBwAHCgAK +CQYKB/////8HBgoHCggICgn/////////BggECwgG/////////////wMGCwMABgAEBv////////8I +BgsIBAYJAAH/////////CQQGCQYDCQMBCwMG/////wYIBAYLCAIKAf////////8BAgoDAAsABgsA +BAb/////BAsIBAYLAAIJAgoJ/////woJAwoDAgkEAwsDBgQGA/8IAgMIBAIEBgL/////////AAQC +BAYC/////////////wEJAAIDBAIEBgQDCP////8BCQQBBAICBAb/////////CAEDCAYBCAQGBgoB +/////woBAAoABgYABP////////8EBgMEAwgGCgMAAwkKCQP/CgkEBgoE/////////////wQJBQcG +C/////////////8ACAMECQULBwb/////////BQABBQQABwYL/////////wsHBggDBAMFBAMBBf// +//8JBQQKAQIHBgv/////////BgsHAQIKAAgDBAkF/////wcGCwUECgQCCgQAAv////8DBAgDBQQD +AgUKBQILBwb/BwIDBwYCBQQJ/////////wkFBAAIBgAGAgYIB/////8DBgIDBwYBBQAFBAD///// +BgIIBggHAgEIBAgFAQUI/wkFBAoBBgEHBgEDB/////8BBgoBBwYBAAcIBwAJBQT/BAAKBAoFAAMK +BgoHAwcK/wcGCgcKCAUECgQICv////8GCQUGCwkLCAn/////////AwYLAAYDAAUGAAkF/////wAL +CAAFCwABBQUGC/////8GCwMGAwUFAwH/////////AQIKCQULCQsICwUG/////wALAwAGCwAJBgUG +CQECCv8LCAULBQYIAAUKBQIAAgX/BgsDBgMFAgoDCgUD/////wUICQUCCAUGAgMIAv////8JBQYJ +BgAABgL/////////AQUIAQgABQYIAwgCBgII/wEFBgIBBv////////////8BAwYBBgoDCAYFBgkI +CQb/CgEACgAGCQUABQYA/////wADCAUGCv////////////8KBQb/////////////////CwUKBwUL +/////////////wsFCgsHBQgDAP////////8FCwcFCgsBCQD/////////CgcFCgsHCQgBCAMB//// +/wsBAgsHAQcFAf////////8ACAMBAgcBBwUHAgv/////CQcFCQIHCQACAgsH/////wcFAgcCCwUJ +AgMCCAkIAv8CBQoCAwUDBwX/////////CAIACAUCCAcFCgIF/////wkAAQUKAwUDBwMKAv////8J +CAIJAgEIBwIKAgUHBQL/AQMFAwcF/////////////wAIBwAHAQEHBf////////8JAAMJAwUFAwf/ +////////CQgHBQkH/////////////wUIBAUKCAoLCP////////8FAAQFCwAFCgsLAwD/////AAEJ +CAQKCAoLCgQF/////woLBAoEBQsDBAkEAQMBBP8CBQECCAUCCwgEBQj/////AAQLAAsDBAULAgsB +BQEL/wACBQAFCQILBQQFCAsIBf8JBAUCCwP/////////////AgUKAwUCAwQFAwgE/////wUKAgUC +BAQCAP////////8DCgIDBQoDCAUEBQgAAQn/BQoCBQIEAQkCCQQC/////wgEBQgFAwMFAf////// +//8ABAUBAAX/////////////CAQFCAUDCQAFAAMF/////wkEBf////////////////8ECwcECQsJ +Cgv/////////AAgDBAkHCQsHCQoL/////wEKCwELBAEEAAcEC/////8DAQQDBAgBCgQHBAsKCwT/ +BAsHCQsECQILCQEC/////wkHBAkLBwkBCwILAQAIA/8LBwQLBAICBAD/////////CwcECwQCCAME +AwIE/////wIJCgIHCQIDBwcECf////8JCgcJBwQKAgcIBwACAAf/AwcKAwoCBwQKAQoABAAK/wEK +AggHBP////////////8ECQEEAQcHAQP/////////BAkBBAEHAAgBCAcB/////wQAAwcEA/////// +//////8ECAf/////////////////CQoICgsI/////////////wMACQMJCwsJCv////////8AAQoA +CggICgv/////////AwEKCwMK/////////////wECCwELCQkLCP////////8DAAkDCQsBAgkCCwn/ +////AAILCAAL/////////////wMCC/////////////////8CAwgCCAoKCAn/////////CQoCAAkC +/////////////wIDCAIICgABCAEKCP////8BCgL/////////////////AQMICQEI//////////// +/wAJAf////////////////8AAwj//////////////////////////////////////w== +""" + +#static const char cases[256][2] +CASES = (256, 2), """ +AP8BAAEBAgABAgMAAgMFAAEDAgEDAwUBAgUFBAUJCAABBAICAwQFAgQCBgIGCQsAAwgFBQcDCQEG +EA4DDAwFGAEFAwECBAUDAwYHAAUKCQAEAwYEBgsOAQYRDAQLBgUZAggFBwUMCAEGEgwFDgcFHAYV +CwQMDwUeCgUGIAYnAgwBBgQAAwUGAAIGBgMFCw4AAwkGBQcEDAEFDgsDCQQFGgMKBgYHBQwCBhMK +AQwNBhgHBwwJDQEHCQwUBiEHDQMMAgoGBwUNCwIFEAwHCAMFHQYWCgIMEQYbDgkGIgUnAg4FFA4F +CQUFIAsKBiMFKQIQDBcGJQcOAxAGLgQGAxUBCAEHAwIEAQYBAwcHAQYKDAACBwUGBgwLAQUPCQIO +BgUbAgkFCAYNDgIGFAwGCgMGGQUSCAIMEAUfCwkFIgYoAg0DCwcCBg4MAwcGDQAMDgcIBhcMCgoE +BhwMFQcKBikDDQUVCQMLCAUhDBYHCwYqAw4OCwUkBiwCEQYvAxIEBwEJAgsGCAYPCgAFEQwICwcG +GgUTDgQMEgYdCAQFIwUoAg8FFgsFDBMGHg4KBiQGKwQECQcFJQcPAxEFLAITAxYBCgUXDAsOCAYf +CQYHDAUqAw8LCwYmBi0EBQUtAxMCFQELCAUFJgUrAhIFLgMUAhYBDAUvAhQDFwENAhcBDgEPAP8= +""" + +#static const char tiling1[16][3] +TILING1 = (16, 3), """ +AAgDAAEJAQIKAwsCBAcICQUECgYFBwYLBwsGCgUGCQQFBAgHAwILAQoCAAkBAAMI +""" + +#static const char tiling2[24][6] +TILING2 = (24, 6), """ +AQgDCQgBAAsCCAsABAMABwMECQIKAAIJAAUEAQUAAwoBCwoDAQYFAgYBBwIDBgIHCQcIBQcJBggE +CwgGCgQJBgQKCwUKBwULCwoFBwsFCgkEBgoEBgQICwYICQgHBQkHBwMCBgcCAQUGAgEGAwEKCwMK +AAQFAQAFCQoCAAkCBAADBwQDAAILCAALAQMICQEI +""" + +#static const char tiling3_1[24][6] +TILING3_1 = (24, 6), """ +AAgDAQIKCQUEAAgDAwAICwcGAQkAAgMLAAEJCAQHCQABBQoGAQIKCQUECgECBgsHCAQHAwsCAgML +CgYFBQoGBAcIBAkFBwYLBQkECwYHBgoFCAcECwMCBQYKBwQIAgsDAgEKBwsGCgIBBAUJAQAJBgoF +CQEABwQIAAkBCwMCCAADBgcLBAUJAwgAAwgACgIB +""" + +#static const char tiling3_2[24][12] +TILING3_2 = (24, 12), """ +CgMCCggDCgEACAoAAwQIAwUEAwAJBQMJBggHBgAIBgsDAAYDCwADCwkACwIBCQsBBwkEBwEJBwgA +AQcABgEKBgABCQAGCQYFBAoFBAIKBAkBAgQBBwILBwECBwYKAQcKAgcLAgQHAgMIBAIIBQsGBQML +BQoCAwUCCAYHCAoGCAQFCggFCwUGCwkFCwcECQsEBgULBQkLBAcLBAsJBwYIBgoIBQQIBQgKBgsF +CwMFAgoFAgUDCwcCBwQCCAMCCAIECwIHAgEHCgYHCgcBBQoECgIEAQkEAQQCCgEGAQAGBgAJBQYJ +BAkHCQEHAAgHAAcBAwALAAkLAQILAQsJBwgGCAAGAwsGAwYACAQDBAUDCQADCQMFAgMKAwgKAAEK +AAoI +""" + +#static const char tiling4_1[8][6] +TILING4_1 = (8, 6), """ +AAgDBQoGAAEJCwcGAQIKCAQHCQUEAgMLBAUJCwMCCgIBBwQICQEABgcLAwgABgoF +""" + +#static const char tiling4_2[8][18] +TILING4_2 = (8, 18), """ +CAUABQgGAwYIBgMKAAoDCgAFCQYBBgkHAAcJBwALAQsACwEGCgcCBwoEAQQKBAEIAggBCAIHCwQD +BAsFAgULBQIJAwkCCQMEAwQLBQsECwUCCQIFAgkDBAMJAgcKBAoHCgQBCAEEAQgCBwIIAQYJBwkG +CQcACwAHAAsBBgELAAUIBggFCAYDCgMGAwoABQAK +""" + +#static const char tiling5[48][9] +TILING5 = (48, 9), """ +AggDAgoICgkIAQsCAQkLCQgLBAEJBAcBBwMBCAUECAMFAwEFAAoBAAgKCAsKCwQHCwIEAgAEBwAI +BwYABgIACQMACQUDBQcDAwYLAwAGAAQGAwkAAwsJCwoJBQIKBQQCBAACCQYFCQAGAAIGAAcIAAEH +AQUHCgABCgYABgQABgMLBgUDBQEDCgcGCgEHAQMHAQQJAQIEAgYECwECCwcBBwUBCAIDCAQCBAYC +AgUKAgMFAwcFBwoGBwgKCAkKBgkFBgsJCwgJBQgEBQoICgsIBAsHBAkLCQoLBAcLBAsJCQsKBQQI +BQgKCggLBgUJBgkLCwkIBwYKBwoICAoJAgoFAgUDAwUHCAMCCAIEBAIGCwIBCwEHBwEFAQkEAQQC +AgQGCgYHCgcBAQcDBgsDBgMFBQMBCgEACgAGBgAEAAgHAAcBAQcFCQUGCQYAAAYCBQoCBQIEBAIA +AwAJAwkLCwkKAwsGAwYAAAYECQADCQMFBQMHBwgABwAGBgACCwcECwQCAgQAAAEKAAoICAoLCAQF +CAUDAwUBBAkBBAEHBwEDAQILAQsJCQsIAgMIAggKCggJ +""" + +#static const char tiling6_1_1[48][9] +TILING6_1_1 = (48, 9), """ +BgUKAwEICQgBCwcGCQMBAwkIAQIKBwAEAAcDAwAIBQIGAgUBBQQJAgALCAsACgYFCAIAAggLCgYF +AAQDBwMEAwAIBgQKCQoECAMACgcFBwoLCAQHCgACAAoJBwYLAAIJCgkCAgMLBAEFAQQAAAEJBgMH +AwYCCQABCwQGBAsICwcGAQUABAAFAAEJBwULCgsFBAcIAQMKCwoDCQUECwEDAQsKCgECCAUHBQgJ +CAQHAgYBBQEGAQIKBAYICwgGAgMLBQcJCAkHCwIDCQYEBgkKCQUEAwcCBgIHBAUJAgcDBwIGAwIL +BAYJCgkGCwMCCQcFBwkICgIBCAYEBggLBwQIAQYCBgEFAgEKBwUICQgFBAUJAwELCgsBCAcECgMB +AwoLCQEACwUHBQsKBgcLAAUBBQAEAQAJBgQLCAsECQEABwMGAgYDCwMCBQEEAAQBCwYHCQIAAgkK +BwQIAgAKCQoAAAMIBQcKCwoHCAADCgQGBAoJBQYKAwQABAMHBQYKAAIICwgCCQQFCwACAAsICAAD +BgIFAQUCCgIBBAAHAwcABgcLAQMJCAkDCgUGCAEDAQgJ +""" + +#static const char tiling6_1_2[48][27] +TILING6_1_2 = (48, 27), """ +AQwDDAoDBgMKAwYIBQgGCAUMDAkIAQkMDAUKAQwDAQsMCwEGCQYBBgkHDAcJCQgMDAgDCwcMBAwA +BAEMAQQKBwoECgcCDAIHBwMMDAMAAQIMBgwCBgMMAwYIBQgGCAUADAAFBQEMDAECAwAMAAwCDAkC +BQIJAgULBAsFCwQMDAgLAAgMDAQJAAwCAAoMCgAFCAUABQgGDAYICAsMDAsCCgYMBAwADAUACgAF +AAoDBgMKAwYMDAcDBAcMDAYFBAwGDAgGAwYIBgMKAAoDCgAMDAkKBAkMDAAIBQwHBQgMCAUACgAF +AAoDDAMKCgsMDAsHCAMMAgwAAggMCAIHCgcCBwoEDAQKCgkMDAkACAQMAgwADAsABwALAAcJBgkH +CQYMDAoJAgoMDAYLBQwBBQIMAgULBAsFCwQDDAMEBAAMDAABAgMMBwwDBwAMAAcJBgkHCQYBDAEG +BgIMDAIDAAEMBgwEBgkMCQYBCwEGAQsADAALCwgMDAgECQAMBQwBDAYBCwEGAQsABwALAAcMDAQA +BQQMDAcGBQwHDAkHAAcJBwALAQsACwEMDAoLBQoMDAEJAwwBDAgBBAEIAQQKBwoECgcMDAsKAwsM +DAcIAwwBAwkMCQMECwQDBAsFDAULCwoMDAoBCQUMBwwFBwoMCgcCCAIHAggBDAEICAkMDAkFCgEM +BgwCDAcCCAIHAggBBAEIAQQMDAUBBgUMDAQHBgwEDAoEAQQKBAEIAggBCAIMDAsIBgsMDAIKBwwF +DAsFAgULBQIJAwkCCQMMDAgJBwgMDAMLBAwGBAsMCwQDCQMEAwkCDAIJCQoMDAoGCwIMBwwDDAQD +CQMEAwkCBQIJAgUMDAYCBwYMDAUEAwwHAwQMBAMJAgkDCQIFDAUCAgYMDAYHBAUMBgwEDAsEAwQL +BAMJAgkDCQIMDAoJBgoMDAILBQwHBQsMCwUCCQIFAgkDDAMJCQgMDAgHCwMMBAwGBAoMCgQBCAEE +AQgCDAIICAsMDAsGCgIMAgwGAgcMBwIIAQgCCAEEDAQBAQUMDAUGBwQMBQwHDAoHAgcKBwIIAQgC +CAEMDAkIBQkMDAEKAQwDDAkDBAMJAwQLBQsECwUMDAoLAQoMDAUJAQwDAQgMCAEECgQBBAoHDAcK +CgsMDAsDCAcMBwwFBwkMCQcACwAHAAsBDAELCwoMDAoFCQEMAQwFAQYMBgELAAsBCwAHDAcAAAQM +DAQFBgcMBAwGDAkGAQYJBgELAAsBCwAMDAgLBAgMDAAJAwwHDAAHCQcABwkGAQYJBgEMDAIGAwIM +DAEAAQwFDAIFCwUCBQsEAwQLBAMMDAAEAQAMDAMCAAwCAAsMCwAHCQcABwkGDAYJCQoMDAoCCwYM +AAwCDAgCBwIIAgcKBAoHCgQMDAkKAAkMDAQIBwwFDAgFAAUIBQAKAwoACgMMDAsKBwsMDAMIBgwE +BggMCAYDCgMGAwoADAAKCgkMDAkECAAMAAwEAAUMBQAKAwoACgMGDAYDAwcMDAcEBQYMAgwADAoA +BQAKAAUIBggFCAYMDAsIAgsMDAYKAgwAAgkMCQIFCwUCBQsEDAQLCwgMDAgACQQMAgwGDAMGCAYD +BggFAAUIBQAMDAEFAgEMDAADAAwEDAEECgQBBAoHAgcKBwIMDAMHAAMMDAIBAwwBDAsBBgELAQYJ +BwkGCQcMDAgJAwgMDAcLAwwBAwoMCgMGCAYDBggFDAUICAkMDAkBCgUM +""" + +#static const char tiling6_2[48][15] +TILING6_2 = (48, 15), """ +AQoDBgMKAwYIBQgGCAUJAQsDCwEGCQYBBgkHCAcJBAEAAQQKBwoECgcCAwIHBgMCAwYIBQgGCAUA +AQAFAAkCBQIJAgULBAsFCwQIAAoCCgAFCAUABQgGCwYIBAUACgAFAAoDBgMKAwYHBAgGAwYIBgMK +AAoDCgAJBQgHCAUACgAFAAoDCwMKAggACAIHCgcCBwoECQQKAgsABwALAAcJBgkHCQYKBQIBAgUL +BAsFCwQDAAMEBwADAAcJBgkHCQYBAgEGBgkECQYBCwEGAQsACAALBQYBCwEGAQsABwALAAcEBQkH +AAcJBwALAQsACwEKAwgBBAEIAQQKBwoECgcLAwkBCQMECwQDBAsFCgULBwoFCgcCCAIHAggBCQEI +BgcCCAIHAggBBAEIAQQFBgoEAQQKBAEIAggBCAILBwsFAgULBQIJAwkCCQMIBAsGCwQDCQMEAwkC +CgIJBwQDCQMEAwkCBQIJAgUGAwQHBAMJAgkDCQIFBgUCBgsEAwQLBAMJAgkDCQIKBQsHCwUCCQIF +AgkDCAMJBAoGCgQBCAEEAQgCCwIIAgcGBwIIAQgCCAEEBQQBBQoHAgcKBwIIAQgCCAEJAQkDBAMJ +AwQLBQsECwUKAQgDCAEECgQBBAoHCwcKBwkFCQcACwAHAAsBCgELAQYFBgELAAsBCwAHBAcABAkG +AQYJBgELAAsBCwAIAwAHCQcABwkGAQYJBgECAQIFCwUCBQsEAwQLBAMAAAsCCwAHCQcABwkGCgYJ +AAgCBwIIAgcKBAoHCgQJBwgFAAUIBQAKAwoACgMLBggECAYDCgMGAwoACQAKAAUEBQAKAwoACgMG +BwYDAgoABQAKAAUIBggFCAYLAgkACQIFCwUCBQsECAQLAgMGCAYDBggFAAUIBQABAAEECgQBBAoH +AgcKBwIDAwsBBgELAQYJBwkGCQcIAwoBCgMGCAYDBggFCQUI +""" + +#static const char tiling7_1[16][9] +TILING7_1 = (16, 9), """ +CQUECgECCAMACwcGCAMACgECAwAIBQQJBwYLCAQHCQABCwIDCgYFCwIDCQABAAEJBgUKBAcIAQIK +BwYLBQQJAgMLBAcIBgUKCwMCCAcECgUGCgIBCwYHCQQFCQEACgUGCAcEBQYKAwILAQAJBwQIAQAJ +AwILCAADCQQFCwYHBgcLAAMIAgEKBAUJAgEKAAMI +""" + +#static const char tiling7_2[16][3][15] +TILING7_2 = (16, 3, 15), """ +AQIKAwQIBAMFAAUDBQAJAwAICQEEAgQBBAIFCgUCCQUEAAoBCgAICggCAwIIAwAIAQYKBgEHAgcB +BwILAQIKCwMGAAYDBgAHCAcACwcGAggDCAIKCAoAAQAKCQUECwMGAAYDBgAHCAcACwcGAwQIBAMF +AAUDBQAJAwAIBAkHCwcJBQsJCwUGAAEJAgcLBwIEAwQCBAMIAgMLCAAHAQcABwEECQQBCAQHAwkA +CQMLCQsBAgELAgMLAAUJBQAGAQYABgEKAAEJCgIFAwUCBQMGCwYDBgUKAQsCCwEJCwkDAAMJBgUK +CAAHAQcABwEECQQBCAQHAAUJBQAGAQYABgEKAAEJBQoECAQKBggKCAYHCwcGCQEEAgQBBAIFCgUC +CQUEAQYKBgEHAgcBBwILAQIKBgsFCQULBwkLCQcECAQHCgIFAwUCBQMGCwYDBgUKAgcLBwIEAwQC +BAMIAgMLBwgGCgYIBAoICgQFBwQIBQIKAgUDBgMFAwYLCgUGCwcCBAIHAgQDCAMECwMCBggHCAYK +CAoEBQQKBgcLBAEJAQQCBQIEAgUKBAUJCgYBBwEGAQcCCwIHCgIBBQsGCwUJCwkHBAcJCgUGBwAI +AAcBBAEHAQQJBwQICQUABgAFAAYBCgEGCQEABAoFCgQICggGBwYICwMCCQUABgAFAAYBCgEGCQEA +BQIKAgUDBgMFAwYLCgUGAgsBCQELAwkLCQMACQEACwcCBAIHAgQDCAMECwMCBwAIAAcBBAEHAQQJ +BwQIAAkDCwMJAQsJCwECBAUJBgMLAwYABwAGAAcIBgcLCAQDBQMEAwUACQAFCAADBwkECQcLCQsF +BgULCAADCgYBBwEGAQcCCwIHCgIBBgMLAwYABwAGAAcIBgcLAwgCCgIIAAoICgABCgIBCAQDBQME +AwUACQAFCAADBAEJAQQCBQIEAgUKBAUJAQoACAAKAggKCAID +""" + +#static const char tiling7_3[16][3][27] +TILING7_3 = (16, 3, 27), """ +DAIKDAoFDAUEDAQIDAgDDAMADAAJDAkBDAECDAUEDAQIDAgDDAMCDAIKDAoBDAEADAAJDAkFBQQM +CgUMAgoMAwIMCAMMAAgMAQAMCQEMBAkMDAAIDAgHDAcGDAYKDAoBDAECDAILDAsDDAMADAcGDAYK +DAoBDAEADAAIDAgDDAMCDAILDAsHBwYMCAcMAAgMAQAMCgEMAgoMAwIMCwMMBgsMCQUMAAkMAwAM +CwMMBgsMBwYMCAcMBAgMBQQMAwAMCwMMBgsMBQYMCQUMBAkMBwQMCAcMAAgMDAMADAAJDAkFDAUG +DAYLDAsHDAcEDAQIDAgDDAEJDAkEDAQHDAcLDAsCDAIDDAMIDAgADAABDAQHDAcLDAsCDAIBDAEJ +DAkADAADDAMIDAgEBAcMCQQMAQkMAgEMCwIMAwsMAAMMCAAMBwgMDAMLDAsGDAYFDAUJDAkADAAB +DAEKDAoCDAIDDAYFDAUJDAkADAADDAMLDAsCDAIBDAEKDAoGBgUMCwYMAwsMAAMMCQAMAQkMAgEM +CgIMBQoMCgYMAQoMAAEMCAAMBwgMBAcMCQQMBQkMBgUMAAEMCAAMBwgMBgcMCgYMBQoMBAUMCQQM +AQkMDAABDAEKDAoGDAYHDAcIDAgEDAQFDAUJDAkACwcMAgsMAQIMCQEMBAkMBQQMCgUMBgoMBwYM +AQIMCQEMBAkMBwQMCwcMBgsMBQYMCgUMAgoMDAECDAILDAsHDAcEDAQJDAkFDAUGDAYKDAoBCAQM +AwgMAgMMCgIMBQoMBgUMCwYMBwsMBAcMAgMMCgIMBQoMBAUMCAQMBwgMBgcMCwYMAwsMDAIDDAMI +DAgEDAQFDAUKDAoGDAYHDAcLDAsCDAQIDAgDDAMCDAIKDAoFDAUGDAYLDAsHDAcEDAMCDAIKDAoF +DAUEDAQIDAgHDAcGDAYLDAsDAwIMCAMMBAgMBQQMCgUMBgoMBwYMCwcMAgsMDAcLDAsCDAIBDAEJ +DAkEDAQFDAUKDAoGDAYHDAIBDAEJDAkEDAQHDAcLDAsGDAYFDAUKDAoCAgEMCwIMBwsMBAcMCQQM +BQkMBgUMCgYMAQoMDAYKDAoBDAEADAAIDAgHDAcEDAQJDAkFDAUGDAEADAAIDAgHDAcGDAYKDAoF +DAUEDAQJDAkBAQAMCgEMBgoMBwYMCAcMBAgMBQQMCQUMAAkMCwMMBgsMBQYMCQUMAAkMAQAMCgEM +AgoMAwIMBQYMCQUMAAkMAwAMCwMMAgsMAQIMCgEMBgoMDAUGDAYLDAsDDAMADAAJDAkBDAECDAIK +DAoFCQEMBAkMBwQMCwcMAgsMAwIMCAMMAAgMAQAMBwQMCwcMAgsMAQIMCQEMAAkMAwAMCAMMBAgM +DAcEDAQJDAkBDAECDAILDAsDDAMADAAIDAgHDAUJDAkADAADDAMLDAsGDAYHDAcIDAgEDAQFDAAD +DAMLDAsGDAYFDAUJDAkEDAQHDAcIDAgAAAMMCQAMBQkMBgUMCwYMBwsMBAcMCAQMAwgMCAAMBwgM +BgcMCgYMAQoMAgEMCwIMAwsMAAMMBgcMCgYMAQoMAAEMCAAMAwgMAgMMCwIMBwsMDAYHDAcIDAgA +DAABDAEKDAoCDAIDDAMLDAsGCgIMBQoMBAUMCAQMAwgMAAMMCQAMAQkMAgEMBAUMCAQMAwgMAgMM +CgIMAQoMAAEMCQAMBQkMDAQFDAUKDAoCDAIDDAMIDAgADAABDAEJDAkE +""" + +#static const char tiling7_4_1[16][15] +TILING7_4_1 = (16, 15), """ +AwQIBAMKAgoDBAoFCQEAAQYKBgEIAAgBBggHCwMCCwMGCQYDBgkFAAkDBwQIAgcLBwIJAQkCBwkE +CAADAAUJBQALAwsABQsGCgIBCAAHCgcABwoGAQoABAUJCQEECwQBBAsHAgsBBQYKCgIFCAUCBQgE +AwgCBgcLBQIKAgUIBAgFAggDCwcGBAEJAQQLBwsEAQsCCgYFBwAIAAcKBgoHAAoBCQUECQUACwAF +AAsDBgsFAQIKCwcCCQIHAgkBBAkHAwAIBgMLAwYJBQkGAwkACAQHCgYBCAEGAQgABwgGAgMLCAQD +CgMEAwoCBQoEAAEJ +""" + +#static const char tiling7_4_2[16][27] +TILING7_4_2 = (16, 27), """ +CQQIBAkFCgUJAQoJCgECAAIBAgADCAMACQgACwYKBgsHCAcLAwgLCAMAAgADAAIBCgECCwoCCwMI +AAgDCAAJCAkEBQQJBAUHBgcFBwYLBwsICAcLBwgECQQIAAkICQABAwEAAQMCCwIDCAsDCgUJBQoG +CwYKAgsKCwIDAQMCAwEACQABCgkBCAAJAQkACQEKCQoFBgUKBQYEBwQGBAcIBAgJCQEKAgoBCgIL +CgsGBwYLBgcFBAUHBQQJBQkKCgILAwsCCwMICwgHBAcIBwQGBQYEBgUKBgoLCwIKAgsDCAMLBwgL +CAcEBgQHBAYFCgUGCwoGCgEJAQoCCwIKBgsKCwYHBQcGBwUECQQFCgkFCQAIAAkBCgEJBQoJCgUG +BAYFBgQHCAcECQgECQUKBgoFCgYLCgsCAwILAgMBAAEDAQAJAQkKCwcIBAgHCAQJCAkAAQAJAAED +AgMBAwILAwsICAMLAwgACQAIBAkICQQFBwUEBQcGCwYHCAsHCgYLBwsGCwcICwgDAAMIAwACAQIA +AgEKAgoLCAQJBQkECQUKCQoBAgEKAQIAAwACAAMIAAgJ +""" + +#static const char tiling8[6][6] +TILING8 = (6, 6), """ +CQgKCggLAQUDAwUHAAQCBAYCAAIEBAIGAQMFAwcFCQoICgsI +""" + +#static const char tiling9[8][12] +TILING9 = (8, 12), """ +AgoFAwIFAwUEAwQIBAcLCQQLCQsCCQIBCgcGAQcKAQgHAQAIAwYLAAYDAAUGAAkFAwsGAAMGAAYF +AAUJCgYHAQoHAQcIAQgABAsHCQsECQILCQECAgUKAwUCAwQFAwgE +""" + +#static const char tiling10_1_1[6][12] +TILING10_1_1 = (6, 12), """ +BQoHCwcKCAEJAQgDAQIFBgUCBAMAAwQHCwAIAAsCBAkGCgYJCQAKAgoABggECAYLBwIDAgcGAAEE +BQQBBwkFCQcICgELAwsB +""" + +#static const char tiling10_1_1_[6][12] +TILING10_1_1_ = (6, 12), """ +BQkHCAcJCwEKAQsDAwIHBgcCBAEAAQQFCgAJAAoCBAgGCwYICAALAgsABgkECQYKBQIBAgUGAAME +BwQDBwoFCgcLCQEIAwgB +""" + +#static const char tiling10_1_2[6][24] +TILING10_1_2 = (6, 24), """ +AwsHAwcICQgHBQkHCQUKCQoBAwEKCwMKBwYFBwUEAAQFAQAFAAECAAIDBwMCBgcCCwIKBgsKCwYE +CwQIAAgECQAEAAkKAAoCCwIKCwoGBAYKCQQKBAkABAAICwgAAgsABwYFBAcFBwQABwADAgMAAQIA +AgEFAgUGBwgDCwcDBwsKBwoFCQUKAQkKCQEDCQMI +""" + +#static const char tiling10_2[6][24] +TILING10_2 = (6, 24), """ +DAUJDAkIDAgDDAMBDAEKDAoLDAsHDAcFDAEADAAEDAQHDAcDDAMCDAIGDAYFDAUBBAgMBgQMCgYM +CQoMAAkMAgAMCwIMCAsMDAkEDAQGDAYLDAsIDAgADAACDAIKDAoJAAMMBAAMBQQMAQUMAgEMBgIM +BwYMAwcMCgUMCwoMAwsMAQMMCQEMCAkMBwgMBQcM +""" + +#static const char tiling10_2_[6][24] +TILING10_2_ = (6, 24), """ +CAcMCQgMAQkMAwEMCwMMCgsMBQoMBwUMBAUMAAQMAwAMBwMMBgcMAgYMAQIMBQEMDAsGDAYEDAQJ +DAkKDAoCDAIADAAIDAgLBgoMBAYMCAQMCwgMAgsMAAIMCQAMCgkMDAcEDAQADAABDAEFDAUGDAYC +DAIDDAMHDAcLDAsKDAoBDAEDDAMIDAgJDAkFDAUH +""" + +#static const char tiling11[12][12] +TILING11 = (12, 12), """ +AgoJAgkHAgcDBwkEAQYCAQgGAQkICAcGCAMBCAEGCAYEBgEKAAgLAAsFAAUBBQsGCQUHCQcCCQIA +AgcLBQAEBQsABQoLCwMABQQABQALBQsKCwADCQcFCQIHCQACAgsHAAsIAAULAAEFBQYLCAEDCAYB +CAQGBgoBAQIGAQYIAQgJCAYHAgkKAgcJAgMHBwQJ +""" + +#static const char tiling12_1_1[24][12] +TILING12_1_1 = (24, 12), """ +BwYLCgMCAwoICQgKBgUKCQIBAgkLCAsJCgYFBwkECQcBAwEHBwYLBAgFAwUIBQMBBQQJCAEAAQgK +CwoIAQIKAAkDBQMJAwUHCgECAAsDCwAGBAYACAMAAgkBCQIEBgQCAwAIAgsBBwELAQcFBgUKBwsE +AgQLBAIACQUEBggHCAYAAgAGCAMABwQLCQsECwkKBAcICwADAAsJCgkLBAcIBQkGAAYJBgACCwcG +BAoFCgQCAAIECwIDAQgACAEHBQcBAAEJAwgCBAIIAgQGAgMLAQoABgAKAAYECQABAwoCCgMFBwUD +CQABBAUICggFCAoLCAQHBQsGCwUDAQMFBQQJBgoHAQcKBwEDCgECBQYJCwkGCQsICwIDBgcKCAoH +CggJ +""" + +#static const char tiling12_1_1_[24][12] +TILING12_1_1_ = (24, 12), """ +AwILCgcGBwoICQgKAgEKCQYFBgkLCAsJCQQFBwoGCgcBAwEHBwQIBgsFAwULBQMBAQAJCAUEBQgK +CwoIAQAJAgoDBQMKAwUHCwMCAAoBCgAGBAYACQEAAggDCAIEBgQCAwILAAgBBwEIAQcFBgcLBQoE +AgQKBAIACAcEBgkFCQYAAgAGCAcEAwALCQsACwkKAAMICwQHBAsJCgkLBAUJBwgGAAYIBgACCgUG +BAsHCwQCAAIECAADAQsCCwEHBQcBAAMIAQkCBAIJAgQGAgEKAwsABgALAAYECgIBAwkACQMFBwUD +CQQFAAEICggBCAoLCwYHBQgECAUDAQMFBQYKBAkHAQcJBwEDCgUGAQIJCwkCCQsICwYHAgMKCAoD +CggJ +""" + +#static const char tiling12_1_2[24][24] +TILING12_1_2 = (24, 24), """ +BwMLAwcICQgHBgkHCQYKAgoGCwIGAgsDBgIKAgYLCAsGBQgGCAUJAQkFCgEFAQoCCgkFCQoBAwEK +BgMKAwYHBAcGBQQGBAUJBwgLAwsICwMBCwEGBQYBBgUEBgQHCAcEBQEJAQUKCwoFBAsFCwQIAAgE +CQAEAAkBAQkKBQoJCgUHCgcCAwIHAgMAAgABCQEACgsCCwoGBAYKAQQKBAEAAwABAgMBAwILCAkA +CQgEBgQIAwYIBgMCAQIDAAEDAQAJAwsIBwgLCAcFCAUAAQAFAAECAAIDCwMCBgsKAgoLCgIACgAF +BAUABQQHBQcGCwYHCQgECAkAAgAJBQIJAgUGBwYFBAcFBwQICAQACQAEAAkKAAoDCwMKAwsHAwcI +BAgHBAAIAAQJCgkEBwoECgcLAwsHCAMHAwgABAkIAAgJCAACCAIHBgcCBwYFBwUECQQFCwoGCgsC +AAILBwALAAcEBQQHBgUHBQYKCwgDCAsHBQcLAgULBQIBAAECAwACAAMIAAgJBAkICQQGCQYBAgEG +AQIDAQMACAADAgoLBgsKCwYECwQDAAMEAwABAwECCgIBCQoBCgkFBwUJAAcJBwADAgMAAQIAAgEK +CQUBCgEFAQoLAQsACAALAAgEAAQJBQkECAsHCwgDAQMIBAEIAQQFBgUEBwYEBgcLBQoJAQkKCQED +CQMEBwQDBAcGBAYFCgUGCgYCCwIGAgsIAggBCQEIAQkFAQUKBgoFCwcDCAMHAwgJAwkCCgIJAgoG +AgYLBwsG +""" + +#static const char tiling12_2[24][24] +TILING12_2 = (24, 24), """ +CQgMCgkMAgoMAwIMCwMMBgsMBwYMCAcMCAsMCQgMAQkMAgEMCgIMBQoMBgUMCwYMAwEMBwMMBAcM +CQQMBQkMBgUMCgYMAQoMDAMBDAEFDAUGDAYLDAsHDAcEDAQIDAgDCwoMCAsMAAgMAQAMCQEMBAkM +BQQMCgUMDAUHDAcDDAMCDAIKDAoBDAEADAAJDAkFBAYMAAQMAQAMCgEMAgoMAwIMCwMMBgsMBgQM +AgYMAwIMCAMMAAgMAQAMCQEMBAkMDAcFDAUBDAEADAAIDAgDDAMCDAILDAsHDAIADAAEDAQFDAUK +DAoGDAYHDAcLDAsCAgAMBgIMBwYMCAcMBAgMBQQMCQUMAAkMDAkKDAoLDAsHDAcEDAQIDAgDDAMA +DAAJCgkMCwoMBwsMBAcMCAQMAwgMAAMMCQAMDAACDAIGDAYHDAcIDAgEDAQFDAUJDAkAAAIMBAAM +BQQMCgUMBgoMBwYMCwcMAgsMBQcMAQUMAAEMCAAMAwgMAgMMCwIMBwsMDAQGDAYCDAIDDAMIDAgA +DAABDAEJDAkEDAYEDAQADAABDAEKDAoCDAIDDAMLDAsGBwUMAwcMAgMMCgIMAQoMAAEMCQAMBQkM +DAoLDAsIDAgADAABDAEJDAkEDAQFDAUKAQMMBQEMBgUMCwYMBwsMBAcMCAQMAwgMDAEDDAMHDAcE +DAQJDAkFDAUGDAYKDAoBDAsIDAgJDAkBDAECDAIKDAoFDAUGDAYLDAgJDAkKDAoCDAIDDAMLDAsG +DAYHDAcI +""" + +#static const char tiling12_2_[24][24] +TILING12_2_ = (24, 24), """ +DAILDAsHDAcGDAYKDAoJDAkIDAgDDAMCDAEKDAoGDAYFDAUJDAkIDAgLDAsCDAIBDAQFDAUKDAoG +DAYHDAcDDAMBDAEJDAkEBwYMCAcMBAgMBQQMAQUMAwEMCwMMBgsMDAAJDAkFDAUEDAQIDAgLDAsK +DAoBDAEAAQIMCQEMAAkMAwAMBwMMBQcMCgUMAgoMDAECDAILDAsDDAMADAAEDAQGDAYKDAoBDAMA +DAAJDAkBDAECDAIGDAYEDAQIDAgDAwAMCwMMAgsMAQIMBQEMBwUMCAcMAAgMBgUMCwYMBwsMBAcM +AAQMAgAMCgIMBQoMDAcEDAQJDAkFDAUGDAYCDAIADAAIDAgHCAcMAAgMAwAMCwMMCgsMCQoMBAkM +BwQMDAcIDAgADAADDAMLDAsKDAoJDAkEDAQHBAcMCQQMBQkMBgUMAgYMAAIMCAAMBwgMDAUGDAYL +DAsHDAcEDAQADAACDAIKDAoFDAADDAMLDAsCDAIBDAEFDAUHDAcIDAgAAAMMCQAMAQkMAgEMBgIM +BAYMCAQMAwgMAgEMCwIMAwsMAAMMBAAMBgQMCgYMAQoMDAIBDAEJDAkADAADDAMHDAcFDAUKDAoC +CQAMBQkMBAUMCAQMCwgMCgsMAQoMAAEMDAYHDAcIDAgEDAQFDAUBDAEDDAMLDAsGBQQMCgUMBgoM +BwYMAwcMAQMMCQEMBAkMCgEMBgoMBQYMCQUMCAkMCwgMAgsMAQIMCwIMBwsMBgcMCgYMCQoMCAkM +AwgMAgMM +""" + +#static const char tiling13_1[2][12] +TILING13_1 = (2, 12), """ +CwcGAQIKCAMACQUECAQHAgMLCQABCgYF +""" + +#static const char tiling13_1_[2][12] +TILING13_1_ = (2, 12), """ +BwQICwMCAQAJBQYKBgcLCgIBAAMIBAUJ +""" + +#static const char tiling13_2[2][6][18] +TILING13_2 = (2, 6, 18), """ +AQIKCwcGAwQIBAMFAAUDBQAJCAMACwcGCQEEAgQBBAIFCgUCCQUECAMAAQYKBgEHAgcBBwILCQUE +AQIKCwMGAAYDBgAHCAcACQUECwcGAAoBCgAICggCAwIIAQIKAwAIBAkHCwcJBQsJCwUGAgMLCAQH +AAUJBQAGAQYABgEKCQABCAQHCgIFAwUCBQMGCwYDBgUKCQABAgcLBwIEAwQCBAMIBgUKAgMLCAAH +AQcABwEECQQBBgUKCAQHAQsCCwEJCwkDAAMJAgMLAAEJBQoECAQKBggKCAYH +""" + +#static const char tiling13_2_[2][6][18] +TILING13_2_ = (2, 6, 18), """ +CgUGCwMCBwAIAAcBBAEHAQQJCwMCBwQICQUABgAFAAYBCgEGAQAJBwQIBQIKAgUDBgMFAwYLCgUG +AQAJCwcCBAIHAgQDCAMECgUGBwQIAgsBCQELAwkLCQMACwMCCQEABAoFCgQICggGBwYIBgcLCAAD +BAEJAQQCBQIEAgUKCAADBAUJCgYBBwEGAQcCCwIHAgEKBAUJBgMLAwYABwAGAAcIBgcLAgEKCAQD +BQMEAwUACQAFBgcLBAUJAwgCCgIIAAoICgABCAADCgIBBQsGCwUJCwkHBAcJ +""" + +#static const char tiling13_3[2][12][30] +TILING13_3 = (2, 12, 30), """ +CwcGDAIKDAoFDAUEDAQIDAgDDAMADAAJDAkBDAECAQIKCQUMAAkMAwAMCwMMBgsMBwYMCAcMBAgM +BQQMCwcGDAUEDAQIDAgDDAMCDAIKDAoBDAEADAAJDAkFAQIKDAMADAAJDAkFDAUGDAYLDAsHDAcE +DAQIDAgDCAMACwcMAgsMAQIMCQEMBAkMBQQMCgUMBgoMBwYMCwcGBQQMCgUMAgoMAwIMCAMMAAgM +AQAMCQEMBAkMCAMAAQIMCQEMBAkMBwQMCwcMBgsMBQYMCgUMAgoMCQUEDAAIDAgHDAcGDAYKDAoB +DAECDAILDAsDDAMACQUEDAcGDAYKDAoBDAEADAAIDAgDDAMCDAILDAsHCAMADAECDAILDAsHDAcE +DAQJDAkFDAUGDAYKDAoBCQUEBwYMCAcMAAgMAQAMCgEMAgoMAwIMCwMMBgsMAQIKAwAMCwMMBgsM +BQYMCQUMBAkMBwQMCAcMAAgMCAQHDAMLDAsGDAYFDAUJDAkADAABDAEKDAoCDAIDAgMLCgYMAQoM +AAEMCAAMBwgMBAcMCQQMBQkMBgUMCAQHDAYFDAUJDAkADAADDAMLDAsCDAIBDAEKDAoGAgMLDAAB +DAEKDAoGDAYHDAcIDAgEDAQFDAUJDAkAAAEJCAQMAwgMAgMMCgIMBQoMBgUMCwYMBwsMBAcMCAQH +BgUMCwYMAwsMAAMMCQAMAQkMAgEMCgIMBQoMCQABAgMMCgIMBQoMBAUMCAQMBwgMBgcMCwYMAwsM +BgUKDAEJDAkEDAQHDAcLDAsCDAIDDAMIDAgADAABBgUKDAQHDAcLDAsCDAIBDAEJDAkADAADDAMI +DAgECQABDAIDDAMIDAgEDAQFDAUKDAoGDAYHDAcLDAsCBgUKBAcMCQQMAQkMAgEMCwIMAwsMAAMM +CAAMBwgMAgMLAAEMCAAMBwgMBgcMCgYMBQoMBAUMCQQMAQkM +""" + +#static const char tiling13_3_[2][12][30] +TILING13_3_ = (2, 12, 30), """ +AwILCAcMAAgMAQAMCgEMBgoMBQYMCQUMBAkMBwQMBQYKDAILDAsHDAcEDAQJDAkBDAEADAAIDAgD +DAMCCgUGDAcEDAQJDAkBDAECDAILDAsDDAMADAAIDAgHCwMCDAEADAAIDAgHDAcGDAYKDAoFDAUE +DAQJDAkBBwQICwMMBgsMBQYMCQUMAAkMAQAMCgEMAgoMAwIMBwQIBQYMCQUMAAkMAwAMCwMMAgsM +AQIMCgEMBgoMCwMCAQAMCgEMBgoMBwYMCAcMBAgMBQQMCQUMAAkMAQAJDAQIDAgDDAMCDAIKDAoF +DAUGDAYLDAsHDAcEBwQIDAUGDAYLDAsDDAMADAAJDAkBDAECDAIKDAoFAQAJDAMCDAIKDAoFDAUE +DAQIDAgHDAcGDAYLDAsDCgUGBwQMCwcMAgsMAQIMCQEMAAkMAwAMCAMMBAgMCQEAAwIMCAMMBAgM +BQQMCgUMBgoMBwYMCwcMAgsMAAMICQQMAQkMAgEMCwIMBwsMBgcMCgYMBQoMBAUMCwYHDAMIDAgE +DAQFDAUKDAoCDAIBDAEJDAkADAADBgcLDAQFDAUKDAoCDAIDDAMIDAgADAABDAEJDAkECAADDAIB +DAEJDAkEDAQHDAcLDAsGDAYFDAUKDAoCBAUJCAAMBwgMBgcMCgYMAQoMAgEMCwIMAwsMAAMMBAUJ +BgcMCgYMAQoMAAEMCAAMAwgMAgMMCwIMBwsMCAADAgEMCwIMBwsMBAcMCQQMBQkMBgUMCgYMAQoM +AgEKDAUJDAkADAADDAMLDAsGDAYHDAcIDAgEDAQFBAUJDAYHDAcIDAgADAABDAEKDAoCDAIDDAML +DAsGAgEKDAADDAMLDAsGDAYFDAUJDAkEDAQHDAcIDAgABgcLBAUMCAQMAwgMAgMMCgIMAQoMAAEM +CQAMBQkMCgIBAAMMCQAMBQkMBgUMCwYMBwsMBAcMCAQMAwgM +""" + +#static const char tiling13_4[2][4][36] +TILING13_4 = (2, 4, 36), """ +DAIKDAoFDAUGDAYLDAsHDAcEDAQIDAgDDAMADAAJDAkBDAECCwMMBgsMBwYMCAcMBAgMBQQMCQUM +AAkMAQAMCgEMAgoMAwIMCQEMBAkMBQQMCgUMBgoMBwYMCwcMAgsMAwIMCAMMAAgMAQAMDAAIDAgH +DAcEDAQJDAkFDAUGDAYKDAoBDAECDAILDAsDDAMADAMLDAsGDAYHDAcIDAgEDAQFDAUJDAkADAAB +DAEKDAoCDAIDCAAMBwgMBAcMCQQMBQkMBgUMCgYMAQoMAgEMCwIMAwsMAAMMCgIMBQoMBgUMCwYM +BwsMBAcMCAQMAwgMAAMMCQAMAQkMAgEMDAEJDAkEDAQFDAUKDAoGDAYHDAcLDAsCDAIDDAMIDAgA +DAAB +""" + +#static const char tiling13_5_1[2][4][18] +TILING13_5_1 = (2, 4, 18), """ +BwYLAQAJCgMCAwoFAwUIBAgFAQIKBwQIAwALBgsACQYABgkFAwAIBQYKAQIJBAkCCwQCBAsHBQQJ +AwILCAEAAQgHAQcKBgoHBAcIAgEKCwADAAsGAAYJBQkGAgMLBAUJAAEIBwgBCgcBBwoGAAEJBgcL +AgMKBQoDCAUDBQgEBgUKAAMICQIBAgkEAgQLBwsE +""" + +#static const char tiling13_5_2[2][4][30] +TILING13_5_2 = (2, 4, 30), """ +AQAJBwQIBwgDBwMLAgsDCwIKCwoGBQYKBgUHBAcFBwQICwMCBgsCCgYCBgoFCQUKAQkKCQEAAgAB +AAIDBQYKCQEABAkACAQABAgHCwcIAwsICwMCAAIDAgABAwILBQYKBQoBBQEJAAkBCQAICQgEBAgH +BAcFBgUHAgEKBAUJBAkABAAIAwgACAMLCAsHBgcLBwYEBQQGBAUJCAADBwgDCwcDBwsGCgYLAgoL +CgIBAwECAQMABgcLCgIBBQoBCQUBBQkECAQJAAgJCAADAQMAAwECAAMIBgcLBgsCBgIKAQoCCgEJ +CgkFBQkEBQQGBwYE +""" + +#static const char tiling14[12][12] +TILING14 = (12, 12), """ +BQkIBQgCBQIGAwIIAgEFAgUIAggLBAgFCQQGCQYDCQMBCwMGAQsKAQQLAQAEBwsECAIACAUCCAcF +CgIFAAcDAAoHAAkKBgcKAAMHAAcKAAoJBgoHCAACCAIFCAUHCgUCAQoLAQsEAQQABwQLCQYECQMG +CQEDCwYDAgUBAggFAgsIBAUIBQgJBQIIBQYCAwgC +""" + +#static const char test3[24] +TEST3 = (24,), """ +BQEEBQECAgMEAwYG+vr9/P3+/v/7/P/7 +""" + +#static const char test4[8] +TEST4 = (8,), """ +BwcHB/n5+fk= +""" + +#static const char test6[48][3] +TEST6 = (48, 3), """ +AgcKBAcLBQcBBQcDAQcJAwcKBgcFAQcIBAcIAQcIAwcLBQcCBQcAAQcJBgcGAgcJBAcIAgcJAgcK +BgcHAwcKBAcLAwcLBgcE+vkE/fkL/PkL/fkK+vkH/vkK/vkJ/PkI/vkJ+vkG//kJ+/kA+/kC/fkL +//kI/PkI//kI+vkF/fkK//kJ+/kD+/kB/PkL/vkK +""" + +#static const char test7[16][5] +TEST7 = (16, 5), """ +AQIFBwEDBAUHAwQBBgcEBAEFBwACAwUHAgECBgcFAgMGBwYDBAYHB/38+vkH/v36+Qb//vr5Bf79 ++/kC/P/7+QD8//r5BP38+/kD//77+QE= +""" + +#static const char test10[6][3] +TEST10 = (6, 3), """ +AgQHBQYHAQMHAQMHBQYHAgQH +""" + +#static const char test12[24][4] +TEST12 = (24, 4), """ +BAMHCwMCBwoCBgcFBgQHBwIBBwkFAgcBBQMHAgUBBwAFBAcDBgMHBgEGBwQBBAcIBAEHCAYBBwQD +BgcGBAUHAwEFBwADBQcCAgUHAQECBwkEBgcHBgIHBQIDBwoDBAcL +""" + +#static const char test13[2][7] +TEST13 = (2, 7), """ +AQIDBAUGBwIDBAEFBgc= +""" + +#static const char subconfig13[64] +SUBCONFIG13 = (64,), """ +AAECBwP/C/8ECP//Dv///wUJDBcP/xUmERT/JBohHiwGCg0TEP8ZJRIY/yMWIB0r////Iv//HCr/ +H/8pGygnLQ== +""" diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index ef92450a..cea657c5 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None): cython(['_find_contours_cy.pyx'], working_path=base_path) cython(['_moments_cy.pyx'], working_path=base_path) cython(['_marching_cubes_cy.pyx'], working_path=base_path) + cython(['_marching_cubes_lewiner_cy.pyx'], working_path=base_path) cython(['_pnpoly.pyx'], working_path=base_path) config.add_extension('_ccomp', sources=['_ccomp.c'], @@ -27,6 +28,9 @@ def configuration(parent_package='', top_path=None): config.add_extension('_marching_cubes_cy', sources=['_marching_cubes_cy.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_marching_cubes_lewiner_cy', + sources=['_marching_cubes_lewiner_cy.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('_pnpoly', sources=['_pnpoly.c'], include_dirs=[get_numpy_include_dirs(), '../_shared'])