4. Signed distance

With the dragon placed in the domain, the first question a grid asks is how far each point is from the wall, and on which side. This chapter queries single points.

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

implicit none
type(surface_stl_object) :: surface
type(vector_R8P)         :: p, closest, normal
real(R8P)                :: d_pn, d_wn, d
integer(I4P)             :: k, facet, edge, vertex

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)
print '(A)', '    y      d (pseudo-n.)  d (winding n.) inside'
do k = 0, 4
   p = 0._R8P * ex_R8P + (-0.6_R8P + 0.3_R8P * k) * ey_R8P + 0._R8P * ez_R8P
   d_pn = surface%distance(point=p, is_signed=.true., is_square_root=.true.)
   d_wn = surface%distance(point=p, is_signed=.true., is_square_root=.true., sign_algorithm=SIGN_FAST_WINDING_NUMBER)
   print '(F7.3,2F15.6,L5)', p%y, d_pn, d_wn, surface%is_point_inside(point=p, sign_algorithm=SIGN_FAST_WINDING_NUMBER)
enddo

• distance(point=) returns the distance from a point to the closest facet, found by a best-first traversal of the AABB tree. The defaults are unsigned and squared: pass is_square_root=.true. for the distance itself and is_signed=.true. for a sign, negative inside.
• sign_algorithm= chooses how the sign is decided:
  • SIGN_PSEUDO_NORMAL (the default): the angle-weighted pseudo-normal of the closest face, edge or vertex, read in the same traversal. Exact on a closed, outward-oriented surface.
  • SIGN_FAST_WINDING_NUMBER: inside when the winding number of the surface around the point exceeds 0.5, evaluated hierarchically on the tree. It tolerates holes and gaps.
  • SIGN_RAY_INTERSECTIONS and SIGN_SOLID_ANGLE: ray parity and the solid-angle sum, O(N) per query.
• is_point_inside(point=, sign_algorithm=) returns the inside test alone, with the same algorithms and the same default.

p = 0._R8P * ex_R8P + 0.6_R8P * ey_R8P + 0._R8P * ez_R8P
call surface%compute_distance(point=p, distance=d, is_signed=.true., is_square_root=.true., &
                              facet_index=facet, edge_index=edge, vertex_index=vertex, closest_point=closest, normal=normal)
print '(A,F10.6,A,I0,A,I0,A,I0)', 'd ', d, '  facet ', facet, '  edge ', edge, '  vertex ', vertex
print '(A,3F10.6)', 'closest point ', closest%x, closest%y, closest%z
print '(A,3F10.6)', 'normal        ', normal%x, normal%y, normal%z

• compute_distance is the subroutine form with the details of the closest feature: facet_index, the local edge_index (1 to 3, 0 when not on an edge) and vertex_index (1 to 3, 0 when not at a vertex), the closest_point and the pseudo-normal normal of that feature. Chapter 8 gets the same for many points at once.

Running it

$ tut04_distance
    y      d (pseudo-n.)  d (winding n.) inside
 -0.600       0.150398       0.150398    F
 -0.300      -0.000289      -0.000289    T
  0.000      -0.015419      -0.015419    T
  0.300      -0.002330      -0.002330    T
  0.600       0.114892       0.114892    F
d   0.114892  facet 1289  edge 0  vertex 3
closest point  -0.017632  0.500000  0.053752
normal         -0.014229  0.824610 -0.565522

The five points on the y axis go from outside (0.150398) through the body to outside again (0.114892). The two signs agree at these points. The point (0, 0.6, 0) is closest to vertex 3 of its facet, so its normal is the vertex pseudo-normal, not the facet normal. On this open dragon the two sign algorithms do not agree everywhere: chapter 5 counts where they differ.

What you learned

distance is unsigned and squared unless asked otherwise; sign_algorithm picks the inside test; compute_distance adds the closest facet, feature, point and normal. Reference: Generalized winding number, Sign algorithms at sharp features.

Next: 5. Batched queries.