Craft tube

Here’s an example that shows a tube created with Craft. It’s a tube that runs along a Bézier curve. Each run creates a different curve. Press the restart icon for a different curve. You can zoom in from an end or just pass thru the skin to see the inside. You can change the diameter or the number of sides. 2 sides creates a ribbon, and 3 sides just looks interesting.

EDIT: Made a change so the red dots line up with the skin. There was 1 dot that stuck out on one end and on the other end it was one dot in.

displayMode(FULLSCREEN)

function setup()
    diameter=25
    sides=20
    assert(craft, "Please include Craft as a dependency")
    assert(OrbitViewer, "Please include Cameras (not Camera) as a dependency")        
    scene = craft.scene()
    v=scene.camera:add(OrbitViewer,vec3(0,0,0), 600, 0, 2000)
    v.camera.farPlane=3000
    pos,ind,col={},{},{}
    x1,y1,z1=math.random(-300,300),math.random(-300,300),math.random(-300,300)
    x2,y2,z2=math.random(-300,300),math.random(-300,300),math.random(-300,300)
    bezier()
    createSkin()
    createTube()
end

function draw()
    update(DeltaTime)
    scene:draw() 
    text("Slide your finger around the screen to rotate the image.",WIDTH/2,HEIGHT-25)
    text("Use two fingers to zoom in, zoom out or to move the image.",WIDTH/2,HEIGHT-50)
end

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

function bezier()
    local tab1={}
    local x,y,z={0,0},{0,0},{-200,200}
    for t=0,1,.02 do
        local xt=(1-t)^3*x[1]+3*t*(1-t)^2*x1+3*t^2*(1-t)*x2+t^3*x[2]
        local yt=(1-t)^3*y[1]+3*t*(1-t)^2*y1+3*t^2*(1-t)*y2+t^3*y[2]
        local zt=(1-t)^3*z[1]+3*t*(1-t)^2*z1+3*t^2*(1-t)*z2+t^3*z[2]           
        table.insert(tab1,vec3(xt,yt,zt))
    end
    for z=1,#tab1-1 do
        local p1=vec3(tab1[z].x,tab1[z].y,tab1[z].z)
        local p2=vec3(tab1[z+1].x,tab1[z+1].y,tab1[z+1].z)
        createLoop(p1,p2)
        createSphere(vec3(tab1[z].x,tab1[z].y,tab1[z].z))
    end    
end

function createLoop(p1,p2)
    local rp=vec3(80,80,80)
    local v1=rp-p1
    local r1=v1:cross(p2-p1)
    local s1=r1:cross(p2-p1)
    r1,s1=r1:normalize(),s1:normalize()    
    local n=vec3(0,0,0)
    for a=0,359,360/sides do
        n.x = r1.x * math.cos(math.rad(a)) + s1.x * math.sin(math.rad(a))
        n.y = r1.y * math.cos(math.rad(a)) + s1.y * math.sin(math.rad(a))
        n.z = r1.z * math.cos(math.rad(a)) + s1.z * math.sin(math.rad(a))
        n=n:normalize()*diameter
        table.insert(pos,vec3(n.x+p1.x,n.y+p1.y,n.z+p1.z))
    end
end

function createSkin()
    o,p={1,2,3,4,5,6,3,2,1,6,5,4},{}
    for z=1,#pos-sides do
        p[1],p[2],p[3],p[4],p[5],p[6]=z,z+1,z+sides+1,z,z+sides+1,z+sides
        if z%sides==0 then
            p[2]=z-sides+1
            p[3]=z+1
            p[5]=z+1
        end
        for t=1,12 do
            table.insert(ind,p[o[t]])
        end
    end 
    for z=1,#pos do
        table.insert(col,color(math.random(255),math.random(255),math.random(255)))
    end
end

function createSphere(p)
    local pt=scene:entity()
    pt.position=p
    pt.model = craft.model.icosphere(2,2)
    pt.material = craft.material("Materials:Specular")
    pt.material.diffuse=color(255,0,0)
end

function createTube()
    local pt=scene:entity()
    pt.model = craft.model.cube()
    pt.model.positions=pos
    pt.model.indices=ind
    pt.model.colors=col
    pt.material = craft.material("Materials:Basic")    
