5. Batched queries

A grid asks for the distance of every one of its points, not of five. distance_many takes the whole batch in one call. This chapter runs it on 16³ points over the domain, on both kinds of AABB tree and with both device-ready sign algorithms.

use fossil, only : surface_stl_object, AABB_TREE_OCTREE, SIGN_FAST_WINDING_NUMBER, STATUS_OK
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
integer(I4P), parameter       :: n = 16                      ! points per direction
type(surface_stl_object)      :: bvh, octree
type(vector_R8P), allocatable :: points(:)
real(R8P),        allocatable :: d_bvh(:), d_oct(:), d_wn(:)
real(R8P)                     :: h
integer(I4P)                  :: i, j, k, p, status

allocate(points(n**3), d_bvh(n**3), d_oct(n**3), d_wn(n**3))
h = 1.2_R8P / (n - 1)                                        ! the cube [-0.6, 0.6]^3
p = 0
do k = 1, n ; do j = 1, n ; do i = 1, n
   p = p + 1
   points(p) = (-0.6_R8P + (i - 1) * h) * ex_R8P + (-0.6_R8P + (j - 1) * h) * ey_R8P + (-0.6_R8P + (k - 1) * h) * ez_R8P
enddo ; enddo ; enddo

call bvh%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)       ! default AABB_TREE_SAH_BVH
call octree%load_from_file(file_name='dragon-placed.stl', guess_format=.true., aabb_tree_kind=AABB_TREE_OCTREE)
call bvh%distance_many(points=points, distances=d_bvh, is_signed=.true., is_square_root=.true., status=status)
if (status /= STATUS_OK) error stop 'distance_many failed'
call octree%distance_many(points=points, distances=d_oct, is_signed=.true., is_square_root=.true.)
print '(A,I0)',     'points                ', size(points)
print '(A,ES10.3)', 'max ||bvh| - |oct||  ', maxval(abs(abs(d_bvh) - abs(d_oct)))
print '(A,I0)',     'bvh, oct signs differ ', count(d_bvh * d_oct < 0._R8P)

• distance_many(points=, distances=) is a loop of compute_distance over points, with the same is_signed, sign_algorithm and is_square_root for every point (unsigned and squared by default, as distance). distances is allocated by the caller, with the size of points; a size mismatch gives status = STATUS_INVALID_INPUT and leaves distances untouched.
• The loop is an OpenMP parallel do: the tree is read-only and each iteration writes its own slot. Built with -fopenmp it runs on all threads, without it serially, with the same results.
• aabb_tree_kind= on load_from_file (also on analyze and adopt_facets) selects the tree: AABB_TREE_SAH_BVH (the default), a binary tree split by the surface-area heuristic, or AABB_TREE_OCTREE, the older 8-way space partition. The octree is kept for comparison; the BVH is the faster one and the only one with a device path (chapter 12).

call bvh%distance_many(points=points, distances=d_wn, is_signed=.true., is_square_root=.true., &
                       sign_algorithm=SIGN_FAST_WINDING_NUMBER)
print '(A,I5,A,F10.6)', 'pseudo-normal:  inside ', count(d_bvh < 0._R8P), '  min d ', minval(d_bvh)
print '(A,I5,A,F10.6)', 'winding number: inside ', count(d_wn  < 0._R8P), '  min d ', minval(d_wn)
print '(A,I0,A,ES10.3)', 'signs differ at ', count(d_bvh * d_wn < 0._R8P), ' points, max ||d_pn| - |d_wn|| ', &
                        maxval(abs(abs(d_bvh) - abs(d_wn)))

Running it

$ tut05_batched
points                4096
max ||bvh| - |oct||   0.000E+00
bvh, oct signs differ 1
pseudo-normal:  inside   152  min d  -0.647299
winding number: inside   108  min d  -0.069887
signs differ at 50 points, max ||d_pn| - |d_wn||  0.000E+00

• The two trees give the same distance magnitudes bit for bit: the traversal never prunes the true closest facet. One sign differs. At that point several facets are equally close, the trees report different ones, and on a non-manifold surface their pseudo-normals disagree.
• The two sign algorithms give the same magnitudes and different signs at 50 of the 4096 points. The pseudo-normal sign puts 152 points inside, one of them at d = -0.647, deeper than the dragon is wide: a sign read from the closest feature is wrong where the surface is not a closed manifold. The fast winding number puts 108 points inside, none deeper than 0.070. On a surface that is not closed, use SIGN_FAST_WINDING_NUMBER; the rest of the tutorial does.

Most ember points of the left figure lie far from the body, near the corners of the domain. Each of the 47 takes its sign from the pseudo-normal of a vertex or an edge of the dragon, which is still neither closed nor manifold after chapter 2; one wrong pseudo-normal flips every point whose closest feature it is.

The timing is not printed: it depends on the machine. The figures are in GPU distance queries (the BVH is about 110 times faster than the octree on a mesh of 24k facets).

What you learned

distance_many answers a whole grid in one OpenMP call; the SAH BVH and the octree agree on the distances; the pseudo-normal sign needs a closed surface, the fast winding number does not. Reference: surface_stl_object.

Next: 6. A distance field.