moebius.py

   1
   2 from __future__ import print_function
   3
   4 import numpy as np
   5 from numpy import cos, sin
   6
   7 from bezier import BezierSegment
   8 from recursivebezier import RecursiveBezierishCurve
   9 from moenp import *
  10
  11 import sys
  12
  13 class DoubleCubicBezier():
  14   def __init__(db, cp):
  15     single = BezierSegment(cp)
  16     midpoint = np.array(single.point_at_t(0.5))
  17     mid_dirn = single.point_at_t(0.5 + 0.001) - midpoint
  18     mid_dirn /= np.linalg.norm(mid_dirn)
  19     ocp_factor = 0.5
  20     mid_scale = ocp_factor * 0.5 * (np.linalg.norm(cp[1] - cp[0]) +
  21                                     np.linalg.norm(cp[3] - cp[2]))
  22     db.b0 = BezierSegment([ cp[0],
  23                             cp[1] * ocp_factor + cp[0] * (1-ocp_factor),
  24                             midpoint - mid_dirn * mid_scale,
  25                             midpoint ])
  26     db.b1 = BezierSegment([ midpoint,
  27                             midpoint + mid_dirn * mid_scale,
  28                             cp[2] * ocp_factor + cp[3] * (1-ocp_factor),
  29                             cp[3] ])
  30   def point_at_t(db, t):
  31     if t < 0.5:
  32       return db.b0.point_at_t(t*2)
  33     else:
  34       return db.b1.point_at_t(t*2 - 1)
  35
  36 class ParametricCircle:
  37   def __init__(pc, c, r0, r1):
  38     ''' circle centred on c
  39         with theta=0 point at c+r0
  40         and with theta=tau/4 point at c+r1 '''
  41     pc._c  = c
  42     pc._r0 = r0
  43     pc._r1 = r1
  44   def radius(pc, theta):
  45     return pc._r0 * cos(theta) + pc._r1 * sin(theta)
  46   def point(pc, theta):
  47     return pc._c + pc.radius(theta)
  48
  49 class Twirler(ParametricCircle):
  50   def __init__(tw, c, r0, r1, cycles, begin_zeta):
  51     ''' circle centred on c, etc.
  52         but with an orientation at each point, orthogonal to
  53           the circle
  54         the orientation circles round cycles times during the
  55           whole cycle (if cycles<1, to get a whole circling of
  56           the dirn around the circle, must pass some theta>tau)
  57         begin_zeta is the angle from outwards at theta==0
  58           positive meaning in the direction of r0 x r1 from r0
  59     '''
  60     ParametricCircle.__init__(tw, c, r0, r1)
  61     tw._cycles = cycles
  62     tw._begin_zeta = begin_zeta
  63     tw._axis = np.cross(r0, r1)
  64   def dirn(tw, theta, extra_zeta=0):
  65     ''' unit vector for dirn at theta,
  66          but rotated extra_theta more around the circle
  67     '''
  68     zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
  69     r = tw.radius(theta)
  70     return cos(zeta) * r + sin(zeta) * tw._axis
  71
  72 class MoebiusHalf:
  73   def __init__(m, nu):
  74     '''
  75      MoebiusHalf().edge is a Twirler for the edge,
  76       expecting theta = u * tau (see MoebiusHalf().point)
  77       with dirn pointing into the surface
  78     '''
  79     m.edge    = Twirler(origin,  unit_z, unit_x, -2, tau/2)
  80     m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
  81     m.nu      = nu
  82     m._thetas = [ u * tau for u in np.linspace(0, 1, nu+1) ]
  83     m._cp2b = BezierSegment([ (c,) for c in [0.33,0.33, 1.50]])
  84     m._beziers = [ m._bezier(theta) for theta in m._thetas ]
  85   def _bezier(m, theta, constructor=RecursiveBezierishCurve):
  86     cp = [None] * 4
  87     cp[0] =               m.edge   .point(theta)
  88     cp[3] =               m.midline.point(theta*2)
  89     ncp3 = np.linalg.norm(cp[3])
  90     cpt = ncp3 * ncp3 / 4
  91     cp2scale = m._cp2b.point_at_t(ncp3/2)
  92     cp1scale = 1.0 * cpt + 0.33 * (1-cpt)
  93     #print('u=%d ncp3=%f cp2scale=%f' % (u, ncp3, cp2scale),
  94     #      file=sys.stderr)
  95     cp[1] = cp[0] + cp1scale * m.edge   .dirn (theta)
  96     cp[2] = cp[3] + cp2scale * m.midline.dirn (theta*2)
  97     return constructor(cp)
  98   def point(m, iu, t):
  99     '''
 100     0 <= iu <= nu     meaning 0 <= u <= 1
 101                       along the extent (well, along the edge)
 102                        0 and 1 are both the top half of the flat traverse
 103                        0.5        is the bottom half of the flat traverse
 104     0 <=  t   <= 1    across the half-traverse
 105                        0 is the edge, 1 is the midline
 106     '''
 107     return np.array(m._beziers[iu].point_at_t(t))
 108
 109   def details(m, iu, t):
 110     '''
 111     returns tuple of 4 vectors:
 112            - point location
 113            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 114            - along extent   (towrds +ve iu)                   unit vector
 115            - along traverse (towards +ve t)                   unit vector
 116     '''
 117     if iu == m.nu:
 118       return m.details(0, t)
 119     p = m.point(iu, t)
 120     vec_t = unit_v( m.point(iu, t + 0.01) - p )
 121     if t == 0:
 122       normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
 123       vec_u = np.cross(vec_t, normal)
 124     else:
 125       vec_u = unit_v( m.point(iu+1, t) - p )
 126       normal = np.cross(vec_u, vec_t)
 127     return p, normal, vec_u, vec_t
 128
 129   def point_offset(m, iu, t, offset):
 130     '''
 131     offset by offset perpendicular to the surface
 132      at the top (iu=t=0), +ve offset is in the +ve y direction
 133     '''
 134     p, normal, dummy, dummy = m.details(iu, t)
 135     return p + offset * normal
 136
 137 class Moebius():
 138   def __init__(m, nv, nw):
 139     '''
 140       0 <= v <= nw    along the extent,    v=0 is the flat traverse
 141       0 <= w <= nv    across the traverse  nw must be even
 142       the top is both   v=0, w=0  v=nv, w=nw
 143     '''
 144     assert(nw % 1 == 0)
 145     m.nv = nv
 146     m.nw = nw
 147     m.nt = nw/2
 148     m._t_vals = np.linspace(0, 1, m.nt+1)
 149     m.h = MoebiusHalf(nu=nv*2)
 150
 151   def _vw2tiu_kw(m, v, w):
 152     if w <= m.nt:
 153       it = w
 154       iu = v
 155     else:
 156       it = m.nw - w
 157       iu = m.nv + v
 158     #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
 159     #      file=sys.stderr)
 160     return { 't': m._t_vals[it], 'iu': iu }
 161
 162   def point(m, v, w):
 163     return m.h.point(**m._vw2tiu_kw(v,w))
 164
 165   def point_offset(m, v, w, offset):
 166     return m.h.point_offset(offset=
 167                             offset if w <= m.nt else -offset,
 168                             **m._vw2tiu_kw(v,w))
 169
 170   def details(m, v, w):
 171     '''
 172     returns tuple of 4 vectors:
 173            - point location
 174            - normal (+ve is in the +ve y direction at iu=t=0) unit vector
 175            - along extent   (towrds +ve v)                    unit vector
 176            - along traverse (towards +ve w)                   unit vector
 177     '''
 178     p, normal, vec_v, vec_w = m.h.details(**m._vw2tiu_kw(v,w))
 179     if w > m.nt:
 180       normal = -normal
 181       vec_w  = -vec_w
 182     return p, normal, vec_v, vec_w