File size: 8,233 Bytes
9d6c005 | 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 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | """Research implementation of exact-evidence certificates for opaque spheres.
The mathematical guarantee is in real arithmetic. Float64 plus conservative
tolerance is tested here, not a formally verified interval-arithmetic kernel.
Only renderer-authoritative geometry and fixed receiver/emitter identities
are accepted. Unreported geometry changes invalidate the contract.
"""
import numpy as np
def visibility_certificate(scene, points, lights, guard=1e-9):
"""One segment query returns visibility and a conservative clearance.
Distance to the trimmed segment is evaluated for every sphere. This costs
more arithmetic than a Boolean early-out shadow query; timings include it.
"""
p, l = np.broadcast_arrays(np.asarray(points, float), np.asarray(lights, float))
if p.shape[-1] != 3 or not np.isfinite(p).all() or not np.isfinite(l).all():
raise ValueError("Finite 3D segment endpoints required")
d = l-p
length2 = (d*d).sum(-1)
if (length2 <= 0).any() or guard < 0:
raise ValueError("Positive segment length and nonnegative guard required")
clearance = np.full(length2.shape, np.inf)
for sphere in scene.spheres:
a = np.clip(((sphere[:3]-p)*d).sum(-1)/length2, 1e-6, 1-1e-6)
signed = np.linalg.norm(p+a[..., None]*d-sphere[:3], axis=-1)-sphere[3]
clearance = np.minimum(clearance, signed)
v = (clearance > 0).astype(float)
margin = np.maximum(np.abs(clearance)-guard, 0)
ambiguous = np.abs(clearance) <= guard
if np.any(ambiguous):
# Preserve the reference oracle's endpoint/tangency convention. Such
# entries have no positive robustness margin and require no rounding
# claim about which side of zero a floating-point clearance lies on.
shape = v.shape
flat = v.reshape(-1).copy(); mask = ambiguous.reshape(-1)
flat[mask] = scene.visibility(p.reshape(-1,3)[mask], l.reshape(-1,3)[mask])
v = flat.reshape(shape)
return v, np.nextafter(margin, 0)
class CertificateMemory:
"""Canonical evidence plus a cumulative geometry-motion budget.
Clearance shrinks by at most max_s(||delta center_s||+|delta radius_s|)
per update. The accumulated bound is conservative, including on return paths.
Light endpoints are fixed. New receiver positions can query old certificates
with an additional Euclidean displacement bound, but cannot commit to old IDs.
"""
def __init__(self, points, lights, spheres, namespace="scene", mode="margin"):
self.points = np.array(points, float, copy=True)
self.lights = np.array(lights, float, copy=True)
self.geometry = np.array(spheres, float, copy=True)
if (self.points.ndim != 2 or self.points.shape[1] != 3 or len(self.points) < 1
or self.lights.ndim != 2 or self.lights.shape[1] != 3 or len(self.lights) < 1
or self.geometry.ndim != 2 or self.geometry.shape[1] != 4 or len(self.geometry) < 1
or not all(np.isfinite(a).all() for a in (self.points, self.lights, self.geometry))
or (self.geometry[:, 3] <= 0).any() or not namespace
or mode not in ("margin", "epoch", "unsafe")):
raise ValueError("Invalid canonical certificate domain")
self.namespace, self.mode = str(namespace), mode
shape = (len(self.points), len(self.lights))
self.values = np.full(shape, np.nan, np.float32)
self.margins = np.zeros(shape, np.float64)
self.stamps = np.zeros(shape, np.float64)
self.budget = 0.0
@property
def nbytes(self):
return sum(a.nbytes for a in (self.points, self.lights, self.geometry,
self.values, self.margins, self.stamps))
def _ids(self, ids):
x = np.asarray(ids)
if x.ndim != 1 or not np.issubdtype(x.dtype, np.integer) or (x < 0).any() or (x >= len(self.points)).any():
raise ValueError("Invalid canonical receiver IDs")
return x
def begin_geometry(self, spheres, namespace=None):
g = np.array(spheres, float, copy=True)
if (g.shape != self.geometry.shape or not np.isfinite(g).all()
or (g[:, 3] <= 0).any()):
raise ValueError("Topology change requires a new canonical memory")
if namespace is not None and str(namespace) != self.namespace:
raise ValueError("Namespace mismatch: create fresh memory")
displacement = np.max(np.linalg.norm(g[:, :3]-self.geometry[:, :3], axis=1)
+ np.abs(g[:, 3]-self.geometry[:, 3]))
if displacement > 0:
self.budget = float(np.nextafter(self.budget+displacement, np.inf))
if self.mode == "epoch":
self.values.fill(np.nan)
self.geometry = g
def lookup(self, ids, query_points=None, extra_motion=0.0):
ids = self._ids(ids)
delta = np.zeros(len(ids))
if query_points is not None:
p = np.asarray(query_points, float)
if p.shape != self.points[ids].shape or not np.isfinite(p).all():
raise ValueError("Query points must match canonical anchor IDs")
delta = np.linalg.norm(p-self.points[ids], axis=1)
if not np.isfinite(extra_motion) or extra_motion < 0:
raise ValueError("Nonnegative future-motion bound required")
movement = self.budget-self.stamps[ids]+delta[:, None]+extra_motion
known = np.isfinite(self.values[ids])
if self.mode != "unsafe":
known &= (self.margins[ids] > movement) | (movement == 0)
return np.nan_to_num(self.values[ids], nan=0).astype(float), known
def commit(self, ids, indices, visibility, margins):
ids = self._ids(ids)
j, v, m = np.asarray(indices), np.asarray(visibility), np.asarray(margins)
if (j.ndim != 1 or j.shape != ids.shape or not np.issubdtype(j.dtype, np.integer)
or v.shape != j.shape or m.shape != j.shape or (j < 0).any()
or (j >= len(self.lights)).any() or not np.isfinite(v).all()
or ((v != 0) & (v != 1)).any() or not np.isfinite(m).all() or (m < 0).any()):
raise ValueError("Expected one exact observation and margin per receiver")
self.values[ids, j] = v
self.margins[ids, j] = m
self.stamps[ids, j] = self.budget
def save(self, path):
np.savez_compressed(path, points=self.points, lights=self.lights, geometry=self.geometry,
values=self.values, margins=self.margins, stamps=self.stamps,
budget=np.float64(self.budget), namespace=np.array(self.namespace), mode=np.array(self.mode))
@classmethod
def load(cls, path, namespace):
with np.load(path, allow_pickle=False) as d:
if str(d['namespace']) != str(namespace):
raise ValueError("Namespace mismatch")
obj = cls(d['points'], d['lights'], d['geometry'], namespace, str(d['mode']))
shape = obj.values.shape
v, m, s = d['values'], d['margins'], d['stamps']
budget = float(d['budget'])
if (v.shape != shape or m.shape != shape or s.shape != shape or v.dtype != np.float32
or not np.isfinite(m).all() or not np.isfinite(s).all() or (m < 0).any()
or not np.isfinite(budget) or budget < 0 or (s < 0).any() or (s > budget).any()
or np.isinf(v).any() or ((v[np.isfinite(v)] != 0) & (v[np.isfinite(v)] != 1)).any()):
raise ValueError("Malformed certificate checkpoint")
obj.values, obj.margins, obj.stamps = v.copy(), m.copy(), s.copy()
obj.budget = budget
return obj
def enclosure(bound, values, known):
"""Deterministic linear RGB interval, conditional on valid certificates."""
b, v, k = np.asarray(bound), np.asarray(values), np.asarray(known, bool)
if b.ndim != 3 or v.shape != b.shape[:2] or k.shape != v.shape or not np.isfinite(b).all() or (b < 0).any():
raise ValueError("Expected nonnegative contribution bounds and known mask")
lo = (b*np.where(k, v, 0)[..., None]).sum(1)
return lo, lo+(b*(~k)[..., None]).sum(1)
|