Hey guys!
I'm trying to write a program that calculates the area of the convex hull of a set of points in a plane. I can find which points construct the convex hull but calculating the area is a little bit difficult for me. I guess the problem is that I sort the points, and then remove duplicates before calculating the area.
Any help will be appreciated :).

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <algorithm>

using namespace std;

class Point
{
private:
	double x;
	double y;
public:
	void get_coords(double a, double b) 
	{
		x = a;
		y = b;
	}
	double get_x() 
	{
		return x;
	}
	double get_y()
	{
		return y;
	}
};

int get_orientation(Point P1, Point P2, Point P);
bool is_internal(Point P1, Point P2, Point P3, Point P);

int main (int argc, char *argv[])
{
    if(argc < 2) {
            cout << "Not enough arguments in the command line." << endl;
            system("PAUSE");
            return 0;
    }
    
    const char* filename = argv[1];
    ifstream inFile(filename);
    
	vector<Point> P;
	vector<int> convex;
	
    int pts, cvx = 1;
	
    if(!inFile) {
       cout << endl << "Can't open input file." << filename;
       system("PAUSE");
       return 1;
    }
    
    pts = 0;
    while(!inFile.eof()) {
          inFile >> a;
          inFile >> b;
          Temp.get_coords(a, b); 
          P.push_back(Temp); 
          pts++;        
    }
    
    for(int i = 0; i < pts; i++) {
            convex.push_back(i);
    }
	
    for(int i = 0; i < pts; i++) {
        for(int j = 0; j < pts; j++) {
                if(j != i) {
                     for(int k = 0; k < pts; k++) {
                             if(k != j && k != i) {
                                  for(int l = 0; l < pts; l++) {
                                          if(l != i && l != j && l != k) {
                                               if(is_internal(P[i], P[j], P[k], P[l])) {
                                                  convex[l] = -1; 
                                               }
                                          }
                                  }
                             }
                     }
                }
        }
    }
    
    cout << endl;
    
    sort(convex.begin(), convex.end());
    convex.erase(unique(convex.begin(), convex.end()), convex.end());  
    cout << endl << convex.size() << endl << endl;
    area = 0;
    cvx = convex.size();
    for(int j, i = 0; i < cvx; i++) {
            if(convex[i] >= 0) {
                  cout << convex[i] << endl;
                  j = (i + 1) % cvx;
                  area += P[convex[i]].get_x()*P[(convex[j])].get_y();
                  area -= P[convex[i]].get_y()*P[(convex[j])].get_x();
            }
    }
    cout << "Area of the convex hull is " << area/2 << endl;                                                                                
         
    	
	system("PAUSE");
	return 0;	
}

int get_orientation(Point P1, Point P2, Point P) 
{
    double orientation = ((P2.get_x() - P1.get_x()) * (P.get_y() - P1.get_y()) - (P2.get_y() - P1.get_y()) * (P.get_x() - P1.get_x()));
    if (orientation > 0) {
        return 1;
    }
    else if (orientation < 0) {
        return -1;
    }
    else {
        return 0;
    }
}

bool is_internal(Point P1, Point P2, Point P3, Point P) 
{    
    int line1 = get_orientation(P1, P2, P);
    int line2 = get_orientation(P2, P3, P);
    int line3 = get_orientation(P3, P1, P);
    return ((line1 == line2) && (line2 == line3));
}

Everything actually looks OK, except that following your use of sort->unique->erase, there's still potentially a single -1 in your convex array (if there were any points not on the convex hull), and while you skip it for i, you include it for j when i is at the end of your list. Perhaps a better option than sort->unique->erase is simply:

convex.remove(-1)

Also, if it isn't already assumed, your points must be specified in counter-clockwise order, and there can be no loops in the order of the resulting convex-hull vertices. (Loops in the polygon are OK as long as their points are eliminated from the convex hull. Mostly.)

While your solution for eliminating points looks like it should work fine, it's grossly inefficient. Imagine if you had 1,000 points: you'd be calling is_internal() on the order of a trillion times! So first of all, once P[l] is eliminated from the convex hull, it doesn't need to be checked again:

...
for(int l = 0; l < pts; l++) {
    if(l != i && l != j && l != k && convex[l] != -1) {
        if(is_internal(P[i], P[j], P[k], P[l])) {
            convex[l] = -1;
        }
    }
}
...

Also, since you brute-force iterate over all of the points with i, j and k, you end up checking each triangle 6 times! (e.g., 012, 021, 102, 120, 201, 210.) Since each triangle is uniquely identified by ordering its point-indices, you only need to check larger values from the ones you have so far:

for(int i = 0; i < pts; i++) {
    for(int j = i+1; j < pts; j++) {
        for(int k = j+1; k < pts; k++) {
            ...
        }
    }
}

I'm sure there are other improvements, such as skipping over any i, j, or k if you've determined that (by looping over l) its corresponding point isn't on the convex hull. And then you're still hitting many more triangles than would be required if you already knew the convex hull, but I'll leave that for now.

Good luck!

Thank you for your answer!
The program runs much faster now but the problem with the area still persists. I think I know what is causing it - the input points are randomly positioned and when I calculate the area in the last loop they are not in the same order any more.

"In a random order" would definitely be a problem! When you calculate the area, they should be in the same order they were in before (with the internal points removed), but that doesn't help if they weren't intelligently ordered previously.

You can order the points by determining the angle of each point from an arbitrary point. A good arbitrary point (if not just starting with the first one in your list) is somewhere in the middle of your convex hull. Write a function that determines that point. Determine the angle to each point on the convex hull as the arctangent of the slope of the line (function atan2() is perfect for this).

Let me know how it goes.

How to find the middle point? I was thinking of sorting the points by their x and y coordinates, put them in 2 vectors and then pick a point that is near the middle of both of the vectors, but I don't know how to sort a class vector. And how should I use atan2()(if you are talking about the cmath) to find the angle?

Edited 5 Years Ago by nightbreed: n/a

hey nightbreed,

you don't need to re-sort the points, a "half-way" point can be found simply by iterating over all your points and keeping track of the min & max value for x & y (hint: initialize min_x and max_x to the x of one of the points, same for y). or you can compute the "centroid" (i believe i've got the terminology right) by taking the mean of x and y values ... just compute a running x_sum and y_sum, and at the end divide each by the number of points.

i was just referring to the plain C function atan2() in math.h ... 'man atan2' if you have a *nix command-line, or see e.g. http://en.wikipedia.org/wiki/Atan2 (Google is your friend). it's nice because it returns a value all the way around the circle (-PI,PI]. this will sort with a start-point on the negative-x axis. If you like non-negative values, with a start-point on the positive-x axis (like I do), then it's just:

theta = atan2(dy, dx);
if (theta < 0.0)
    theta += 2*M_PI;

Do dy and dx make sense to you here? The signs (+/-) of each are significant if you use atan2(), so think about what the angle "to" another point means.

I added to the class 1 more variable to store the slope angle. For arbitrary point I took the first one (later I will try with a middle one) and I calculate the slope in the following way:

class Point
{
private:
	double x;
	double y;
	double angle;
public:
...
    void set_angle(double ang)
    {
        angle = ang;
    }

};

...
P[i].set_angle(atan2(P[i].get_y() - P[0].get_y(), P[i].get_x() - P[0].get_x()));

Also I store all the angles in another vector so I can compare them with each convex hull vertex when I calculate the area, but with no success :/

Edited 5 Years Ago by nightbreed: n/a

So what does "no success" mean? Do you now sort your convex hull points by angle?

If you're going to use one of the existing points as the start-point, then (if it's also a convex-hull point) its angle will be undefined. Since you're using the output of atan2() directly, your angle values start and end at the -x-axis (as I mentioned previously), so you should start with the point farthest to the west (most negative X). Then angles to all other points will be in the range [-PI/2, PI/2], and if you can set the angle for the starting point to be either -PI or PI, your entire set will be correctly sortable.

Stay in touch!

I think I got it right! For arbitrary point I used the most x-negative point and set its slope angle to Pi. Then I put all convex hull vertices in another vector and sort it by the slope angle of each point in ascending order (the arbitrary point is always the last in the vector) and then I calculate the area. Everything looks good except with this certain set of coordinates:

48 45
0 34
41 54
50 34
41 57
38 45
37 64
50 47
43 94
50 43
4 91
53 2
82 92
48 46
95 18
41 48
49 58
40 56
99 67
53 46
11 3
47 47
40 51
50 38
45 52
49 42
41 50
44 50
41 46
42 42
47 42
40 72
35 68
42 40
43 48
37 39
36 66
48 51
1 6
44 70
37 49
42 47
40 56
41 65
39 36
39 44
45 51
38 37
40 47

For some reason the point with coordinates (1, 6) is flagged as internal point, while it's actually a convex hull vertex.

I check my results here.

Here is the source code:

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <algorithm>

#define PI 3.1415926535897

using namespace std;

class Point
{
private:
	double x;
	double y;
	double angle;
public:
	void get_coords(double a, double b) 
	{
		x = a;
		y = b;
	}
	void set_angle(double ang)
	{
        angle = ang;
    }
	double get_x() 
	{
		return x;
	}
	double get_y()
	{
		return y;
	}
	double get_angle() const 
	{
           return angle;
    }
};

int get_orientation(Point P1, Point P2, Point P);
bool is_internal(Point P1, Point P2, Point P3, Point P);
int arbitrary_point(vector<Point> P);
bool sort_by_angle(Point &P1, Point &P2);

