from numpy import cos, sin
from bezier import BezierSegment
+from moenp import *
+from moebez import *
+from curveopt 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 ParametricCircle:
def __init__(pc, c, r0, r1):
''' circle centred on c
def radius(pc, theta):
return pc._r0 * cos(theta) + pc._r1 * sin(theta)
def point(pc, theta):
- print('point %s %f' % (pc, theta), file=sys.stderr)
return pc._c + pc.radius(theta)
class Twirler(ParametricCircle):
but with an orientation at each point, orthogonal to
the circle
the orientation circles round cycles times during the
- whole cycle
+ whole cycle (if cycles<1, to get a whole circling of
+ the dirn around the circle, must pass some theta>tau)
begin_zeta is the angle from outwards at theta==0
- positive meaning in the direction of r0 x r1
+ positive meaning in the direction of r0 x r1 from r0
'''
ParametricCircle.__init__(tw, c, r0, r1)
tw._cycles = cycles
tw._begin_zeta = begin_zeta
tw._axis = np.cross(r0, r1)
- def dirn(tw, theta):
- zeta = tw._begin_zeta + theta * tw._cycles
+ def dirn(tw, theta, extra_zeta=0):
+ ''' unit vector for dirn at theta,
+ but rotated extra_theta more around the circle
+ '''
+ zeta = tw._begin_zeta + theta * tw._cycles + extra_zeta
r = tw.radius(theta)
return cos(zeta) * r + sin(zeta) * tw._axis
-class Moebius:
- def __init__(m, n_u): # ix_u will be in [0, n_u>
- m.edge = Twirler(origin, unit_z, unit_x, 2, 0)
- m.midline = Twirler(-unit_z, unit_z, unit_y, 1, 0)
- m._beziers = [ m._bezier(u) for u in np.linspace(0, 1, n_u) ]
- def _bezier(m,u):
- theta = u * tau
+class MoebiusHalf:
+ def __init__(m, nu, nt):
+ '''
+ MoebiusHalf().edge is a Twirler for the edge,
+ expecting theta = u * tau (see MoebiusHalf().point)
+ with dirn pointing into the surface
+ '''
+ m.edge = Twirler(origin, unit_z, unit_x, -2, tau/2)
+ m.midline = Twirler(-unit_z, unit_z, unit_y, -0.5, 0)
+ m.nu = nu
+ m.nt = nt
+ 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._dbeziers = [ m._dbezier(theta) for theta in m._thetas ]
+ #check = int(nu/3)-1
+ #checks = (
+ # m._dbezier(m._thetas[check], OptimisedCurve),
+ # m._dbezier(m._thetas[check+1], OptimisedCurve),
+ # m._dbezier(m._thetas[check+2], OptimisedCurve),
+ #)
+ #for c in checks: c.point_at_it(0)
+ def _dbezier(m, theta, dconstructor=OptimisedCurve):
cp = [None] * 4
cp[0] = m.edge .point(theta)
- cp[1] = cp[0] + 0.5 * m.edge .dirn (theta)
cp[3] = m.midline.point(theta*2)
- cp[2] = cp[3] + 0.5 * m.midline.dirn (theta*2)
- return BezierSegment(cp)
- def point(m, ix_u, t):
- return m._beziers[ix_u].point_at_t(t)
+ ncp3 = np.linalg.norm(cp[3])
+ cpt = ncp3 * ncp3 / 4
+ cp2scale = m._cp2b.point_at_t(ncp3/2)
+ 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 dconstructor(cp, m.nt)
+ def point(m, iu, it):
+ '''
+ 0 <= iu <= nu meaning 0 <= u <= 1
+ along the extent (well, along the edge)
+ 0 and 1 are both the top half of the flat traverse
+ 0.5 is the bottom half of the flat traverse
+ 0 <= it <= nt across the half-traverse
+ 0 is the edge, 1 is the midline
+ '''
+ return np.array(m._dbeziers[iu].point_at_it(it))
+
+ def details(m, iu, t):
+ '''
+ 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:
+ normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
+ vec_u = np.cross(vec_t, normal)
+ else:
+ 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):
+ '''
+ 0 <= v <= nw along the extent, v=0 is the flat traverse
+ 0 <= w <= nv across the traverse nw must be even
+ the top is both v=0, w=0 v=nv, w=nw
+ '''
+ assert(nw % 1 == 0)
+ m.nv = nv
+ m.nw = nw
+ m.nt = nw/2
+ m.h = MoebiusHalf(nu=nv*2, nt=m.nt)
+
+ def _vw2tiu_kw(m, v, w):
+ if w <= m.nt:
+ it = w
+ iu = v
+ else:
+ it = m.nw - w
+ iu = m.nv + v
+ #print('v,w=%d,%d => it,iu=%d,%d' % (v,w,it,iu),
+ # file=sys.stderr)
+ return { 'it': it, 'iu': iu }
+
+ def point(m, v, w):
+ return m.h.point(**m._vw2tiu_kw(v,w))
+
+ def point_offset(m, v, w, offset):
+ 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
~ian / moebius3.git / blobdiff