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