File size: 4,500 Bytes
9cc92da | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 | """Developable sheet geometry: a generalised cylinder.
The sheet bends along one direction: a unit-speed planar profile curve, given by its turning angle
theta(s), is swept along straight rulings. Paper lengths are preserved exactly. theta is a sum of
components: cylinder, wave, ridge, crease, curl (an edge roll) and fold (an accordion step).
"""
import math
import numpy as np
S_MAX = 200.0
DS = 0.02
THETA_MAX = 1.10
def bent_components(rng, s_view=0.0):
"""Random components of a smoothly bent sheet around s_view."""
comps = []
radius = math.exp(rng.uniform(math.log(50.0), math.log(900.0)))
comps.append(dict(type="cylinder", s0=float(s_view + rng.uniform(-30, 30)),
R=float(radius * rng.choice([-1.0, 1.0]))))
for _ in range(int(rng.integers(0, 2))):
comps.append(dict(type="wave", s0=float(s_view + rng.uniform(-40, 40)),
wavelength=float(rng.uniform(40.0, 150.0)),
amp=float(rng.uniform(0.08, 0.45))))
for _ in range(int(rng.integers(0, 3))):
comps.append(dict(type="ridge", s0=float(s_view + rng.uniform(-35, 35)),
sigma=float(rng.uniform(3.0, 16.0)),
amp=float(rng.uniform(0.15, 0.80) * rng.choice([-1.0, 1.0]))))
for _ in range(int(rng.integers(0, 3))):
comps.append(dict(type="crease", s0=float(s_view + rng.uniform(-35, 35)),
width=float(rng.uniform(0.5, 2.6)), half=float(rng.uniform(0.4, 2.5)),
amp=float(rng.uniform(0.15, 0.80) * rng.choice([-1.0, 1.0]))))
if rng.random() < 0.4:
comps.append(dict(type="curl", s0=float(s_view + rng.uniform(-60, 60)),
R=float(rng.uniform(15.0, 60.0) * rng.choice([-1.0, 1.0])),
side=float(rng.choice([-1.0, 1.0])), width=float(rng.uniform(2.0, 8.0))))
return comps
def eval_theta(s, comps):
"""Turning angle at s. Summation order (smooth components, curls, folds) is part of the
reproducibility of the published datasets."""
th = np.zeros_like(s)
for c in comps:
d = s - c["s0"]
if c["type"] == "cylinder":
th = th + d / c["R"]
elif c["type"] == "wave":
th = th + c["amp"] * np.sin(2.0 * math.pi * d / c["wavelength"])
elif c["type"] == "ridge":
q = d / c["sigma"]
th = th - c["amp"] * q * np.exp(-0.5 * q * q)
elif c["type"] == "crease":
th = th + 0.5 * c["amp"] * (np.tanh((d + c["half"]) / c["width"]) - np.tanh((d - c["half"]) / c["width"]))
for c in comps:
if c["type"] == "curl":
d = c["side"] * (s - c["s0"]) / c["width"]
th = th + c["width"] * np.logaddexp(0.0, d) / c["R"]
for c in comps:
if c["type"] == "fold":
th = th + 0.5 * c["amp"] * (1.0 + np.tanh((s - c["s0"]) / c["width"]))
return th
def _cumtrapz(y, dx):
out = np.zeros_like(y)
out[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * dx)
return out
class Profile:
"""Profile curve tabulated on s in [-S_MAX, S_MAX]; theta limited per point to +-THETA_MAX."""
def __init__(self, comps):
self.s = np.arange(-S_MAX, S_MAX + DS, DS)
th = eval_theta(self.s, comps)
self.theta = THETA_MAX * np.tanh(th / THETA_MAX)
self.f = _cumtrapz(np.cos(self.theta), DS)
self.zeta = _cumtrapz(np.sin(self.theta), DS)
self.f -= np.interp(0.0, self.s, self.f)
self.zeta -= np.interp(0.0, self.s, self.zeta)
def at(self, s):
return np.interp(s, self.s, self.f), np.interp(s, self.s, self.zeta)
def sheet_points(shape, U, V, centre_mm):
"""Paper (U, V) mm -> 3-D points; the sheet is near z = 0, +zeta toward the camera = -z."""
prof = Profile(shape["comps"])
phi = math.radians(shape["ruling_deg"])
N = np.array([math.cos(phi), math.sin(phi), 0.0])
T = np.array([-math.sin(phi), math.cos(phi), 0.0])
X, Y = U - centre_mm[0], V - centre_mm[1]
s = X * N[0] + Y * N[1]
w = X * T[0] + Y * T[1]
f, zeta = prof.at(s)
return f[..., None] * N + w[..., None] * T - zeta[..., None] * np.array([0, 0, 1.0])
def rot(axis, a):
c, s = math.cos(a), math.sin(a)
if axis == "x":
return np.array([[1, 0, 0], [0, c, -s], [0, s, c]], float)
if axis == "y":
return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]], float)
return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]], float)
|