end

@dave1707 very neat, given me a few ideas. Thanks.

I’m pleased to see that you are beginning to work with vectors! I’ve cleaned up your code a little bit to simplify the vector work. In a few places, you basically break a vector into its components, do some arithmetic with those, and then put it back together again. In those cases, you can just do the arithmetic with the vectors themselves.

My programming style has fewer globals, but I left yours alone for the main part. The only ones I shifted were the control points for the bezier. It seemed odd that you define the control points in setup but the end points in bezier.

The function createTube had me baffled for a bit as you use a cube model. It took me a few moments to realise that you throw away all the cube information, in which case why not just use a blank model?

Your createLoop function will fail if the curve aligns with your reference vector (rp - incidentally, there’s no need for this to have components 80, simply 1 will do as all you care about is its direction not its magnitude). The cross products will then vanish. One of the reasons why my code (in this gist ) is complicated is because it is robust against that situation. The normal directions in my code adjust as they go around the path rather than being set from outside.

Quite by chance, in one run I produced a tube showing the issues that this can produce. The issue here is that the loops get out of sync so one loop starts from a radically different place to the next one. I’ve attached an image to this post (I changed the colours to white to better see the structure).

I also added normals, and changed the material to Standard instead of Basic. The difference isn’t so obvious with the random colours, but if you set them to white (or some other fixed colour) then you get a nice surface effect.

-- DaveBezierTubes

displayMode(FULLSCREEN)

function setup()
    diameter=25
    sides=20
    assert(craft, "Please include Craft as a dependency")
    assert(OrbitViewer, "Please include Cameras (not Camera) as a dependency")        
    scene = craft.scene()
    v=scene.camera:add(OrbitViewer,vec3(0,0,0), 600, 0, 2000)
    v.camera.farPlane=3000
    pos,ind,col,nor={},{},{},{}
    
    bezier()
    createSkin()
    createTube()
end

function draw()
    update(DeltaTime)
    scene:draw() 
    text("Slide your finger around the screen to rotate the image.",WIDTH/2,HEIGHT-25)
    text("Use two fingers to zoom in, zoom out or to move the image.",WIDTH/2,HEIGHT-50)
end

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

function bezier()
    local tab1={}
    local x0 = vec3(0,0,-200)
    local x1 = vec3(math.random(-300,300),math.random(-300,300),math.random(-300,300))
    local x2 = vec3(math.random(-300,300),math.random(-300,300),math.random(-300,300))
    local x3 = vec3(0,0,200)
    local xt
    for t=0,1,.02 do
        xt = (1-t)^3*x0 + 3*t*(1-t)^2*x1 + 3*t^2*(1-t)*x2 + t^3*x3         
        table.insert(tab1,xt)
    end
    for z=1,#tab1-1 do
        createLoop(tab1[z],tab1[z+1])
        createSphere(tab1[z])
    end    
end

function createLoop(p1,p2)
    local rp=vec3(1,1,1)
    local v1=rp-p1
    local r1=v1:cross(p2-p1)
    local s1=r1:cross(p2-p1)
    r1,s1=r1:normalize(),s1:normalize()    
    local n
    for a=0,359,360/sides do
        n = r1 * math.cos(math.rad(a)) + s1 * math.sin(math.rad(a))
        n=n*diameter
        table.insert(pos,n + p1)
        table.insert(nor,n)
    end
end

function createSkin()
    o,p={1,2,3,4,5,6,3,2,1,6,5,4},{}
    for z=1,#pos-sides do
        p[1],p[2],p[3],p[4],p[5],p[6]=z,z+1,z+sides+1,z,z+sides+1,z+sides
        if z%sides==0 then
            p[2]=z-sides+1
            p[3]=z+1
            p[5]=z+1
        end
        for t=1,12 do
            table.insert(ind,p[o[t]])
        end
    end 
    for z=1,#pos do
        table.insert(col,color(math.random(255),math.random(255),math.random(255)))
    end
end

function createSphere(p)
    local pt=scene:entity()
    pt.position=p
    pt.model = craft.model.icosphere(2,2)
    pt.material = craft.material("Materials:Specular")
    pt.material.diffuse=color(255,0,0)
end