int main (int argc, char *argv[])
{
    if(argc < 2) {
            cout << "Not enough arguments in the command line." << endl;
            system("PAUSE");
            return 0;
    }
    
    const char* filename = argv[1];
    ifstream inFile(filename);
    
    Point Temp;
	vector<Point> P, c_hull;
	vector<int> external;
	
    int pts, cvx = 1;
    int a, b;
    int count = 0;
    int a_point;
    double area = 0;
	
    if(!inFile) {
       cout << endl << "Can't open input file - " << filename;
       system("PAUSE");
       return 1;
    }
    
    pts = 0;
    while(!inFile.eof()) {
          inFile >> a;
          inFile >> b;
          Temp.get_coords(a, b); 
          P.push_back(Temp); 
          pts++;        
    }
    
    for(int i = 0; i < pts; i++) {
            external.push_back(i);
    }
    
    a_point = arbitrary_point(P);
    P[a_point].set_angle(PI);
    cout << "Arbitrary point -> " << a_point << endl;
    
    for(int i = 0; i < P.size(); i++) {
            if(P[i].get_angle() != P[a_point].get_angle()) {
                    P[i].set_angle(atan2(P[i].get_y() - P[a_point].get_y(), P[i].get_x() - P[a_point].get_x())); //atan2(y2 - y1, x2 - x1) * 180 / PI;
            }
    }
    
    cout << "There are " << pts << " points." << endl;
	cout << "Internal points: " << endl;
    for(int i = 0; i < pts; i++) {
        for(int j = i+1; j < pts; j++) {
                if(j != i) {
                     for(int k = j+1; k < pts; k++) {
                             if(k != j && k != i) {
                                  for(int l = 0; l < pts; l++) {
                                          if(l != i && l != j && l != k && external[l] != -1) { //&& external[l] != -1
                                               if(is_internal(P[i], P[j], P[k], P[l])) {
												  cout << external[l] << ": " << P[external[l]].get_x() << " " << P[external[l]].get_y() << " " << P[external[l]].get_angle() << endl;
												  count++;
                                                  external[l] = -1;
                                               }
                                          }
                                  }
                             }
                     }
                }
        }
	}
    
    external.erase(remove(external.begin(), external.end(), -1), external.end());

    cvx = external.size();

    area = 0;

	for(int i = 0; i < cvx; i++ ) {
		if(external[i] != -1) {
			c_hull.push_back(P[external[i]]);
		}
	}

	cout << endl << endl;

	sort(c_hull.begin(), c_hull.end(), sort_by_angle);

	cout << "Coordinates  / slope angle" << endl;
	for(int j, i = 0; i < c_hull.size(); i++) {
		cout << "(" << c_hull[i].get_x() << ", " << c_hull[i].get_y() << ")  " << c_hull[i].get_angle() << endl;
        j = (i + 1) % cvx;
        area += c_hull[i].get_x()*c_hull[j].get_y();
        area -= c_hull[i].get_y()*c_hull[j].get_x();
	}

	cout << "Area of the convex hull is " << abs(area) << endl;
    	
	system("PAUSE");
	return 0;	
}

int get_orientation(Point P1, Point P2, Point P) 
{
    double orientation = ((P2.get_x() - P1.get_x()) * (P.get_y() - P1.get_y()) - (P2.get_y() - P1.get_y()) * (P.get_x() - P1.get_x()));
    if (orientation > 0) {
        return 1;
    }
    else if (orientation < 0) {
        return -1;
    }
    else {
        return 0;
    }
}

bool is_internal(Point P1, Point P2, Point P3, Point P) 
{    
    int line1 = get_orientation(P1, P2, P);
    int line2 = get_orientation(P2, P3, P);
    int line3 = get_orientation(P3, P1, P);
    return ((line1 == line2) && (line2 == line3));
}

int arbitrary_point(vector<Point> P) 
{
    double min_x;
    int min_index;
	min_index = 0;
    min_x = P[0].get_x();
    for(int i = 1; i < P.size(); i++) {
            if(P[i].get_x() < min_x) {
                            min_x = P[i].get_x();
                            min_index = i;
            }
    }
    return min_index;
}

bool sort_by_angle(Point &P1, Point &P2) 
{
     return (P1.get_angle() < P2.get_angle());
}

Hmmm, not sure right off. To see what's going on with the mystery point, add some debug printing. When it's incorrectly removed from the convex hull, also print out the three points supposedly around it so you can better determine what's going on.

The triangle, which "surrounds" the point, appears to be somewhere in the middle of the set. I tried a lot of combinations and I couldn't replicate the same output. I will try again tomorrow.
Anyways thanks for your time!

This article has been dead for over six months. Start a new discussion instead.