This function simply takes all points $\{A, B, C, D, \ldots\}$ connected to a point $P$, and constructs the triangles $PAB, PBC, PCD, \ldots$. Since $A,B,C,...$ are in counterclockwise order, I assumed all these would be valid triangles. But take the following case:

$PAB$ will not be a valid triangle, even though $ABC$ are in counterclockwise order.

I've moved Needs["ComputationalGeometry`"] outside of the function definition of triangles. As Szabolcs correctly remarked in his comments, putting Needs["ComputationalGeometry`"] inside the definition will cause shadowing problems because of the creation of the symbol Global`DelaunayTriangulation

Edit 2

Apparently ListDensityPlot uses a Delaunay triangulation as well, and is much faster than TriangularSurfacePlot, so the first part of this answer could be made much more efficient by rewriting it according to

Note that since ListDensityPlot returns a GraphicsComplex and it keeps the order of the points the same, the index lists for the polygons can extracted directly from the plot without having to lookup the indices of the vertices in points.

A comment on putting the Needs inside the function definition: this will only work if you explicitly include the context names in the symbols, i.e. use ComputationalGeometry`DelaunayTriangulation. This is because the context of symbols is determined at parse time. When the function definition is evaluated in a fresh kernel, the package is not loaded yet, so there's no DelaunayTriangulation symbol yet. But as soon as Mathematica sees this symbol in the definition, it creates it---in the Global` context. When the function is called for the first time, ...
–
SzabolcsJan 19 '12 at 19:29

... the package will be loaded, but ComputationalGeometry`DelaunayTriangulation will be shadowed by Global`DelaunayTriangulation.
–
SzabolcsJan 19 '12 at 19:30

The idea is as follows: firstly, for an every edge {a,b} we find additional vertices attached to it (they form triangle with the edge), see Intersection in tri. Then we gather these vertices according to parts of the plane divided by the edge linear continuation, see GatherBy in findPartners. Finally, we find the triangle of the smallest area in each part of the plane, see SortBy in findSmallest.

This is the same what I did (with minor variations) and has the same problem. Not all of the points that are connected to a given point need to be connected together, not even if they are adjacent in the counterclockwise order.
–
SzabolcsJan 19 '12 at 17:15

Just test on points = Tuples[{0, 1}, 2]; tri = {{1, {2, 3, 4}}, {2, {1, 4}}, {3, {4, 1}}, {4, {2, 1, 3}}}; DelaunayTriangulationQ[points, tri]. This is a manually constructed example, but it is a valid triangulation, and I think there's no guarantee something like this won't be returned by DelaunayTriangulation[].
–
SzabolcsJan 19 '12 at 17:15

I should have explained a bit better what was wrong with my original (flawed) implementation. I have now updated the question to avoid confusion.
–
SzabolcsJan 19 '12 at 17:32

I tried the DelaunayTriangulation directly on your points, and I get a different triangulation than your example: same points in each list, but the last list {{4,{2,1,3}} is a cyclic permutation in the numbers of the first list {{1,{3,4,2}}). And this one works with my code, while yours does not. I wonder if it possible to prove whether DT will always return lists whose repeated elements are cyclic permutations only. (If that makes any sense!)
–
JxBJan 19 '12 at 21:27

Yes, I constructed the example by hand to show that your function is not robust. This data is also a valid triangulation, and the points are still in counterclockwise order. I cannot trust that DelaunayTriangulation won't ever return such a result because according to the documentation any cyclic permutation of the points is a possible output (i.e. it mentions no restrictions).
–
SzabolcsJan 19 '12 at 21:32

In Mathematica 10, you can use DelaunayMesh on a set of points. This returns a MeshRegion. You can use MeshCoordinates to return a list of coordinates of the points (should be the same as the initial set of points) and then MeshCells to return the triangles.

Mathematica is a registered trademark of Wolfram Research, Inc. While the mark is used herein with the limited permission of Wolfram Research, Stack Exchange and this site disclaim all affiliation therewith.