Tube plots?

I use a lot of 3D line plots in my teaching, and I’ve been trying different visualization packages to make them . Is there anything in the Julia graphics ecosystem that could show the line as a ‘tube’ with a lighting model, to make it easier to see the line’s path in 3D space?

you can plot cylinder meshes in Makie.jl for example, it would look a bit like the 3d arrows on this page, just without the arrow tips arrows · Makie Plotting Ecosystem

2 Likes

Reusing some thing from arrows would be one option. You could use Makie._mantle to get open cylinders/tubes and use meshscatter with fitting rotations and scaling to trace some path. The problem with that is that you end up with gaps (or edges if you elongate the cylinders).

positions = map(range(0, 10pi, length=200)) do t
    r = (1.0 + 0.1t)
    Point3f0(r * cos(t), r * sin(t), 0.5t)
end

# origin, extremity, radius1, radius2, segments
m = Makie._mantle(Point3f0(0), Point3f0(0,0,1), 1f0, 1f0, 16)
scales = map(s -> Vec3f0(0.2f0, 0.2f0, s), norm.(positions[2:end] .- positions[1:end-1]))
rots = normalize(positions[2:end] .- positions[1:end-1])

fig, ax, p = meshscatter(positions[1:end], rotation = rots, markersize=scales, color = :red, marker=m)

A quick fix for that would be adding some spheres, but that also doesn’t look quite right.

# ...
meshscatter!(ax, positions[2:end-1], markersize=0.2f0, color = :red)

To really get clean joints you would have to generate a full mesh from your path. That’s more difficult and I don’t know if my solution is always going to work, but it does look good on my test path:

using LinearAlgebra, GeometryBasics, GLMakie


function generate_circle!(
        vertices, i0, origin, normal, radius, segments, 
        starting_unit_vec = abs(dot(Vec3f0(0,0,1), normal)) > 0.1 ? Vec3f0(0,0,1) : Vec3f0(0,1,0)
    )
    u1 = normalize(cross(starting_unit_vec, normal))
    u2 = normalize(cross(normal, u1))
    dphi = 2pi / segments
    for (i, phi) in enumerate(0:dphi:2pi-0.5dphi)
        vertices[i0 + i - 1] = radius * (u1 * sin(phi) + u2 * cos(phi)) + origin
    end
    return i0 + segments, u2
end


function tube(positions; radius = 0.2, segments = 16)
    vertices = Vector{Point3f0}(undef, length(positions) * segments)
    # This call fills vertices[1:segments] with points on a circle centered
    # at positions[1] with surface normal normalize(positions[2] - positions[1])
    # The output is the next index to write to (segments+1) and one of the
    # generated unit vectors. The unit vector is passed to the next iteration 
    # to keep vertices aligned (I hope).
    i, u = generate_circle!(
        vertices, 1, positions[1], 
        normalize(positions[2] - positions[1]), radius, segments
    )
    for k in 2:length(positions)-1
        p = positions[k]
        n = normalize(positions[k+1] - positions[k-1])
        i, u = generate_circle!(vertices, i, p, n, radius, segments, u)
    end
    i, u = generate_circle!(
        vertices, i, positions[end], 
        normalize(positions[end] - positions[end-1]), radius, segments, u
    )
    # finish mesh
    faces = Vector{GLTriangleFace}(undef, 2 * segments * (length(positions)-1))
    idx = 1
    for k in 0:length(positions)-2
        for j in 1:segments
            # creates faces:
            # circle k    circle k+1
            #     i  ------------  i + segments
            #       | `-.        |
            #       |     `-.    |
            #       |         `-.|
            # ishift ------------ ishift + segments
            i = segments * k + j
            ishift = segments * k + mod1(j+1, segments)
            faces[idx] = GLTriangleFace(i, i+segments, ishift + segments)
            faces[idx+1] = GLTriangleFace(i, ishift + segments, ishift)
            idx += 2
        end
    end
    normal_mesh(vertices, faces)
end

positions = map(range(0, 10pi, length=200)) do t
    r = (1.0 + 0.1t)
    Point3f0(r * cos(t), r * sin(t), 0.5t)
end
fig, ax, p = mesh(tube(positions), color = :blue)

7 Likes

An alternative way is to generate a parametric tube surface along the line path.

using GLMakie, LinearAlgebra
set_theme!(backgroundcolor = :black)

L(s)  = [R*cos(s), R*sin(s), H*s]
∇L(s) = [-R*sin(s), R*cos(s), H]

function C(s,u)
    n1, n2 = eachcol(nullspace(adjoint(∇L(s))))
    return L(s) + r*cos(u) * n1 + r*sin(u) * n2
end

