Skip to content

squared_distance

squared_distance(P, V, F=None, use_cpp=False, use_aabb=False, C=None, W=None, CH=None, tri_ind=None, split_dir=None)

Squared distances from a set of points in space.

General-purpose function which computes the squared distance from a set of points to a mesh, point cloud or polyline, in two or three dimensions. Optionally, uses an aabb tree for efficient computation.

Parameters:

Name Type Description Default
P (p,dim) numpy double array

Matrix of query point positions

required
V (v,dim) numpy double array

Matrix of mesh/polyline/pointcloud coordinates

required
F (f,s) numpy int array (optional

Matrix of mesh/polyline/pointcloud indices into V. If None, input is assumed to be point cloud.

None)
use_cpp bool, optional (default False)

Whether to use the C++ libigl implementation of the AABB tree (much faster). If True, the following parameters are ignored.

False
use_aabb bool, optional (default False)

Whether to use an AABB tree for logarithmic computation

False
C numpy double array, optional (default None)

Matrix of AABB box centers (if None, will be computed)

None
W numpy double array, optional (default None)

Matrix of AABB box widths (if None, will be computed)

None
CH numpy int array, optional (default None)

Matrix of child indeces (-1 if leaf node). If None, will be computed

None
tri_indices numpy int array, optional (default None)

Vector of AABB element indices (-1 if not leaf node). If None, will be computed

required
split_dir numpy double array, optional (default None)

Vector of AABB split directions (if None, will be computed)

None

Returns:

Name Type Description
squared_distances (p,) numpy double array

Vector of minimum squared distances

indices (p,) numpy int array

Indices into F (or V, if F is None) of closest elements to each query point

lmbs (p,s) numpy double array

Barycentric coordinates into the closest element of each closest mesh point to each query point

See Also

squared_distance_to_element, initialize_aabb.

Examples:

v,f = gpytoolbox.read_mesh("bunny.obj") # Read a mesh
v = gpytoolbox.normalize_points(v) # Normalize mesh
# Generate query points
P = 2*np.random.rand(num_samples,3)-4
# Compute distances
sqrD_gt,ind,b = gpytoolbox.squared_distance(P,v,F=f,use_aabb=True)
Source code in src/gpytoolbox/squared_distance.py
 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
137
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
165
166
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
192
193
194
195
196
197
198
def squared_distance(P,V,F=None,use_cpp=False,use_aabb=False,C=None,W=None,CH=None,tri_ind=None,split_dir=None):
    """Squared distances from a set of points in space.

    General-purpose function which computes the squared distance from a set of points to a mesh, point cloud or polyline, in two or three dimensions. Optionally, uses an aabb tree for efficient computation.

    Parameters
    ----------
    P : (p,dim) numpy double array
        Matrix of query point positions
    V : (v,dim) numpy double array
        Matrix of mesh/polyline/pointcloud coordinates
    F : (f,s) numpy int array (optional, default None)
        Matrix of mesh/polyline/pointcloud indices into V. If None, input is assumed to be point cloud.
    use_cpp : bool, optional (default False)
        Whether to use the C++ libigl implementation of the AABB tree (much faster). If True, the following parameters are ignored.
    use_aabb : bool, optional (default False)
        Whether to use an AABB tree for logarithmic computation 
    C : numpy double array, optional (default None)
        Matrix of AABB box centers (if None, will be computed)
    W : numpy double array, optional (default None)
        Matrix of AABB box widths (if None, will be computed)
    CH : numpy int array, optional (default None)
        Matrix of child indeces (-1 if leaf node). If None, will be computed
    tri_indices : numpy int array, optional (default None)
        Vector of AABB element indices (-1 if *not* leaf node). If None, will be computed
    split_dir : numpy double array, optional (default None)
        Vector of AABB split directions (if None, will be computed)

    Returns
    -------
    squared_distances : (p,) numpy double array
        Vector of minimum squared distances
    indices : (p,) numpy int array
        Indices into F (or V, if F is None) of closest elements to each query point
    lmbs : (p,s) numpy double array
        Barycentric coordinates into the closest element of each closest mesh point to each query point


    See Also
    --------
    squared_distance_to_element, initialize_aabb.

    Examples
    --------
    ```python
    v,f = gpytoolbox.read_mesh("bunny.obj") # Read a mesh
    v = gpytoolbox.normalize_points(v) # Normalize mesh
    # Generate query points
    P = 2*np.random.rand(num_samples,3)-4
    # Compute distances
    sqrD_gt,ind,b = gpytoolbox.squared_distance(P,v,F=f,use_aabb=True)
    ```
    """
    if (F is None):
        F = np.linspace(0,V.shape[0]-1,V.shape[0],dtype=int)[:,None]
    dim = V.shape[1]
    P = np.reshape(P,(-1,dim),order='F')


    if use_cpp:
        try:
            from gpytoolbox_bindings import _point_mesh_squared_distance_cpp_impl
        except:
            raise ImportError("Gpytoolbox cannot import its C++ point_mesh_squared_distance binding.")
        squared_distances, indices, closest_points = _point_mesh_squared_distance_cpp_impl(V.astype(np.float64),F.astype(np.int32),P.astype(np.float64))
        lmbs = np.zeros((P.shape[0],F.shape[1]))
        if (F.shape[1] == 1):
            lmbs = np.ones((P.shape[0],1))
        elif(F.shape[1] == 2):
            lmbs = np.zeros((P.shape[0],2))
            lmbs[:,0] = np.linalg.norm(V[F[indices,0],:]-closest_points,axis=1)
            lmbs[:,1] = np.linalg.norm(V[F[indices,1],:]-closest_points,axis=1)
            lmbs = lmbs/np.sum(lmbs,axis=1)[:,None]
            lmbs = np.fliplr(lmbs)
        elif(F.shape[1] == 3):
            lmbs = barycentric_coordinates(closest_points,V[F[indices,0],:],V[F[indices,1],:],V[F[indices,2],:])
            # for j in range(P.shape[0]):
            #     lmbs[j,:] = barycentric_coordinates(closest_points[j,:],V[F[indices[j],0],:],V[F[indices[j],1],:],V[F[indices[j],2],:])
        return squared_distances, indices, lmbs


    squared_distances = -np.ones(P.shape[0])
    indices = -np.ones(P.shape[0],dtype=int)
    lmbs = np.zeros((P.shape[0],F.shape[1]))
    if use_aabb:
        # Build tree once
        if ((C is None) or (W is None) or (tri_ind is None) or (CH is None) or (split_dir is None)):
            C,W,CH,_,_,tri_ind,split_dir = initialize_aabbtree(V,F=F)
            # V,Q,H = bad_quad_mesh_from_quadtree(C,W,CH)
            # ps.init()
            # ps.register_surface_mesh("test",V,Q)
            # ps.show()
        for j in range(P.shape[0]):
            t = closest_point_traversal(V,F,P[j,:])
            traverse_fun = t.traversal_function
            add_to_queue_fun = t.add_to_queue
            # print(C)
            # print(W)
            # print(CH)
            # print(tri_ind)
            _ = traverse_aabbtree(C,W,CH,tri_ind,split_dir,traverse_fun,add_to_queue=add_to_queue_fun)
            # print(t.num_traversal)
            indices[j] = np.squeeze(t.current_best_element)
            squared_distances[j] = np.squeeze(t.current_best_guess)
            lmbs[j,:] = t.current_best_lmb
    else:
        # Loop over every element
        t = None
        for j in range(P.shape[0]):
            min_sqrd_dist = np.Inf
            ind = -1
            best_lmb = []
            for i in range(F.shape[0]):
                this_sqrd_dist,lmb = squared_distance_to_element(P[j,:],V,F[i,:])
                if this_sqrd_dist<min_sqrd_dist:
                    ind = i
                    min_sqrd_dist = this_sqrd_dist
                    best_lmb = lmb
            squared_distances[j] = min_sqrd_dist
            indices[j] = ind
            lmbs[j,:] = best_lmb
    return squared_distances, indices, lmbs