helixish: attempt at the whole thing

[moebius3.git] / moebius.py

diff --git a/moebius.py b/moebius.py
index d9888d2cea0a896b511dd417ab9d898e73120db0..e5a81acad9d0d458bbc823769b5cea55176e0801 100644 (file)
--- a/moebius.py
+++ b/moebius.py
@@ -5,15 +5,33 @@ import numpy as np
 from numpy import cos, sin
 
 from bezier import BezierSegment
+from helixish import HelixishCurve
+from moenp import *
 
 import sys
 
-tau = np.pi * 2
-
-origin = np.array((0,0,0))
-unit_x = np.array((1,0,0))
-unit_y = np.array((0,1,0))
-unit_z = np.array((0,0,1))
+class DoubleCubicBezier():
+  def __init__(db, cp):
+    single = BezierSegment(cp)
+    midpoint = np.array(single.point_at_t(0.5))
+    mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint
+    mid_dirn /= np.linalg.norm(mid_dirn)
+    ocp_factor = 0.5
+    mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) +
+                                    np.linalg.norm(cp[3] - cp[2]))
+    db.b0 = BezierSegment([ cp[0],
+                            cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
+                            midpoint - mid_dirn * mid_scale,
+                            midpoint ])
+    db.b1 = BezierSegment([ midpoint,
+                            midpoint + mid_dirn * mid_scale,
+                            cp[2] * ocp_factor + cp[3] * (1-ocp_factor),
+                            cp[3] ])
+  def point_at_t(db, t):
+    if t < 0.5:
+      return db.b0.point_at_t(t*2)
+    else:
+      return db.b1.point_at_t(t*2 - 1)
 
 class ParametricCircle:
   def __init__(pc, c, r0, r1):
@@ -63,20 +81,22 @@ class MoebiusHalf:
     m.nu      = nu
     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
-    m._beziers = [ m._bezier(theta) for theta in m._thetas ]
-  def _bezier(m, theta):
+    m._beziers = [ m._bezier(theta, HelixishCurve) for theta in m._thetas ]
+    #check = int(nu/3)
+    #m._beziers[check] = m._bezier(m._thetas[check], HelixishCurve)
+  def _bezier(m, theta, constructor=DoubleCubicBezier):
     cp = [None] * 4
     cp[0] =               m.edge   .point(theta)
     cp[3] =               m.midline.point(theta*2)
     ncp3 = np.linalg.norm(cp[3])
     cpt = ncp3 * ncp3 / 4
     cp2scale = m._cp2b.point_at_t(ncp3/2)
-    cp1scale = 0.75 * cpt + 0.33 * (1-cpt)
+    cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
     #      file=sys.stderr)
     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
-    return BezierSegment(cp)
+    return constructor(cp)
   def point(m, iu, t):
     '''
     0 <= iu <= nu     meaning 0 <= u <= 1
@@ -88,21 +108,33 @@ class MoebiusHalf:
     '''
     return np.array(m._beziers[iu].point_at_t(t))
 
-  def point_offset(m, iu, t, offset):
+  def details(m, iu, t):
     '''
-    offset by offset perpendicular to the surface
-     at the top (iu=t=0), in the y direction
+    returns tuple of 4 vectors:
+           - point location
+           - normal (+ve is in the +ve y direction at iu=t=0) unit vector
+           - along extent   (towrds +ve iu)                   unit vector
+           - along traverse (towards +ve t)                   unit vector
     '''
+    if iu == m.nu:
+      return m.details(0, t)
     p = m.point(iu, t)
+    vec_t = unit_v( m.point(iu, t + 0.01) - p )
     if t == 0:
-      dirn = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
-    elif iu == m.nu:
-      return m.point_offset(0, t, offset)
+      normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
+      vec_u = np.cross(vec_t, normal)
     else:
-      vec_t = m.point(iu,   t + 0.01) - p
-      vec_u = m.point(iu+1, t)        - p
-      dirn =  np.cross(vec_u, vec_t)
-    return p + offset * dirn / np.linalg.norm(dirn)
+      vec_u = unit_v( m.point(iu+1, t) - p )
+      normal = np.cross(vec_u, vec_t)
+    return p, normal, vec_u, vec_t
+
+  def point_offset(m, iu, t, offset):
+    '''
+    offset by offset perpendicular to the surface
+     at the top (iu=t=0), +ve offset is in the +ve y direction
+    '''
+    p, normal, dummy, dummy = m.details(iu, t)
+    return p + offset * normal
 
 class Moebius():
   def __init__(m, nv, nw):
@@ -136,3 +168,17 @@ class Moebius():
     return m.h.point_offset(offset=
                             offset if w <= m.nt else -offset,
                             **m._vw2tiu_kw(v,w))
+
+  def details(m, v, w):
+    '''
+    returns tuple of 4 vectors:
+           - point location
+           - normal (+ve is in the +ve y direction at iu=t=0) unit vector
+           - along extent   (towrds +ve v)                    unit vector
+           - along traverse (towards +ve w)                   unit vector
+    '''
+    p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
+    if w > m.nt:
+      normal = -normal
+      vec_w  = -vec_w
+    return p, normal, vec_v, vec_w