r = 0.5; R = 5.0; H = 1.0   # define constants
s = LinRange(0, 8π, 1000);  u = LinRange(0, 2π, 36)
xs, ys, zs = ([p[i] for p in C.(s, u')] for i in 1:3)

scene = surface(xs, ys, zs; lightposition = Vec3f0(0,0,-100))

NB:
_ Written more compactly with simpler coordinates
_ For more complex functions, use ForwardDiff to compute the gradient

12 Likes

This looks gorgeous! I love the colors!

1 Like

The correct way of finding the gradient’s normal plane (unit vectors) is to compute the nullspace of its adjoint. Edited the code above accordingly, as original was not general enough.

3 Likes

Thanks @rafael.guerra, that looks like just what I want, and the prospect of using AD is nice too. Unfortunately when I run your code I get a method error for svd!, am I missing a package? I’m on v1.6.2 and the full stacktrace is

ERROR: MethodError: no method matching svd!(::Vector{Float64}; full=true, alg=LinearAlgebra.DivideAndConquer())
Closest candidates are:
  svd!(::LinearAlgebra.AbstractTriangular; kwargs...) at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\triangular.jl:2212
  svd!(::StridedMatrix{T}; full, alg) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\svd.jl:92
  svd!(::StridedMatrix{T}, ::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\svd.jl:361 got unsupported keyword arguments "full", "alg"
  ...
Stacktrace:
  [1] svd(A::Vector{Float64}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\svd.jl:157  
  [2] svd(A::Adjoint{Float64, Vector{Float64}}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\svd.jl:166  
  [3] #nullspace#30
    @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\dense.jl:1446 [inlined]   
  [4] nullspace(A::Adjoint{Float64, Vector{Float64}})
    @ LinearAlgebra C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\dense.jl:1444
  [5] C(s::Float64, u::Float64)
    @ Main .\REPL[9]:2
  [6] _broadcast_getindex_evalf
    @ .\broadcast.jl:648 [inlined]
  [7] _broadcast_getindex
    @ .\broadcast.jl:621 [inlined]
  [8] getindex
    @ .\broadcast.jl:575 [inlined]
  [9] copy
    @ .\broadcast.jl:922 [inlined]
 [10] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(C), Tuple{LinRange{Float64}, Adjoint{Float64, LinRange{Float64}}}})
    @ Base.Broadcast .\broadcast.jl:883
 [11] (::var"#29#31")(i::Int64)
    @ Main .\none:0
 [12] iterate
    @ .\generator.jl:47 [inlined]
 [13] collect(itr::Base.Generator{UnitRange{Int64}, var"#29#31"})
    @ Base .\array.jl:678
 [14] top-level scope
    @ REPL[12]:1
```

@mcmwright, you are right, the code works as is in Julia 1.7.1. For Julia 1.6.2, you may need to collect and permutedims at:

n1, n2 = collect(eachcol(nullspace(permutedims(∇L(s)))))

NB:
Hope @DNF is away while recommending doing this…

2 Likes

Why do you need collect here? Won’t

n1, n2 = eachcol(nullspace(permutedims(∇L(s))))

work as well, with less allocations?

1 Like

Here you can try

xs, ys, zs = ([p[i] for p in C.(s, u')] for i in 1:3)
2 Likes

That’s great, thank you!

@DNF, we can always count on you to save on collect calls :slight_smile:

NB: have VS Code running Julia 1.7.1, and a separate REPL running 1.6.2 with different package versions installed and often get such quick tests all messed up

2 Likes

The @rafael.guerra’s solution seems to be great, but from mathematical point of view it has a drawback, because it does not ensure
the continuity of the vector fields n1, n2, along the line L(s). I searched for an example to illustrate it visually.

Consider a closed line and its tangent vector field, v:

f(t) = [a*cos(t)+a/2, a*sin(t)+a*sqrt(3)/2, b*cos(3*t)]
v(t)  =  = [-a*sin(t), a*cos(t), -3*b*sin(3*t)]

function tube(f, v,  t, s)
    n1, n2 = eachcol(nullspace(permutedims(v(t))))  #v(t) n1(t), n2(t) - orthogonal basis at each point of the curve f(t)   
    return f(t) + r*cos(2π*s) * n1 + r*sin(2π*s)* n2 #parameterization of a tubular surface of directrice f(t), and radius r  
end

r = 0.1; a = 1; b = 0.1   
t = LinRange(0, 2π, 200)
s = LinRange(0, 1, 50);    
x, y, z = [[p[i] for p in tube.(f, v, t, s')] for i in 1:3]    

The corresponding tube has a gap along the small circle, corresponding to t=π:

At t=π, and t=0 the vector n1 (the green one in this animation https://github.com/empet/Datasets/blob/master/Images/basis.gif’)
suddenly change its direction.
This particularity can be confirmed by computing the determinant det(cat(v(t), n1, n2, dims=2)) at each t[k] k=1:length(t) (n1, n2 are orthogonal to v(t[k])).
We can notice a change of the determinant sign at t=π, and at t=0, i.e. the basis orientation changes and this property leads to a gap at t=π.
Why not at t=0, as well? Due to periodicity of the function f,
f(0)=f(2*pi), the points on the tube circle corresponding to t=2π are just the points on the circle
corresponding to t=0, but slightly rotated.
That’s why here there isn’t a discontinuity (due to interpolation).

After seeing the animation I’ve just changed the direction of the vectors n1 at the points
where I got a negative determinant and the corresponding tube is this one:

For an arbitrary curve we cannot control nullspace(permutedims(v(t))) to return a basis that ensure tube continuity.
Hence we must use the theoretical definition of a perfect tube which involves not an arbitrary orthogonal basis at each point, but the Frenet basis consisting in vectors:
(v(t), (v(t) x f''(t)) x v(t), v(t) x f''(t)) =(tangent, principal normal, binormal)(t)

4 Likes

@empet, what happens if the curve has straight line portions where the curvature and torsion are zero? Will the Frenet principal normal and binormal vectors still be well defined?

@rafael.guerra Obviously you cannot define the Frenet frame along any curve. Its parameterization, f, must have linear independent tangent and acceleration at each point., i.e. f must be of class C^2, and
f’(t) x f’’(t)\neq 0 for any t (<=> f’(t)\neq 0, f’’(t)\neq 0, and the two vectors are not colinear; in fact if their cross product is different from 0 implies that both are non-null vectors)).