Tubular surfaces

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$.

Tubular surface defined by a spine curve parameterized by arc length

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]:

Tubular surface with a spine curve of given parameterization

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]:

Canal (Channels) surfaces

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]: