From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Wed, 22 Nov 2017 00:24:19 +0000 (+0000)
Subject: helixish: wip
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=24f59fed463189b12d123ee2d71461e402a37e13;p=moebius3.git

helixish: wip

Signed-off-by: Ian Jackson <ijackson@chiark.greenend.org.uk>
---

diff --git a/helixish.py b/helixish.py
new file mode 100644
index 0000000..541bb79
--- /dev/null
+++ b/helixish.py
@@ -0,0 +1,110 @@
+
+from __future__ import print_function
+
+import numpy as np
+from numpy import cos, sin
+
+import sys
+
+def augment(v): return np.append(v, 1)
+def augment0(v): return np.append(v, 0)
+def unaugment(v): return v[0:3]
+
+class HelixishCurve():
+  def __init__(hc, cp):
+    p = cp[0]
+    q = cp[3]
+    dp = unit_v(cp[1]-cp[0])
+    dq = unit_v(cp[3]-cp[2])
+
+    # the initial attempt
+    #   - solve in the plane containing dP and dQ
+    #   - total distance normal to that plane gives mu
+    #   - now resulting curve is not parallel to dP at P
+    #     nor dQ at Q, so tilt it
+    #   - [[ pick as the hinge point the half of the curve
+    #     with the larger s or t ]] not yet implemented
+    #   - increase the other distance {t,s} by a bodge factor
+    #     approx distance between {Q,P} and {Q,P}' due to hinging
+    #     but minimum is 10% of (wlog) {s,t} [[ not quite like this ]]
+
+    dPQplane_normal = np.cross(dp, dq)
+    if (np.norm(dPQplane_normal) < 1E6):
+      dPQplane_normal += [0, 0, 1E5]
+    dPQplane_normal = unit_v(dPQplane_normal)
+
+    dPQplane_basis = np.column_stack(np.cross(dp, dPQplane_normal),
+                                     dp,
+                                     dPQplane_normal,
+                                     p);
+    dPQplane_basis = np.vstack(dPQplane_basis, [0,0,0,1])
+    dPQplane_into = np.linalg.inv(dPQplane_basis)
+
+    dp_plane = unaugment(dPQplane_into * augment0(dp))
+    dq_plane = unaugment(dPQplane_into * augment0(dq))
+    q_plane  = unaugment(dPQplane_into * augment(q))
+    dist_pq_plane = np.linalg.norm(q_plane)
+
+    # two circular arcs of equal maximum possible radius
+    # algorithm courtesy of Simon Tatham (`Railway problem',
+    # pers.comm. to ijackson@chiark 23.1.2004)
+    railway_angleoffset = atan2(*q_plane[0:1])
+    railway_theta =                      tau/4 - railway_angleoffset
+    railway_phi   = atan2(*dq_plane[0:1]) - railway_angleoffset
+    railway_cos_theta = cos(railway_theta)
+    railway_cos_phi   = cos(railway_phi)
+    if railway_cos_theta**2 + railway_cos_phi**2 > 1E6:
+      railway_roots = np.roots([
+        2 * (1 + cos(railway_theta - railway_phi)),
+        2 * (railway_cos_theta - railway_cos_phi),
+        -1
+        ])
+      for railway_r in railway_roots:
+        def railway_CPQ(pq, dpq):
+          nonlocal railway_r
+          return pq + railway_r * [-dpq[1], dpq[0]]
+
+        railway_CP = railway_CPQ([0,0,0],       dp_plane)
+        railway_QP = railway_CPQ(q_plane[0:2], -dq_plane)
+        railway_midpt = 0.5 * (railway_CP + railway_QP)
+
+        best_st = None
+        def railway_ST(C, start, end):
+          nonlocal railway_r
+          delta = atan2(*(end - C)[0:2]) - atan2(start - C)[0:2]
+          s = delta * railway_r
+
+        try_s = railway_ST(railway_CP, [0,0], midpt)
+        try_t = railway_ST(railway_CP, midpt, q_plane)
+        try_st = try_s + try_t
+        if best_st is None or try_st < best_st:
+          start_la = 1/r
+          start_s = try_s
+          start_t = try_t
+          best_st = try_st
+          start_mu = q_plane[2] / (start_s + start_t)
+
+    else: # twoarcs algorithm is not well defined
+      start_la = 0.1
+      start_s = dist_pq_plane * .65
+      start_t = dist_pq_plane * .35
+      start_mu = 0.05
+
+    bodge = max( q_plane[2] * mu,
+                 (start_s + start_t) * 0.1 )
+    start_s += 0.5 * bodge
+    start_t += 0.5 * bodge
+    start_kappa = 0
+    start_gamma = 1
+
+    
+
+    # we work in two additional coordinate systems:
+    # for both these:
+    #   P is at the origin
+    #   |PQ| = 1
+    # for findcurve:
+    #   dP is the +ve y axis
+    #   Q lies in the x/y plane
+    # for calculating the initial attempt:
+    #   P is at the origin