A tubular surface (or tube surface) is generated by a 3D curve, called spine, and a moving circle of radius r, with center on the spine and included in planes orthogonal to curve.
Tubular surfaces are associated to spines that are biregular, that is, they have a $C^2$ parameterization, $c:[a,b]\to \mathbb{R}^3$, with velocity, $\dot{c}(t)$, and acceleration, $\ddot{c}(t)$, that are non-null and non-colinear vectors: $\dot{c}(t)\times \ddot{c}(t)\neq 0$.
A tube of prescribed curvature and torsion is defined by a spine parameterized by the arc length, i.e. by $c(s)$, with constant speed, $||\dot{c}(s)||=1$, and non-null acceleration, $\ddot{c}(s)\neq 0$, for all $s$.
The given curvature and torsion, $\kappa(s)$, $\tau(s)$, define the Frenet-Serre equations: $$\begin{array}{lll} \dot{e}_1(s)&=&\kappa(t)e_2(s)\\ \dot{e}_2(s)&=&-\kappa(s)e_1(s)+\tau(s)e_3(s)\\ \dot{e}_3(s)&=&-\tau(s)e_2(s),\\ \end{array} $$
where $e_1(s), e_2(s), e_3(s)$ are respectively the unit vectors of tangent, principal normal and binormal along the curve.
Frenet-Serre equations completed with the equation $ \dot{c}(s)=e_1(s)$ define a system of ordinary differential equations, with 12 equations and 12 unknown functions. The last three coordinates of a solution represent the discretized curve, $c(s)$, starting from an initial point, with a prescribed Frenet frame at that point.
We define below a tubular surface with highly oscillating curvature and constant torsion of the spine.
In [1]:
import numpy as np
from scipy import integrate
In [2]:
def curv(s):#curvature
return 3*np.sin(s/10.)*np.sin(s/10.)
def tors(s):#torsion is constant
return 0.35
def Frenet_eqns(x, s):# right side vector field of the system of ODE
return [ curv(s)*x[3],
curv(s)*x[4],
curv(s)*x[5],
-curv(s)*x[0]+tors(s)*x[6],
-curv(s)*x[1]+tors(s)*x[7],
-curv(s)*x[2]+tors(s)*x[8],
-tors(s)*x[3],
-tors(s)*x[4],
-tors(s)*x[5],
x[0], x[1], x[2]]
Integrate the system, with an initial point consisting in the initial Frenet frame (of three orthonormal vectors) and the initial position of the curve, $c(0)$:
In [3]:
x_init=np.array([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0])
s_final=150# [0, s_final] is the interval of integration
N=1000
s_div=np.linspace(0, s_final, N)
X=integrate.odeint(Frenet_eqns, x_init, s_div)
normal=X[:, 3:6].T
binormal=X[:, 6:9].T
curve=X[:, 9:].T
xc, yc, zc=curve# lists of coordinates of the spine points
Now we define a tubular surface that has as spine the above curve.
A tubular surface having as spine a curve, $c(s)$, parameterized by the arclength, is defined as follows: $r(s,u)=c(s)+\varepsilon(e_2(s)cos(u)+e_3(s)sin(u))$, $0<\varepsilon <<1$, $u\in[0, 2\pi]$. $\varepsilon$ is the radius of circles orthogonal to the spine.
In [4]:
import plotly.plotly as py
from plotly.graph_objs import *
Define a function that sets the plot layout:
In [9]:
axis = dict(
showbackground=True,
backgroundcolor="rgb(230, 230,230)",
gridcolor="rgb(255, 255, 255)",
zerolinecolor="rgb(255, 255, 255)",
)
noaxis=dict(showbackground=False,
showgrid=False,
showline=False,
showticklabels=False,
ticks='',
title='',
zeroline=False
)
def set_layout(title='', width=800, height=800, axis_type=axis, aspect=(1, 1, 1)):
return Layout(
title=title,
autosize=False,
width=width,
height=height,
showlegend=False,
scene=Scene(xaxis=XAxis(axis_type),
yaxis=YAxis(axis_type),
zaxis=ZAxis(axis_type),
aspectratio=dict(x=aspect[0],
y=aspect[1],
z=aspect[2]
)
)
)
The colorscale for the tubular surface:
In [10]:
my_colorscale=[[0.0, 'rgb(46, 107, 142)'],
[0.1, 'rgb(41, 121, 142)'],
[0.2, 'rgb(36, 134, 141)'],
[0.3, 'rgb(31, 147, 139)'],
[0.4, 'rgb(30, 160, 135)'],
[0.5, 'rgb(40, 174, 127)'],
[0.6, 'rgb(59, 186, 117)'],
[0.7, 'rgb(85, 198, 102)'],
[0.8, 'rgb(116, 208, 84)'],
[0.9, 'rgb(151, 216, 62)'],
[1.0, 'rgb(189, 222, 38)']]
Define a function that evaluates the tube parameterization, $r(s,u)=(x, y, z)$, at the meshgrid np.meshgrid(s_div, u):
In [13]:
def create_tube(spine_points, normal, binormal,
epsilon=0.2, colorscale=my_colorscale, zmin=None, zmax=None):
#returns an instance of the Plotly Surface, representing a tube
u=np.linspace(0, 2*np.pi, 100)
x,y,z=[np.outer(spine_points[k,:], np.ones(u.shape))+
epsilon*(np.outer(normal[k, :], np.cos(u))+np.outer(binormal[k,:], np.sin(u)))
for k in range(3)]
if zmin is not None and zmax is not None:
return Surface(x=x, y=y, z=z, zmin=zmin, zmax=zmax,
colorscale=colorscale,
colorbar=dict(thickness=25, lenmode='fraction', len=0.75))
else:
return Surface(x=x, y=y, z=z,
colorscale=colorscale,
colorbar=dict(thickness=25, lenmode='fraction', len=0.75))
The keywords zmin, zmax are set when we connect at least two tubular surfaces. They define the color bounds for
the tubular structure.
In [14]:
tube=create_tube(curve, normal, binormal, epsilon=0.1)
In [15]:
data1=Data([tube])
layout1=set_layout(title='Tubular surface', aspect=(1,1,1.05))
fig1 = Figure(data=data1, layout=layout1)
py.sign_in('empet', '')
py.iplot(fig1, filename='tubular-cst-torsion')
Out[15]:
If a general biregular parameterization, $c(t)$, of the spine is given, then we have to do some analytical computations by hand, in order to get the directions $\dot{c}(t)$, $\ddot{c}(t)$, $\dot{c}(t)\times \ddot{c}(t)$, of the velocity (tangent), acceleration, and binormals along the curve.
Then we define Python functions, tangent, acceleration, curve_normals, that compute the unit vectors of these directions.
Finally the unit vector of the principal normal is computed as $n(t)=b(t)\times tg(t)$, where $b(t), tg(t)$ are the unit vectors of binormals and tangents.
The tube parameterization, $$r(t,u)=c(t)+\varepsilon(n(t)\cos(u)+b(t)\sin(u)), t\in[tm, tM], u\in[0,2\pi],$$ is evaluated at a meshgrid.
We illustrate a tubular structure, called Hopf link, defined by two tubes, having the spines parameterized by: $$c(t)=(\pm a+\cos(t), \sin(t), \pm b\sin(t)), t\in[0, 2\pi]$$ The first spine corresponds to $a=0.5, b=0.2$, and the second one, to $a=-0.5, b=-0.2$.
In [16]:
from numpy import sin, cos, pi
In [17]:
def spine_param( a, b, tm, tM, nr):
#spine parameterization c:[tm, tM]-->R^3
# a, b are parameters on which the spine parameterization depends
# nr is the number of points to be evaluated on spine
t=np.linspace(tm, tM, nr )# nr is the number of points to ve evaluated on spine
return t, a+cos(t), sin(t), b*sin(t)
In [18]:
def tangent( a, b, t):
# returns the unit tangent vectors along the spine curve
v=np.vstack((-sin(t), cos(t), b*cos(t)))
return v/np.vstack((np.linalg.norm(v, axis=0),)*3)
def acceleration( a, b, t):
# returns the unit acceleration vectors along the spine
v=np.array([ -cos(t), -sin(t), -b*sin(t)])
return v/np.vstack((np.linalg.norm(v, axis=0),)*3)
In [19]:
def curve_normals(a, b):
# computes and returns the point coordinates on spine, and the unit normal vectors
t, xc, yc, zc=spine_param(a,b, 0.0, 2*pi, 100)
tang=tangent(a,b, t)
binormal=np.cross(tang, acceleration(a, b, t), axis=0)
binormal=binormal/np.vstack((np.linalg.norm(binormal, axis=0),)*3)
normal=np.cross(binormal, tang, axis=0)
return np.vstack((xc, yc, zc)), normal, binormal
In [20]:
epsilon=0.025 # the radius of each tube
zm=[]# list of min z-values on both tubes
zM=[]# list of max z-values on both tubes
spine1, normal1, binormal1=curve_normals(0.5, 0.2)
zm.append(min(spine1[2,:]))
zM.append(max(spine1[2,:]))
spine2, normal2, binormal2=curve_normals(-0.5, -0.2)
zm.append(min(spine2[2,:]))
zM.append(max(spine2[2,:]))
zmin=min(zm)
zmax=max(zM)
tube1=create_tube(spine1, normal1, binormal1, epsilon=epsilon, zmin=zmin, zmax=zmax)
tube2=create_tube(spine2, normal2, binormal2, epsilon=epsilon, zmin=zmin, zmax=zmax)
layout2=set_layout(title='Hopf link', aspect=(1, 0.75, 0.35))
In [21]:
data2=Data([tube1,tube2])
fig2 = Figure(data=data2, layout=layout2)
In [22]:
py.sign_in('empet', '')
py.iplot(fig2, filename='Hopf-link')
Out[22]:
If we take all combinations of signs for the parameters, a, b, we get an interesting configuration of tubes communicating with each other:
In [23]:
from IPython.display import HTML
HTML('<iframe src=https://plot.ly/~empet/13930/comunicating-rings/ width=900 height=700></iframe>')
Out[23]:
Tubular surfaces are particular surfaces in the class of canal surfaces. A canal surface is again defined by a biregular spine, $c(t)$, but the circles ortogonal to spine have variable radii, gigen by a $C^1$-function, $r(t)$, with $|r'(t)|<||\dot{c}(t)||$.
The parameterization of a canal surface is:
$$r(t,u)=c(t)-\displaystyle\frac{r(t)r'(t)}{||\dot{c}(t)||^2}\dot{c}(t)+ \displaystyle\frac{r(t)\sqrt{||\dot{c}(t) ||^2-r'(t)^2}}{||\dot{c}(t) ||}(n(t)\cos(u)+b(t)\sin(u))$$We plot the canal surface of spine, $c(t)=(10\cos(t), 10\sin(t), 0)$, and radius function $r(t)=2+\cos(2t)$, $t\in[0,2\pi]$.
In [24]:
def radius_deriv(t):
return 2+cos(2*t), -2*sin(2*t)
In [25]:
def create_canal(spine, normal, binormal, term,
colorscale=my_colorscale, zmin=None, zmax=None):
#returns an instance of the Plotly Surface, representing a canal surface
#term is the second term in the parameterization
u=np.linspace(0, 2*np.pi, 100)
x,y,z=[np.outer(spine[k,:]-term[k, :], np.ones(u.shape))+\
np.outer(normal[k, :], np.cos(u))+np.outer(binorm[k,:], np.sin(u)) for k in range(3)]
if zmin is not None and zmax is not None:
return Surface(x=x, y=y, z=z, zmin=zmin, zmax=zmax,
colorscale=colorscale,
colorbar=dict(thickness=25, lenmode='fraction', len=0.75))
else:
return Surface(x=x, y=y, z=z,
colorscale=colorscale,
colorbar=dict(thickness=25, lenmode='fraction', len=0.75))
In [26]:
t=np.linspace(0, 3*pi/2, 50)
xc, yc, zc= 10*cos(t), 10*sin(t), np.zeros(t.shape)
spine=np.vstack((xc,yc, zc))
rt,rdt=radius_deriv(t)# rt is the variable radius r(t), and rdt its derivative
In [27]:
tang=np.vstack((-10*sin(t), 10*cos(t), np.zeros(t.shape))) #c'(t)
cdot_norm=np.vstack((np.linalg.norm(tang, axis=0),)*3)# ||c'(t)||
factor=rt*rdt/cdot_norm**2
term=factor*tang#term.shape=(3, t.shape[0])# second term in canal surface parameterization
R=rt*np.sqrt(cdot_norm**2-rdt**2)/cdot_norm # R.shape (3, t.shape[0]) is the scalar factor in the third term
tangu= (tang/cdot_norm) #unit tangent vector
acceler=np.vstack((-10*cos(t), -10*sin(t), np.zeros(t.shape)))
acceler= acceler/np.vstack((np.linalg.norm(acceler, axis=0),)*3)#unit acceleration vector
binorm=np.cross(tangu, acceler, axis=0)
binorm=binorm/np.vstack((np.linalg.norm(binorm, axis=0),)*3)#unit binormal vector
normal=np.cross(binorm, tangu, axis=0)# unit normal vector
binorm=R*binorm
normal=R*normal
In [28]:
canal=create_canal(spine, normal, binorm, term, colorscale=my_colorscale)
In [29]:
layout3=set_layout(title='Canal surface', axis_type=axis, aspect=(1, 1, 0.25))
data3=Data([canal])
fig3 = Figure(data=data3, layout=layout3)
py.sign_in('empet', '')
py.iplot(fig3, filename='Canal-surf')
Out[29]:
Finally, we stress that in order to get a tubular looking surface, we have to set the aspect ratio of the plot that respects the real ratios between axes lengths. Otherwise the tube is deformed.
In [31]:
from IPython.core.display import HTML
def css_styling():
styles = open("./custom.css", "r").read()
return HTML(styles)
css_styling()
Out[31]: