# Tube plots?
**URL:** https://discourse.julialang.org/t/tube-plots/65999
**Category:** Visualization
**Tags:** makie
**Created:** [August 7, 2021, 8:30pm UTC](https://discourse.julialang.org/t/tube-plots/65999 "2021-08-07T20:30:52Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@mcmwright](https://discourse.julialang.org/u/mcmwright)
#### Post date: [August 7, 2021, 8:30pm UTC](https://discourse.julialang.org/t/tube-plots/65999/1 "2021-08-07T20:30:52Z")
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?
---
### Author:  [@jules](https://discourse.julialang.org/u/jules)
#### Post date: [August 7, 2021, 9:57pm UTC](https://discourse.julialang.org/t/tube-plots/65999/2 "2021-08-07T21:57:20Z")
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](http://makie.juliaplots.org/stable/plotting_functions/arrows.html)
---
### Author:  [@ffreyer](https://discourse.julialang.org/u/ffreyer)
#### Post date: [August 8, 2021, 1:32am UTC](https://discourse.julialang.org/t/tube-plots/65999/3 "2021-08-08T01:32:36Z")
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).
```julia
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.
```julia
# ...
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:
```julia
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)
```
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [August 8, 2021, 11:27am UTC](https://discourse.julialang.org/t/tube-plots/65999/4 "2021-08-08T11:27:27Z")
An alternative way is to generate a parametric tube surface along the line path.
```julia
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_
---
### Author:  [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)
#### Post date: [August 8, 2021, 1:12pm UTC](https://discourse.julialang.org/t/tube-plots/65999/5 "2021-08-08T13:12:03Z")
This looks gorgeous! I love the colors!
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [August 8, 2021, 1:50pm UTC](https://discourse.julialang.org/t/tube-plots/65999/6 "2021-08-08T13:50:12Z")
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.
---
### Author:  [@mcmwright](https://discourse.julialang.org/u/mcmwright)
#### Post date: [August 8, 2021, 4:58pm UTC](https://discourse.julialang.org/t/tube-plots/65999/7 "2021-08-08T16:58:01Z")
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
````julia
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
```
````
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [August 8, 2021, 5:14pm UTC](https://discourse.julialang.org/t/tube-plots/65999/8 "2021-08-08T17:14:40Z")
@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:
```julia
n1, n2 = collect(eachcol(nullspace(permutedims(∇L(s)))))
```
_ **NB:** _
_Hope @DNF is away while recommending doing this…_
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 8, 2021, 5:38pm UTC](https://discourse.julialang.org/t/tube-plots/65999/9 "2021-08-08T17:38:51Z")
> [@rafael.guerra](#):
>
> `n1, n2 = collect(eachcol(nullspace(permutedims(∇L(s)))))`
Why do you need `collect` here? Won’t
```julia
n1, n2 = eachcol(nullspace(permutedims(∇L(s))))
```
work as well, with less allocations?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [August 8, 2021, 5:40pm UTC](https://discourse.julialang.org/t/tube-plots/65999/10 "2021-08-08T17:40:36Z")
> [@rafael.guerra](#):
>
> `xs, ys, zs = [[p[i] for p in C.(s, u')] for i in 1:3]`
Here you can try
```julia
xs, ys, zs = ([p[i] for p in C.(s, u')] for i in 1:3)
```
---
### Author:  [@mcmwright](https://discourse.julialang.org/u/mcmwright)
#### Post date: [August 8, 2021, 7:15pm UTC](https://discourse.julialang.org/t/tube-plots/65999/11 "2021-08-08T19:15:33Z")
That’s great, thank you!
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [August 8, 2021, 8:07pm UTC](https://discourse.julialang.org/t/tube-plots/65999/12 "2021-08-08T20:07:13Z")
@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_
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [August 9, 2021, 10:54pm UTC](https://discourse.julialang.org/t/tube-plots/65999/13 "2021-08-09T22:54:32Z")
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:
```julia
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’](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)`
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [August 9, 2021, 11:40pm UTC](https://discourse.julialang.org/t/tube-plots/65999/14 "2021-08-09T23:40:34Z")
@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?
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [August 10, 2021, 5:14am UTC](https://discourse.julialang.org/t/tube-plots/65999/15 "2021-08-10T05:14:06Z")
@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)).
---
### Author:  [@lazarusA](https://discourse.julialang.org/u/lazarusA)
#### Post date: [October 6, 2023, 10:32am UTC](https://discourse.julialang.org/t/tube-plots/65999/16 "2023-10-06T10:32:00Z")
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?
---
### Author:  [@ffreyer](https://discourse.julialang.org/u/ffreyer)
#### Post date: [October 6, 2023, 11:56am UTC](https://discourse.julialang.org/t/tube-plots/65999/17 "2023-10-06T11:56:51Z")
You could just push `positions[1]` and `positions[end]` to the vertices and add two sets of faces connecting to those points. Something like
```julia
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
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [October 6, 2023, 8:00pm UTC](https://discourse.julialang.org/t/tube-plots/65999/18 "2023-10-06T20:00:05Z")
**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:
```julia
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](https://discourse.julialang.org/t/drawing-a-tubular-path-with-meshes-jl-blog-post/104203/2), as well as a link to examples.
---
### Author:  [@lazarusA](https://discourse.julialang.org/u/lazarusA)
#### Post date: [October 6, 2023, 9:36pm UTC](https://discourse.julialang.org/t/tube-plots/65999/19 "2023-10-06T21:36:15Z")
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.
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [October 7, 2023, 7:15am UTC](https://discourse.julialang.org/t/tube-plots/65999/20 "2023-10-07T07:15:59Z")
Use DelaunayTriangulations.jl to triangulate a tube.
Example:
```julia
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:
[Next page](https://discourse.julialang.org/t/tube-plots/65999.md?page=2)