Skip to content

Commit 6dbb136

Browse files
treegardeltreegardel
committed
Fix toroid phantom collisions in non-sequential tracing (#467)
The Newton solver in getToroidCollision converges to within NEW_TOLERANCE = 1e-4 mm. After a hit, the ray is nudged forward by COLLISION_EPSILON = 1e-6 mm — much less than the Newton residual. The forward check tested only the sign, so a solution within the residual on the wrong side of the new origin still passed and each toroid bounce was counted ~4x in non-sequential tracing. Reject hits inside Newton's own uncertainty by checking parametric forward distance, local to the toroid solver: const double t = dot(rayToHitpoint, rayDirection) / dot(rayDirection, rayDirection); if (t <= NEW_TOLERANCE) return std::nullopt; COLLISION_EPSILON stays at 1e-6 mm, all other surfaces keep their precision. Tests: - toroid_iteration integration test (grazing toroid vs. RAY-UI reference, regenerated at proper 1e-4 tolerance — old reference was at 1e-1 and hid the agreement). Now matches to ~1e-10 mm. - Three unit tests on getToroidCollision: forward and backward grazing hits at z ~ +/-600 mm walking a tolerance series from 1e-2 down to 1e-6 mm (1 nm), plus a phantom-rejection regression test. Newton's quadratic convergence delivers sub-nm precision in practice. - Re-enable ReflectionZonePlateDefault200Toroid with Sequential::Yes. Bonus: getCubicCollision had 29 pow(float(expr), N) calls truncating doubles to single precision on every Newton iteration, capping cubic-surface accuracy at ~1e-7 relative regardless of tolerance. Removed; all tests still pass. Also export getToroidCollision via RAYX_API so the test binary can call it directly.
1 parent 91978c1 commit 6dbb136

6 files changed

Lines changed: 720 additions & 35 deletions

File tree

‎Intern/rayx-core/src/Shader/Collision.cpp‎

Lines changed: 34 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -184,20 +184,20 @@ OptCollisionPoint getCubicCollision(const glm::dvec3& __restrict rayPosition, co
184184
y1 = y - aml * (x - xx);
185185
z1 = z - anl * (x - xx);
186186

187-
double func = (2 * ((x1 - xx) * an - al * z1) * cu.m_a23 - (2 * cu.m_a24 + cu.m_b12 * pow((xx), 2.0)) * al +
187+
double func = (2 * ((x1 - xx) * an - al * z1) * cu.m_a23 - (2 * cu.m_a24 + cu.m_b12 * pow(xx, 2.0)) * al +
188188
((x1 - xx) * am - al * y1) * (cu.m_a22 + cu.m_b21 * xx)) *
189189
((x1 - xx) * am - al * y1);
190-
func = func + pow((((x1 - xx) * an - al * z1)), 2.0) * cu.m_a33;
191-
func = func - ((x1 - xx) * an - al * z1) * (2 * cu.m_a34 + cu.m_b13 * pow((xx), 2.0) * al + cu.m_a44 * pow((al), 2.0));
190+
func = func + pow(((x1 - xx) * an - al * z1), 2.0) * cu.m_a33;
191+
func = func - ((x1 - xx) * an - al * z1) * (2 * cu.m_a34 + cu.m_b13 * pow(xx, 2.0) * al + cu.m_a44 * pow(al, 2.0));
192192
func = (func -
193193
(2 * ((x1 - xx) * an - al * z1) * cu.m_a13 - (cu.m_a11 * xx + 2 * cu.m_a14) * al + 2 * ((x1 - xx) * am - al * y1) * cu.m_a12) *
194194
al * xx) *
195195
al;
196-
func = (func - (pow((((x1 - xx) * am - al * y1)), 2.0) * cu.m_b23 +
196+
func = (func - (pow(((x1 - xx) * am - al * y1), 2.0) * cu.m_b23 +
197197
((x1 - xx) * am - al * y1) * ((x1 - xx) * an - al * z1) * cu.m_b32 - ((x1 - xx) * an - al * z1) * al * cu.m_b31 * xx) *
198-
((x1 - xx) * an - al * z1) / pow((al), 3));
198+
((x1 - xx) * an - al * z1) / pow(al, 3));
199199

200-
double dfunc = (2 * ((x1 - xx) * an - al * z1) * cu.m_a23 - (2 * cu.m_a24 + cu.m_b12 * pow((xx), 2)) * al +
200+
double dfunc = (2 * ((x1 - xx) * an - al * z1) * cu.m_a23 - (2 * cu.m_a24 + cu.m_b12 * pow(xx, 2)) * al +
201201
((x1 - xx) * am - al * y1) * (cu.m_a22 + cu.m_b21 * xx)) *
202202
am;
203203
dfunc = dfunc - (2 * (cu.m_a12 * am + cu.m_a13 * an) + cu.m_a11 * al) * al * xx;
@@ -207,15 +207,15 @@ OptCollisionPoint getCubicCollision(const glm::dvec3& __restrict rayPosition, co
207207
dfunc * al + ((cu.m_a22 + cu.m_b21 * xx) * am + 2 * (cu.m_a23 * an + al * cu.m_b12 * xx) - ((x1 - xx) * am - al * y1) * cu.m_b21) *
208208
((x1 - xx) * am - al * y1);
209209
dfunc = (dfunc + 2 * ((x1 - xx) * an - al * z1) * (cu.m_a33 * an + al * cu.m_b13 * xx) -
210-
(2 * cu.m_a34 + cu.m_b13 * pow((xx), 2)) * al * an) *
210+
(2 * cu.m_a34 + cu.m_b13 * pow(xx, 2)) * al * an) *
211211
al;
212212
dfunc = (dfunc - ((((x1 - xx) * an - al * z1) * (al * cu.m_b31 + am * cu.m_b32) - al * an * cu.m_b31 * xx +
213213
((x1 - xx) * am - al * y1) * (2 * am * cu.m_b23 + an * cu.m_b32)) *
214214
((x1 - xx) * an - al * z1) +
215-
(pow((((x1 - xx) * am - al * y1)), 2) * cu.m_b23 +
215+
(pow(((x1 - xx) * am - al * y1), 2) * cu.m_b23 +
216216
((x1 - xx) * am - al * y1) * ((x1 - xx) * an - al * z1) * cu.m_b32 - ((x1 - xx) * an - al * z1) * al * cu.m_b31 * xx) *
217217
an));
218-
dfunc = dfunc / pow((al), 3);
218+
dfunc = dfunc / pow(al, 3);
219219

220220
if (glm::abs(dfunc) < 0.001) { dfunc = 0.001; }
221221

@@ -246,33 +246,33 @@ OptCollisionPoint getCubicCollision(const glm::dvec3& __restrict rayPosition, co
246246
((y1 - yy) * al - am * x1) * (cu.m_a11 + cu.m_b12 * yy)) *
247247
((y1 - yy) * al - am * x1);
248248
func = func + (((y1 - yy) * an - am * z1) * cu.m_a33 - 2 * (cu.m_a23 * yy + cu.m_a34) * am) * ((y1 - yy) * an - am * z1) +
249-
(2 * cu.m_a24 * yy + cu.m_a44 + cu.m_a22 * pow((yy), 2) * pow((am), 2));
249+
(2 * cu.m_a24 * yy + cu.m_a44 + cu.m_a22 * pow(yy, 2) * pow(am, 2));
250250
func = func * am +
251-
((((y1 - yy) * an - am * z1) * cu.m_b32 - am * cu.m_b23 * yy) * am * yy - pow(((y1 - yy) * al - am * x1), 2) * cu.m_b13) *
251+
((((y1 - yy) * an - am * z1) * cu.m_b32 - am * cu.m_b23 * yy) * am * yy - pow((y1 - yy) * al - am * x1, 2) * cu.m_b13) *
252252
((y1 - yy) * an - am * z1);
253-
func = func - (pow((((y1 - yy) * an - am * z1)), 2) * cu.m_b31 + pow((am), 2) * cu.m_b21 * pow((yy), 2)) *
253+
func = func - (pow(((y1 - yy) * an - am * z1), 2) * cu.m_b31 + pow(am, 2) * cu.m_b21 * pow(yy, 2)) *
254254
((y1 - yy) * al - am * x1);
255-
func = func / pow((am), 3);
255+
func = func / pow(am, 3);
256256

257-
double dfunc = (pow(((y1 - yy) * an - am * z1), 2) * cu.m_b31 + pow((am), 2) * cu.m_b21 * pow((yy), 2) * al +
258-
2 * (((y1 - yy) * an - am * z1) * an * cu.m_b31 - pow((am), 2)) * cu.m_b21 * yy) *
257+
double dfunc = (pow((y1 - yy) * an - am * z1, 2) * cu.m_b31 + pow(am, 2) * cu.m_b21 * pow(yy, 2) * al +
258+
2 * (((y1 - yy) * an - am * z1) * an * cu.m_b31 - pow(am, 2)) * cu.m_b21 * yy) *
259259
((y1 - yy) * al - am * x1);
260260
dfunc =
261261
dfunc -
262-
((((y1 - yy) * an - am * z1) * cu.m_b32 - am * cu.m_b23 * yy) * am * yy - pow(((y1 - yy) * al - am * x1), 2) * cu.m_b13) * an;
262+
((((y1 - yy) * an - am * z1) * cu.m_b32 - am * cu.m_b23 * yy) * am * yy - pow((y1 - yy) * al - am * x1, 2) * cu.m_b13) * an;
263263
dfunc = dfunc + (2 * ((y1 - yy) * al - am * x1) * al * cu.m_b13 - (am * cu.m_b23 + an * cu.m_b32) * am * yy +
264264
(((y1 - yy) * an - am * z1) * cu.m_b32 - am * cu.m_b23 * yy) * am) *
265265
((y1 - yy) * an - am * z1);
266266
dfunc = dfunc - (((cu.m_a11 + cu.m_b12 * yy) * al + 2 * (cu.m_a12 * am + cu.m_a13 * an) - ((y1 - yy) * al - am * x1) * cu.m_b12) *
267267
((y1 - yy) * al - am * x1) -
268-
2 * (cu.m_a22 * pow((am), 2) * yy + cu.m_a23 * pow((am), 2) * z1 - cu.m_a23 * am * an * y1 +
269-
2 * cu.m_a23 * am * an * yy + cu.m_a24 * pow((am), 2) + cu.m_a33 * am * an * z1 -
270-
cu.m_a33 * pow((an), 2) * y1 + cu.m_a33 * pow((an), 2) * yy + cu.m_a34 * am * an) +
268+
2 * (cu.m_a22 * pow(am, 2) * yy + cu.m_a23 * pow(am, 2) * z1 - cu.m_a23 * am * an * y1 +
269+
2 * cu.m_a23 * am * an * yy + cu.m_a24 * pow(am, 2) + cu.m_a33 * am * an * z1 -
270+
cu.m_a33 * pow(an, 2) * y1 + cu.m_a33 * pow(an, 2) * yy + cu.m_a34 * am * an) +
271271
(2 * (((y1 - yy) * an - am * z1) * cu.m_a13 - (cu.m_a12 * yy + cu.m_a14) * am) +
272272
((y1 - yy) * al - am * x1) * (cu.m_a11 + cu.m_b12 * yy)) *
273273
al) *
274274
am;
275-
dfunc = dfunc / pow((am), 3);
275+
dfunc = dfunc / pow(am, 3);
276276

277277
if (glm::abs(dfunc) < 0.001) { dfunc = 0.001; }
278278

@@ -302,29 +302,29 @@ OptCollisionPoint getCubicCollision(const glm::dvec3& __restrict rayPosition, co
302302
double func = ((2 * (((z1 - zz) * am - an * y1) * cu.m_a12 - (cu.m_a13 * zz + cu.m_a14) * an) + ((z1 - zz) * al - an * x1) * cu.m_a11) *
303303
((z1 - zz) * al - an * x1) +
304304
(((z1 - zz) * am - an * y1) * cu.m_a22 - 2 * (cu.m_a23 * zz + cu.m_a24) * an) * ((z1 - zz) * am - an * y1) +
305-
(2 * cu.m_a34 * zz + cu.m_a44 + cu.m_a33 * pow((zz), 2)) * pow((an), 2)) *
305+
(2 * cu.m_a34 * zz + cu.m_a44 + cu.m_a33 * pow(zz, 2)) * pow(an, 2)) *
306306
an;
307-
func = func - ((((z1 - zz) * am - an * y1) * cu.m_b12 - an * cu.m_b13 * zz) * pow((((z1 - zz) * al - an * x1)), 2) -
307+
func = func - ((((z1 - zz) * am - an * y1) * cu.m_b12 - an * cu.m_b13 * zz) * pow(((z1 - zz) * al - an * x1), 2) -
308308
(((z1 - zz) * am - an * y1) * cu.m_b23 - an * cu.m_b32 * zz) * ((z1 - zz) * am - an * y1) * an * zz +
309-
(pow((((z1 - zz) * am - an * y1)), 2) * cu.m_b21 + pow((an), 2) * cu.m_b31 * pow((zz), 2)) *
309+
(pow(((z1 - zz) * am - an * y1), 2) * cu.m_b21 + pow(an, 2) * cu.m_b31 * pow(zz, 2)) *
310310
((z1 - zz) * al - an * x1));
311-
func = func / pow((an), 3);
311+
func = func / pow(an, 3);
312312

313313
double dfunc = (((z1 - zz) * am - an * y1) * cu.m_a22 - 2 * (cu.m_a23 * zz + cu.m_a24) * an) * am +
314314
(2 * (cu.m_a12 * am + cu.m_a13 * an) + cu.m_a11 * al) * ((z1 - zz) * al - an * x1);
315-
dfunc = dfunc + ((z1 - zz) * am - an * y1) * (cu.m_a22 * am + 2 * cu.m_a23 * an) - 2 * (cu.m_a33 * zz + cu.m_a34) * pow((an), 2);
315+
dfunc = dfunc + ((z1 - zz) * am - an * y1) * (cu.m_a22 * am + 2 * cu.m_a23 * an) - 2 * (cu.m_a33 * zz + cu.m_a34) * pow(an, 2);
316316
dfunc = (dfunc +
317317
(2 * (((z1 - zz) * am - an * y1) * cu.m_a12 - (cu.m_a13 * zz + cu.m_a14) * an) + ((z1 - zz) * al - an * x1) * cu.m_a11) * al) *
318318
an;
319319
dfunc = dfunc - (2 *
320-
(((z1 - zz) * am - an * y1) * am * cu.m_b21 - pow((an), 2) * cu.m_b31 * zz +
320+
(((z1 - zz) * am - an * y1) * am * cu.m_b21 - pow(an, 2) * cu.m_b31 * zz +
321321
(((z1 - zz) * am - an * y1) * cu.m_b12 - an * cu.m_b13 * zz) * al) *
322322
((z1 - zz) * al - an * x1) +
323-
(pow((((z1 - zz) * am - an * y1)), 2) * cu.m_b21 + pow((an), 2) * cu.m_b31 * pow((zz), 2)) * al +
324-
pow((((z1 - zz) * al - an * x1)), 2) * (am * cu.m_b12 + an * cu.m_b13) -
323+
(pow(((z1 - zz) * am - an * y1), 2) * cu.m_b21 + pow(an, 2) * cu.m_b31 * pow(zz, 2)) * al +
324+
pow(((z1 - zz) * al - an * x1), 2) * (am * cu.m_b12 + an * cu.m_b13) -
325325
((z1 - zz) * am - an * y1) * (am * cu.m_b23 + an * cu.m_b32) * an * zz +
326326
(((z1 - zz) * am - an * y1) * cu.m_b23 - an * cu.m_b32 * zz) * (am * z1 - 2 * am * zz - an * y1) * an);
327-
dfunc = (-dfunc) / pow((an), 3);
327+
dfunc = (-dfunc) / pow(an, 3);
328328

329329
if (glm::abs(dfunc) < 0.001) { dfunc = 0.001; }
330330

@@ -426,7 +426,11 @@ OptCollisionPoint getToroidCollision(const glm::dvec3& __restrict rayPosition, c
426426
// Note that multiplying the rays direction with -1 SHOULD totally have an effect on the collision detection - in most cases this 180° rotation
427427
// will make the ray point away from the toroid, and hence preventing a Collision completely. The above code however, is unaffected when
428428
// multiplying the ray direction with -1. Due to it having no effect on `glm::dvec3 normalized_dir = glm::dvec3(rayDirection) / rayDirection.z;`
429-
if (dot(rayToHitpoint, rayDirection) <= 0.0) return std::nullopt;
429+
//
430+
// Use NEW_TOLERANCE as minimum forward distance: a sign-only check (> 0) accepts solutions that fall within the Newton residual on
431+
// the wrong side of the ray origin, producing phantom self-intersections. Requiring t > NEW_TOLERANCE rejects exactly those cases.
432+
const double t = dot(rayToHitpoint, rayDirection) / dot(rayDirection, rayDirection);
433+
if (t <= NEW_TOLERANCE) return std::nullopt;
430434

431435
return col;
432436
}

‎Intern/rayx-core/src/Shader/Collision.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,8 @@ RAYX_FN_ACC OptCollisionPoint getQuadricCollision(const glm::dvec3& __restrict r
3030
RAYX_FN_ACC OptCollisionPoint getCubicCollision(const glm::dvec3& __restrict rayPosition, const glm::dvec3& __restrict rayDirection,
3131
const Surface::Cubic& __restrict cu);
3232

33-
RAYX_FN_ACC OptCollisionPoint getToroidCollision(const glm::dvec3& __restrict rayPosition, const glm::dvec3& __restrict rayDirection,
34-
const Surface::Toroid& __restrict toroid, bool isTriangul);
33+
RAYX_FN_ACC OptCollisionPoint RAYX_API getToroidCollision(const glm::dvec3& __restrict rayPosition, const glm::dvec3& __restrict rayDirection,
34+
const Surface::Toroid& __restrict toroid, bool isTriangul);
3535

3636
RAYX_FN_ACC OptCollisionPoint RAYX_API findCollisionInElementCoordsWithoutSlopeError(const glm::dvec3& __restrict rayPosition,
3737
const glm::dvec3& __restrict rayDirection,

0 commit comments

Comments
 (0)