Spatial Searching

Introduction

MATLAB provides you with the necessary functions to perform a spatial search using a Delaunay triangulation. The Delaunay-based search queries that MATLAB supports are:

• Nearest-neighbor search (sometimes called closest-point search or proximity search),

• Point-location search (sometimes called point-in-triangle search or point-in-simplex search, where a simplex is a triangle, tetrahedron or higher dimensional equivalent).

Given a set of points X and a query point q in Euclidean space, the nearest-neighbor search locates a point p in X that is closer to q than to any other point in X. Given a Delaunay triangulation of X, the point location search locates the simplex that contains the query point q.

While MATLAB supports these search schemes in N dimensions, exact spatial searches usually become prohibitive in 6 dimensions. You should consider approximate alternatives for large problems in up to 10 dimensions.

Back to Top

Nearest-Neighbor Search

There are two Delaunay-based approaches to computing nearest neighbors in MATLAB, depending on the dimensionality of the problem:

• For 2-D and 3-D search, use the nearestNeighbor method provided by the DelaunayTri class.

• For 4-D and higher, use the delaunayn function to construct the triangulation and the complementary dsearchn function to perform the search. While these N-D functions support 2-D and 3-D, they are not as efficient as the search methods provided by DelaunayTri.

The following example demonstrates how to perform a nearest-neighbor search in 2-D with DelaunayTri. Begin by creating a random set of 15 points and plotting the points showing ID labels:

X = [3.5 8.2; 6.8 8.3; 1.3 6.5; 3.5 6.3; 5.8 6.2; 8.3 6.5;...
    1 4; 2.7 4.3; 5 4.5; 7 3.5; 8.7 4.2; 1.5 2.1; 4.1 1.1; ...
    7 1.5; 8.5 2.75];

Plot the points and add annotations:

plot(X(:,1),X(:,2),'ob')
hold on
vxlabels = arrayfun(@(n) {sprintf('X%d', n)}, (1:15)');
Hpl = text(X(:,1)+0.2, X(:,2)+0.2, vxlabels, 'FontWeight', ...
  'bold', 'HorizontalAlignment','center', 'BackgroundColor', ...
  'none');
hold off

Create a Delaunay triangulation from the points:

dt = DelaunayTri(X);

Create some query points and for each query point find the index of its corresponding nearest neighbor in X using the nearestNeighbor method:

numq = 10;
q = 2+rand(numq,2)*6;
xi = nearestNeighbor(dt, q);

Add the query points to the plot and add line segments joining the query points to their nearest neighbors:

xnn = X(xi,:);

hold on
plot(q(:,1),q(:,2),'or');
plot([xnn(:,1) q(:,1)]',[xnn(:,2) q(:,2)]','-r');

vxlabels = arrayfun(@(n) {sprintf('q%d', n)}, (1:numq)');
Hpl = text(q(:,1)+0.2, q(:,2)+0.2, vxlabels, 'FontWeight', ...
     'bold', 'HorizontalAlignment','center', ...
     'BackgroundColor','none');

hold off

Performing a nearest-neighbor search in 3-D is a direct extension of the 2-D example based on DelaunayTri. For 4-D and higher use the delaunayn and dsearchn functions as illustrated in the following example:

Create a random sample of points in 4-D and triangulate the points using delaunayn:

X = 20*rand(50,4) -10;
tri = delaunayn(X);

Create some query points and for each query point find the index of its corresponding nearest neighbor in X using the dsearchn function:

q = rand(5,4);
xi = dsearchn(X,tri, q)

The nearestNeighbor method and the dsearchn function allow the Euclidean distance between the query point and its nearest neighbor to be returned as an optional argument. In the 4-D example, you can compute the distances (dnn) as follows:

[xi dnn] = dsearchn(X,tri, q)

Back to Top

Point Location

A point location search is a triangulation-search algorithm that locates the simplex (triangle, tetrahedron, and so on) enclosing a query point. As in the case of the nearest-neighbor search, there are two approaches to performing a point-location search in MATLAB, depending on the dimensionality of the problem:

• For 2-D and 3-D, use the class-based approach with the pointLocation method provided by the DelaunayTri class.

• For 4-D and higher, use the delaunayn function to construct the triangulation and the complementary tsearchn function performs the point location search. Although supporting 2-D and 3-D, these N-D functions are not as efficient as the search methods provided by DelaunayTri.

The following example demonstrates how to use the DelaunayTri class to perform a point location search in 2-D.

Begin with a set of 2-D points. Create the triangulation and plot it showing the triangle ID labels at the incenters of the triangles:

X = [3.5 8.2; 6.8 8.3; 1.3 6.5; 3.5 6.3; 5.8 6.2; ...
     8.3 6.5; 1 4; 2.7 4.3; 5 4.5; 7 3.5; 8.7 4.2; ...
     1.5 2.1; 4.1 1.1; 7 1.5; 8.5 2.75];
dt = DelaunayTri(X);
triplot(dt);

hold on
ic = incenters(dt);
numtri = size(dt,1);
trilabels = arrayfun(@(x) {sprintf('T%d', x)}, (1:numtri)');
Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', ...
      'bold', 'HorizontalAlignment', 'center', 'Color', ...
      'blue');
hold off

Now create some query points and add them to the plot. Then find the index of the corresponding enclosing triangles using the pointLocation method:

numq = 10;
q = 2+rand(numq,2)*6;
hold on; 
plot(q(:,1),q(:,2),'*r'); 
vxlabels = arrayfun(@(n) {sprintf('q%d', n)}, (1:numq)');
Hpl = text(q(:,1)+0.2, q(:,2)+0.2, vxlabels, 'FontWeight', ...
      'bold', 'HorizontalAlignment','center', ... 
      'BackgroundColor', 'none');
hold off

ti = pointLocation(dt, q)

Performing a point-location search in 3-D is a direct extension of performing a point-location search in 2-D with DelaunayTri. For 4-D and higher, use the delaunayn and tsearchn functions as illustrated in the following example:

Create a random sample of points in 4-D and triangulate them using delaunayn:

X = 20*rand(50,4) -10;
tri = delaunayn(X);

Create some query points and find the index of the corresponding enclosing simplices using the tsearchn function:

q = rand(5,4);
ti = tsearchn(X,tri, q)

The pointLocation method and the tsearchn function allow the corresponding barycentric coordinates to be returned as an optional argument. In the 4-D example, you can compute the barycentric coordinates as follows:

[ti bc] = tsearchn(X,tri, q)

The barycentric coordinates are useful for performing linear interpolation. These coordinates provide you with weights that you can use to scale the values at each vertex of the enclosing simplex. See Interpolating Scattered Data for further details.

Back to Top

Searching Non-Delaunay Triangulations

You only can use the search functions tsearchn and dsearchn to search convex, fully connected Delaunay triangulations that were created in MATLAB. If you create a Delaunay triangulation using delaunayn and subsequently modify the location of the coordinates, then the Delaunay criterion could be invalid. Therefore, you should not use the search functions as they may not converge. The search functions do not check the validity of the Delaunay triangulation as this would introduce a significant performance penalty. Imported triangulations, such as triangulations used for finite element analysis, tend to be nonconvex. They could have holes, or they might not respect the Delaunay criterion. Hence, you cannot reliably use the search functions tsearchn and dsearchn on triangulations like these.

The class DelaunayTri is more secure in this respect because you only can use the search methods to search a DelaunayTri. However, you cannot perform a nearestNeighbor search on a DelaunayTri that has constrained edges as it is locally non-Delaunay.

To perform a nearest-neighbor search on a constrained DelaunayTri, do one of the following:

• Delete the constraints

• Create a new unconstrained DelaunayTri using the same point set, and use the unconstrained triangulation to perform the search.

Similarly, to perform a nearest-neighbor search on a triangulation that is represented in face-vertex format, create a DelaunayTri from the vertex coordinates and use it to perform the search.

MATLAB does not provide dedicated functions for performing a point location search on a 2-D (triangle) or 3-D (tetrahedral) triangulations that are in face vertex format and do not meet the search function requirements. However, for relatively small triangulations a brute-force search might be viable. For a 2-D triangulation, you might can recreate the triangulation topology using a constrained DelaunayTri

You can perform the brute-force search by checking the query point against every triangle (or tetrahedron) and selecting the one that contains the point. For example, suppose you have the following triangulation in face-vertex format and you want to perform a brute-force point location search that finds a triangle enclosing a point.

X = [3.5 8.2; 6.8 8.3; 1.3 6.5; 3.5 6.3; 5.8 6.2; 8.3 6.5; ...
     1 4; 2.7 4.3; 5 4.5; 7 3.5; 8.7 4.2; 1.5 2.1; 4.1 1.1; ...
     7 1.5; 8.5 2.75];

tri = delaunay(X);
% Remove first triangle to make the triangulation nonconvex
tri(1,:) = []; 

The first thing you can do is create a TriRep to represent this triangulation:

tr = TriRep(tri,X)

The TriRep has many useful query methods, but it does not have a pointLocation or a nearestNeighbor method. A TriRep can represent a nonconvex triangulation with holes and non-Delaunay triangulations, but the search algorithms preclude these.

TriRep provides a useful method that you can use to check if a query point lies within a specified triangle (or tetrahedron). The method is cartToBary. Given a point expressed in Cartesian coordinates and a specified triangle in the triangulation, this method returns the barycentric coordinates of the point with respect to the triangle. If all three barycentric coordinates are positive, the point is within the triangle . Otherwise, the point is outside the triangle. Plot the triangulation to explore this method:

triplot(tr)
axis equal

hold on
ic = incenters(tr);
numtri = size(tr,1);
trilabels = arrayfun(@(x) {sprintf('T%d', x)}, (1:numtri)');
Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', ...
      'bold', 'HorizontalAlignment', 'center', 'Color', ...
      'blue');
hold off

Search for the triangle enclosing (6,3):

numtri = size(tr,1);
B = cartToBary(tr, (1:numtri)', repmat([6 3], numtri,1));
posB = (B >= 0);
enclosingTri = find(posB(:,1) & posB(:,2) & posB(:,3))

You can now use a constrained DelaunayTri to perform the search. You should start with the TriRep you created in the brute-force search and use it to find all edges in the triangulation so you can constrain them.:

Cedges = edges(tr);

Now create the constrained DelaunayTri and plot it in a separate figure so you can compare both triangulations.:

dt = DelaunayTri(X, Cedges)
figure
triplot(dt)
axis equal
hold on
ic = incenters(dt);
numtri = size(dt,1);
trilabels = arrayfun(@(x) {sprintf('T%d', x)}, (1:numtri)');
Htl = text(ic(:,1), ic(:,2), trilabels, 'FontWeight', ...
      'bold', 'HorizontalAlignment', 'center', 'Color', ...
      'blue');
hold off

Since the triangle numbering is different in each triangle, you need to build an array of indices to map triangle IDs from the DelaunayTri to the TriRep. To do this mapping you will compute the incenters of the TriRep and find the corresponding enclosing triangles in the DelaunayTri:

dtTotrIndex = nan(size(dt,1),1);
ic = incenters(tr);

dtid = pointLocation(dt,ic)
dtTotrIndex(dtid) = 1:size(tr,1);

Now you use the DelaunayTri to perform the point location search. Next, use the map to find the corresponding triangle in the TriRep:

trid = pointLocation(dt, 6,3);
if isfinite(trid)
   trid = dtTotrIndex(trid);
end
trid

Back to Top

Delaunay Triangulation

Voronoi Diagrams

Recommended Products

• MATLAB
• Signal Processing Toolbox
• Optimization Toolbox
• Statistics Toolbox
• Control System Toolbox

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit