# approximate_hausdorff_distance

## `approximate_hausdorff_distance(v1, f1, v2, f2, use_cpp=True)`

Approximate the Hausdorff distance between two triangle meshes in 3D, i.e. the maximum distance between any point on one mesh and the closest point on the other mesh.

d = max { d(pA,B), d(pB,A) }

where pA is a point on mesh A and pB is a point on mesh B, and d(pA,B) is the distance between pA and the closest point on B. Our approximation will instead compute

d = max { d(vA,B), d(vB,A) }

where vA is a vertex on mesh A and vB is a vertex on mesh B.

Parameters:

Name Type Description Default
`v1` `(n1,3) array`

Vertices of first mesh.

required
`f1` `(m1,3) array`

Faces of first mesh.

required
`v2` `(n2,3) array`

Vertices of second mesh.

required
`f2` `(m2,3) array`

Faces of second mesh.

required
`use_cpp` `bool, optional (default`

Whether to use the C++ implementation of triangle_triangle_distance.

`True`

Returns:

Name Type Description
`d` `float`

Minimum distance value.

#### Notes

If you want an exact Hausdorff distance, you can heavily upsample both meshes using gpytoolbox.subdivide(method='upsample'). The approximated Hausdorff distance returned by this function will converge to the exact Hausdorff distance as the number of iterations of subdivision goes to infinity.

Examples:

``````# meshes in v,f and u,g
# Minimum distance value
d = gpytoolbox.minimum_distance(v,f,u,g)
``````
Source code in `src/gpytoolbox/approximate_hausdorff_distance.py`
 ``` 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136``` ``````def approximate_hausdorff_distance(v1,f1,v2,f2,use_cpp=True): """ Approximate the Hausdorff distance between two triangle meshes in 3D, i.e. the maximum distance between any point on one mesh and the closest point on the other mesh. d = max { d(pA,B), d(pB,A) } where pA is a point on mesh A and pB is a point on mesh B, and d(pA,B) is the distance between pA and the closest point on B. Our approximation will instead compute d = max { d(vA,B), d(vB,A) } where vA is a vertex on mesh A and vB is a vertex on mesh B. Parameters ---------- v1 : (n1,3) array Vertices of first mesh. f1 : (m1,3) array Faces of first mesh. v2 : (n2,3) array Vertices of second mesh. f2 : (m2,3) array Faces of second mesh. use_cpp : bool, optional (default: True) Whether to use the C++ implementation of triangle_triangle_distance. Returns ------- d : float Minimum distance value. Notes ----- If you want an **exact** Hausdorff distance, you can heavily upsample both meshes using gpytoolbox.subdivide(method='upsample'). The approximated Hausdorff distance returned by this function will converge to the exact Hausdorff distance as the number of iterations of subdivision goes to infinity. Examples -------- ```python # meshes in v,f and u,g # Minimum distance value d = gpytoolbox.minimum_distance(v,f,u,g) ``` """ # cpp implementation if use_cpp: try: from gpytoolbox_bindings import _hausdorff_distance_cpp_impl except: raise ImportError("Gpytoolbox cannot import its C++ fast winding number binding.") return _hausdorff_distance_cpp_impl(v1,f1,v2,f2) # Let's start by computing the one-sided distance, i.e., max(d(vA,B)). We will do this with a loop, but with pre-computed trees for efficient queriying. dim = v1.shape[1] # Initialize AABB tree for mesh 1 C1,W1,CH1,PAR1,D1,tri_ind1,_ = initialize_aabbtree(v1,f1) # Initialize AABB tree for mesh 2 C2,W2,CH2,PAR2,D2,tri_ind2,_ = initialize_aabbtree(v2,f2) # print("Computing one-sided distance...") current_best_guess_hd = 0.0 for i in range(v1.shape[0]): # print("Vertex %d of %d" % (i+1,v1.shape[0])) current_best_guess_dviB = np.Inf # current best guess for d(vi,B) queue = [0] while (len(queue)>0 and current_best_guess_dviB>current_best_guess_hd): q2 = queue.pop() # print("-----------") # print("Queue length : {}".format(len(queue))) # print("q1: ",q1) # print("q2: ",q2) # print("CH1[q1,]: ",CH1[q1,:]) # print("CH2[q2,]: ",CH2[q2,:]) # print("current_best_guess: ",current_best_guess) is_leaf2 = (CH2[q2,1]==-1) if (is_leaf2): # Compute distance between vi and triangle q2 d = np.sqrt(squared_distance_to_element(v1[i,:],v2,f2[tri_ind2[q2],:])[0]) if d < current_best_guess_dviB: current_best_guess_dviB = d else: # Compute distance between vi and bounding box of q2 # Distance from vi to the bounding box of q2 d = point_to_box_distance(v1[i,:],C2[q2,:],W2[q2,:]) d_max = point_to_box_max_distance(v1[i,:],C2[q2,:],W2[q2,:]) # print(d_max) # If d_max is smaller than the current best guess for HD, then we don't need to pursue this part of the tree if ((d < current_best_guess_dviB) and (d_max > current_best_guess_hd)): # Add children to queue queue.append(CH2[q2,0]) queue.append(CH2[q2,1]) # We have computed d(vi,B). Does it change our guess for hausdorff? if current_best_guess_dviB > current_best_guess_hd: current_best_guess_hd = current_best_guess_dviB # print("One-sided distance: {}".format(current_best_guess_hd)) # print("Computing two-sided distance...") # Now we do the other side, i.e., max(d(vB,A)) for i in range(v2.shape[0]): # print("Vertex %d of %d" % (i+1,v2.shape[0])) current_best_guess_dviA = np.Inf queue = [0] while (len(queue)>0 and current_best_guess_dviA>current_best_guess_hd): q2 = queue.pop() is_leaf2 = (CH1[q2,1]==-1) if (is_leaf2): # Compute distance between vi and triangle q2 d = np.sqrt(squared_distance_to_element(v2[i,:],v1,f1[tri_ind1[q2],:])[0]) if d < current_best_guess_dviA: current_best_guess_dviA = d else: # Compute distance between vi and bounding box of q2 # Distance from vi to the bounding box of q2 d = point_to_box_distance(v2[i,:],C1[q2,:],W1[q2,:]) d_max = point_to_box_max_distance(v2[i,:],C1[q2,:],W1[q2,:]) # If d_max is smaller than the current best guess for HD, then we don't need to pursue this part of the tree if (d < current_best_guess_dviA and d_max > current_best_guess_hd): # Add children to queue queue.append(CH1[q2,0]) queue.append(CH1[q2,1]) # We have computed d(vi,A). Does it change our guess for hausdorff? if current_best_guess_dviA > current_best_guess_hd: current_best_guess_hd = current_best_guess_dviA # We are done! return current_best_guess_hd ``````

## `point_to_box_distance(p, C, W)`

Compute the distance between a point and a box.

Parameters:

Name Type Description Default
`p` `(3,) array`

Point.

required
`C` `(3,) array`

Center of box.

required
`W` `(3,) array`

Width of box.

required

Returns:

Name Type Description
`d` `float`

Distance between point and box.

Source code in `src/gpytoolbox/approximate_hausdorff_distance.py`
 ```167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191``` ``````def point_to_box_distance(p,C,W): """ Compute the distance between a point and a box. Parameters ---------- p : (3,) array Point. C : (3,) array Center of box. W : (3,) array Width of box. Returns ------- d : float Distance between point and box. """ d = 0.0 for i in range(3): if p[i] < C[i]-0.5*W[i]: d += (p[i]-C[i]+0.5*W[i])**2 elif p[i] > C[i]+0.5*W[i]: d += (p[i]-C[i]-0.5*W[i])**2 return np.sqrt(d) ``````

## `point_to_box_max_distance(p, C, W)`

Compute the maximum distance between a point and a box.

Parameters:

Name Type Description Default
`p` `(3,) array`

Point.

required
`C` `(3,) array`

Center of box.

required
`W` `(3,) array`

Width of box.

required

Returns:

Name Type Description
`d` `float`

Maximum distance between point and box.

Source code in `src/gpytoolbox/approximate_hausdorff_distance.py`
 ```138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164``` ``````def point_to_box_max_distance(p,C,W): """ Compute the maximum distance between a point and a box. Parameters ---------- p : (3,) array Point. C : (3,) array Center of box. W : (3,) array Width of box. Returns ------- d : float Maximum distance between point and box. """ d = 0.0 for i in range(3): if p[i] < C[i]: # Distance to oppoisite corner d += (p[i]-C[i]-0.5*W[i])**2 elif p[i] > C[i]: # Distance to oppoisite corner d += (p[i]-C[i]+0.5*W[i])**2 return np.sqrt(d) ``````