Codea Tube

-- BezierTubes

function setup()
    assert(craft, "Please include Craft as a dependency")
    assert(OrbitViewer, "Please include Cameras (not Camera) as a dependency")

    scene = craft.scene()
    scene.camera:add(OrbitViewer)
    local track = scene:entity()
    local p = f(0)
    local np,nt,nq,ax,ang

    local t = (f(.001) - f(-.001))/0.002
    local it = t
    local nml = vec3(1,0,0)
    local pts = {{p,it,nml}}
    for k=1,100 do
        np = f(k/100*2*math.pi)
        nt = (f(k/100*2*math.pi+.001) - f(k/100*2*math.pi-.001))/0.002
        ax = t:cross(nt)
        ang = math.asin(ax:len()/(t:len()*nt:len()))
        nml = nml:rotate(ang,ax)

        table.insert(pts,{np,nt,nml})
        t = nt
        p = np
    end
    local __track = PseudoMesh()
    local ends = 1
    for k=1,100 do
        if k == 100 then
            ends = 2
        end
        __track:addCylinder({
            startCentre = pts[k][1],
            startWidth = pts[k][3],
            startHeight = pts[k][2]:cross(pts[k][3]):normalize(),
            endCentre = pts[k+1][1],
            endWidth = pts[k+1][3],
            endHeight = pts[k+1][2]:cross(pts[k+1][3]):normalize(),
            ends = ends,
            faceted = false,
            size=50
        })
        ends = 0
    end
    track.model = __track:toModel()
    track.material = craft.material("Materials:Standard")
end

function draw()
    background(40,40,50)
    update(DeltaTime)
    scene:draw()
end

function update(dt)
    scene:update(dt)
end

local mt = getmetatable(vec3())
mt.rotate = function(v,a,u)
    u = u or vec3(1,0,0)
    u = u:normalize()
    local x = v - v:dot(u)*u
    local y = u:cross(x)
    return v - x + math.cos(a)*x + math.sin(a)*y
end

function f(t)
    local r,h=5,5
    return vec3(r*math.cos(t),r*math.sin(t),h*t)
end

PseudoMesh = class()

function PseudoMesh:init()
    self.vertices = {}
    self.texCoords = {}
    self.normals = {}
    self.colors = {}
    self.size = 0
end

function PseudoMesh:vertex(k,v)
    if v then
        self.vertices[k] = v
    else
        return self.vertices[k]
    end
end

function PseudoMesh:normal(k,v)
    if v then
        self.normals[k] = v
    else
        return self.normals[k]
    end
end

function PseudoMesh:texCoord(k,v)
    if v then
        self.texCoords[k] = v
    else
        return self.texCoords[k]
    end
end

function PseudoMesh:color(k,v)
    if v then
        self.colors[k] = v
    else
        return self.colors[k]
    end
end

function PseudoMesh:resize(k)
    self.size = k
end

function PseudoMesh:addCylinder(t)
    t = t or {}
    local nt = {}
    for k,v in pairs(t) do
        nt[k] = v
    end
    nt.mesh = self
    return addCylinder(nt)
end

function PseudoMesh:invertNormals()
    for k,v in ipairs(self.normals) do
        self.normals[k] = -v
    end
end

function PseudoMesh:toModel()
    local m = craft.model()
    local i = {}
    local n = #self.vertices
    for k=1,n,3 do
        if (self.vertices[k+1] - self.vertices[k]):cross(self.vertices[k+2] - self.vertices[k]):dot(self.normals[k+1] + self.normals[k+2] + self.normals[k]) < 0 then
            table.insert(i,k)
            table.insert(i,k+1)
            table.insert(i,k+2)
        else
            table.insert(i,k)
            table.insert(i,k+2)
            table.insert(i,k+1)
        end
    end
    m.positions = self.vertices
    m.normals = self.normals
    m.uvs = self.texCoords
    m.colors = self.colors
    m.indices = i
    return m
end

-- MeshExt

local __doJewel, __doSuspension, __doPyramid, __doBlock, __addTriangle, __doSphere, __threeFromTwo, __orthogonalTo, __doCylinder, __discreteNormal, __doCone, __doPoly, __doFacetedClosedCone, __doFacetedOpenCone, __doSmoothClosedCone, __doSmoothOpenCone, __doFacetedClosedCylinder, __doFacetedOpenCylinder, __doSmoothClosedCylinder, __doSmoothOpenCylinder, __initmesh



--[[
A "cone" is formed by taking a curve in space and joining each of its points to an apex.
If the original curve is made from line segments, the resulting cone has a natural triangulation which can be used to construct it as a mesh.

m mesh
p position in mesh
n number of points
o "internal" point (to ensure that normals point outwards)
a apex of cone
v table of base points
t table of texture points
col colour
to texture offset
ts not used
f faceted or not
cl closed curve or not
l light vector
--]]
function __doCone(m,p,n,o,a,v,t,col,to,ts,f,cl,l,am)
    if f then
        if cl then
            return __doFacetedClosedCone(m,p,n,o,a,v,t,col,to,ts,l,am)
        else
            return __doFacetedOpenCone(m,p,n,o,a,v,t,col,to,ts,l,am)
        end
    else
        if cl then
            return __doSmoothClosedCone(m,p,n,o,a,v,t,col,to,ts,l,am)
        else
            return __doSmoothOpenCone(m,p,n,o,a,v,t,col,to,ts,l,am)
        end
    end
end

function __doFacetedClosedCone(m,p,n,o,a,v,t,col,to,ts,l,am)
    local j,nml,c
    for k=1,n do
        j = k%n + 1
        nml = (v[k] - a):cross(v[j] - a)
        if nml:dot(a - o) < 0 then
            nml = -nml
        end
        nml = nml:normalize()
        c = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nml)))
        __addTriangle(m,p,v[j],v[k],a,c,c,c,nml,nml,nml,to+t[j],to+t[k],to)
        p = p + 3
    end
    return p
end

function __doFacetedOpenCone(m,p,n,o,a,v,t,col,to,ts,l,am)
    local j,nml,c
    for k=1,n-1 do
        j = k + 1
        nml = (v[k] - a):cross(v[j] - a)
        if nml:dot(a - o) < 0 then
            nml = -nml
        end
        nml = nml:normalize()
        c = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nml)))
        __addTriangle(m,p,v[j],v[k],a,c,c,c,nml,nml,nml,to+t[j],to+t[k],to)
        p = p + 3
    end
    return p
end

function __doSmoothClosedCone(m,p,n,o,a,v,t,col,to,ts,l,am)
    local j,nmla,nmlb,nmlc,cc,ca,nb,cb
    nmlb = vec3(0,0,0)
    nmlc = __discreteNormal(v[1],o,v[n],a,v[2])
    cc = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nmlc)))
    nb = vec3(0,0,0)
    for k=1,n do
        j = k%n + 1
        nb = nb + __discreteNormal(v[j],o,v[k],a,v[j%n+1])
    end
    nb = nb:normalize()
    cb = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nb)))
    for k=1,n do
        j = k%n + 1
        nmla = nmlc
        ca = cc
        nmlc = __discreteNormal(v[j],o,v[k],a,v[j%n+1])
        cc = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nmlc)))
        __addTriangle(m,p,v[j],v[k],a,cc,ca,cb,nmlc,nmla,nmlb,to+t[j],to+t[k],to)
        p = p + 3
    end
    return p
end

function __doSmoothOpenCone(m,p,n,o,a,v,t,col,to,ts,l,am)
    local j,nmla,nmlb,nmlc,cc,ca,ll,nb,cb
    ll = l:len()
    nmlb = vec3(0,0,0)
    nmlc = __discreteNormal(v[1],o,a,v[2])
    cc = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nmlc)))
    nb = vec3(0,0,0)
    for k=1,n-2 do
        j = k + 1
        nb = nb + __discreteNormal(v[j],o,v[k],a,v[j%n+1])
    end
    nb = nb + __discreteNormal(v[n],o,v[n-1],a)
    nb = nb:normalize()
    cb = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nb)))
    for k=1,n-2 do
        j = k + 1
        nmla = nmlc
        ca = cc
        nmlc = __discreteNormal(v[j],o,v[k],a,v[j%n+1])
        cc = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nmlc)))
        __addTriangle(m,p,v[j],v[k],a,cc,ca,cb,nmlc,nmla,nmlb,to+t[j],to+t[k],to)
        p = p + 3
    end
    nmla = nmlc
    ca = cc
    nmlc = __discreteNormal(v[n],o,v[n-1],a)
    cc = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nmlc)))
    __addTriangle(m,p,v[n],v[n-1],a,cc,ca,cb,nmlc,nmla,nmlb,to+t[n],to+t[n-1],to)
    return p + 3
