Sort points by angle from given axis?

Most straightforward, but possibly not the optimal way is to shift the cartesian coordinates to be relative to center point and then convert them to polar coordinates. Then just subtract the angle of the "starting vector" modulo 360, and finally sort by angle.

Or, you could make a custom comparator for just handling all the possible slopes and configurations, but I think the polar coordinates are little more transparent.

Yes, you can do it with a custom comparator based on the cross-product. The only problem is that a naive comparator won't have the transitivity property. So an extra step is needed, to prevent angles either side of the reference from being considered close.

This will be MUCH faster than anything involving trig. There's not even any need to normalize first.

Here's the comparator:

class angle_sort
{
    point m_origin;
    point m_dreference;

    // z-coordinate of cross-product, aka determinant
    static double xp(point a, point b) { return a.x * b.y - a.y * b.x; }
public:
    angle_sort(const point origin, const point reference) : m_origin(origin), m_dreference(reference - origin) {}
    bool operator()(const point a, const point b) const
    {
        const point da = a - m_origin, db = b - m_origin;
        const double detb = xp(m_dreference, db);

        // nothing is less than zero degrees
        if (detb == 0 && db.x * m_dreference.x + db.y * m_dreference.y >= 0) return false;

        const double deta = xp(m_dreference, da);

        // zero degrees is less than anything else
        if (deta == 0 && da.x * m_dreference.x + da.y * m_dreference.y >= 0) return true;

        if (deta * detb >= 0) {
            // both on same side of reference, compare to each other
            return xp(da, db) > 0;
        }

        // vectors "less than" zero degrees are actually large, near 2 pi
        return deta > 0;
    }
};

Demo: http://ideone.com/YjmaN

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

struct Point {
    static double base_angle;
    static void set_base_angle(double angle){
        base_angle = angle;
    }
    double x;
    double y;
    Point(double x, double y):x(x),y(y){}
    double Angle(Point o = Point(0.0, 0.0)){
        double dx = x - o.x;
        double dy = y - o.y;
        double r = sqrt(dx * dx + dy * dy);
        double angle = atan2(dy , dx);
        angle -= base_angle;
        if(angle < 0) angle += M_PI * 2;
        return angle;
    }
};
double Point::base_angle = 0;

ostream& operator<<(ostream& os, Point& p){
    return os << "Point(" << p.x << "," << p.y << ")";
}

bool comp(Point a, Point b){
    return a.Angle() < b.Angle();
}

int main(){
    Point p[] = { Point(-4., -4.), Point(-6., 3.), Point(2., -4.), Point(1., 5.) };
    Point::set_base_angle(p[0].Angle());
    sort(p, p + 4, comp);
    Point::set_base_angle(0.0);
    for(int i = 0;i< 4;++i){
        cout << p[i] << " angle:" << p[i].Angle() << endl;
    }
}

DEMO

Point(-4,-4) angle:3.92699
Point(2,-4) angle:5.17604
Point(1,5) angle:1.3734
Point(-6,3) angle:2.67795