{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Parallel Transport Frames scratch work\n", "\n", "This is based on a couple things.\n", "\n", "It follows the paper, [\"Parallel Transport Approach to Curve Framing\" by Andrew J. Hanson & Hui Ma](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.42.8103), which I first encountered via the [geom-types/src/ptf.org](https://github.com/thi-ng/geom/blob/master/geom-types/src/ptf.org) examples from Karsten Schmidt's work with [thi.ng/geom](https://github.com/thi-ng/geom). If I refer to any page number or section, I'm talking about this paper.\n", "\n", "It was also partly an experiment in learning a bit of [JupyterLab](https://github.com/jupyterlab/jupyterlab), which I only tried out after running into some [Julia](https://julialang.org/) problems on Jupyter. I ordinarily also work in [Maxima](http://maxima.sourceforge.net/), but decided to try out [SymPy](https://www.sympy.org/en/index.html). I was rather impressed with both: JupyterLab addresses many of my complaints with Jupyter, and I found SymPy quite capable and part of a more efficient workflow compared to manually transferring results from Maxima into Python.\n", "\n", "## TODO items\n", "\n", "- Don't require template to be a single closed path\n", "- Figure out why I am leaving one lone vertex at the origin\n", "- Figure out \"face_normals didn't match triangles, ignoring!\" (if it matters)\n", "- Find a way to start out the trimesh view with things visible\n", "- Figure out why this knot is not actually a [cinquefoil](https://en.wikipedia.org/wiki/Cinquefoil_knot)..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sympy\n", "from sympy import symbols, cos, sin, diff, lambdify\n", "from sympy.vector import CoordSys3D\n", "import numpy\n", "import quaternion\n", "import stl.mesh\n", "import trimesh\n", "\n", "import parallel_transport\n", "import quat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(p,q) [Torus knot](https://en.wikipedia.org/wiki/Torus_knot) as a parametric formula (torus lies in XY plane):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "r,p,q,t = symbols(\"r,p,q,t\")\n", "r = cos(q*t)+2\n", "N = CoordSys3D('N')\n", "knot_vec = (r * cos(p*t))*N.i + (r * sin(p*t))*N.j + (-sin(q*t))*N.k" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def vec2python(v):\n", " return sympy.python(N.origin.locate_new('na', v).express_coordinates(N))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "q = Symbol('q')\n", "t = Symbol('t')\n", "p = Symbol('p')\n", "e = ((cos(q*t) + 2)*cos(p*t), (cos(q*t) + 2)*sin(p*t), -sin(q*t))\n" ] } ], "source": [ "print(vec2python(knot_vec))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = Symbol('p')\n", "q = Symbol('q')\n", "t = Symbol('t')\n", "e = (-p*(cos(q*t) + 2)*sin(p*t) - q*sin(q*t)*cos(p*t), p*(cos(q*t) + 2)*cos(p*t) - q*sin(p*t)*sin(q*t), -q*cos(q*t))\n" ] } ], "source": [ "print(vec2python(diff(knot_vec,t)))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 3. , 0. , -0. ],\n", " [ 2.99960977, 0.04714521, -0.00628629],\n", " [ 2.9984392 , 0.09427692, -0.01257233],\n", " [ 2.99648866, 0.14138161, -0.01885787],\n", " [ 2.99375876, 0.18844581, -0.02514266],\n", " [ 2.99025034, 0.23545602, -0.03142647],\n", " [ 2.9859645 , 0.28239879, -0.03770903],\n", " [ 2.98090257, 0.32926067, -0.0439901 ],\n", " [ 2.97506613, 0.37602826, -0.05026943],\n", " [ 2.968457 , 0.42268818, -0.05654678]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# can I do this without having to make it a Point?\n", "torus52_expr = knot_vec.subs(p, 5).subs(q, 2)\n", "torus52 = sympy.lambdify([t], N.origin.locate_new('na', torus52_expr).express_coordinates(N))\n", "torus52_tan = sympy.lambdify([t], N.origin.locate_new('na', diff(torus52_expr, t)).express_coordinates(N))\n", "#ts = numpy.linspace(0, 2*numpy.pi, 50000)\n", "ts = numpy.linspace(0, 2*numpy.pi, 2000)\n", "points = numpy.array(torus52(ts)).T\n", "points[:10,:]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-0. , 0.9912279 , -0.13216372],\n", " [-0.01640804, 0.99109232, -0.13216196],\n", " [-0.0328116 , 0.99068562, -0.1321567 ],\n", " [-0.04920622, 0.99000792, -0.13214791],\n", " [-0.06558742, 0.98905939, -0.13213562],\n", " [-0.08195074, 0.98784029, -0.13211981],\n", " [-0.09829171, 0.98635095, -0.13210049],\n", " [-0.11460587, 0.98459179, -0.13207766],\n", " [-0.13088879, 0.98256327, -0.1320513 ],\n", " [-0.14713601, 0.98026595, -0.13202143]])" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# analytical tangents:\n", "tans = numpy.array(torus52_tan(ts)).T\n", "tans = tans / numpy.linalg.norm(tans, axis=1)[:,numpy.newaxis]\n", "tans[:10,:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "SymPy is failing at the analytical integration, but also [numerical](https://docs.sympy.org/latest/modules/evalf.html#sums-and-integrals). I tried oscillatory quadrature, but it seems that despite the cos/sin everywhere, they're not in the right form." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# tor = parallel_transport.torsion_integral(torus52_expr, t, 0, 2*sympy.pi)\n", "# Works (returns 1.408) but is horrendously slow\n", "#tor_num = tor.evalf(4)\n", "#tor_num" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "My own numerical approach seems to work:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-1.4071908841197978" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "parallel_transport.torsion_integral_approx(torus52_expr, t, 0, 2*numpy.pi, 0.001)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do I want analytical tangents, or approximations? approximations may, oddly enough, be 'more correct' given that I am already approximating the curve." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "points = numpy.array(torus52(ts)).T\n", "tan_approx = parallel_transport.approx_tangents(points)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are rather close regardless:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9.38346943765395e-06" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "numpy.mean(numpy.abs(tans - tan_approx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, make some sort of template I can duplicate. Use a vertex ordering that can go directly into the STL library." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "template = numpy.array([\n", " [0.0, 0.0, 0.0],\n", " [0.5, 0.0, 0.0],\n", " [0.5, 0.0, 0.5],\n", "])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and attempt the algorithm from \"Parallel Transport to Curve Framing\":" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Angular change: -1.4071908841197978, -80.62609862934731 deg\n", "Approximate arc length: 64.25469846213409\n" ] } ], "source": [ "# Precompute cross-products, angles, arc length, and angular difference:\n", "t1 = tans\n", "t2 = numpy.roll(tans, -1, axis=0)\n", "cps = numpy.cross(t1, t2)\n", "cps /= numpy.linalg.norm(cps, axis=1)[:,numpy.newaxis]\n", "thetas = numpy.arccos(numpy.sum(t1*t2, axis=1))\n", "ang = parallel_transport.torsion_integral_approx(torus52_expr, t, 0, 2*numpy.pi, 0.001)\n", "print(\"Angular change: {}, {} deg\".format(ang, ang/numpy.pi*180))\n", "total_length = parallel_transport.approx_arc_length(points)\n", "print(\"Approximate arc length: {}\".format(total_length))\n", "# One can add additional twists here if a multiple of a whole rotation:\n", "# ang += 40 * 2*numpy.pi" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "done\n" ] } ], "source": [ "# Build up the actual mesh, correcting spin as we go:\n", "count = len(points)\n", "data = numpy.zeros(count*2*3*len(template), dtype=stl.mesh.Mesh.dtype)\n", "# Each step on the curve involves len(template) vertices which\n", "# must connect to respective vertices of the next step. This\n", "# requires two triangles for each vertex of the template, and\n", "# each triangle requires 3 vertices.\n", "vectors = data[\"vectors\"]\n", "length = 0.0\n", "# Use q for a running transformation of our template:\n", "q = numpy.quaternion(1, 0, 0, 0)\n", "for i in range(count):\n", " t1 = tans[i,:]\n", " b = cps[i]\n", " theta = thetas[i]\n", " # Add another transform (note we post-multiply):\n", " q = quat.rotation_quaternion(b, theta) * q\n", " length += numpy.linalg.norm(points[i,:] - points[(i+1) % count,:])\n", " # Correction by arc length (see section 3.1):\n", " theta2 = -ang * length / total_length\n", " # Add a correction for twist. This doesn't update our\n", " # running transform - it's only for final geometry:\n", " q_corrected = quat.rotation_quaternion(t1, theta2) * q\n", " # Connect every vertex in the frame (the template now rotated &\n", " # shifted) to the 'last' frame location via 2 triangles\n", " p_adj = quat.conjugate_by(template.copy(), q_corrected) + points[i,:]\n", " if i == 0:\n", " first_p_adj = p_adj\n", " else:\n", " for j,_ in enumerate(p_adj):\n", " idx = len(template)*2*3*(i-1) + j*2\n", " f1,f2 = parallel_transport.gen_faces(\n", " p_adj[(j+1) % len(template),:],\n", " p_adj[j,:],\n", " last_p_adj[(j+1) % len(template),:],\n", " last_p_adj[j,:],\n", " )\n", " vectors[idx] = f1\n", " vectors[idx+1] = f2\n", " last_p_adj = p_adj\n", "\n", "# Patch up first & last:\n", "for j,_ in enumerate(last_p_adj):\n", " idx = len(template)*2*3*(count-1) + j*2\n", " f1,f2 = parallel_transport.gen_faces(\n", " last_p_adj[j,:],\n", " last_p_adj[(j+1) % len(template),:],\n", " first_p_adj[j,:],\n", " first_p_adj[(j+1) % len(template),:])\n", " vectors[idx] = f1\n", " vectors[idx+1] = f2\n", "\n", "mesh = stl.mesh.Mesh(data)\n", "mesh_fname = \"test_2k.stl\"\n", "mesh.save(mesh_fname)\n", "print(\"done\")" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "face_normals didn't match triangles, ignoring!\n", "/home/hodapp/.local/lib/python3.6/site-packages/IPython/core/display.py:694: UserWarning: Consider using IPython.display.IFrame instead\n", " warnings.warn(\"Consider using IPython.display.IFrame instead\")\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = trimesh.load_mesh(mesh_fname)\n", "m.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.8" } }, "nbformat": 4, "nbformat_minor": 4 }