function createTube()
    local pt=scene:entity()
    pt.model = craft.model()
    pt.model.positions=pos
    pt.model.indices=ind
    pt.model.colors=col
    pt.model.normals=nor
    pt.material = craft.material("Materials:Standard")    
end

@LoopSpace I’ve seen the twists that you show and figured they were from the rp value. I used to have rp created with random values which created a lot more twists so I set it to 80to get it out of the way. I wasn’t sure what it’s purpose was for since it didn’t change. I can’t take all the credit for the code in createLoop. I still don’t understand why it works. I kept searching on google for code to create a plane at right angles to a line. I found a line by line explanation and coded what they said. I was surprised it actually worked after a lot of trial and error coding. I used cube in createTube because I did that in another program and it seemed to work, so I kept it. A lot of my coding is by habit so I tend to use n.x, n.y, n.z instead of just n. While you’re giving me coding hints, what are normals and uvs and what are they used for. I see you added normal to createTube.

@LoopSpace, @dave1707 much appreciated. i have an interesting app idea that will make use of this-but it will take some time to implement. will report back if it ever works!

@LoopSpace I looked up normals and uvs on google and I see what they are. The uvs are texture mapping and I used them with meshes/texture, but never referred to them as uvs. Normals are for lighting/shadows and I never included them in anything.

Okay, here’s a simplified version of my code with no external dependencies and a bit closer to @dave1707’s code above.

The key part is the answer to this part of @dave1707’s post above:

I kept searching on google for code to create a plane at right angles to a line.

The code below contains such code. It’s easy to define something that almost always works, but getting something that is guaranteed to work needs a little more care.

