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][1].)

Visualize them:

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

![enter image description here][2]

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 :)


 [1]: http://mathematica.stackexchange.com/questions/13437/how-to-speed-up-the-function-delaunaytriangulation
 [2]: https://i.sstatic.net/1RQ67.png