end


--[[
This forms a surface which has boundary a given curve by forming a cone with the barycentre of the curve as its apex.

m mesh
p position in mesh
n number of points
o "internal" point (for normals)
v table of base points
t table of texture points
col colour
to texture offset
ts not used
f faceted or not
cl closed curve or not
l light vector
--]]
function __doPoly(m,p,n,o,v,t,col,to,ts,f,cl,l,am)
    local a,b,r = vec3(0,0,0),vec2(0,0),0
    for k,u in ipairs(v) do
        a = a + u
        r = r + 1
    end
    a = a / r
    for k,u in ipairs(t) do
        b = b + u
    end
    b = b / r
    for k=1,r do
        t[k] = t[k] - b
    end
    return __doCone(m,p,n,o,a,v,t,col,to+b,ts,f,cl,l,am)
end

--[[
| Option | Default | Description |
|:-------|:--------|:------------|
| `mesh` | new mesh | The mesh to add the shape to. |
| `position` | end of mesh | The place in the mesh at which to add the shape. |
| `colour`/`color` | white | The colour of the shape. |
| `faceted` | true | Whether to make it faceted or smooth. |
| `ends` | `0` | Which ends to fill in (`0` for none, `1` for start, `2` for end, `3` for both) |
| `texOrigin`    | `vec2(0,0)`              | If using a sprite sheet, this is the lower left corner of the rectangle associated with this shape. |
| `texSize`      | `vec2(1,1)`              | This is the width and height of the rectangle of the texture associated to this shape. |

There are various ways to specify the dimensions of the cylinder.
If given together, the more specific overrides the more general.

`radius` and `height` (`number`s) can be combined with `axes` (table of three `vec3`s) to specify the dimensions, where the first axis vector lies along the cylinder.  The vector `origin` then locates the cylinder in space.

`startCentre`/`startCenter` (a `vec3`), `startWidth` (`number` or `vec3`), `startHeight` (`number` or `vec3`), `startRadius` (`number`) specify the dimensions at the start of the cylinder (if numbers, they are taken with respect to certain axes).

Similarly named options control the other end.

If axes are needed, these can be supplied via the `axes` option.
If just the `axis` option is given (a single `vec3`), this is the direction along the cylinder.
Other directions (if needed) are found by taking orthogonal vectors to this axis.
--]]
										
function addCylinder(t)
    t = t or {}
    local m = t.mesh
    local p = t.position or (m.size + 1)
    p = p + (1-p)%3
    local ip = p
    local col = t.colour or t.color or color(255, 255, 255, 255)
    local f = true
    local ends
    local solid = t.solid
    if solid then
        ends = t.ends or 3
    else
        ends = t.ends or 0
    end
    if t.faceted ~= nil then
        f = t.faceted
    end
    local l,am
    if rl then
        l = vec3(0,0,0)
        am = 1
    else
        l = t.light or vec3(0,0,0)
        if t.intensity then
            l = l:normalize()*t.intensity
        elseif l:lenSqr() > 1 then
            l = l:normalize()
        end
        am = t.ambience or (1 - l:len())
    end
    local r = t.radius or 1
    local h = t.height or 1
    local to = t.texOrigin or vec2(0,0)
    local ts = t.texSize or vec2(1,1)
    local sc,si,sj,ec,ei,ej,a,o
    
    if t.axis or t.axes or t.origin or t.centre or t.center then
        if t.axis then
            a = t.axis
        elseif t.axes then
            a = t.axes[1]
        else
            a = vec3(0,1,0)
        end
        if t.height then
            a = h*a:normalize()
        end
        if t.origin or t.centre or t.center then
            local o = t.origin or t.centre or t.center
            sc,ec = o - a/2,o + a/2
        end
    end
    sc = t.startCentre or t.startCenter or sc
    ec = t.endCentre or t.endCenter or ec
    sc,ec,a = __threeFromTwo(sc,ec,a,vec3(0,-h/2,0),vec3(0,h/2,0),vec3(0,h,0))
    si = t.startWidth or t.startRadius or t.radius or 1
    sj = t.startHeight or t.startRadius or t.radius or 1
    ei = t.endWidth or t.endRadius or t.radius or 1
    ej = t.endHeight or t.endRadius or t.radius or 1
    local c,d
    if t.axes then
        a,c,d = unpack(t.axes)
    end
    if type(si) == "number" then
        if type(sj) == "number" then
            if not c then
                c = __orthogonalTo(a)
            end
            si = si*c:normalize()
            if not d then
                sj = sj*a:cross(c):normalize()
            else
                sj = sj*d:normalize()
            end
        else
            si = si*sj:cross(a):normalize()
        end
    elseif type(sj) == "number" then
        sj = sj*a:cross(si):normalize()
    end
    if type(ei) == "number" then
        if type(ej) == "number" then
            if not c then
                c = __orthogonalTo(a)
            end
            ei = ei*c:normalize()
            if not d then
                ej = ej*a:cross(c):normalize()
            else
                ej = ej*d:normalize()
            end
        else
            ei = ei*ej:cross(a):normalize()
        end
    elseif type(ej) == "number" then
        ej = ej*a:cross(ei):normalize()
    end

    local n = t.size or 12
    local sa,ea,da = __threeFromTwo(t.startAngle,t.endAngle,t.deltaAngle,0,360,360)
    local closed
    if da == 360 then
        closed = true
        solid = false
    else
        closed = false
    end
    sa = math.rad(sa)
    ea = math.rad(ea)
    da = math.rad(da)/n
    o = (sc + math.cos((sa+ea)/2)*si/2 + math.sin((sa+ea)/2)*sj/2 + ec + math.cos((sa+ea)/2)*ei/2 + math.sin((sa+ea)/2)*ej/2)/2
    local ss = 1 + math.floor((ends+1)/2)
    if solid then
        ss = ss + 2 
    end
    ts.x = ts.x / ss
    local cs,sn,ti,tj
    ti,tj = vec2(ts.x/2,0),vec2(0,ts.y/2)
    cs = math.cos(sa)
    sn = math.sin(sa)
    si,sj = cs*si + sn*sj, -sn*si + cs*sj
    ei,ej = cs*ei + sn*ej, -sn*ei + cs*ej
    ti,tj = cs*ti + sn*tj, -sn*ti + cs*tj
    local u,v,tu,tv,tw,cnrs = {},{},{},{},{},{}
    cnrs[1] = {sc,sc+si,ec+ei,ec}
    cs = math.cos(da)
    sn = math.sin(da)
    u[0] = sc+cs*si - sn*sj
    v[0] = ec+cs*ei - sn*ej
    for k=0,n do
        table.insert(u,sc+si)
        table.insert(v,ec+ei)
        table.insert(tu,to + vec2(ts.x*k/n,0))
        table.insert(tv,to + vec2(ts.x*k/n,ts.y))
        table.insert(tw,ti)
        si,sj = cs*si + sn*sj, -sn*si + cs*sj
        ei,ej = cs*ei + sn*ej, -sn*ei + cs*ej
        ti,tj = cs*ti + sn*tj, -sn*ti + cs*tj
    end
    u[n+2] = sc+cs*si + sn*sj
    v[n+2] = ec+cs*ei + sn*ej
    cnrs[2] = {sc, sc+cs*si - sn*sj, ec+cs*ei - sn*ej, ec}
    local size = 6*n + math.floor((ends+1)/2)*3*n
    if closed then
        size = size + math.floor((ends+1)/2)*3
    elseif solid then
        size = size + 24
    end
    if p - 1 + size > m.size then
        m:resize(p-1+size)
    end
    n = n + 1
    p = __doCylinder(m,p,n,o,u,v,tu,tv,col,f,closed,l,am)
    to = to + ts/2
    if solid and not closed then
        local tex = {-ts/2,vec2(ts.x/2,-ts.y/2),ts/2,vec2(-ts.x/2,ts.y/2)}
        for i=1,2 do
            to.x = to.x + ts.x
            p = __doPoly(m,p,4,o,cnrs[i],tex,col,to,ts,f,true,l,am)
        end
    end
    if ends%2 == 1 then
        to.x = to.x + ts.x
        p = __doCone(m,p,n,o,sc,u,tw,col,to,ts,f,closed,l,am)
    end
    if ends >= 2 then
        to.x = to.x + ts.x
        p = __doCone(m,p,n,o,ec,v,tw,col,to,ts,f,closed,l,am)
    end
    if ret then
        return m,ip, p
    else
        return m
    end
end

--[[
This adds a cylinder to the mesh.

m mesh to add shape to
p position of first vertex of shape
n number of points
o centre of shape (vec3)
a apexes (table of 2 vec3s)
v vertices (table of vec3s in cyclic order)
t texture coordinates corresponding to vertices (relative to centre)
col colour
to offset of texture region (vec2)
ts size of texture region (vec2)
f faceted
cl closed
l light vector
--]]
function __doCylinder(m,p,n,o,u,v,ut,vt,col,f,cl,l,am)
    if f then
        if cl then
            return __doFacetedClosedCylinder(m,p,n-1,o,u,v,ut,vt,col,l,am)
        else
            return __doFacetedOpenCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
        end
    else
        if cl then
            return __doSmoothClosedCylinder(m,p,n-1,o,u,v,ut,vt,col,l,am)
        else
            return __doSmoothOpenCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
        end
    end
end

function __doFacetedClosedCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
    local i,j,nv,nu,cu,cv,ll
    ll = l:len()
    for k=1,n do
        j = k%n + 1
        nu = (u[j] - u[k]):cross(v[k] - u[k]):normalize()
        nv = (v[j] - v[k]):cross(u[k] - v[k]):normalize()
        if nu:dot(u[k]-o) < 0 then
            nu = -nu
        end
        if nv:dot(v[k]-o) < 0 then
            nv = -nv
        end
        cv = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv)))
        cu = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu)))
        __addTriangle(m,p,v[j],v[k],u[j],cv,cv,cu,nv,nv,nu,vt[j],vt[k],ut[j])
        p = p + 3
        __addTriangle(m,p,v[k],u[j],u[k],cv,cu,cu,nv,nu,nu,vt[k],ut[j],ut[k])
        p = p + 3
    end
    return p
end

function __doFacetedOpenCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
    local i,j,nv,nu,cu,cv,ll
    ll = l:len()
    for k=1,n-1 do
        j = k + 1
        nu = (u[j] - u[k]):cross(v[k] - u[k]):normalize()
        nv = (v[j] - v[k]):cross(u[k] - v[k]):normalize()
        if nu:dot(u[k]-o) < 0 then
            nu = -nu
        end
        if nv:dot(v[k]-o) < 0 then
            nv = -nv
        end
        cv = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv)))
        cu = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu)))
        __addTriangle(m,p,v[j],v[k],u[j],cv,cv,cu,nv,nv,nu,vt[j],vt[k],ut[j])
        p = p + 3
        __addTriangle(m,p,v[k],u[j],u[k],cv,cu,cu,nv,nu,nu,vt[k],ut[j],ut[k])
        p = p + 3
    end
    return p
end

function __doSmoothClosedCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
    local i,j,nv,nu,cu,cv
    nv,nu,cv,cu = {},{},{},{}
    nv[1] = __discreteNormal(v[1],o,v[n],u[1],v[2])
    nu[1] = __discreteNormal(u[1],o,u[n],v[1],u[2])
    cv[1] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv[1])))
    cu[1] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu[1])))
    for k=1,n do
        j = k%n + 1
        i = j%n + 1
        nv[j] = __discreteNormal(v[j],o,v[k],u[j],v[i])
        nu[j] = __discreteNormal(u[j],o,u[k],v[j],u[i])
        cv[j] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv[j])))
        cu[j] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu[j])))
        __addTriangle(m,p,v[j],v[k],u[j],cv[j],cv[k],cu[j],nv[j],nv[k],nu[j],vt[j],vt[k],ut[j])
        p = p + 3
        __addTriangle(m,p,v[k],u[j],u[k],cv[k],cu[j],cu[k],nv[k],nu[j],nu[k],vt[k],ut[j],ut[k])
        p = p + 3
    end
    return p
end

function __doSmoothOpenCylinder(m,p,n,o,u,v,ut,vt,col,l,am)
    local i,j,nv,nu,cu,cv
    nv,nu,cv,cu = {},{},{},{}
    nv[1] = __discreteNormal(v[1],o,v[0],u[1],v[2])
    nu[1] = __discreteNormal(u[1],o,v[0],v[1],u[2])
    cv[1] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv[1])))
    cu[1] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu[1])))
    for k=1,n-2 do
        j = k + 1
        i = j + 1
        nv[j] = __discreteNormal(v[j],o,v[k],u[j],v[i])
        nu[j] = __discreteNormal(u[j],o,u[k],v[j],u[i])
        cv[j] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv[j])))
        cu[j] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu[j])))
        __addTriangle(m,p,v[j],v[k],u[j],cv[j],cv[k],cu[j],nv[j],nv[k],nu[j],vt[j],vt[k],ut[j])
        p = p + 3
        __addTriangle(m,p,v[k],u[j],u[k],cv[k],cu[j],cu[k],nv[k],nu[j],nu[k],vt[k],ut[j],ut[k])
        p = p + 3
    end
    nv[n] = __discreteNormal(v[n],o,v[n-1],u[n],v[n+1])
    nu[n] = __discreteNormal(u[n],o,u[n-1],v[n],u[n+1])
    cv[n] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nv[n])))
    cu[n] = col:mix(color(0,0,0,col.a),am+(1-am)*math.max(0,l:dot(nu[n])))
    __addTriangle(m,p,v[n],v[n-1],u[n],cv[n],cv[n-1],cu[n],nv[n],nv[n-1],nu[n],vt[n],vt[n-1],ut[n])
    p = p + 3
    __addTriangle(m,p,v[n-1],u[n],u[n-1],cv[n-1],cu[n],cu[n-1],nv[n-1],nu[n],nu[n-1],vt[n-1],ut[n],ut[n-1])
    return p+3
end


--[[
This works out a normal vector for a vertex in a triangulated surface by taking an average of the triangles in which it appears.
The normals are weighted by the reciprocal of the size of the corresponding triangle.

a vertex under consideration
o a point to determine which side the normals lie
... a cyclic list of vertices, successive pairs of which make up the triangles
--]]
function __discreteNormal(a,o,...)
    local arg = {}
    local n = 0
    for k,v in ipairs({...}) do
        if v then
            table.insert(arg,v)
            n = n + 1
        end
    end
    local na,nb
    na = vec3(0,0,0)
    for k=2,n do
        nb = (arg[k] - a):cross(arg[k-1] - a)
        na = na + nb/nb:lenSqr()
    end
    na = na:normalize()
    if na:dot(a-o) < 0 then
        na = -na
    end
    return na
end






--[[
Adds a triangle to a mesh, with specific vertices, colours, normals, and texture coordinates.
--]]
function __addTriangle(m,p,a,b,c,d,e,f,g,h,i,j,k,l)
    m:vertex(p,a)
    m:color(p,d)
    m:normal(p,g)
    m:texCoord(p,j)
    p = p + 1
    m:vertex(p,b)
    m:color(p,e)
    m:normal(p,h)
    m:texCoord(p,k)
    p = p + 1
    m:vertex(p,c)
    m:color(p,f)
    m:normal(p,i)
    m:texCoord(p,l)
end

--[[
Returns three things, u,v,w, with the property that u + w = v.  The input can be any number of the three together with three defaults to be used if not enough information is given (so if any two of the first three are given then that is enough information to determine the third).
--]]
function __threeFromTwo(a,b,c,d,e,f)
    local u,v,w = a or d or 0, b or e or 1, c or f or 1
    if not a then
        u = v - w
    end
    if not b then
        v = u + w
    end
    if not c then
        w = v - u
    end
    v = u + w
    return u,v,w
end

--[[
Returns a vector orthogonal to the given `vec3`.
--]]
function __orthogonalTo(v)
    local a,b,c = math.abs(v.x), math.abs(v.y), math.abs(v.z)
    if a < b and a < c then
        return vec3(0,v.z,-v.y)
    end
    if b < c then
        return vec3(-v.z,0,v.x)
    end
    return vec3(v.y,-v.x,0)
end