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