_bitree_marching_cubes_cy.pyx

import numpy as np
cimport numpy as cnp
from stl import mesh

# Enable low level memory management
from libc.stdlib cimport malloc, free

# Define tiny winy number
cdef cnp.float64_t FLT_EPSILON = np.spacing(1.0)

'''
Define main volume type: 
float32 => 2^24 => small size)
float64 => 2^1023 => used for high-resolution volumetric data)
'''

cdef int DirectionX = 0
cdef int DirectionY = 1
cdef int DirectionZ = 2


cdef int[:,:] EDGE_DELTA = np.array([
    [0, 0, 0],[1, 0, 0],[0, 1, 0],[0, 0, 0],[0, 0, 1],[1, 0, 1],
    [0, 1, 1],[0, 0, 1],[0, 0, 0],[1, 0, 0],[1, 1, 0],[0, 1, 0],
])

cdef int[:] EDGE_DIRECTION = np.array([0,1,0,1,0,1,0,1,2,2,2,2])

cdef int[:,:] GEOMETRY_LOOKUP = np.array([
        [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 0 = 00000000
        [0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 1 = 00000001
        [0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 2 = 00000010
        [8, 3, 9, 1, 3, 9, -1, -1, -1, -1, -1, -1],  # 3 = 00000011
        [1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 4 = 00000100
        [0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1],  # 5 = 00000101
        [9, 10, 0, 2, 10, 0, -1, -1, -1, -1, -1, -1],  # 6 = 00000110
        [8, 3, 10, 2, 3, 10, 10, 9, 8, -1, -1, -1],  # 7 = 00000111
        [3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 8 = 00001000
        [0, 2, 8, 11, 2, 8, -1, -1, -1, -1, -1, -1],  # 9 = 00001001
        [1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1],  # 10 = 00001010
        [1, 2, 9, 11, 2, 9, 9, 8, 11, -1, -1, -1],  # 11 = 00001011
        [10, 1, 11, 3, 1, 11, -1, -1, -1, -1, -1, -1],  # 12 = 00001100
        [0, 1, 8, 10, 1, 8, 8, 11, 10, -1, -1, -1],  # 13 = 00001101
        [9, 0, 11, 3, 0, 11, 11, 10, 9, -1, -1, -1],  # 14 = 00001110
        [8, 9, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1],  # 15 = 00001111
        [4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 16 = 00010000
        [3, 0, 7, 4, 0, 7, -1, -1, -1, -1, -1, -1],  # 17 = 00010001
        [0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1],  # 18 = 00010010
        [1, 9, 7, 4, 9, 7, 7, 3, 1, -1, -1, -1],  # 19 = 00010011
        [1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1],  # 20 = 00010100
        [1, 2, 10, 0, 4, 3, 7, 4, 3, -1, -1, -1],  # 21 = 00010101
        [8, 4, 7, 0, 2, 9, 10, 2, 9, -1, -1, -1],  # 22 = 00010110
        [3, 2, 10, 3, 10, 4, 3, 4, 7, 4, 10, 9],  # 23 = 00010111
        [8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1],  # 24 = 00011000
        [11, 7, 2, 4, 7, 2, 2, 0, 4, -1, -1, -1],  # 25 = 00011001
        [8, 4, 7, 9, 0, 1, 11, 2, 3, -1, -1, -1],  # 26 = 00011010
        [4, 7, 11, 9, 4, 2, 11, 4, 2, 9, 2, 1],  # 27 = 00011011
        [4, 7, 8, 1, 3, 10, 11, 10, 3, -1, -1, -1],  # 28 = 00011100
        [10, 7, 11, 10, 0, 7, 10, 1, 0, 4, 7, 0],  # 29 = 00011101
        [4, 7, 8, 11, 0, 3, 0, 11, 9, 10, 9, 11],  # 30 = 00011110
        [11, 7, 9, 4, 7, 9, 9, 11, 10, -1, -1, -1],  # 31 = 00011111
        [9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 32 = 00100000
        [9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1],  # 33 = 00100001
        [0, 4, 1, 5, 4, 1, -1, -1, -1, -1, -1, -1],  # 34 = 00100010
        [8, 4, 3, 5, 4, 3, 3, 1, 5, -1, -1, -1],  # 35 = 00100011
        [1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1],  # 36 = 00100100
        [9, 5, 4, 10, 1, 2, 8, 3, 0, -1, -1, -1],  # 37 = 00100101
        [2, 10, 4, 5, 10, 4, 4, 0, 2, -1, -1, -1],  # 38 = 00100110
        [3, 2, 10, 8, 3, 10, 8, 10, 5, 8, 5, 4],  # 39 = 00100111
        [9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1],  # 40 = 00101000
        [5, 4, 9, 0, 2, 8, 11, 2, 8, -1, -1, -1],  # 41 = 00101001
        [2, 3, 11, 4, 1, 5, 1, 4, 0, -1, -1, -1],  # 42 = 00101010
        [11, 2, 1, 11, 1, 4, 11, 4, 8, 5, 4, 1],  # 43 = 00101011
        [9, 5, 4, 11, 1, 3, 1, 11, 10, -1, -1, -1],  # 44 = 00101100
        [5, 4, 9, 8, 0, 10, 1, 0, 10, 11, 10, 8],  # 45 = 00101101
        [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3],  # 46 = 00101110
        [8, 4, 10, 5, 4, 10, 10, 8, 11, -1, -1, -1],  # 47 = 00101111
        [9, 8, 5, 7, 8, 5, -1, -1, -1, -1, -1, -1],  # 48 = 00110000
        [9, 0, 5, 3, 0, 5, 5, 7, 3, -1, -1, -1],  # 49 = 00110001
        [0, 8, 1, 7, 8, 1, 1, 5, 7, -1, -1, -1],  # 50 = 00110010
        [1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1],  # 51 = 00110011
        [10, 1, 2, 8, 5, 7, 5, 8, 9, -1, -1, -1],  # 52 = 00110100
        [1, 2, 10, 9, 0, 5, 3, 0, 5, 3, 5, 7],  # 53 = 00110101
        [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2],  # 54 = 00110110
        [2, 10, 3, 5, 10, 3, 3, 5, 7, -1, -1, -1],  # 55 = 00110111
        [2, 3, 11, 5, 7, 9, 8, 9, 7, -1, -1, -1],  # 56 = 00111000
        [0, 9, 5, 0, 5, 11, 0, 11, 2, 11, 5, 7],  # 57 = 00111001
        [11, 2, 3, 8, 0, 7, 1, 0, 7, 5, 7, 1],  # 58 = 00111010
        [1, 2, 7, 11, 2, 7, 7, 1, 5, -1, -1, -1],  # 59 = 00111011
        [7, 9, 5, 9, 7, 8, 10, 1, 11, 3, 11, 1],  # 60 = 00111100
        [9, 1, 0, 11, 7, 10, 5, 7, 10, -1, -1, -1],  # 61 = 00111101
        [0, 3, 8, 10, 5, 11, 7, 5, 11, -1, -1, -1],  # 62 = 00111110
        [11, 10, 7, 5, 10, 7, -1, -1, -1, -1, -1, -1],  # 63 = 00111111
        [10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 64 = 01000000
        [0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1],  # 65 = 01000001
        [9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1],  # 66 = 01000010
        [6, 5, 10, 8, 3, 9, 1, 3, 9, -1, -1, -1],  # 67 = 01000011
        [1, 5, 2, 6, 5, 2, -1, -1, -1, -1, -1, -1],  # 68 = 01000100
        [3, 0, 8, 2, 6, 1, 5, 6, 1, -1, -1, -1],  # 69 = 01000101
        [9, 5, 0, 6, 5, 0, 0, 2, 6, -1, -1, -1],  # 70 = 01000110
        [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8],  # 71 = 01000111
        [2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1],  # 72 = 01001000
        [10, 6, 5, 8, 0, 11, 2, 0, 11, -1, -1, -1],  # 73 = 01001001
        [10, 6, 5, 11, 2, 3, 9, 0, 1, -1, -1, -1],  # 74 = 01001010
        [6, 5, 10, 9, 2, 1, 2, 9, 11, 8, 11, 9],  # 75 = 01001011
        [3, 11, 5, 6, 11, 5, 5, 1, 3, -1, -1, -1],  # 76 = 01001100
        [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6],  # 77 = 01001101
        [11, 6, 5, 3, 11, 5, 3, 5, 9, 3, 9, 0],  # 78 = 01001110
        [9, 5, 11, 6, 5, 11, 11, 9, 8, -1, -1, -1],  # 79 = 01001111
        [5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1],  # 80 = 01010000
        [10, 6, 5, 3, 0, 7, 4, 0, 7, -1, -1, -1],  # 81 = 01010001
        [0, 1, 9, 6, 5, 10, 4, 7, 8, -1, -1, -1],  # 82 = 01010010
        [10, 6, 5, 9, 4, 1, 7, 4, 1, 3, 1, 7],  # 83 = 01010011
        [8, 4, 7, 1, 2, 5, 6, 2, 5, -1, -1, -1],  # 84 = 01010100
        [2, 1, 6, 5, 1, 6, 3, 0, 7, 4, 0, 7],  # 85 = 01010101
        [4, 7, 8, 9, 5, 0, 6, 5, 0, 6, 0, 2],  # 86 = 01010110
        [4, 5, 9, 2, 7, 3, 7, 2, 6, -1, -1, -1],  # 87 = 01010111
        [2, 3, 11, 4, 7, 8, 6, 5, 10, -1, -1, -1],  # 88 = 01011000
        [6, 5, 10, 7, 11, 4, 2, 4, 11, 4, 2, 0],  # 89 = 01011001
        [8, 4, 7, 2, 3, 11, 9, 0, 1, 10, 6, 5],  # 90 = 01011010
        [10, 2, 1, 11, 6, 7, 9, 4, 5, -1, -1, -1],  # 91 = 01011011
        [8, 4, 7, 11, 6, 3, 5, 6, 3, 1, 3, 5],  # 92 = 01011100
        [6, 7, 11, 0, 5, 1, 5, 0, 4, -1, -1, -1],  # 93 = 01011101
        [8, 0, 3, 9, 4, 5, 11, 6, 7, -1, -1, -1],  # 94 = 01011110
        [5, 9, 4, 11, 6, 7, -1, -1, -1, -1, -1, -1],  # 95 = 01011111
        [10, 9, 6, 4, 9, 6, -1, -1, -1, -1, -1, -1],  # 96 = 01100000
        [3, 0, 8, 10, 6, 9, 4, 6, 9, -1, -1, -1],  # 97 = 01100001
        [0, 1, 6, 10, 1, 6, 6, 4, 0, -1, -1, -1],  # 98 = 01100010
        [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10],  # 99 = 01100011
        [1, 9, 2, 4, 9, 2, 2, 6, 4, -1, -1, -1],  # 100 = 01100100
        [8, 3, 0, 2, 9, 1, 9, 2, 4, 6, 4, 2],  # 101 = 01100101
        [2, 0, 6, 4, 0, 6, -1, -1, -1, -1, -1, -1],  # 102 = 01100110
        [8, 3, 4, 2, 3, 4, 4, 2, 6, -1, -1, -1],  # 103 = 01100111
        [11, 2, 3, 9, 4, 10, 6, 4, 10, -1, -1, -1],  # 104 = 01101000
        [0, 8, 2, 11, 8, 2, 9, 4, 10, 6, 4, 10],  # 105 = 01101001
        [2, 3, 11, 0, 1, 6, 10, 1, 6, 0, 6, 4],  # 106 = 01101010
        [10, 2, 1, 8, 4, 11, 6, 4, 11, -1, -1, -1],  # 107 = 01101011
        [1, 4, 9, 1, 11, 4, 1, 3, 11, 6, 4, 11],  # 108 = 01101100
        [1, 0, 9, 11, 6, 8, 4, 6, 8, -1, -1, -1],  # 109 = 01101101
        [3, 11, 0, 6, 11, 0, 0, 6, 4, -1, -1, -1],  # 110 = 01101110
        [8, 4, 11, 6, 4, 11, -1, -1, -1, -1, -1, -1],  # 111 = 01101111
        [10, 6, 8, 7, 6, 8, 8, 9, 10, -1, -1, -1],  # 112 = 01110000
        [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10],  # 113 = 01110001
        [1, 10, 6, 0, 1, 7, 6, 1, 7, 0, 7, 8],  # 114 = 01110010
        [10, 6, 1, 7, 6, 1, 1, 7, 3, -1, -1, -1],  # 115 = 01110011
        [9, 1, 2, 9, 2, 7, 9, 7, 8, 7, 2, 6],  # 116 = 01110100
        [9, 1, 0, 7, 3, 6, 2, 6, 3, -1, -1, -1],  # 117 = 01110101
        [0, 8, 6, 7, 8, 6, 6, 0, 2, -1, -1, -1],  # 118 = 01110110
        [2, 3, 6, 7, 3, 6, -1, -1, -1, -1, -1, -1],  # 119 = 01110111
        [11, 2, 3, 10, 6, 8, 7, 6, 8, 10, 8, 9],  # 120 = 01111000
        [11, 6, 7, 9, 0, 10, 2, 0, 10, -1, -1, -1],  # 121 = 01111001
        [6, 7, 11, 0, 3, 8, 2, 1, 10, -1, -1, -1],  # 122 = 01111010
        [2, 1, 10, 7, 11, 6, -1, -1, -1, -1, -1, -1],  # 123 = 01111011
        [6, 7, 11, 9, 1, 8, 3, 1, 8, -1, -1, -1],  # 124 = 01111100
        [9, 1, 0, 6, 7, 11, -1, -1, -1, -1, -1, -1],  # 125 = 01111101
        [8, 0, 3, 6, 7, 11, -1, -1, -1, -1, -1, -1],  # 126 = 01111110
        [7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 127 = 01111111
        [7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 128 = 10000000
        [3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1],  # 129 = 10000001
        [0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1],  # 130 = 10000010
        [11, 7, 6, 9, 1, 8, 3, 1, 8, -1, -1, -1],  # 131 = 10000011
        [10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1],  # 132 = 10000100
        [11, 7, 6, 8, 3, 0, 10, 1, 2, -1, -1, -1],  # 133 = 10000101
        [7, 6, 11, 9, 0, 10, 2, 0, 10, -1, -1, -1],  # 134 = 10000110
        [7, 6, 11, 10, 2, 8, 3, 2, 8, 9, 8, 10],  # 135 = 10000111
        [2, 3, 6, 7, 3, 6, -1, -1, -1, -1, -1, -1],  # 136 = 10001000
        [0, 8, 6, 7, 8, 6, 6, 2, 0, -1, -1, -1],  # 137 = 10001001
        [0, 1, 9, 6, 3, 7, 3, 6, 2, -1, -1, -1],  # 138 = 10001010
        [9, 2, 1, 9, 7, 2, 9, 8, 7, 7, 6, 2],  # 139 = 10001011
        [10, 6, 1, 7, 6, 1, 1, 3, 7, -1, -1, -1],  # 140 = 10001100
        [1, 6, 10, 0, 1, 7, 6, 1, 7, 0, 8, 7],  # 141 = 10001101
        [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7],  # 142 = 10001110
        [10, 6, 8, 7, 6, 8, 8, 10, 9, -1, -1, -1],  # 143 = 10001111
        [8, 4, 11, 6, 4, 11, -1, -1, -1, -1, -1, -1],  # 144 = 10010000
        [3, 11, 0, 6, 11, 0, 0, 4, 6, -1, -1, -1],  # 145 = 10010001
        [9, 0, 1, 11, 6, 8, 4, 6, 8, -1, -1, -1],  # 146 = 10010010
        [1, 9, 4, 1, 4, 11, 1, 11, 3, 6, 11, 4],  # 147 = 10010011
        [1, 2, 10, 8, 4, 11, 6, 4, 11, -1, -1, -1],  # 148 = 10010100
        [10, 1, 2, 0, 3, 6, 11, 3, 6, 4, 6, 0],  # 149 = 10010101
        [9, 0, 10, 2, 10, 0, 6, 8, 4, 8, 6, 11],  # 150 = 10010110
        [3, 2, 11, 9, 4, 10, 6, 4, 10, -1, -1, -1],  # 151 = 10010111
        [8, 3, 4, 2, 3, 4, 4, 6, 2, -1, -1, -1],  # 152 = 10011000
        [2, 0, 6, 4, 0, 6, -1, -1, -1, -1, -1, -1],  # 153 = 10011001
        [0, 1, 9, 3, 8, 2, 4, 2, 8, 2, 4, 6],  # 154 = 10011010
        [1, 9, 2, 4, 9, 2, 2, 4, 6, -1, -1, -1],  # 155 = 10011011
        [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1],  # 156 = 10011100
        [0, 1, 6, 10, 1, 6, 6, 0, 4, -1, -1, -1],  # 157 = 10011101
        [8, 0, 3, 10, 6, 9, 4, 6, 9, -1, -1, -1],  # 158 = 10011110
        [10, 9, 6, 4, 9, 6, -1, -1, -1, -1, -1, -1],  # 159 = 10011111
        [4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1],  # 160 = 10100000
        [3, 0, 8, 5, 4, 9, 7, 6, 11, -1, -1, -1],  # 161 = 10100001
        [11, 7, 6, 1, 5, 0, 4, 0, 5, -1, -1, -1],  # 162 = 10100010
        [7, 6, 11, 8, 4, 3, 5, 4, 3, 5, 3, 1],  # 163 = 10100011
        [1, 2, 10, 7, 6, 11, 5, 4, 9, -1, -1, -1],  # 164 = 10100100
        [11, 7, 6, 1, 2, 10, 8, 3, 0, 9, 5, 4],  # 165 = 10100101
        [11, 7, 6, 4, 10, 5, 10, 4, 2, 0, 2, 4],  # 166 = 10100110
        [11, 3, 2, 8, 7, 4, 10, 5, 6, -1, -1, -1],  # 167 = 10100111
        [9, 5, 4, 3, 7, 2, 6, 2, 7, -1, -1, -1],  # 168 = 10101000
        [9, 5, 4, 8, 7, 0, 6, 7, 0, 2, 0, 6],  # 169 = 10101001
        [7, 2, 3, 2, 7, 6, 0, 1, 4, 5, 4, 1],  # 170 = 10101010
        [7, 4, 8, 1, 2, 5, 6, 2, 5, -1, -1, -1],  # 171 = 10101011
        [5, 4, 9, 10, 6, 1, 7, 6, 1, 7, 1, 3],  # 172 = 10101100
        [9, 1, 0, 10, 5, 6, 8, 7, 4, -1, -1, -1],  # 173 = 10101101
        [5, 6, 10, 3, 0, 7, 4, 0, 7, -1, -1, -1],  # 174 = 10101110
        [6, 10, 5, 8, 7, 4, -1, -1, -1, -1, -1, -1],  # 175 = 10101111
        [9, 5, 11, 6, 5, 11, 11, 8, 9, -1, -1, -1],  # 176 = 10110000
        [11, 5, 6, 3, 5, 11, 3, 9, 5, 3, 0, 9],  # 177 = 10110001
        [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11],  # 178 = 10110010
        [3, 11, 5, 6, 11, 5, 5, 3, 1, -1, -1, -1],  # 179 = 10110011
        [10, 1, 2, 5, 6, 9, 11, 9, 6, 9, 11, 8],  # 180 = 10110100
        [5, 6, 10, 3, 2, 11, 1, 0, 9, -1, -1, -1],  # 181 = 10110101
        [5, 6, 10, 8, 0, 11, 2, 0, 11, -1, -1, -1],  # 182 = 10110110
        [11, 3, 2, 5, 6, 10, -1, -1, -1, -1, -1, -1],  # 183 = 10110111
        [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2],  # 184 = 10111000
        [9, 5, 0, 6, 5, 0, 0, 6, 2, -1, -1, -1],  # 185 = 10111001
        [8, 0, 3, 2, 6, 1, 5, 6, 1, -1, -1, -1],  # 186 = 10111010
        [1, 5, 2, 6, 5, 2, -1, -1, -1, -1, -1, -1],  # 187 = 10111011
        [10, 5, 6, 8, 3, 9, 1, 3, 9, -1, -1, -1],  # 188 = 10111100
        [1, 0, 9, 6, 10, 5, -1, -1, -1, -1, -1, -1],  # 189 = 10111101
        [3, 8, 0, 6, 10, 5, -1, -1, -1, -1, -1, -1],  # 190 = 10111110
        [10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 191 = 10111111
        [11, 10, 7, 5, 10, 7, -1, -1, -1, -1, -1, -1],  # 192 = 11000000
        [8, 3, 0, 10, 5, 11, 7, 5, 11, -1, -1, -1],  # 193 = 11000001
        [0, 1, 9, 11, 7, 10, 5, 7, 10, -1, -1, -1],  # 194 = 11000010
        [10, 5, 11, 7, 5, 11, 8, 9, 3, 1, 9, 3],  # 195 = 11000011
        [1, 2, 7, 11, 2, 7, 7, 5, 1, -1, -1, -1],  # 196 = 11000100
        [3, 0, 8, 1, 2, 7, 11, 2, 7, 1, 7, 5],  # 197 = 11000101
        [0, 5, 9, 0, 11, 5, 0, 2, 11, 11, 7, 5],  # 198 = 11000110
        [11, 3, 2, 9, 7, 5, 7, 9, 8, -1, -1, -1],  # 199 = 11000111
        [2, 10, 3, 5, 10, 3, 3, 7, 5, -1, -1, -1],  # 200 = 11001000
        [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5],  # 201 = 11001001
        [9, 0, 1, 10, 2, 5, 3, 2, 5, 7, 5, 3],  # 202 = 11001010
        [2, 1, 10, 7, 5, 8, 9, 8, 5, -1, -1, -1],  # 203 = 11001011
        [1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1],  # 204 = 11001100
        [0, 8, 1, 7, 8, 1, 1, 7, 5, -1, -1, -1],  # 205 = 11001101
        [9, 0, 5, 3, 0, 5, 5, 3, 7, -1, -1, -1],  # 206 = 11001110
        [9, 8, 5, 7, 8, 5, -1, -1, -1, -1, -1, -1],  # 207 = 11001111
        [8, 4, 10, 5, 4, 10, 10, 11, 8, -1, -1, -1],  # 208 = 11010000
        [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0],  # 209 = 11010001
        [9, 0, 1, 8, 4, 10, 5, 4, 10, 8, 10, 11],  # 210 = 11010010
        [4, 5, 9, 3, 1, 11, 10, 11, 1, -1, -1, -1],  # 211 = 11010011
        [11, 1, 2, 11, 4, 1, 11, 8, 4, 5, 1, 4],  # 212 = 11010100
        [11, 3, 2, 5, 1, 4, 0, 4, 1, -1, -1, -1],  # 213 = 11010101
        [9, 4, 5, 0, 2, 8, 11, 2, 8, -1, -1, -1],  # 214 = 11010110
        [4, 5, 9, 11, 3, 2, -1, -1, -1, -1, -1, -1],  # 215 = 11010111
        [3, 10, 2, 8, 10, 3, 8, 5, 10, 8, 4, 5],  # 216 = 11011000
        [2, 10, 4, 5, 10, 4, 4, 2, 0, -1, -1, -1],  # 217 = 11011001
        [4, 5, 9, 2, 1, 10, 0, 3, 8, -1, -1, -1],  # 218 = 11011010
        [10, 2, 1, 4, 5, 9, -1, -1, -1, -1, -1, -1],  # 219 = 11011011
        [8, 4, 3, 5, 4, 3, 3, 5, 1, -1, -1, -1],  # 220 = 11011100
        [0, 4, 1, 5, 4, 1, -1, -1, -1, -1, -1, -1],  # 221 = 11011101
        [4, 5, 9, 3, 8, 0, -1, -1, -1, -1, -1, -1],  # 222 = 11011110
        [9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 223 = 11011111
        [11, 7, 9, 4, 7, 9, 9, 10, 11, -1, -1, -1],  # 224 = 11100000
        [8, 3, 0, 7, 4, 11, 9, 11, 4, 11, 9, 10],  # 225 = 11100001
        [10, 11, 7, 10, 7, 0, 10, 0, 1, 4, 0, 7],  # 226 = 11100010
        [8, 7, 4, 10, 3, 1, 3, 10, 11, -1, -1, -1],  # 227 = 11100011
        [4, 11, 7, 9, 4, 2, 11, 4, 2, 9, 1, 2],  # 228 = 11100100
        [7, 4, 8, 1, 0, 9, 3, 2, 11, -1, -1, -1],  # 229 = 11100101
        [11, 7, 2, 4, 7, 2, 2, 4, 0, -1, -1, -1],  # 230 = 11100110
        [7, 4, 8, 2, 11, 3, -1, -1, -1, -1, -1, -1],  # 231 = 11100111
        [3, 10, 2, 3, 4, 10, 3, 7, 4, 4, 9, 10],  # 232 = 11101000
        [7, 4, 8, 0, 2, 9, 10, 2, 9, -1, -1, -1],  # 233 = 11101001
        [10, 2, 1, 0, 4, 3, 7, 4, 3, -1, -1, -1],  # 234 = 11101010
        [10, 2, 1, 7, 4, 8, -1, -1, -1, -1, -1, -1],  # 235 = 11101011
        [1, 9, 7, 4, 9, 7, 7, 1, 3, -1, -1, -1],  # 236 = 11101100
        [9, 1, 0, 7, 4, 8, -1, -1, -1, -1, -1, -1],  # 237 = 11101101
        [3, 0, 7, 4, 0, 7, -1, -1, -1, -1, -1, -1],  # 238 = 11101110
        [4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 239 = 11101111
        [8, 9, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1],  # 240 = 11110000
        [9, 0, 11, 3, 0, 11, 11, 9, 10, -1, -1, -1],  # 241 = 11110001
        [0, 1, 8, 10, 1, 8, 8, 10, 11, -1, -1, -1],  # 242 = 11110010
        [10, 1, 11, 3, 1, 11, -1, -1, -1, -1, -1, -1],  # 243 = 11110011
        [1, 2, 9, 11, 2, 9, 9, 11, 8, -1, -1, -1],  # 244 = 11110100
        [0, 9, 1, 11, 3, 2, -1, -1, -1, -1, -1, -1],  # 245 = 11110101
        [0, 2, 8, 11, 2, 8, -1, -1, -1, -1, -1, -1],  # 246 = 11110110
        [3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 247 = 11110111
        [8, 3, 10, 2, 3, 10, 10, 8, 9, -1, -1, -1],  # 248 = 11111000
        [9, 10, 0, 2, 10, 0, -1, -1, -1, -1, -1, -1],  # 249 = 11111001
        [3, 8, 0, 10, 2, 1, -1, -1, -1, -1, -1, -1],  # 250 = 11111010
        [1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 251 = 11111011
        [8, 3, 9, 1, 3, 9, -1, -1, -1, -1, -1, -1],  # 252 = 11111100
        [0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 253 = 11111101
        [0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 254 = 11111110
        [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1],  # 255 = 11111111
])

cdef cnp.float64_t [:] VOLUME_CASE_LOOKUP = np.array([
    0,
    1,         # case 1 
    1/48,      # case 2
    1/8,       # case 3
    1/24,      # case 4
    1/24,      # case 5
    17/48,     # case 6
    7/48,      # case 7
    1/16,      # case 8
    1/2,       # case 9
    1/2,       # case 10
    1/4,       # case 11
    1/2,       # case 12
    139/144,   # case 13
    1/12,      # case 14
    5/12,      # case 15
    1-1,       # case 16 
    1-1/48,    # case 17
    1-1/8,     # case 18
    1-1/24,    # case 19
    1-1/24,    # case 20
    1-17/48,   # case 21
    1-7/48,    # case 22
    1-1/16,    # case 23
    1-1/2,     # case 24
    1-1/2,     # case 25
    1-1/4,     # case 26
    1 - 1/2,   # case 27
    1-139/144, # case 28
    1-1/12,    # case 29
    1-5/12,    # case 30
]).astype(np.float64)

cdef int [:] VOLUME_LOOKUP = np.array([
        1,   # 0 = 00000000
        2,   # 1 = 00000001
        2,   # 2 = 00000010
        3,   # 3 = 00000011
        1,   # 4 = 00000100
        2,   # 5 = 00000101
        3,   # 6 = 00000110
        6,   # 7 = 00000111
        2,   # 8 = 00001000
        2,   # 9 = 00001001
        2,   # 10 = 00001010
        7,   # 11 = 00001011
        3,   # 12 = 00001100
        7,   # 13 = 00001101
        6,   # 14 = 00001110
        9,   # 15 = 00001111
        2,   # 16 = 00010000
        3,   # 17 = 00010001
        4,   # 18 = 00010010
        6,   # 19 = 00010011
        4,   # 20 = 00010100
        7,   # 21 = 00010101
        7,   # 22 = 00010110
        12,  # 23 = 00010111
        4,   # 24 = 00011000
        6,   # 25 = 00011001
        8,   # 26 = 00011010
        15,  # 27 = 00011011
        7,   # 28 = 00011100
        15,  # 29 = 00011101
        13,  # 30 = 00011110
        21,  # 31 = 00011111
        2,   # 32 = 00100000
        4,   # 33 = 00100001
        3,   # 34 = 00100010
        6,   # 35 = 00100011
        4,   # 36 = 00100100
        8,   # 37 = 00100101
        6,   # 38 = 00100110
        10,  # 39 = 00100111
        5,   # 40 = 00101000
        7,   # 41 = 00101001
        7,   # 42 = 00101010
        15,  # 43 = 00101011
        7,   # 44 = 00101100
        13,  # 45 = 00101101
        12,  # 46 = 00101110
        21,  # 47 = 00101111
        3,   # 48 = 00110000
        6,   # 49 = 00110001
        6,   # 50 = 00110010
        9,   # 51 = 00110011
        7,   # 52 = 00110100
        13,  # 53 = 00110101
        15,  # 54 = 00110110
        21,  # 55 = 00110111
        7,   # 56 = 00111000
        12,  # 57 = 00111001
        13,  # 58 = 00111010
        21,  # 59 = 00111011
        11,  # 60 = 00111100
        22,  # 61 = 00111101
        22,  # 62 = 00111110
        18,  # 63 = 00111111
        2,   # 64 = 01000000
        5,   # 65 = 01000001
        4,   # 66 = 01000010
        7,   # 67 = 01000011
        3,   # 68 = 01000100
        7,   # 69 = 01000101
        6,   # 70 = 01000110
        15,  # 71 = 01000111
        4,   # 72 = 01001000
        7,   # 73 = 01001001
        7,   # 74 = 01001010
        15,  # 75 = 01001011
        7,   # 76 = 01001100
        13,  # 77 = 01001101
        12,  # 78 = 01001110
        21,  # 79 = 01001111
        4,   # 80 = 01010000
        7,   # 81 = 01010001
        8,   # 82 = 01010010
        13,  # 83 = 01010011
        3,   # 84 = 01010100
        11,  # 85 = 01010101
        13,  # 86 = 01010110
        7,   # 87 = 01010111
        8,   # 88 = 01011000
        13,  # 89 = 01011001
        14,  # 90 = 01011010
        23,  # 91 = 01011011
        22,  # 92 = 01011100
        20,  # 93 = 01011101
        23,  # 94 = 01011110
        19,  # 95 = 01011111
        3,   # 96 = 01100000
        7,   # 97 = 01100001
        6,   # 98 = 01100010
        12,  # 99 = 01100011
        6,   # 100 = 01100100
        13,  # 101 = 01100101
        9,   # 102 = 01100110
        21,  # 103 = 01100111
        7,   # 104 = 01101000
        11,  # 105 = 01101001
        13,  # 106 = 01101010
        22,  # 107 = 01101011
        15,  # 108 = 01101100
        22,  # 109 = 01101101
        21,  # 110 = 01101110
        18,  # 111 = 01101111
        6,   # 112 = 01110000
        15,  # 113 = 01110001
        10,  # 114 = 01110010
        21,  # 115 = 01110011
        12,  # 116 = 01110100
        20,  # 117 = 01110101
        21,  # 118 = 01110110
        19,  # 119 = 01110111
        13,  # 120 = 01111000
        20,  # 121 = 01111001
        19,  # 122 = 01111010
        19,  # 123 = 01111011
        22,  # 124 = 01111100
        20,  # 125 = 01111101
        19,  # 126 = 01111110
        17,  # 127 = 01111111
        2,   # 128 = 10000000
        4,   # 129 = 10000001
        5,   # 130 = 10000010
        7,   # 131 = 10000011
        4,   # 132 = 10000100
        8,   # 133 = 10000101
        7,  # 134 = 10000110
        13, # 135 = 10000111
        7,  # 136 = 10001000
        6,  # 137 = 10001001
        7,  # 138 = 10001010
        12, # 139 = 10001011
        6,  # 140 = 10001100
        10, # 141 = 10001101
        15, # 142 = 10001110
        21, # 143 = 10001111
        7,  # 144 = 10010000
        6,  # 145 = 10010001
        7,  # 146 = 10010010
        15, # 147 = 10010011
        7,  # 148 = 10010100
        13, # 149 = 10010101
        11, # 150 = 10010110
        22, # 151 = 10010111
        6,  # 152 = 10011000
        9,  # 153 = 10011001
        13, # 154 = 10011010
        21, # 155 = 10011011
        12, # 156 = 10011100
        21, # 157 = 10011101
        22, # 158 = 10011110
        18, # 159 = 10011111
        4,  # 160 = 10100000
        8,  # 161 = 10100001
        7,  # 162 = 10100010
        13, # 163 = 10100011
        8,  # 164 = 10100100
        14, # 165 = 10100101
        13, # 166 = 10100110
        20, # 167 = 10100111
        7,  # 168 = 10101000
        13, # 169 = 10101001
        11, # 170 = 10101010
        7,  # 171 = 10101011
        13, # 172 = 10101100
        23, # 173 = 10101101
        22, # 174 = 10101110
        19, # 175 = 10101111
        6,  # 176 = 10110000
        10, # 177 = 10110001
        12, # 178 = 10110010
        18, # 179 = 10110011
        13, # 180 = 10110100
        19, # 181 = 10110101
        20, # 182 = 10110110
        19, # 183 = 10110111
        15, # 184 = 10111000
        21, # 185 = 10111001
        22, # 186 = 10111010
        18, # 187 = 10111011
        22, # 188 = 10111100
        19, # 189 = 10111101
        20, # 190 = 10111110
        17, # 191 = 10111111
        3,  # 192 = 11000000
        7,  # 193 = 11000001
        7,  # 194 = 11000010
        11, # 195 = 11000011
        6,  # 196 = 11000100
        13, # 197 = 11000101
        12, # 198 = 11000110
        3,  # 199 = 11000111
        6,  # 200 = 11001000
        15, # 201 = 11001001
        13, # 202 = 11001010
        22, # 203 = 11001011
        9,  # 204 = 11001100
        21, # 205 = 11001101
        21, # 206 = 11001110
        18, # 207 = 11001111
        6,  # 208 = 11010000
        12, # 209 = 11010001
        13, # 210 = 11010010
        19, # 211 = 11010011
        15, # 212 = 11010100
        18, # 213 = 11010101
        19, # 214 = 11010110
        20, # 215 = 11010111
        10, # 216 = 11011000
        21, # 217 = 11011001
        23, # 218 = 11011010
        19, # 219 = 11011011
        21, # 220 = 11011100
        18, # 221 = 11011101
        19, # 222 = 11011110
        17, # 223 = 11011111
        6,  # 224 = 11100000
        13, # 225 = 11100001
        15, # 226 = 11100010
        22, # 227 = 11100011
        10, # 228 = 11100100
        23, # 229 = 11100101
        21, # 230 = 11100110
        19, # 231 = 11100111
        12, # 232 = 11101000
        22, # 233 = 11101001
        22, # 234 = 11101010
        20, # 235 = 11101011
        21, # 236 = 11101100
        19, # 237 = 11101101
        18, # 238 = 11101110
        17, # 239 = 11101111
        9,  # 240 = 11110000
        21, # 241 = 11110001
        21, # 242 = 11110010
        18, # 243 = 11110011
        21, # 244 = 11110100
        19, # 245 = 11110101
        18, # 246 = 11110110
        17, # 247 = 11110111
        21, # 248 = 11111000
        18, # 249 = 11111001
        19, # 250 = 11111010
        17, # 251 = 11111011
        18, # 252 = 11111100
        17, # 253 = 11111101
        17, # 254 = 11111110
        16, # 255 = 11111111
])

cdef cnp.float64_t interpolation(cnp.float64_t a, cnp.float64_t b, cnp.float64_t level) nogil: 
    return (a - level) / (a - b + 1e-20)

cdef int calculate_vertex_id (int i, int j, int k, int N, int M, int direction):
    return (i + j * N + k * N * M) * 3 + direction

cdef int edge_to_vertex_id(int i, int j, int k, int N, int M, int edge_number):
    cdef int dx, dy, dz, direction
    dx, dy, dz = EDGE_DELTA[edge_number]
    direction = EDGE_DIRECTION[edge_number]
    return calculate_vertex_id(i + dx, j + dy, k + dz, N, M, direction)


cdef class FenwickTree3D:
    cdef cnp.float64_t [:,:,:] bit
    cdef int n, m, p

    def __cinit__(self, int n, int m, int p):
        self.n = n
        self.m = m
        self.p = p
        self.bit = np.array([[[0.0] * (p + 1) for _ in range(m + 1)] for _ in range(n + 1)]).astype(np.float64)

    cpdef update(self, int x, int y, int z, cnp.float64_t val):
        cdef int x1, y1, z1
        x += 1
        y += 1
        z += 1
        while x <= self.n:
            y1 = y
            while y1 <= self.m:
                z1 = z
                while z1 <= self.p:
                    self.bit[x][y1][z1] += val
                    z1 += z1 & -z1
                y1 += y1 & -y1
            x += x & -x
        
    cpdef cnp.float64_t getSum(self, int x, int y, int z):
        cdef int x1, y1, z1
        x += 1
        y += 1
        z += 1
        cdef cnp.float64_t res = 0
        while x > 0:
            y1 = y
            while y1 > 0:
                z1 = z
                while z1 > 0:
                    res += self.bit[x][y1][z1]
                    z1 -= z1 & -z1
                y1 -= y1 & -y1
            x -= x & -x
        return res
        
    cpdef cnp.float64_t queryByRange(self, int x1, int y1, int z1, int x2, int y2, int z2):
        return self.getSum(x2, y2, z2) - self.getSum(x2, y2, z1 - 1) - self.getSum(x2, y1 - 1, z2) + self.getSum(x2, y1 - 1, z1 - 1) \
        - self.getSum(x1 - 1, y2, z2) + self.getSum(x1 - 1, y2, z1 - 1) + self.getSum(x1 - 1, y1 - 1, z2) - self.getSum(x1 - 1, y1 - 1, z1 - 1)
        
    
def MarchingCubesLorensen(cnp.float64_t [:, :, :] volume not None, cnp.int32_t [:, :, :] mask not None, cnp.float64_t level):

    # Initialize variables
    cdef cnp.float64_t Sum=0  # use this only for checking 
    cdef cnp.float64_t delta, CurVol
    cdef unsigned int i, j, k, t
    cdef unsigned int N, M, P
    cdef list vertices = []
    cdef list vertex_ids = []
    cdef list triangles = []
    cdef list triangle_ids = []
    cdef int volume_type
    cdef int vertex_id0, vertex_id1, vertex_id2
    cdef int edge0, edge1, edge2

    # Intialize 3D Fenwick Tree
    N, M, P = volume.shape[0], volume.shape[1], volume.shape[2]
    fenwick_tree = FenwickTree3D(N, M, P)

    i,j,k = 0,0,0
    for k in range(P):
        for j in range (M):
            for i in range (N):
                if i < (N-1) and mask[i, j, k] != mask[i + 1, j, k]:
                    delta = interpolation(volume[i, j, k], volume[i + 1, j, k], level)
                    vertices.append([i + delta, j, k])
                    vertex_ids.append(calculate_vertex_id(i, j, k, N, M, DirectionX))
                
                if j < (M-1) and mask[i, j, k] != mask[i, j + 1, k]:
                    delta = interpolation(volume[i, j, k], volume[i, j + 1, k], level)
                    vertices.append([i, j + delta, k])
                    vertex_ids.append(calculate_vertex_id(i, j, k, N, M, DirectionY))

                if k < (P-1) and mask[i, j, k] != mask[i, j, k + 1]:
                    delta = interpolation(volume[i, j, k], volume[i, j, k + 1], level)
                    vertices.append([i, j, k + delta])
                    vertex_ids.append(calculate_vertex_id(i, j, k, N, M, DirectionZ))
            
                if i == (N - 1) or j == (M - 1) or k == (P - 1):
                    continue

                volume_type = 0
                if mask[i, j, k]:
                    volume_type |= 1<<0
                if mask[i + 1, j, k]:
                    volume_type |= 1<<1
                if mask[i + 1, j + 1, k]:
                    volume_type |= 1<<2
                if mask[i, j + 1, k]:
                    volume_type |= 1<<3
                if mask[i, j, k + 1]:
                    volume_type |= 1<<4
                if mask[i + 1, j, k + 1]:
                    volume_type |= 1<<5
                if mask[i + 1, j + 1, k + 1]:
                    volume_type |= 1<<6
                if mask[i, j + 1, k + 1]:
                    volume_type |= 1<<7

                CurVol = VOLUME_CASE_LOOKUP[VOLUME_LOOKUP[volume_type]]
                Sum += CurVol  # --> use this only for checking
                fenwick_tree.update(i, j, k, CurVol)

                lookup = GEOMETRY_LOOKUP[volume_type]
                t = 0
                while t < len(lookup):
                    if lookup[t] < 0:
                        break
                    edge0, edge1, edge2 = lookup[t : t + 3]
                    vertex_id0 = edge_to_vertex_id(i, j, k, N, M, edge0)
                    vertex_id1 = edge_to_vertex_id(i, j, k, N, M, edge1)
                    vertex_id2 = edge_to_vertex_id(i, j, k, N, M, edge2)
                    triangle_ids.append([vertex_id0, vertex_id1, vertex_id2])
                    t += 3
    
    order_of_ids = {id: order for order, id in enumerate(vertex_ids)}
    for triangle_corners in triangle_ids:
        triangles.append([order_of_ids[c] for c in triangle_corners])
    
    return vertices, triangles, fenwick_tree, Sum