-- SimpleBezierTubes

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,vec3(0,0,0), 600, 0, 2000)

    -- This entity will contain our tube
    local track = scene:entity()
    -- Some local variables for the calculations
    local np,nt,nq,ax,ang
    -- Start off with the initial position on the bézier curve
    local p = f(0)
    -- This calculates the tangent vector at the initial position
    -- That is, a vector that points along the curve
    local t = (f(.001) - f(-.001)):normalize()
    -- We need a vector orthogonal to the tangent vector
    -- There isn't a canonical way to get such a vector, so the following
    -- picks one of the axis directions and constructs an orthogonal
    -- vector using that.  It picks the best choice for numerical stability
    -- of the following calculation.
    local nml
    if math.abs(t.x) < math.abs(t.y) and math.abs(t.x) < math.abs(t.z) then
        nml = vec3(1,0,0)
    elseif math.abs(t.y) < math.abs(t.z) then
        nml = vec3(0,1,0)
    else
        nml = vec3(0,0,1)
    end
    -- Once we've picked an appropriate axis direction, we use the cross
    -- product to create a vector orthogonal to the tangent vector.
    -- This will be our initial normal vector.
    nml = t:cross(nml):normalize()

    -- We store the position, tangent vector, and normal vector
    local pts = {{p,t,nml}}
    -- We set the number of segments along the tub, the number of sides
    -- around the tube, and radius of the tube.
    local nseg = 50
    local sides = 30
    local radius = 30
    -- Now we add the points, tangent vectors, and normal vectors along
    -- the curve.
    for k=1,nseg do
        -- This is the next point on the curve
        np = f(k/nseg)
        -- This is the (normalised) tangent vector at that point
        nt = (f(k/nseg+.001) - f(k/nseg-.001)):normalize()
        -- Now we compare the old tangent vector and the new one.
        -- We can rotate space to get the old one to the new one, with
        -- axis orthogonal to both of them, this is the axis:
        ax = t:cross(nt)
        -- The length of the axis is the sine of the rotation angle, so
        -- this gets us that angle
        ang = math.asin(ax:len())
        -- Now we apply that rotation to our current normal vector which
        -- gives us the new normal vector
        nml = nml:rotate(ang,ax)
        -- Store the point, tangent vector, and normal vector in the table
        table.insert(pts,{np,nt,nml})
        -- Remember the tangent vector for the next point on the curve
        t = nt
    end

    -- Now that we have our points, we set our model
    local pos, nor, col, ind = {}, {}, {}, {}

    -- Each point is a point along the curve, we need to add a disc of
    -- points for each one.
    for k,v in ipairs(pts) do
        -- Our disc is actually a polygon with "sides" number of sides.
        -- To get its vertices, we rotate the normal vector around the
        -- tangent vector by the appropriate angle.
        for l=1,sides do
            -- To get the actual vertex, we then add that to the position
            table.insert(pos, v[1] + radius*v[3]:rotate(l*2*math.pi/sides,v[2]))
            -- The rotated vector is the normal vector to the tube at that point
            table.insert(nor, v[3]:rotate(l*2*math.pi/sides,v[2]))
            -- In case we want to mess with the colours ...
            table.insert(col,color(255,255,255))
        end
    end
    -- One we have our points in place, we need to knit them together
    -- into triangles.  Each pair of points on a polygonal disc gets joined
    -- in a quadrilateral with the corresponding pair of points on the next
    -- polygonal disc
    for k=1,nseg do
        for l=1,sides do
            table.insert(ind,(k-1)*sides+l)
            table.insert(ind,k*sides+l)
            table.insert(ind,(k-1)*sides+l%sides+1)
            
            table.insert(ind,(k-1)*sides+l%sides+1)
            table.insert(ind,k*sides+l)
            table.insert(ind,k*sides+l%sides+1)

        end
    end
    -- Just for the fun of it, we add another tube of a very slightly smaller
    -- radius and with the triangles oriented the other way around.
    -- This means that when we look inside the tube, we can see it
    -- (Otherwise face culling would mean the far side of the tube would
    -- vanish)
    radius = radius - 1
    for k,v in ipairs(pts) do
        for l=1,sides do
            table.insert(pos, v[1] + radius*v[3]:rotate(l*2*math.pi/sides,v[2]))
            table.insert(nor, -v[3]:rotate(l*2*math.pi/sides,v[2]))
            table.insert(col,color(255,255,255))
        end
    end
    for k=1,nseg do
        for l=1,sides do
            table.insert(ind,(k-1)*sides+l)
            table.insert(ind,(k-1)*sides+l%sides+1)
            table.insert(ind,k*sides+l)
            
            table.insert(ind,(k-1)*sides+l%sides+1)
            table.insert(ind,k*sides+l%sides+1)
            table.insert(ind,k*sides+l)
        end
    end
    -- Now that we have all our arrays, add them to our entity.
    track.model = craft.model()
    track.model.positions = pos
    track.model.normals = nor
    track.model.colors = col
    track.model.indices = ind
    track.material = craft.material("Materials:Standard")
    track.material.diffuse = color(0, 255, 223, 255)
    -- Since @dave1707's code added spheres along the path, so shall we
    local sph
    for k,v in ipairs(pts) do
        sph = scene:entity()
        sph.position = v[1]
        sph.model = craft.model.icosphere(2,2)
        sph.material = craft.material("Materials:Specular")
        sph.material.diffuse = color(255, 0, 0, 255)
    end
end

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

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

-- This is a small extension to the vec3 code which defines a rotation of
-- a vec3 about a given axis.
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

-- This is our bézier curve function.  The end points are fixed but the
-- control points are chosen to be in random directions.
local a,b,c,d
a = vec3(0,0,0)
b = vec3(0,0,100):rotate(math.random()*math.pi/2,vec3(1,0,0):rotate(math.random()*2*math.pi,vec3(0,0,1)))
d = vec3(0,0,300)
c =d + vec3(0,0,-100):rotate(math.random()*math.pi/2,vec3(1,0,0):rotate(math.random()*2*math.pi,vec3(0,0,1)))

function f(t)
    return (1-t)^3*a + 3*(1-t)^2*t*b + 3*(1-t)*t^2*c + t^3*d
end

@LoopSpace I like all the comments you gave, helps me understand somewhat as to what’s happening. Without going back and commenting my code, I’ll give an explanation of what I’m doing in my code. The bezier() function just calculates points along a bezier curve. I put the x,y,z values of each calculated point into tab1 and also plot the red spheres along the curve. I then iterate thru tab1 using the current point and the next point and call createLoop() which calculates a plane perpendicular to those two points. After getting the plane, I then calculate the x,y,z points around a circle in that plane based on the number of sides. I save all those points in the table pos. So table pos contains the points that make up all the circles perpendicular to the Bézier line at each red circle. I then call createSkin() which iterates thru the table pos which creates triangles that connect the points of the first circle with the points of the second circle. Then it connects the second circle with the third and so on until I connected all the points. I then call createTube() which draws the tube based on all the triangles that were calculated. I don’t fully understand how the dot product stuff in createLoop() actually calculates a plane perpendicular to two points, but the more I play around with them, the more I’ll understand it.

@LoopSpace Here’s the link I found that showed me how to calculate a plane perpendicular to a line. Going thru the Procedure section took a lot of trial and error before I got my program to work since I didn’t understand what was going on. Actually, I still don’t, but I will eventually. Looking at the Example program there also helped.

http://homer.com.au/webdoc/geometry/disk.htm

@LoopSpace Something I forgot to mention when I was reading your comments above. You said you created another tube of radius - 1 for your reverse culling. On mine, I did the reverse using the same tube coordinates. That seems to work. Maybe if the outside and the inside were different colors or textures, it wouldn’t.

@LoopSpace i have the suspicion that if the tangent of the previous segment and next segment are the same (i.e. straight line) then your algorithm for the normal fails.

@piinthesky @LoopSpace I’ve noticed that in my tubes where there’s a twist here and there. Since I don’t fully understand what’s happening, I’m wondering if the point used for the calculation of the next tangent has to be at a certain location relative to the previous tangent.

@piinthesky It fails on a technicality: if you pass a zero vector to the mt:rotate as the axis then because of the normalisation step you get (nan, nan, nan) back out again. To guard against this, just modify the function to:

local mt = getmetatable(vec3())
mt.rotate = function(v,a,u)
    if u:len() == 0 then
        return v
    end
    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

So there’s nothing wrong with the algorithm per se.

In my original code, I used quaternions here (because the VecExt makes it possible to easily apply a quaternion to a vec3, which isn’t easy in vanilla Codea) but in making a self-contained version for this thread, I took that out and put in the rotate method instead. I missed the fact that the cross product can be zero and so mess up the rotate function, but what happens then is that the rotation should be by 0 degrees so it’s easy to code for that case.

@dave1707 (Replying to an earlier comment than the most recent) I tried with a cylinder with reversed triangles and it works just fine, so the radius - 1 stuff is overkill. Thanks for pointing that out.

@dave1707 (Replying to the most recent comment) No, it’s not about how far along the curve you calculate. It’s to do with the fact that it’s not enough to calculate the normal plane, you also need a starting point on that plane to start your cylinder from. So you need to keep track of that as well, and that’s what the nml vector does.

@LoopSpace, perfect-that fixed it, thanks.

This might be a useful function for this type of thing. It produces a vector guaranteed to be at 90 degrees to its input. It can be mathematically shown that any such assignment has to have “jumps”, and so this does have those.

It works by picking the axis vector which is furthest from the given vector (the jumps occur when there are two equi-distant) and then taking the cross product with the original vector. By taking the furthest, we reduce the potential for numerical inaccuracy with small numbers.

local floor = math.floor
local abs = math.abs
function Heaviside(t)
    local s = floor(t)+.5
    return s/abs(s)/2+.5
end

function orthogonal(v)
    local x,y,z = abs(v.x),abs(v.y),abs(v.z)
    local xy = Heaviside(x-y) -- 1 if x >= y
    local yz = Heaviside(y-z) -- 1 if y >= z
    local zx = Heaviside(z-x) -- 1 if z >= x
    return v:cross(vec3(
        zx*(1-xy) + xy*yz*zx,
        xy*(1-yz),
        yz*(1-zx)
    ))
end

It also avoids explicit conditionals by using a numerical method: the Heaviside function is 0 for negative numbers and 1 for positive (including 0). So Heaviside(x-y) is 1 if x-y is positive, which is the same as saying that x >= y. Thus xy*(1-yz) is 1 if x >= y and y < z (the 1- flips the inequality); ie, if y is the smallest component (if both x and y are smallest, then y “wins” - the extra term in the first component means that if all three are the same then x wins).

Couple of videos based on stuff from this thread:

  1. https://youtu.be/NNpcx6-HBgI
  2. https://youtu.be/A15Hu6FlF5c