Create non-intersecting polygon passing through all given points

Suppose I have an array of points in random order, and I need to find a polygon (by sorting them, such that every adjacent pair represents a side) which passes through all of the points, and its sides are non-intersecting of course.

I tried to do it by selecting a point, and adding all points to the final array which are below it, sorted left to right. Then, adding all points which are above it, sorted right to left.

I've been told that I can add an additional point and sort naturally to avoid self-intersections.. I am unable to figure out that though. What's a simple way to do this?

As someone said, the minimal length solution is exactly the traveling salesman problem. Here's a non-optimal but feasible approach:

Compute a Delauney triangulation of your points. Successively remove boundary segments until you are left with a boundary that interpolates all points or no more segments can be removed. Don't remove boundary segments if all points of the triangle using that segment are on the boundary. Take this boundary as your path.

I implemented this in Mathematica using 40 random points. Here is a typical result:

The obvious objection is that you might get to a point where not all your points are boundary points, but you can't remove a boundary segment without making the boundary self intersecting. This is a valid objection. It took me dozens of runs to see a case where this happened, but finally got this case:

You can probably see some obvious ways of fixing this using the local topology, but I'll leave the details to you! One thing that might help is "edge flipping" where you take two triangles that share a side, say triangle (p,q,r) and (q,p,s) and replace them with (r,p,s) and (r,s,q) (all coordinates counterclockwise around the triangle). This can be done as long as the resulting triangles in this transformation are also counterclockwise ordered.

To reduce the need for fix-ups, you will want to make good choices of the segments to remove at each step. I used the ratio of the length of the boundary segment to the sum of the lengths of the other side of the candidate triangle (the triangle formed by the potentially incoming point with the segment).