voxelPolygonizer.lua

local bit32 = require("bit32"--[[
   Linearly interpolate the position where an isosurface cuts
   an edge between two vertices, each with their own scalar value
]]
local VertexInterp = function(isolevel,p1,p2,valp1,valp2)
    --double isolevel;
    --XYZ p1,p2;
    --double valp1,valp2;
    local mu --double mu;
    local p = {x = 0, y = 0, z = 0}--XYZ p;

    if (abs(isolevel-valp1) < 0.00001) then
        return p1
    end
    if (abs(isolevel-valp2) < 0.00001) then
        return p2
    end
    if (abs(valp1-valp2) < 0.00001) then
        return p1
    end
    mu = (isolevel - valp1) / (valp2 - valp1);
    p.x = p1.x + mu * (p2.x - p1.x);
    p.y = p1.y + mu * (p2.y - p1.y);
    p.z = p1.z + mu * (p2.z - p1.z);

    return p
end

--[[
    Given a grid cell and an isolevel, calculate the triangular
    facets required to represent the isosurface through the cell.
    Return the number of triangular facets, the array "triangles"
    will be loaded up with the vertices at most 5 triangular facets.
    0 will be returned if the grid cell is either totally above
    of totally below the isolevel.
]]
-- Gridcell, number, triangle
function voxelPolygonizer.Polygonise(grid, isolevel, triangles)
    local i, ntriang
    local cubeindex = 0
    local vertlist = {}

    --[[
        Determine the index into the edge table which
        tells us which vertices are inside of the surface
    ]]
    if (grid.val[0] < isolevel) then cubeindex = bor(cubeindex, 1) end
    if (grid.val[1] < isolevel) then cubeindex = bor(cubeindex, 2) end
    if (grid.val[2] < isolevel) then cubeindex = bor(cubeindex, 4) end
    if (grid.val[3] < isolevel) then cubeindex = bor(cubeindex, 8) end
    if (grid.val[4] < isolevel) then cubeindex = bor(cubeindex, 16) end
    if (grid.val[5] < isolevel) then cubeindex = bor(cubeindex, 32) end
    if (grid.val[6] < isolevel) then cubeindex = bor(cubeindex, 64) end
    if (grid.val[7] < isolevel) then cubeindex = bor(cubeindex, 128) end

    cubeindex = cubeindex + 1

    -- Cube is entirely in/out of the surface
    if (edgeTable[cubeindex] == 0) then
        return 0
    end

    -- Find the vertices where the surface intersects the cube
    if (band(edgeTable[cubeindex], 1) > 0) then
        vertlist[0] = VertexInterp(isolevel,grid.p[0],grid.p[1],grid.val[0],grid.val[1]);
    end
    if (band(edgeTable[cubeindex],  2) > 0) then
        vertlist[1] = VertexInterp(isolevel,grid.p[1],grid.p[2],grid.val[1],grid.val[2]);
    end
    if (band(edgeTable[cubeindex],  4) > 0) then
        vertlist[2] = VertexInterp(isolevel,grid.p[2],grid.p[3],grid.val[2],grid.val[3]);
    end
    if (band(edgeTable[cubeindex],  8) > 0) then
        vertlist[3] = VertexInterp(isolevel,grid.p[3],grid.p[0],grid.val[3],grid.val[0]);
    end
    if (band(edgeTable[cubeindex],  16) > 0) then
        vertlist[4] = VertexInterp(isolevel,grid.p[4],grid.p[5],grid.val[4],grid.val[5]);
    end
    if (band(edgeTable[cubeindex],  32) > 0) then
        vertlist[5] = VertexInterp(isolevel,grid.p[5],grid.p[6],grid.val[5],grid.val[6]);
    end
    if (band(edgeTable[cubeindex],  64) > 0) then
        vertlist[6] = VertexInterp(isolevel,grid.p[6],grid.p[7],grid.val[6],grid.val[7]);
    end
    if (band(edgeTable[cubeindex],  128) > 0) then
        vertlist[7] = VertexInterp(isolevel,grid.p[7],grid.p[4],grid.val[7],grid.val[4]);
    end
    if (band(edgeTable[cubeindex],  256) > 0) then
        vertlist[8] = VertexInterp(isolevel,grid.p[0],grid.p[4],grid.val[0],grid.val[4]);
    end
    if (band(edgeTable[cubeindex],  512) > 0) then
        vertlist[9] = VertexInterp(isolevel,grid.p[1],grid.p[5],grid.val[1],grid.val[5]);
    end
    if (band(edgeTable[cubeindex],  1024) > 0) then
        vertlist[10] = VertexInterp(isolevel,grid.p[2],grid.p[6],grid.val[2],grid.val[6]);
    end
    if (band(edgeTable[cubeindex],  2048) > 0) then
        vertlist[11] = VertexInterp(isolevel,grid.p[3],grid.p[7],grid.val[3],grid.val[7]);
    end

    -- Create the triangle
    ntriang = 0
    i = 0
    --for (i=0, triTable[cubeindex][i]!=-1, i+=3) do
    while triTable[cubeindex][i] ~= -1 do
        triangles[ntriang] = {p = {}}
        triangles[ntriang].p[0] = vertlist[triTable[cubeindex][i+1]];
        triangles[ntriang].p[1] = vertlist[triTable[cubeindex][i+2]];
        triangles[ntriang].p[2] = vertlist[triTable[cubeindex][i+3]];
        ntriang = ntriang + 1;
        i = i + 3
        if i > 90 then break end
    end

    return ntriang
end

return voxelPolygonizer