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

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)
```

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

This looks gorgeous! I love the colors!

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.

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…*

Why do you need `collect`

here? Won’t

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

work as well, with less allocations?

Here you can try

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

That’s great, thank you!

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

**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*

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)`

@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)).