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 http://makie.juliaplots.org/stable/plotting_functions/arrows.html
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)
8 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 
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)
6 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)).
I’m trying to use this code, is there any easy way to close the bottom and upper endings, the first and last circles? Namely, do some faces also there, just there?
You could just push positions[1] and positions[end] to the vertices and add two sets of faces connecting to those points. Something like
for i in 1:segments
push!(faces, GLTriangleFace(end-1, i, mod1(i+1, segments))
end
for i in 1:segments
push!(faces, GLTriangleFace(end, end - 1 - mod1(i+1, segments), end - 1 - i)
end
I didn’t run this, so you may still need to adjust the order of indices to get correct normals and the indices might be off by a bit. But conceptually that should work I think
Hanson algorithm to represent a tubular surface, given its directrice line
This algorithm ensures the continuity of a moving othonormal frame along the directrice.
Let \gamma(t) = (x(t), y(t), z(t)) be the directrice parameterization with respect to the standard orthonormal frame, \mathcal{F}=(O; (e_1, e_2, e_3)).
Denote by q(θ; w) the unit quaternion that defines the rotation about the unit vector w, with an angle θ. It is represented by the 4-vector (\cos(\theta/2), \sin(\theta/2)w).
To this unit quaternion corresponds the rotation, of matrix:
M= Rotations.UnitQuaternion(cos(θ/2, (sin(θ/2)*w)...))
Let τ(t) be the unit tangent vector field along the curve. The vector field of unit quaternions
Q(t) =q(\theta(t), w(t)), with \theta(t)=\arccos(e_1 \cdot \tau(t)), and w=(e_1 \times \tau(t))/||e_1 \times \tau(t)||, associated to the directrice, represents a path in the group of rotations, i.e. a field of orthonormal bases consisting in the columns of rotation matrices.
The first column of such a matrix is the unit tangent vector, \tau(t), and the last two, the orthonormal vectors denoted n_1(t), n_2(t), which define the tube parametrization:
S(u,v)=γ(u) + r\cos(v) * n₁(u) + r\sin(v) * n₂(u), \:\: (u, v)\in[a,b]\times[0, 2\pi]
r is the tube radius.
The function implementing this algorithm is posted on this thread:
https://discourse.julialang.org/t/drawing-a-tubular-path-with-meshes-jl-blog-post/104203/2, as well as a link to examples.
Oh, I have seen those before. They look great, however I want the meshes, hence some extra steps will be needed in order to get there. Do you happened to have a version getting there? Otherwise, it shouldn’t be that hard to get there looking at how Makie does a surface.
Use DelaunayTriangulations.jl to triangulate a tube.
Example:
using DelaunayTriangulation
import ForwardDiff.derivative
import LinearAlgebra: norm, cross
import Rotations.UnitQuaternion
dγ(t)=derivative(γ, t) # t->̇γ(t)
ddγ(t)=derivative(dγ, t) #t->̈γ(t)
function tube(γ::T, u::S, v::S; radius=0.1) where {S<:Real, T<:Function}
.
.
.
#Plot a knot from the following directrice
γ = t->[cos(t)−3*cos(4*t), sin(t)+3*sin(4*t), sin(2*t)] #more knots here: https://www.scirp.org/html/4-7404441_100804.htm
u = range(0, 2π, length=800)
v = range(0, 2π, length=100);
uv=vec(collect(Iterators.product(u,v)))# vec from a meshgrid over the rectangle [0,2π]x[0, 2π]
tri = triangulate(uv)
#tri.triangles
vertices=tube.(γ, first.(uv), last.(uv); radius=0.2)
From tri.triangles and vertices you can get the tube defined as a mesh:
