lib/geometry/src/primitive.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 //! Shapes over f32 points. Every shape is a plain value with no storage of its
2 //! own, so a caller can build one on the stack for each query.
3 const std = @import("std");
4 const linear = @import("linear");
5 const tolerance = @import("tolerance.zig");
6
7 const assert = std.debug.assert;
8 const Mat3 = linear.Mat3;
9 const Vec3 = linear.Vec3;
10
11 /// A half-line from `origin` along `direction`. A query parameter `t` names
12 /// the point `origin + direction * t`, so with a unit direction `t` is a
13 /// distance.
14 pub const Ray = struct {
15 origin: Vec3,
16 direction: Vec3,
17
18 pub fn at(ray: Ray, t: f32) Vec3 {
19 return ray.origin.add(ray.direction.scale(t));
20 }
21 };
22
23 /// The points `a + (b - a) * t` for `t` in [0, 1].
24 pub const Segment = struct {
25 a: Vec3,
26 b: Vec3,
27
28 pub fn at(segment: Segment, t: f32) Vec3 {
29 return segment.a.add(segment.b.sub(segment.a).scale(t));
30 }
31
32 pub fn bounds(segment: Segment) Aabb {
33 return .{ .min = segment.a.min(segment.b), .max = segment.a.max(segment.b) };
34 }
35 };
36
37 pub const Sphere = struct {
38 center: Vec3,
39 radius: f32,
40
41 pub fn bounds(sphere: Sphere) Aabb {
42 assert(sphere.radius >= 0);
43 return Aabb.fromPoint(sphere.center).inflate(sphere.radius);
44 }
45 };
46
47 /// The points within `radius` of the segment from `a` to `b`.
48 pub const Capsule = struct {
49 a: Vec3,
50 b: Vec3,
51 radius: f32,
52
53 pub fn axis(capsule: Capsule) Segment {
54 return .{ .a = capsule.a, .b = capsule.b };
55 }
56
57 pub fn bounds(capsule: Capsule) Aabb {
58 assert(capsule.radius >= 0);
59 return capsule.axis().bounds().inflate(capsule.radius);
60 }
61 };
62
63 /// An axis-aligned box. A box with `min` above `max` on any axis is empty.
64 pub const Aabb = extern struct {
65 min: Vec3,
66 max: Vec3,
67
68 /// The box that contains nothing and grows to fit whatever joins it.
69 pub const empty: Aabb = .{
70 .min = Vec3.splat(std.math.inf(f32)),
71 .max = Vec3.splat(-std.math.inf(f32)),
72 };
73
74 pub fn fromPoint(p: Vec3) Aabb {
75 return .{ .min = p, .max = p };
76 }
77
78 pub fn join(a: Aabb, b: Aabb) Aabb {
79 return .{ .min = a.min.min(b.min), .max = a.max.max(b.max) };
80 }
81
82 pub fn joinPoint(box: Aabb, p: Vec3) Aabb {
83 return .{ .min = box.min.min(p), .max = box.max.max(p) };
84 }
85
86 pub fn inflate(box: Aabb, margin: f32) Aabb {
87 return .{ .min = box.min.sub(Vec3.splat(margin)), .max = box.max.add(Vec3.splat(margin)) };
88 }
89
90 pub fn isEmpty(box: Aabb) bool {
91 return box.min.x > box.max.x or box.min.y > box.max.y or box.min.z > box.max.z;
92 }
93
94 pub fn center(box: Aabb) Vec3 {
95 return box.min.add(box.max).scale(0.5);
96 }
97
98 pub fn extent(box: Aabb) Vec3 {
99 return box.max.sub(box.min);
100 }
101
102 /// Half the surface area, the cost weight of a surface area heuristic.
103 pub fn halfArea(box: Aabb) f32 {
104 if (box.isEmpty()) return 0;
105 const e = box.extent();
106 return e.x * e.y + e.y * e.z + e.z * e.x;
107 }
108
109 /// The largest coordinate magnitude of either corner.
110 pub fn magnitude(box: Aabb) f32 {
111 return @max(box.min.abs().maxComponent(), box.max.abs().maxComponent());
112 }
113
114 pub fn containsPoint(box: Aabb, p: Vec3) bool {
115 return p.x >= box.min.x and p.x <= box.max.x and
116 p.y >= box.min.y and p.y <= box.max.y and
117 p.z >= box.min.z and p.z <= box.max.z;
118 }
119
120 pub fn overlaps(a: Aabb, b: Aabb) bool {
121 return a.min.x <= b.max.x and b.min.x <= a.max.x and
122 a.min.y <= b.max.y and b.min.y <= a.max.y and
123 a.min.z <= b.max.z and b.min.z <= a.max.z;
124 }
125 };
126
127 /// A box rotated by `axes`, whose columns are orthonormal, around `center`.
128 pub const Obb = struct {
129 center: Vec3,
130 axes: Mat3 = .{},
131 half_extents: Vec3,
132
133 /// `p` in the box frame, where the box is the axis-aligned
134 /// `[-half_extents, half_extents]`.
135 pub fn toLocal(box: Obb, p: Vec3) Vec3 {
136 const d = p.sub(box.center);
137 return .{ .x = d.dot(box.axes.cols[0]), .y = d.dot(box.axes.cols[1]), .z = d.dot(box.axes.cols[2]) };
138 }
139
140 pub fn fromLocal(box: Obb, p: Vec3) Vec3 {
141 return box.center.add(box.axes.mulVec(p));
142 }
143 };
144
145 /// A triangle with counter-clockwise winding around its normal.
146 pub const Triangle = extern struct {
147 a: Vec3,
148 b: Vec3,
149 c: Vec3,
150
151 /// The unnormalized normal, twice the area in length.
152 pub fn areaNormal(triangle: Triangle) Vec3 {
153 return triangle.b.sub(triangle.a).cross(triangle.c.sub(triangle.a));
154 }
155
156 /// Whether the corners are collinear, or so close to it that the area
157 /// normal is rounding noise.
158 pub fn isDegenerate(triangle: Triangle) bool {
159 const edge1 = triangle.b.sub(triangle.a);
160 const edge2 = triangle.c.sub(triangle.a);
161 return tolerance.isParallel(edge1.cross(edge2).lengthSq(), edge1.lengthSq() * edge2.lengthSq());
162 }
163
164 /// The unit normal, or null for a degenerate triangle.
165 pub fn normal(triangle: Triangle) ?Vec3 {
166 if (triangle.isDegenerate()) return null;
167 return triangle.areaNormal().normalized(0);
168 }
169
170 pub fn centroid(triangle: Triangle) Vec3 {
171 return triangle.a.add(triangle.b).add(triangle.c).scale(1.0 / 3.0);
172 }
173
174 pub fn bounds(triangle: Triangle) Aabb {
175 return Aabb.fromPoint(triangle.a).joinPoint(triangle.b).joinPoint(triangle.c);
176 }
177
178 pub fn vertex(triangle: Triangle, index: usize) Vec3 {
179 assert(index < 3);
180 return switch (index) {
181 0 => triangle.a,
182 1 => triangle.b,
183 else => triangle.c,
184 };
185 }
186 };
187
188 /// The points `p` with `normal ยท p = offset`. The normal is unit length, so
189 /// `signedDistance` is a distance.
190 pub const Plane = struct {
191 normal: Vec3,
192 offset: f32,
193
194 pub fn fromPointNormal(point: Vec3, unit_normal: Vec3) Plane {
195 return .{ .normal = unit_normal, .offset = unit_normal.dot(point) };
196 }
197
198 /// The plane through a triangle, or null for a degenerate triangle.
199 pub fn fromTriangle(triangle: Triangle) ?Plane {
200 const n = triangle.normal() orelse return null;
201 return fromPointNormal(triangle.a, n);
202 }
203
204 pub fn signedDistance(plane: Plane, p: Vec3) f32 {
205 return plane.normal.dot(p) - plane.offset;
206 }
207 };
208
209 comptime {
210 assert(@sizeOf(Triangle) == @sizeOf([9]f32));
211 assert(@sizeOf(Aabb) == @sizeOf([6]f32));
212 }
213
214 const testing = std.testing;
215
216 test "an empty box absorbs whatever joins it" {
217 const box = Aabb.empty.joinPoint(Vec3.init(1, 2, 3)).joinPoint(Vec3.init(-1, 0, 5));
218 try testing.expect(Aabb.empty.isEmpty());
219 try testing.expectEqual(Vec3.init(-1, 0, 3), box.min);
220 try testing.expectEqual(Vec3.init(1, 2, 5), box.max);
221 try testing.expectEqual(@as(f32, 0), Aabb.empty.halfArea());
222 try testing.expectEqual(@as(f32, 2 * 2 + 2 * 2 + 2 * 2), (Aabb{ .min = .{}, .max = Vec3.splat(2) }).halfArea());
223 }
224
225 test "shape bounds enclose the shape" {
226 const capsule = Capsule{ .a = Vec3.init(0, 0, 0), .b = Vec3.init(0, 2, 0), .radius = 0.5 };
227 try testing.expectEqual(Vec3.init(-0.5, -0.5, -0.5), capsule.bounds().min);
228 try testing.expectEqual(Vec3.init(0.5, 2.5, 0.5), capsule.bounds().max);
229 const sphere = Sphere{ .center = Vec3.init(1, 1, 1), .radius = 1 };
230 try testing.expectEqual(Vec3.init(2, 2, 2), sphere.bounds().max);
231 }
232
233 test "a triangle reports its unit normal or none" {
234 const flat = Triangle{ .a = .{}, .b = Vec3.init(1, 0, 0), .c = Vec3.init(0, 1, 0) };
235 try testing.expectEqual(Vec3.init(0, 0, 1), flat.normal().?);
236 const sliver = Triangle{ .a = .{}, .b = Vec3.init(1, 0, 0), .c = Vec3.init(2, 0, 0) };
237 try testing.expectEqual(@as(?Vec3, null), sliver.normal());
238 const far_pair = Triangle{ .a = Vec3.init(300, 1500, -600), .b = Vec3.init(-0.0, 0, 3300), .c = Vec3.init(0, 0, 3300) };
239 try testing.expect(far_pair.isDegenerate());
240 try testing.expect(!flat.isDegenerate());
241 const plane = Plane.fromTriangle(flat).?;
242 try testing.expectEqual(@as(f32, 3), plane.signedDistance(Vec3.init(5, 5, 3)));
243 }
244
245 test "an oriented box maps points into its frame and back" {
246 const quarter = Mat3.fromCols(Vec3.init(0, 1, 0), Vec3.init(-1, 0, 0), Vec3.init(0, 0, 1));
247 const box = Obb{ .center = Vec3.init(1, 0, 0), .axes = quarter, .half_extents = Vec3.init(2, 1, 1) };
248 const p = Vec3.init(1, 2, 0);
249 try testing.expectEqual(Vec3.init(2, 0, 0), box.toLocal(p));
250 try testing.expectEqual(p, box.fromLocal(box.toLocal(p)));
251 }