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