| """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) |
|
|
