Optimization of a clockwise sort algorithm

Question

So this is an optimization question, which I could not find answered.

I have written some code for computing the convex hull of a given set of random points. For comparative purposes, I have made my own, slow O(n³) algorithm, using some old OpenGL stuff to visualize it. Because of the nature of the slow convex hull algorithm, the points are not sorted at all. So the sort takes place just after the computation of the CH.

My problem is that up till 1000 points, I have a visual result in less than a second. But for 10000 points, it takes over 15-20 minutes (I have not waited beyond that). However if I skip the sorting code and let opengl display the convex hull vertices unsorted, it takes less than a minute. So it's the ClockWise sorting that's eats up all that time. Check the code (some names don't make sense as they are defined elsewhere) :

// This code actually compares every pair iteratively with respect to the center point
// Consider a given vector "convex", which contains all the convex hull points unsorted

.
..
...
....

int avgx,avgy,sumx=0,sumy=0;

for (int i=0;i<convex.size();i++){
    sumx+=convex.at(i).at(0);
    sumy+=convex.at(i).at(1);
}

avgx=sumx/(convex.size()/2.0); //something like an internal point
avgy=sumy/(convex.size()/2.0);

sort(convex.begin(),convex.end()); //x-sort 
int det,tempx,tempy;
for (int i=0;i<convex.size()-1;i++){
    x1=convex.at(i).at(0);
    y1=convex.at(i).at(1);
    x2=convex.at(i+1).at(0);
    y2=convex.at(i+1).at(1);
    det=(x1-avgx)*(y2-avgy)-(x2-avgx)*(y1-avgy); 
    if (det<0){ //on which side of O-X1 lies X2?
        tempx=convex.at(i).at(0); //swapping points 
        tempy=convex.at(i).at(1);
        convex.at(i).at(0)=convex.at(i+1).at(0);
        convex.at(i).at(1)=convex.at(i+1).at(1);
        convex.at(i+1).at(0)=tempx;
        convex.at(i+1).at(1)=tempy;
        i=-1; //check again the new vector from the beginning
    }
    }
return convex;

The display part:

glColor3f(0.8, 0.2, 0.2);
glPointSize(3);
glBegin(GL_LINE_LOOP);
    for (int i=0;i<count;i++){
        glVertex2f(convexHull.at(i).at(0),convexHull.at(i).at(1));
    }
glEnd();

From what I've seen, this approach (by comparing cross products) is the most efficient. However, before this I wrote some filthy code which actually was faster, as it gave me a visual result within 8 minutes. I don't want to keep that because it's too dirty and long, this one is cleaner but very slow (if it actually works at all....). So, am I bound to wait so much for a CW sort of 10000 convex points, or there is something I am missing?

I'd appreciate any thoughts, let me know if I need to include anything else.

Answer 1

Generally speaking, this problem is a little strange. Most 2D convex-hull algorithms (all that I know of) give out the list of points (vertices) in clockwise or counterclockwise order, or they can be modified trivially to do so.

In any case, since there are several good 2D convex-hull determination methods that run in O(N^2) or faster, you can use one of those to "sort" your data into clockwise order. What I'm saying is that you can run a CH algorithm on your data and get the result in the order you want.

Here's a sample code I had lying around that I think would do what you want:

#define TURN_DIR(p1,p2,p3)  (p1.x * p2.y - p1.y * p2.x + \
                             p2.x * p3.y - p2.y * p3.x + \
                             p3.x * p1.y - p3.y * p1.x)
#define LAST(cntnr)         (cntnr).back()
#define BEFORE_LAST(cntnr)  (cntnr)[(cntnr).size() - 2]

void ConvexHull (std::vector<Point> & pts)
{
    std::sort (pts.begin(), pts.end());

    std::vector<Point> lower, upper;
    for (unsigned i = 0; i < pts.size(); ++i)
    {
        while (lower.size() > 1 && TURN_DIR(BEFORE_LAST(lower), LAST(lower), pts[i]) <= 0)
            lower.pop_back ();
        while (upper.size() > 1 && TURN_DIR(BEFORE_LAST(upper), LAST(upper), pts[i]) >= 0)
            upper.pop_back ();

        lower.push_back (pts[i]);
        upper.push_back (pts[i]);
    }

    upper.insert (upper.end(), lower.rbegin() + 1, lower.rend() - 1);
    pts.swap (upper);
}

There are a few points to consider:

1. The above code receives its input and returns its output in the same argument: pts.
2. The Point structure is a simple struct (or class) with two public members x and y, and a less-than operator (an operator <) that compares them very simply based on x first and then y.
3. I believe the running time of the above code is O(N*log(N)), but it certainly is no worse than O(N^2).
4. The returned points will be in clockwise order. If you want counterclockwise, you only need to change that last two lines.
5. This code will not handle the case that all points have the same X coordinates (I think!)
6. Other than that, this is a functional and fast and simple 2D convex-hull implementation.
7. If there are consecutive points in the input that lie on the same line, this implementation removes them. If you don't want that, you can replace the <= 0 and >= 0 tests in the while loops with < 0 and > 0 respectively.

Let me emphasize this: although the above code is a CH implementation, you can use it just to sort your points in a clockwise winding order (if they already form a convex hull.)

Answer 2

You have implemented Bubble Sort, which is O(n²). You can get O(n log(n)) by using the STL sort with an appropriate comparison functor.

Here's a first effort; it uses a non-transitive comparitor, but it seems to work in my test cases and I think it's generally correct:

struct clockwise : public binary_function<vector<int>, vector<int>, bool>
{
  bool operator()(vector<int> A, vector<int> B)
  {
    int det=(A[0]-avgx)*(B[1]-avgy)-(B[0]-avgx)*(A[1]-avgy);
    return(det<0);
  }

  static int avgx, avgy;
};

int clockwise::avgx = 0;
int clockwise::avgy = 0;

...

int clockwise::avgx=sumx/(convex.size()/2.0);
int clockwise::avgy=sumy/(convex.size()/2.0);
sort(convex.begin(),convex.end(), clockwise()); //clockwise-sort

EDIT:

The non-transitive comparitor wasn't a good idea. This is more reliable:

#include <math.h>

...

struct clockwise : public binary_function<vector<int>, vector<int>, bool>
{
  bool operator()(vector<int> A, vector<int> B)
  {
    return(atan2(A[0]-avgx, A[1]-avgy) < atan2(B[0]-avgx, B[1]-avgy));
  }
}