Find local maximum in 3D list

I have a set of data of the kind {{$x_i$, $y_i$, $z_i$}, ...} at randomly chosen points $x_i, y_i$.

The $z_i$ are supposed to be a smooth function of the $x_i$ and $y_i$. Unfortunately they turn out not to be. A few points stick out. These are local extrema - I want to find them and remove them from the data set.

The problem is that the $x_i$ and $y_i$ are chosen randomly, so I cannot choose nearest neighbors for comparison easily. Are there any methods to find the nearest neighbors in a given list, or even directly to find the local maxima in a random grid?

• One possibility would be to construct the Delaunay triangulation ( = find nearest neighbours), then pick those $(x,y)$ points which have a higher $z$ value than all of their neighbours. Check the ComputationalGeometry package for Delaunay triangulation. – Szabolcs Jun 24 '13 at 7:46
• Other way is to create Median filter. en.wikipedia.org/wiki/Median_filter – Kuba Jun 24 '13 at 7:54

The idea is to construct the Delaunay triangulation and get the points that have higher values than all of their neighbours (according to this triangulation).

Here's a worked example:

Generate the sample data:

f[x_, y_] := Sin[x] Sin[y]
points = RandomReal[{-3, 3}, {300, 2}];

zval = f @@@ points;


Find the local maxima:

<< ComputationalGeometry
tri = DelaunayTriangulation[points];
tests = And @@ Thread[zval[[#1]] > zval[[#2]]] & @@@ tri
maxima = Pick[points, tests]

(* ==> {{1.3437, 1.68432}, {2.75095, -2.69678}, {-2.97601, 1.98261}, {-1.57433, -1.57115}} *)


(For your noisy data you might want to use a stricter criterion than just "larger than all neighbours", e.g. larger by a threshold. Also, for simplicity I used the slow ComputationalGeometry package but you'll probably want something faster.)

Visualize them:

ListDensityPlot[ArrayFlatten[{{points, List /@ zval}}],
Epilog -> {Red, PointSize[Large], Point[maxima]}]


Two of them are on the convex hull. You can filter those out:

Complement[maxima, points[[ConvexHull[points]]]]

(* ==> {{-1.57433, -1.57115}, {1.3437, 1.68432}} *)


(Note: it would have been better to filter these based on the indexes of points rather than the coordinates to avoid unnecessary floating point comparisons, but I was lazy :)

• Thank you for the excellent answer. – Neuneck Jun 24 '13 at 8:02

An alternative way (notice usage of Nearest in case of random distributed points):

pts = Table[RandomReal[{0, 2 Pi}], {1000}, {2}];
data = {##, Sin[#1] Cos[#2]} & @@@ pts;
(*and some noise*)
dev = Table[{RandomReal[{0, 2 Pi}], RandomReal[{0, 2 Pi}], RandomReal[{2, 5}]}, {20}];
data = data~Join~dev;
ListDensityPlot[data]


set = Thread[data[[ ;; , ;; 2]] -> data[[ ;; , 3]]];
dane = With[{k = 1, med = Median@Nearest[set, {#1, #2}, 10]},
If[Abs[(#3 - med)] > k, {#1, #2, med}, {##}]] & @@@ data;

ListDensityPlot[dane]


You can change If[Abs[(#3 - med)] > k and k`, it depends what do You want, and what do You know about deviations.