Surface extraction algorithms in F#, currently contains only Surface Nets

SurfaceExtraction.fs

/// <summary>
/// Precalculated values which are power of two. Exponent begins with zero!
/// </summary>
let PowerOfTwo = [| for i in 0..20 -> int (System.Math.Pow(2., float i)) |]
let toBoolean<'T> (v:'T) = System.Convert.ToBoolean(v)

// size = 12
// This maps a edge index to a vertex pair
// 12 edges = 12 vertex pairs
let cubeEdgesSurfacenets = 
    [| 
        for i in 0..7 do
            for j in [|1; 2; 4|] do
                let p = i ^^^ j
                if i <= p then yield i, p 
    |]
    
let cubeEdgesSurfacenetsStatic =
    [|  
        (0, 1) // 0
        (0, 2) // 1
        (0, 4) // 2
        (1, 3) // 3
        (1, 5) // 4
        (2, 3) // 5
        (2, 6) // 6
        (3, 7) // 7
        (4, 5) // 8
        (4, 6) // 9
        (5, 7) // 10
        (6, 7) // 11
    |]

let edge_table_surfacenets = 
    let convert i (a, b) = 
        i &&& (1 <<< a) |> toBoolean,
        i &&& (1 <<< b)|> toBoolean

    let calcEm i =
        let rec calc j em =
            if j < 12 then
                let a, b = cubeEdgesSurfacenets.[j] |> convert i
                let res = em ||| (if a <> b then (1 <<< j) else 0)
                calc (j + 1) res 
            else
                em
        calc 0 0
        
    Array.init 256 calcEm

module Interpolate=
    /// <summary>
    /// Generates vertieces and faces for a given voxel input.
    /// </summary>
    /// <param name="scalars">Voxel already ordered in 3D</param>
    /// <param name="isovalue">The value for identifying a surface crossing</param>
    /// <param name="gridSize">The resolution of the grid. Warning, currently ignored!</param>
    let generateSurfaceNet (scalars:float[,,]) isovalue gridSize = 
        let float a = float(a)

        let vertices = ResizeArray()

        let indicies = System.Collections.Generic.Dictionary()
        let xMax = (Array3D.length1 scalars) - 2
        let yMax = (Array3D.length2 scalars) - 2
        let zMax = (Array3D.length3 scalars) - 2
        let xMaxHalf = float xMax / 2.
        let yMaxHalf = float yMax / 2.
        let zMaxHalf = float zMax / 2.
        
        let faces = ResizeArray()
        let mutable buf_no = 0
        for x1 in 0..xMax do
            for y1 in 0..yMax do
                for z1 in 0..zMax do
                        // create a temporary array containing the scalar values of our current block.
                        // Keep in mind that the order must respect / be the same as the order of the pre calculated table
                        let block = 
                            [|  scalars.[x1, y1, z1]
                                scalars.[x1, y1, z1 + 1]
                                scalars.[x1, y1 + 1, z1]
                                scalars.[x1, y1 + 1, z1 + 1]
                                scalars.[x1 + 1, y1, z1]
                                scalars.[x1 + 1, y1, z1 + 1]
                                scalars.[x1 + 1, y1 + 1, z1]
                                scalars.[x1 + 1, y1 + 1, z1 + 1] |]
                            
                        // 8 bit mask representing which vertex isn't in the surface
                        let mask =
                            let rec calcMask i mask =
                                if i < block.Length then
                                    let res =
                                        if block.[i] < isovalue then 
                                            mask ||| (1 <<< i)
                                        else
                                            mask
                                    calcMask (i + 1) res
                                else
                                    mask
                            calcMask 0 0
                            
                        // Skip if there aren't any intersections
                        if not ((mask = 0x00) || (mask = 0xff)) then 
                            let edge_mask = edge_table_surfacenets.[mask]

                            // Fast interpolation algorithm, multiplication is done in the end
                            let interpolateVertices () =
                                let rec iterate i e_count (vertices:float[]) =
                                    if i < 12 then
                                        if hasFlag PowerOfTwo.[i] edge_mask then
                                            let e0, e1 = cubeEdgesSurfacenets.[i]
                                            let g0, g1 = block.[e0], block.[e1]
                                            let diff = g0 - g1
                                            // Here we can specify whether an interpolation should be discarded if the values are too near
                                            if abs diff > 1e-6 then 
                                                let t = g0 / diff
                                                // See blog post to understand what's happening here!
                                                let res = 
                                                    [|  for j in 0..2 do
                                                            let k = PowerOfTwo.[j]
                                                            let a = hasFlag k e0 //(e0 &&& k)
                                                            let b = hasFlag k e1 //(e1 &&& k)
                                                            let res = 
                                                                if (a <> b) then 
                                                                    t
                                                                elif a then
                                                                    1.
                                                                else
                                                                    0.0
                                                            yield vertices.[j] + res |]
                                                iterate (i + 1) (e_count + 1) res
                                            else
                                                iterate (i + 1) (e_count + 1) vertices
                                        else
                                            iterate (i + 1) e_count vertices
                                    else
                                        vertices, e_count
                                let vertices, e_count = iterate 0 0 [| 0.0; 0.0; 0.0 |]
                                let s = 1.0 / float e_count

                                // finally normalize the values again and afterwards translate them centric to the grid.
                                float z1 + s * vertices.[0] - zMaxHalf,
                                float y1 + s * vertices.[1] - yMaxHalf,
                                float x1 + s * vertices.[2] - xMaxHalf
                                
                            // Add vertex index to the indicies array.
                            indicies.[(x1, y1, z1)] <- vertices.Count
                            // Add the interpolated vertex
                            vertices.Add(interpolateVertices ())

                            // if we're the last point of a combination towards x we create a face
                            if hasFlag 1 edge_mask && not((y1 = 0) || (x1 = 0)) then
                                let a, b, c, d =
                                    if hasFlag 0x1 mask then 
                                        (x1, y1, z1), (x1, y1 - 1, z1), (x1 - 1, y1 - 1, z1), (x1 - 1, y1, z1) 
                                    else 
                                        (x1, y1, z1), (x1 - 1, y1, z1), (x1 - 1, y1 - 1, z1), (x1, y1 - 1, z1)
                                faces.Add((indicies.[a], indicies.[b], indicies.[c], indicies.[d]))
                                
                            // if we're the last point of a combination towards y we create a face
                            if hasFlag 2 edge_mask && not((x1 = 0) || (z1 = 0)) then
                                let a, b, c, d =
                                    if hasFlag 0x1 mask then 
                                        (x1, y1, z1), (x1 - 1, y1, z1), (x1 - 1, y1, z1 - 1), (x1, y1, z1 - 1) 
                                    else 
                                        (x1, y1, z1), (x1, y1, z1 - 1), (x1 - 1, y1, z1 - 1), (x1 - 1, y1, z1)
                                faces.Add((indicies.[a], indicies.[b], indicies.[c], indicies.[d]))
                                
                            // if we're the last point of a combination towards z we create a face
                            if hasFlag 4 edge_mask && not((z1 = 0) || (y1 = 0)) then
                                let a, b, c, d =
                                    if hasFlag 0x1 mask then 
                                        (x1, y1, z1), (x1, y1, z1 - 1), (x1, y1 - 1, z1 - 1), (x1, y1 - 1, z1) 
                                    else 
                                        (x1, y1, z1), (x1, y1 - 1, z1), (x1, y1 - 1, z1 - 1), (x1, y1, z1 - 1)
                                faces.Add((indicies.[a], indicies.[b], indicies.[c], indicies.[d]))
        (vertices, faces)