Skip to content

Prune shadowed neighbours in alpha-shape search and add FastInSphereTesterSoS - #6566

Merged
Fedr merged 24 commits into
masterfrom
alpha-shape-triangle-shadow
Aug 18, 2026
Merged

Prune shadowed neighbours in alpha-shape search and add FastInSphereTesterSoS#6566
Fedr merged 24 commits into
masterfrom
alpha-shape-triangle-shadow

Conversation

@Fedr

@Fedr Fedr commented Aug 8, 2026

Copy link
Copy Markdown
Contributor

A follow-up to #6551. There the filter was pairwise: a neighbour p shadows a farther neighbour x
when p is strictly inside every ball of the given radius via #v that contains x. Here the same
idea is applied to a pair: as soon as the two balls of the given radius touching #v, p and q are
known to exist, they shadow a whole wedge, and every neighbour inside it is dropped from the not yet
visited tail of neis.

The output is bit-identical to master, the win is fewer pairs in the quadratic loops and shorter
ballEmpty scans.

The conditions

All relative to the origin at #v. Let p, q be the two other points of the pair, c and
cbar the centers of the two radius-r balls touching 0, p, q (mirror images of each other in the
triangle's plane), and n the doubled normal of the triangle. With pp = p*p, qq = q*q, pq = p*q
and, for a candidate x, pd = p*x, qd = q*x, nd = n*x, xx = x*x:

W = n*n
M = qq*(pp-pq)*p + pp*(qq-pq)*q = 2*W*( circumcenter - 0 )
E = 4*rSq*W^2 - M*M
S = sqrt( E*W ) >= 0

W, M, E are exactly the quantities InSphereTester<int>::reset already computed (E does not
depend on which point of the triangle is taken as the origin, since that point is on the circumcircle,
so |M| = 2*W*circumradius either way), and the two centers satisfy exactly

2*W^2*c = W*M + S*n     2*W^2*cbar = W*M - S*n

x is dropped iff all three hold, strictly:

  1. on q's side of both planes (0, p, c) and (0, p, cbar), which is S*bp > |cp*nd| with
    bp = pp*qd - pq*pd and cp = W*pp*(qq-pq), i.e. bp > 0 && bp^2*E*W > (cp*nd)^2;
  2. symmetrically on p's side of both planes (0, q, c) and (0, q, cbar):
    bq > 0 && bq^2*E*W > (cq*nd)^2 with bq = qq*pd - pq*qd, cq = W*qq*(pp-pq);
  3. strictly outside both balls: W*(W*xx - M*x) > S*|nd|, which in the notation of
    InSphereTester<int>::operator() is A*W > S*|t|, so it reuses the tester's own arithmetic.

Substituting the center identity into cross(p, c)*x gives 1 for the c plane, and S -> -S gives the
cbar plane; the two differ only in the sign of the right-hand side, so requiring both is the single
comparison with |cp*nd| above — no sign analysis and no case split anywhere. Writing
x = alpha*p + beta*q + gamma*c, conditions 1 and 2 say alpha, beta > 0, so the shadow is the union of
the cones spanned by (p, q, c) and (p, q, cbar) together with the in-plane sector between p and q,
minus the two balls. Since it is symmetric in c and cbar, it does not depend on which side's ball was
found empty, and one pass per pair of neighbours is enough.

All arithmetic is exact integer: bp, bq fit Int256, cp, cq fit Int512, and the squared
comparisons stay below 2^705, so Int1024 is enough — the same tier the SoS path uses. No division and
no square root anywhere.

Why dropping is safe

For x satisfying 1-3, every radius-r ball via #v containing x inside or on its boundary contains
p or q strictly inside. Hence x can neither make a triangle with #v (both balls of any
triangle (#v, x, y) are blocked) nor be the last point blocking a ball of others, and the strictness
makes it immune to the simulation-of-simplicity tie resolution — exactly as in #6551. Note that p and
q are on both balls, so condition 3 never removes them, and the strictly-inside point the claim
provides is never one of the two points ballEmpty skips by id.

The centers of the radius-r balls via #v that contain neither p nor q strictly form a region on
the sphere |c'| = r whose two corners are exactly c and cbar; conditions 1 and 2 place x in the
normal cone at a corner, and condition 3 says the corner value is negative. The cbar half of the region
also follows from the c half by the reflection in the triangle's plane, which fixes 0, p, q, swaps the
two balls and maps balls via #v to balls via #v.

Condition 3 must be strict for the same reason condition 2 of #6551 must be: for x exactly on a ball
the touching ball itself leaves p and q on the boundary, and 0, p, q, x are cospherical, so
simulation-of-simplicity decides the tie by ids and x may well be needed.

Dropping composes: a point that shadowed x may itself be dropped later, but only by points present at
that later moment, and the removal times strictly increase, so the chain always ends at a point that is
still in neis when the ball is tested.

Verified numerically (scripts are re-derivable from the formulas above):

  • 183k random configurations: the integer conditions agree with the floating-point geometry
    (sides of the four planes and distances to the two explicitly constructed centers) in every
    non-degenerate case;
  • 3000 configurations passing 1-3, each swept over 300k balls of radius r via #v containing x:
    no violation, worst relative margin 1.4e-4; half of them are not caught by the one-sided form;
  • 1500 fully coplanar configurations (the crossing-grids case) swept the same way: no violation,
    worst relative margin 1.4e-4;
  • controls with exactly one condition violated: 98-100% of the configurations violate the claim, so
    every one of the three is load-bearing.

Benchmarks

VS2026 x64 Release, Ryzen 3900X, master and this branch interleaved, 18 samples each (6 for the
cloud), median. The hash of the whole triangle list is identical to master everywhere.

test master this PR consideredTris inBallTests
Fibonacci sphere 20k, r = 1.5*sqrt(4*pi/n) 33.0 ms 17.7 ms (-46%) 265080 -> 117965 (-56%) 1581652 -> 894344 (-43%)
Fibonacci sphere 50k, same r 71.7 ms 38.8 ms (-46%) 561398 -> 296901 (-47%) 3453163 -> 2185402 (-37%)
two crossing 61x61 grids, step 0.05, r = 0.1 49.4 ms 37.7 ms (-24%) 260811 -> 109771 (-58%) 1378089 -> 727343 (-47%)
a scanned cloud, 18820 points, r = 0.5 10024 ms 4460 ms (-56%) 57380489 -> 7754243 (-86%) 261200891 -> 41810772 (-84%)

The scanned cloud is the case the search is actually used on, and the one with resolution to spare: the
three synthetic tests now run in tens of milliseconds, where the run-to-run spread is wider than most of
the effects below, so every isolated measurement here is taken on the cloud.

Both floating-point pre-filters still earn their place after the exact arithmetic moved to FastInt,
which is worth re-checking since that made the integers several times faster. Removing either, on the
cloud, median of 6:

ms cost
both filters 4430 --
without the wedge filter 5748 +30%
without the query filter of the tester 6240 +41%

Casting the shadow from every touchable pair rather than only from the pairs that yield a triangle is
also still right: 5006 ms against 7801 ms for the latter, which does 5.6x more ball tests (235M against
41.8M) and 6.6x more pairs. That balance depends on the ball tests being cheap, so it moves with the
arithmetic -- it was the other way round before the pre-filters were restored.

The commits are worth comparing on the crossing grids: with the one-sided wedge of the first commit that
case dropped nothing at all (consideredTris -0.6%) and cost 2.4%, because the cloud is piecewise planar
and every candidate lies exactly in the triangle's plane, on the boundary of the one-sided wedge. The
second commit is both simpler and strictly stronger.

The remaining three commits are independent of the geometry:

  • swapping the two loops so that the farther point of the pair is taken outside. neis is sorted by
    distance, so every pair of the closest neighbours is then processed before any farther point is
    touched; those pairs have the smallest touching balls and hence the largest shadows, and what those
    shadows drop is exactly what the outer loop has not reached yet. Free, and worth another -11% / -5% /
    -3% of consideredTris.
  • rejecting most candidates by the signs of bp and bq computed in double, falling back to the exact
    values only within the tolerance of zero. The operands are dot products of point differences below
    2^31, which Vector3d represents exactly, so they are recomputed in double rather than converted
    from the bignums -- 0.7 bits worse and no conversion needed. Their dot products with d stay below
    2^129 and the roundings keep the error below 2^79, so a double below -2^79 is negative exactly as
    well; the tolerance takes 32 times that margin. shadowedNeis is identical with and without the
    rejection on every test above, and the rejection skips 84-90% of the exact evaluations.
  • casting the shadow from every touchable pair instead of only from the pairs that give a triangle. The
    claim needs only the existence of the touching balls, never their emptiness, and dropping before the
    emptiness tests shortens their scans too: consideredTris -29% and inBallTests -22% on the crossing
    grids. On the spheres almost every touchable pair yields a triangle anyway, so there it only saves the
    scans (-6%).
  • FastInSphereTesterSoS, a new class next to InSphereTesterSoS giving exactly its answers but
    filtering the queries in floating point first. reset() also computes the sphere's center as
    M/(2W) in the plane of the triangle plus the height S*w/(2W^2) orthogonal to it -- the two cannot
    cancel each other, and flip() only negates the second -- and a query farther from the surface than
    rSq*2^-44 is then decided by one squared distance. That is sixteen times the error of the distance,
    which also means a point exactly on the sphere always reaches the tie resolution: on the integer grid
    one unit of movement changes the squared distance by about 2r, far more than the tolerance, so
    nothing but an exactly cospherical point can land inside it. The existing classes and their tests are
    untouched, and a new test checks the two agree over the cospherical square and thousands of
    pseudo-random spheres and queries at three magnitudes. Alone it is worth 34% / 36% / 16% of the whole
    search -- more than everything above together, since the emptiness tests are what dominates once the
    pairs have been cut down.

The shadow is close to the largest one possible for a pair of blockers, but not exactly it. The maximum
of 2c'x - xx over the centers c' of the radius-r balls via #v having neither p nor q strictly
inside can sit only in three places: the free maximizer r*x/|x| (there it is non-negative for any x
within the search radius, so x is not redundant), one of the two corners c, cbar (what the
conditions above test), or the relative interior of one of the two bounding arcs (there redundancy is
decided by a single blocker, which is #6551's case). Enumerating those exactly on 20000 random
configurations that are redundant:

share of the redundant configurations
caught by the conditions above 62.4%
caught by #6551 with p or with q 59.6%
caught by neither 0.09%

so 99.91% of the redundancy available from a pair of blockers is already taken. The remainder is 13 arc
cases, where #6551's predicate is sufficient but not necessary, and 4 corner cases, where the wedge is
slightly conservative because the normal cone at a corner of a spherical region is not exactly
alpha, beta > 0. Closing either would need more predicates per candidate than the 0.09% is worth.

Finally, AlphaShapeStats gained six counters, none of them measurable in the timings. shadowTests,
exactShadowTests and shadowedNeis, together with inBallTests, show what the shadow filter really
does -- it trades the expensive point-in-ball tests for cheap shadow tests, 84-90% of which end at two
double dot products:

test inBallTests, master inBallTests, here shadowTests of them exact shadowed point tests in total
Fibonacci sphere 20k 1581652 894344 (-43%) 381631 40211 33425 1581652 -> 1275975 (-19%)
Fibonacci sphere 50k 3453163 2185402 (-37%) 840270 84267 62742 3453163 -> 3025672 (-12%)
two crossing 61x61 grids 1378089 727343 (-47%) 604657 98659 22007 1378089 -> 1332000 (-3%)

On the crossing grids the total number of point tests is almost unchanged, so the whole gain there comes
from their cost, not from their number.

collectedNeis, redundancyTests and redundantNeis cover the one stage that had no counters at all,
the distance filter of #6551. On the scanned cloud it turns out to be 383M tests dropping 8.42M of the
10.73M neighbours found -- the same scale as the shadow filter's 387M tests, and invisible until now.
collectedNeis is also the size the quadratic loops start from, which is what makes the time per point
differ by two orders of magnitude between one cloud and another.

The API this adds

MRInSphere.h gains FastInSphereTesterSoS, with reset, flip, operator() and
outsideBothSpheres; InSphereTesterSoS gains the accessors normal(), normalSq() and heightSq(),
which describe the prepared sphere exactly and were protected before. All of them carry
MR_BIND_IGNORE, as everything mentioning FastInt does. InSphereTester<int>, InSphereTesterSoS
and their tests are otherwise untouched, and AlphaShapeData and the signature of
findAlphaShapeNeiTriangles do not change -- only AlphaShapeStats grows the counters above.

Fedr added 13 commits August 8, 2026 18:10
After a triangle (#v, p, q) is appended, every neighbour strictly outside the ball
just proven empty and strictly inside the wedge of the planes (#v, p, center) and
(#v, q, center) is dropped: each ball of the given radius via #v containing such a
point has p or q strictly inside, so it can neither make a triangle with #v nor be
the last blocker of a ball of others. The output is bit-identical.
The two touching balls of the triangle are mirror images of each other in its
plane, so the shadow is the wedge of four planes via #v, one of p and q, and one
of the two centers, minus both balls. Each pair of mirrored planes collapses into
a single squared comparison without any sign analysis, the filter no longer depends
on which side's ball was found empty, and one pass per pair replaces two.

The region is twice as large: on a piecewise planar cloud the previous one-sided
form dropped nothing at all, since every candidate is in the triangle's plane.
neis is sorted by distance from #v, so with the farther point outside every pair
of the closest neighbours is processed before any farther point is touched. Their
triangles have the smallest touching balls and thus the largest shadows, and the
points those shadows drop are exactly the ones the outer loop has not reached yet.

consideredTris -11% / -5% / -3% and inBallTests -7% / -3% / -2% on the spheres and
the crossing grids, with the found triangles unchanged.
…oint

The two dot products deciding the dihedral angle are computed in double first, and
the exact Int256 values only when the approximation is within the tolerance of zero,
which cannot change any verdict: the error of the dot products is below 2^78 here.

The scan of the neighbours is what the filter spends its time on, so this is worth
6..8% of the whole search on the spheres and the crossing grids.
…iangles

The claim behind the filter needs only the existence of the two balls touching #v,
p and q, never their emptiness, so the pair does not have to give a triangle. Doing
it before the emptiness tests shortens their scans as well.

consideredTris -29% and inBallTests -22% on the crossing grids, -0.4% and -6% on
the spheres, where almost every touchable pair yields a triangle anyway. The found
triangles are unchanged.
shadowTests are the neighbours tested for being shadowed, exactShadowTests those
of them the floating-point rejection could not decide, and shadowedNeis those
found shadowed and excluded from the search. The increments cost nothing
measurable, and they show the rejection stopping 84..90% of the tests:

  sphere 20k  381631 tested / 40211 exact / 33425 shadowed
  sphere 50k  840270 / 84267 / 62742
  grids 61x61 604657 / 98659 / 22007
pp and qq are neis[i].distSq and neis[j].distSq, which #6551 introduced exactly to
avoid recomputing them, and the tail is checked for being empty before anything is
prepared. Timing-neutral, the found triangles and all the counters are unchanged.
The chain now ends at the constant itself: coordinate differences below 2^31 make
pp and pq below 2^64, the components of up and uq below 2^96 and the dot products
below 2^129, whose roundings keep the error below 2^78, and the tolerance takes 64
times that. Measured over 2M configurations at the full coordinate range the worst
error is 2^74, so the constant has three decimal orders of margin in practice.
The new class answers exactly as InSphereTesterSoS, only sooner: reset() also computes
the sphere's center in double, and a query farther from the surface than the tolerance
is decided by one squared distance, leaving the exact predicates and the tie resolution
to the points near it. The center splits in the in-plane part M/(2W) and the height
S*w/(2W^2), which are orthogonal, so their sum loses nothing to cancellation and flip()
only negates the second. The tolerance is rSq*2^-44, sixteen times the error of that
distance, which also guarantees that a point exactly on the sphere always reaches the
simulation-of-simplicity - on the integer grid nothing else can come that close.

The existing classes and their tests are untouched; findAlphaShapeNeiTriangles takes
the new one, which is 34% / 36% / 16% of the whole search on the two Fibonacci spheres
and the crossing grids, with the found triangles unchanged.
…erSoS

The shadow filter had its own copy of the sphere arithmetic; as a method of the
tester it is filtered in floating point like the queries, and the same tolerance
argument covers it. Both centers are cc_ +- hn_, so the two squared distances come
from one subtraction each, and the exact fallback is the A * W > S * |t| of
operator() without any tie to resolve: a point on either sphere is not outside it.

Worth 5% on the Fibonacci spheres and 1.5% on the crossing grids, and the shadow
filter no longer duplicates the predicate.
master replaced boost multiprecision with FastInt<N>/VarBigInt, so the shadow filter and outsideBothSpheres no longer compiled on this branch. Exact path only; the floating-point pre-filters need a bignum->double that master dropped and are re-added separately once measured.
meshinspector-agent Bot and others added 9 commits August 16, 2026 20:47
Branch was 14 commits behind; master's FastInt rework (toDouble, MRFastInt.cpp,
doubleFromWords) is needed by the pending double pre-filters. File sets are
disjoint, so the merge is textual only.
The port to FastInt dropped the floating-point pre-filter because the family
has no conversion to double, so every candidate was evaluated exactly. The
filter's operands are dot products of point differences below 2^31, which
Vector3d represents exactly, so they can be recomputed in double directly
instead of converting the exact bignums - 0.7 bits worse, and no dependency
on a bignum-to-double conversion.

Co-authored-by: meshinspector-agent[bot] <meshinspector-agent[bot]@users.noreply.github.com>
* Add toDouble for the FastInt family

The pre-filters that reject a case in double before evaluating it exactly
need a bignum-to-double conversion, which the family lacked: MSVC's
std::_Signed128 has no such operator at all. It is correctly rounded, and
deliberately platform-independent rather than deferring to the built-in
__int128_t conversion, since the stated error bound is what makes those
rejections safe.

* Do not generate a C binding for toDouble( FastInt128 )

mrbind emitted MR_toDouble(__int128 v), which no build accepts: gcc rejects
__int128 in a header compiled with -pedantic-errors, and MSVC has no such
keyword at all. FastInt128 has no C spelling, and the rest of the family is
already ignored for the same reason.

* Restore the double query filter in FastInSphereTesterSoS

The center was computed exactly only, because the port to FastInt lost the
bignum-to-double conversion the floating-point pre-filter needs; toDouble brings
it back. W, M and E are converted rather than recomputed in double: E is a
difference of two values below 2^322 that cancels completely on exactly the
near-degenerate triangles the filter has to decide.

---------

Co-authored-by: meshinspector-agent[bot] <meshinspector-agent[bot]@users.noreply.github.com>
Both comments asserted a bound; now they show where it comes from, step by step,
including which term dominates - for the wedge filter it is the error of the double
dot products propagated through up and uq, not the roundings of the final dot
product, and for the tester it is the displacement of the center, of which the
conversion of w is a part since its components reach 2^63.

The tester comment also records why the tolerance is proportional to rSq instead of
absolute: that is what bounds the error below it on both sides of the surface, and
an absolute one fitted to the largest coordinates would decide nothing at small ones.

Both derivations were re-checked against the exact values: worst error 2^74.6 over
400k configurations for the wedge filter (bound 2^79) and rSq * 2^-50 over 60k
near-surface queries for the tester (bound rSq * 2^-48), no wrong decision.
The shadow filter reached w, W and E through a local class deriving from the tester
only to widen their access; as accessors of InSphereTesterSoS they serve any caller
needing its own exact predicate about the prepared sphere, and findAlphaShapeNeiTriangles
declares a plain FastInSphereTesterSoS. Their doc comments give the exact center
identity they belong to, and heightSq() states that it is sqr( 2 * h * normalSq() )
rather than h^2, which the member name E did not.

MR_BIND_IGNORE on all three, as on everything else mentioning FastInt.
collectedNeis, redundancyTests and redundantNeis were the only stage of the search
with no counters at all, although the filter of #6551 runs on every point and is
quadratic in its neighbours. On a scanned cloud of 18820 points at radius 0.5 it
turns out to be 383M tests dropping 8.42M of the 10.73M neighbours found - the same
scale as the shadow filter's 387M tests, and invisible until now.

collectedNeis also gives the size the quadratic loops start from, which is what
explains the time per point differing by two orders between clouds.
dropShadowed had grown to eighty lines of body and half as much comment inside
findAlphaShapeNeiTriangles, reached through a [&] capture of four locals. As a
function of an anonymous namespace it states what it needs - the neighbours it
prunes, the pair shadowing them, the central point, the tester already reset on
their two balls, and the counters - and the derivations read as its documentation
rather than as an aside in the middle of another function. insideWedge follows it,
having captured nothing in the first place.

Pure move: the body is identical after dedenting and renaming the two captures,
and on a scanned cloud every counter and the triangle hash are unchanged, with
the timings within noise of each other.
@Fedr Fedr changed the title Filter out the neighbours shadowed by a found alpha-shape triangle Prune shadowed neighbours in alpha-shape search and add FastInSphereTesterSoS Aug 18, 2026
The product is below 2^353 and so fits in six words, while its type is seven wide as
the sum of its factors'. mulWords skips the zero words of its left operand but runs
over every word of its right one, so narrowing first makes the square cost six words
by six. The left side stays 1024 bits, which is what its factors' widths add up to,
and the right one is widened to it for the comparison as before.

Six limb products out of the hundred a full evaluation takes, in the branch reached
only by a candidate passing both sign tests: not measurable on the benchmarks, the
counters and the triangle hashes are unchanged. It is the narrowest width the bound
allows, which is how the rest of the FastInt code is written.
@Fedr
Fedr merged commit 91a54df into master Aug 18, 2026
50 checks passed
@Fedr
Fedr deleted the alpha-shape-triangle-shadow branch August 18, 2026 10:02
Fedr added a commit that referenced this pull request Aug 18, 2026
Besides the trivial include conflict, FastInSphereTesterSoS from #6566 meets
the degenerate-triangle mode of this branch: there is no exact sphere center
to compute in floating point when W == 0, so reset() skips it, the query
delegates to the symbolic evaluation of the base class, and outsideBothSpheres
conservatively never prunes (an exact answer for the perturbed spheres would
need the id of the query point).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant