zig/mach
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

math: Have Ray return fitting subtype, scale precision automatically

Browse source
This commit is contained in:
Joel D. Schüller 2023-10-03 23:03:06 +02:00 • committed by Stephen Gutekanst
parent 0273e12902
commit c03c627780
1 changed files with 157 additions and 127 deletions
Download patch file Download diff file Expand all files Collapse all files

284
src/math/ray.zig
View file

@ -32,141 +32,171 @@ fn maxDim(v: math.Vec3) u8 {
} }
} }
pub const RayHit = packed struct { u: f32, v: f32, w: f32, t: f32 }; // A Ray in three-dimensional space
pub fn Ray(comptime Vec3P: type) type {
// Floating precision, will be either f16, f32, or f64
const P: type = Vec3P.T;
pub const Ray = struct { // Fallback floating point precision to scale fallback according to
origin: math.Vec3, // input precision
direction: math.Vec3, const PP: type = floatFallbackPrecision(P);
// Algorithm based on: return extern struct {
// https://www.jcgt.org/published/0002/01/05/ origin: Vec3P,
/// Check for collision of a ray and a triangle in 3D space. direction: Vec3P,
/// Triangle winding, which determines front- and backface of
/// the given triangle, matters if backface culling is to be
/// enabled. Without backface culling it does not matter.
/// On hit, will return a RayHit which contains distance t
/// and barycentric coordinates.
pub fn triangleIntersect(
ray: *const Ray,
va: *const math.Vec3,
vb: *const math.Vec3,
vc: *const math.Vec3,
backface_culling: bool,
) ?RayHit {
var kz: u32 = maxDim(math.vec3(
@abs(ray.direction.v[0]),
@abs(ray.direction.v[1]),
@abs(ray.direction.v[2]),
));
var kx: u32 = kz + 1;
if (kx == 3)
kx = 0;
var ky: u32 = kx + 1;
if (ky == 3)
ky = 0;
if (ray.direction.v[kz] < 0.0) { /// A ray hit for which xyz represent the barycentric coordinates
const tmp = kx; /// and w represents hit distance t
kx = ky; pub const Hit = math.Vec4;
ky = tmp;
}
const sx: f32 = ray.direction.v[kx] / ray.direction.v[kz]; pub usingnamespace switch (Vec3P) {
const sy: f32 = ray.direction.v[ky] / ray.direction.v[kz]; math.Vec3, math.Vec3h, math.Vec3d => struct {
const sz: f32 = 1.0 / ray.direction.v[kz]; // Algorithm based on:
// https://www.jcgt.org/published/0002/01/05/
/// Check for collision of a ray and a triangle in 3D space.
/// Triangle winding, which determines front- and backface of
/// the given triangle, matters if backface culling is to be
/// enabled. Without backface culling it does not matter.
/// On hit, will return a RayHit which contains distance t
/// and barycentric coordinates.
pub inline fn triangleIntersect(
ray: *const math.Ray,
va: *const Vec3P,
vb: *const Vec3P,
vc: *const Vec3P,
backface_culling: bool,
) ?Hit {
var kz: u8 = maxDim(math.vec3(
@abs(ray.direction.v[0]),
@abs(ray.direction.v[1]),
@abs(ray.direction.v[2]),
));
var kx: u8 = kz + 1;
if (kx == 3)
kx = 0;
var ky: u8 = kx + 1;
if (ky == 3)
ky = 0;
const a: @Vector(3, f32) = va.v - ray.origin.v; if (ray.direction.v[kz] < 0.0) {
const b: @Vector(3, f32) = vb.v - ray.origin.v; const tmp = kx;
const c: @Vector(3, f32) = vc.v - ray.origin.v; kx = ky;
ky = tmp;
}
const ax: f32 = a[kx] - sx * a[kz]; const sx: P = ray.direction.v[kx] / ray.direction.v[kz];
const ay: f32 = a[ky] - sy * a[kz]; const sy: P = ray.direction.v[ky] / ray.direction.v[kz];
const bx: f32 = b[kx] - sx * b[kz]; const sz: P = 1.0 / ray.direction.v[kz];
const by: f32 = b[ky] - sy * b[kz];
const cx: f32 = c[kx] - sx * c[kz];
const cy: f32 = c[ky] - sy * c[kz];
var u: f32 = cx * by - cy * bx; const a: @Vector(3, P) = va.v - ray.origin.v;
var v: f32 = ax * cy - ay * cx; const b: @Vector(3, P) = vb.v - ray.origin.v;
var w: f32 = bx * ay - by * ax; const c: @Vector(3, P) = vc.v - ray.origin.v;
// Double precision fallback //const a: Vec3P = va.sub(&ray.origin);
if (u == 0.0 or v == 0.0 or w == 0.0) { //const b: Vec3P = vb.sub(&ray.origin);
var cxby: f64 = @as(f64, @floatCast(cx)) * @as(f64, @floatCast(by)); //const c: Vec3P = vc.sub(&ray.origin);
var cybx: f64 = @as(f64, @floatCast(cy)) * @as(f64, @floatCast(bx));
u = @floatCast(cxby - cybx);
var axcy: f64 = @as(f64, @floatCast(ax)) * @as(f64, @floatCast(cy)); const ax: P = a[kx] - sx * a[kz];
var aycx: f64 = @as(f64, @floatCast(ay)) * @as(f64, @floatCast(cx)); const ay: P = a[ky] - sy * a[kz];
v = @floatCast(axcy - aycx); const bx: P = b[kx] - sx * b[kz];
const by: P = b[ky] - sy * b[kz];
const cx: P = c[kx] - sx * c[kz];
const cy: P = c[ky] - sy * c[kz];
var bxay: f64 = @as(f64, @floatCast(bx)) * @as(f64, @floatCast(ay)); var u: P = cx * by - cy * bx;
var byax: f64 = @as(f64, @floatCast(by)) * @as(f64, @floatCast(ax)); var v: P = ax * cy - ay * cx;
v = @floatCast(bxay - byax); var w: P = bx * ay - by * ax;
}
if (backface_culling) { // Double precision fallback
if (u < 0.0 or v < 0.0 or w < 0.0) if (u == 0.0 or v == 0.0 or w == 0.0) {
return null; // no hit const cxby: PP = @as(PP, @floatCast(cx)) *
} else { @as(PP, @floatCast(by));
if ((u < 0.0 or v < 0.0 or w < 0.0) and var cybx: PP = @as(PP, @floatCast(cy)) *
(u > 0.0 or v > 0.0 or w > 0.0)) @as(PP, @floatCast(bx));
return null; // no hit u = @floatCast(cxby - cybx);
}
var det: f32 = u + v + w; var axcy: PP = @as(PP, @floatCast(ax)) *
if (det == 0.0) @as(PP, @floatCast(cy));
return null; // no hit var aycx: PP = @as(PP, @floatCast(ay)) *
@as(PP, @floatCast(cx));
v = @floatCast(axcy - aycx);
// Calculate scaled z-coordinates of vertices and use them to calculate var bxay: PP = @as(PP, @floatCast(bx)) *
// the hit distance @as(PP, @floatCast(ay));
const az: f32 = sz * a[kz]; var byax: PP = @as(PP, @floatCast(by)) *
const bz: f32 = sz * b[kz]; @as(PP, @floatCast(ax));
const cz: f32 = sz * c[kz]; v = @floatCast(bxay - byax);
var t: f32 = u * az + v * bz + w * cz; }
// hit.t counts as a previous hit for backface culling, in which if (backface_culling) {
// case triangle behind will no longer be considered a hit if (u < 0.0 or v < 0.0 or w < 0.0)
var hit: RayHit = RayHit{ return null; // no hit
.u = undefined, } else {
.v = undefined, if ((u < 0.0 or v < 0.0 or w < 0.0) and
.w = undefined, (u > 0.0 or v > 0.0 or w > 0.0))
.t = std.math.inf(f32), return null; // no hit
}
var det: P = u + v + w;
if (det == 0.0)
return null; // no hit
// Calculate scaled z-coordinates of vertices and use them to calculate
// the hit distance
const az: P = sz * a[kz];
const bz: P = sz * b[kz];
const cz: P = sz * c[kz];
var t: P = u * az + v * bz + w * cz;
// hit.t counts as a previous hit for backface culling, in which
// case triangle behind will no longer be considered a hit
// Since Ray.Hit is represented by a Vec4, t is the last element
// of that vector
var hit: Hit = math.vec4(
undefined,
undefined,
undefined,
std.math.inf(f32),
);
if (backface_culling) {
if ((t < 0.0) or (t > hit.v[3] * det))
return null; // no hit
} else {
if (det < 0) {
t = -t;
det = -det;
}
if ((t < 0.0) or (t > hit.v[3] * det))
return null; // no hit
}
// Normalize u, v, w and t
const rcp_det = 1.0 / det;
hit.v[0] = u * rcp_det;
hit.v[1] = v * rcp_det;
hit.v[2] = w * rcp_det;
hit.v[3] = t * rcp_det;
return hit;
}
},
else => @compileError("Expected Vec3, Vec3h, or Vec3d, found '" ++
@typeName(Vec3P) ++ "'"),
}; };
};
}
if (backface_culling) { test "triangleIntersect_basic_frontface_bc_hit" {
if ((t < 0.0) or (t > hit.t * det))
return null; // no hit
} else {
if (det < 0) {
t = -t;
det = -det;
}
if ((t < 0.0) or (t > hit.t * det))
return null; // no hit
}
// Normalize u, v, w and t
const rcp_det = 1.0 / det;
hit.u = u * rcp_det;
hit.v = v * rcp_det;
hit.w = w * rcp_det;
hit.t = t * rcp_det;
return hit;
}
};
test "triIntersect_basic_frontface_bc_hit" {
const a: math.Vec3 = math.vec3(0, 0, 0); const a: math.Vec3 = math.vec3(0, 0, 0);
const b: math.Vec3 = math.vec3(1, 0, 0); const b: math.Vec3 = math.vec3(1, 0, 0);
const c: math.Vec3 = math.vec3(0, 1, 0); const c: math.Vec3 = math.vec3(0, 1, 0);
const ray0: Ray = Ray{ const ray0: math.Ray = math.Ray{
.origin = math.vec3(0.1, 0.1, 1), .origin = math.vec3(0.1, 0.1, 1),
.direction = math.vec3(0.1, 0.1, -1), .direction = math.vec3(0.1, 0.1, -1),
}; };
const result: RayHit = ray0.triangleIntersect( const result: math.Ray.Hit = ray0.triangleIntersect(
&a, &a,
&b, &b,
&c, &c,
@ -177,23 +207,23 @@ test "triIntersect_basic_frontface_bc_hit" {
const expected_u: f32 = 0.6; const expected_u: f32 = 0.6;
const expected_v: f32 = 0.2; const expected_v: f32 = 0.2;
const expected_w: f32 = 0.2; const expected_w: f32 = 0.2;
try testing.expect(f32, expected_t).eql(result.t); try testing.expect(f32, expected_u).eql(result.v[0]);
try testing.expect(f32, expected_u).eql(result.u); try testing.expect(f32, expected_v).eql(result.v[1]);
try testing.expect(f32, expected_v).eql(result.v); try testing.expect(f32, expected_w).eql(result.v[2]);
try testing.expect(f32, expected_w).eql(result.w); try testing.expect(f32, expected_t).eql(result.v[3]);
} }
test "triIntersect_basic_backface_no_bc_hit" { test "triangleIntersect_basic_backface_no_bc_hit" {
const a: math.Vec3 = math.vec3(0, 0, 0); const a: math.Vec3 = math.vec3(0, 0, 0);
const b: math.Vec3 = math.vec3(1, 0, 0); const b: math.Vec3 = math.vec3(1, 0, 0);
const c: math.Vec3 = math.vec3(0, 1, 0); const c: math.Vec3 = math.vec3(0, 1, 0);
const ray0: Ray = Ray{ const ray0: math.Ray = math.Ray{
.origin = math.vec3(0.1, 0.1, 1), .origin = math.vec3(0.1, 0.1, 1),
.direction = math.vec3(0.1, 0.1, -1), .direction = math.vec3(0.1, 0.1, -1),
}; };
// Reverse winding from previous test // Reverse winding from previous test
const result: RayHit = ray0.triangleIntersect( const result: math.Ray.Hit = ray0.triangleIntersect(
&a, &a,
&c, &c,
&b, &b,
@ -204,28 +234,28 @@ test "triIntersect_basic_backface_no_bc_hit" {
const expected_u: f32 = -0.6; const expected_u: f32 = -0.6;
const expected_v: f32 = -0.2; const expected_v: f32 = -0.2;
const expected_w: f32 = -0.2; const expected_w: f32 = -0.2;
try testing.expect(f32, expected_t).eql(result.t); try testing.expect(f32, expected_u).eql(result.v[0]);
try testing.expect(f32, expected_u).eql(result.u); try testing.expect(f32, expected_v).eql(result.v[1]);
try testing.expect(f32, expected_v).eql(result.v); try testing.expect(f32, expected_w).eql(result.v[2]);
try testing.expect(f32, expected_w).eql(result.w); try testing.expect(f32, expected_t).eql(result.v[3]);
} }
test "triIntersect_basic_backface_bc_miss" { test "triangleIntersect_basic_backface_bc_miss" {
const a: math.Vec3 = math.vec3(0, 0, 0); const a: math.Vec3 = math.vec3(0, 0, 0);
const b: math.Vec3 = math.vec3(1, 0, 0); const b: math.Vec3 = math.vec3(1, 0, 0);
const c: math.Vec3 = math.vec3(0, 1, 0); const c: math.Vec3 = math.vec3(0, 1, 0);
const ray0: Ray = Ray{ const ray0: math.Ray = math.Ray{
.origin = math.vec3(0.1, 0.1, 1), .origin = math.vec3(0.1, 0.1, 1),
.direction = math.vec3(0.1, 0.1, -1), .direction = math.vec3(0.1, 0.1, -1),
}; };
// Reverse winding from previous test // Reverse winding from previous test
const result: ?RayHit = ray0.triangleIntersect( const result: ?math.Ray.Hit = ray0.triangleIntersect(
&a, &a,
&c, &c,
&b, &b,
true, true,
); );
try testing.expect(?RayHit, null).eql(result); try testing.expect(?math.Ray.Hit, null).eql(result);
} }