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):
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
'''
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
'''
- p = MoebiusHalf.point(m, iu, t)
+ 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 MoebiusHalf.point_offset(m, 0, t, offset)
+ normal = m.edge.dirn(m._thetas[iu], extra_zeta=-tau/4)
+ vec_u = np.cross(vec_t, normal)
else:
- vec_t = MoebiusHalf.point(m, iu, t + 0.01) - p
- vec_u = MoebiusHalf.point(m, 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(MoebiusHalf):
+class Moebius():
def __init__(m, nv, nw):
'''
0 <= v <= nw along the extent, v=0 is the flat traverse
m.nw = nw
m.nt = nw/2
m._t_vals = np.linspace(0, 1, m.nt+1)
- MoebiusHalf.__init__(m, nu=nv*2)
+ m.h = MoebiusHalf(nu=nv*2)
def _vw2tiu_kw(m, v, w):
if w <= m.nt:
return { 't': m._t_vals[it], 'iu': iu }
def point(m, v, w):
- return MoebiusHalf.point(m, **m._vw2tiu_kw(v,w))
+ return m.h.point(**m._vw2tiu_kw(v,w))
def point_offset(m, v, w, offset):
- return MoebiusHalf.point_offset(m,
- offset=
- offset if w <= m.nt else -offset,
- **m._vw2tiu_kw(v,w))
+ 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