Iso-Suface Generator

The iso-surface of an implicit function as

    f(x, y, z) = 0

can be created using the iso-suface function of zeMake. The Lua iso-suface function below standarizes the procedure:

require("register")

function iso_surface(iso, xarr, yarr, zarr, func)
    local nx, ny, nz = xarr:size(), yarr:size(), zarr:size()

    local xyz, nor = zeGrf.new("vertex", "vertex")
    
    local v1, v2, vol, tri = zeUtl.new("double", "double", "double", "double")
    v1:resize(ny, nx)
    v2:resize(ny, nx)
    vol:resize(8, 4)
    
    local function get_volume(v, z)
        for i = 0, ny-1 do
            local y = yarr:getele(i, 0)
            for j = 0, nx-1 do
                v:setele(i, j, func(xarr:getele(j, 0), y, z))
            end
        end
    end
    
    get_volume(v1, zarr:getele(0, 0))
        
    for k = 1, nz-1 do
        get_volume(v2, zarr:getele(k, 0))
        for i = 1, ny-1 do
            for j = 1, nx-1 do
                vol:setele(0, -1,
                           xarr:getele(j-1,0),
                           yarr:getele(i-1,0),
                           zarr:getele(k-1,0),
                           v1:getele(i-1,j-1))
                vol:setele(1, -1,
                           xarr:getele(j,0),
                           yarr:getele(i-1,0),
                           zarr:getele(k-1,0),
                           v1:getele(i-1,j))
                vol:setele(2, -1,
                           xarr:getele(j,0),
                           yarr:getele(i,0),
                           zarr:getele(k-1,0),
                           v1:getele(i, j))
                vol:setele(3, -1,
                           xarr:getele(j-1,0),
                           yarr:getele(i,0),
                           zarr:getele(k-1,0),
                           v1:getele(i, j-1))
                vol:setele(4, -1,
                           xarr:getele(j-1,0),
                           yarr:getele(i-1,0),
                           zarr:getele(k,0),
                           v2:getele(i-1,j-1))
                vol:setele(5, -1,
                           xarr:getele(j,0),
                           yarr:getele(i-1,0),
                           zarr:getele(k,0),
                           v2:getele(i-1,j))
                vol:setele(6, -1,
                           xarr:getele(j,0),
                           yarr:getele(i,0),
                           zarr:getele(k,0),
                           v2:getele(i, j))
                vol:setele(7, -1,
                           xarr:getele(j-1,0),
                           yarr:getele(i,0),
                           zarr:getele(k,0),
                           v2:getele(i, j-1))
                zeMake.isosurf(tri, iso, vol)
                local n = tri:size()
                if n >= 3 then
                    xyz:add(tri)
                    tri:shift(3)
                    nor:add(tri)
                end
            end
        end
        v2:copy(v1)
    end
        
    return xyz, nor
end

Inputs to the function are the iso-value; x-, y-, and z-range as zeArray object of double; and the callback function that returns a value for given x, y, and z. The generator return two zeVertex objects containing surface vertices and their normals.

The signs of normals may need to be reversed with the scaling function of zeVertex so that the surface shows properly. If the partial derivatives of the iso-surface function is available, you can use the generator to produce normal and the surface will looks smooth. The partial derivatives may be estimated approximatelly with small x-, y-, and z-step and the iso-surface function. The following example demonstrate the application of the generator to the Gumdrop torus function:

        4 * [x^4 + (y^2 + z^2)^2] +
        17 * x^2 * (y^2 + z^2) -
        20 * (x^2 + y^2 + z^2) + 17 = 0 

The following figure is produced by this testing code for the gumdrop_torus() function

        require("plot_simple")
        plot = plot_simple.new()
        plot:add_static(zeGrf.new("light"))
        require("gumdrop_torus")
        xyz, nor = gumdrop_torus(0.05)
        xyz:scale(80, 80, 80)
        nor:scale(-1, -1, -1)
        shape = zeGrf.new("polygon")
        shape:set{type = "triangles", vertex = xyz,
                  vertex_normal = nor, color = {0, .5, .5, 1}}
        plot:add(shape)
        require("custom_materials")
        mat = custom_materials.new()
        mat:pewter()
        plot:add_static(mat.material)
        plot:animate()

require("iso_surface")

function gumdrop_torus(step)
    assert(step > 0)
    assert(step < 1)

    local function sfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        return 4*(x2*x2 + (y2 + z2)*(y2 + z2)) +
               17*x2*(y2 + z2) -
               20*(x2 + y2 + z2) + 17
    end
    
    local function nfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        local nx = 16*x*x2 + 34*x*(y2 + z2) - 40*x
        local ny = 16*y*(y2 + z2) + 34*x2*y - 40*y
        local nz = 16*z*(y2 + z2) + 34*x2*z - 40*z
        local d = math.sqrt(nx*nx + ny*ny + nz*nz)
        return -nx/d, -ny/d, -nz/d
    end
    
    local xarr, yarr, zarr = zeUtl.new("double", "double", "double")
    local n = 2 / step
    xarr:range(-2, step, 2*n+1)
    yarr:range(-2, step, 2*n+1)
    zarr:range(-2, step, 2*n+1)
    
    local xyz, nor = iso_surface(0, xarr, yarr, zarr, sfunc)
    n = xyz:size()-1
    nor:clear()
    for k = 0, n do
        nor:add(nfunc(xyz:get(k)))
    end
    
    return xyz, nor
end