Saturday, August 20, 2011

Maximum Triangle Area (Part 3 and final)

Part 1
Part 2

After second part, our algorithm became something like this -
for each 0 <= i < n, initialize A = Pi, B = P(i + 1) % n, C = P(i + 2) % n, then move C counter clockwise until area non-decreasing, move B to its adjacent point counter clockwise, if area non-decreasing, move C again until area non decreasing and so on. You will get the maximum area for each A and output the maximum of these areas. For a given point starting point Pi, if the maximum area is ABC. We can say ABC is 'stable' triangle for A, in other word s(pi, Pi + 1, Pi + 2) = ABC, that means we cannot move B or C any further by not decreasing the area.

Now, after getting a maximum ABC for a fixed first point Pi, we will start with Pi + 1, Pi + 2 and Pi + 3. We get A'B'C' as a stable triangle i.e. s(pi, pi + 1, pi + 2) = A'B'C'. Hence our algorithm is O(n ^ 2). But it can be proved that, rather than, start from Pi + 1 as a second point and Pi + 2 as a third point, we can start with B and C as a second and third point and get the same result, i.e. s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C), where s(Pi, Pi + 1, Pi + 2) = ABC.

Lets assume that we are wrong, and s(Pi + 1, Pi + 2, Pi + 3) = A'B'C' > s(pi, Pi + 1, Pi + 2).

Case 1:
B' < B and C' = C


We know that,
ABC >= AB'C (otherwise largest triangle for A would be AB'C)
=> ABC + BB'A >= AB'C + BB'A
=> AB'BC >= AB'C + BB'A
=> AB'C + BB'C >= AB'C + BB'A
=> BB'C >= BB'A
=> BB'C >= BB'A' (for base BB' area is non decreasing from BB'C to BB'A, so it
will be non decreasing from BB'A to BB'C)
=> BB'C + A'B'C >= BB'A' + A'B'C
=> A'B'BC >= BB'A' + A'B'C
=> A'BC + BB'A' >= BB'A' + A'B'C
=> A'BC >= A'B'C
=> A'BC >= A'B'C' ...... [1] (as C = C')

So, for case 1, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

Case 2:
Let's say, A'B'C' is the largest triangle for first point A'.
Where, B' < B and C' > C. See the following image.



For first point A, ABC is the largest trianlge,
So, ABC >= AB'C
=> ABC + BB'C >= ABC + BB'C
=> ABC + BB'C >= AB'BC
=> ABC + BB'C >= ABC + BB'A
=> BB'C >= BB'A
i.e. BB'A <= BB'C Hence, BB'A <= BB'C' (otherwise, BB'A > BB'C and BB'A <= BB'C, which is not possible for a convex polygon) So, BB'A' <= BB'C' ... [2] (as the triangle area non increasing by moving C' to A, it will be non increasing from A to A') Now as we are considering, for first point A', the largest area triangle is A'B'C'. So, A'B'C' > A'BC'
=> A'B'C' + BB'A' > A'BC' + BB'A
=> A'B'C' + BB'A' > A'B'BC'
=> A'B'C' + BB'A' > A'B'C' + BB'C'
=> BB'A' > BB'C', this contradicts to [2].

So it is not possible to have, A'B'C' > A'BC' i.e. for case 2, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

Case 3:
Let's say, A'B'C' is the largest triangle for fixed first point A', where,
B' < B and C' < C

As ABC is the largest triangle having first point A,
either ABC' >= AB'C' (by moving B' to B, area increases) .. (3A)
or AB'C >= AB'C' (by moving C' to C, area increases) .. (3B)

Case 3A
ABC' >= AB'C'
=> ABC' + BB'A >= AB'C' + BB'A
=> AB'BC' >= AB'C' + BB'A
=> AB'C' + BB'C' >= AB'C' + BB'A
=> BB'C' >= BB'A
=> BB'C' >= BB'A' (As by moving C' to A area non increasing, moving A to A' will area be non increasing)
=> BB'C + A'B'C' >= BB'A' + A'B'C'
=> BB'A'C' >= BB'A' + A'B'C'
=> BB'A' + A'BC' >= BB'A' + A'B'C'
=> A'BC >= A'B'C'

So it is not possible to have, A'B'C' > A'BC' i.e. for case 3A, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

Case 3B
AB'C >= AB'C'
=> AB'C + CC'B' >= AB'C' + CC'B'
=> AB'C'C >= AB'C' + CC'B'
=> AB'C' + CC'A >= AB'C' + CC'B'
=> CC'A >= CC'B'
=> CC'A' >= CC'B' (As moving A to B', area decreases, it's not possible moving A' to B' area increases)
=> CC'A' + A'B'C' >= CC'B' + A'B'C'
=> A'B'C'C >= CC'B' + A'B'C'
=> CC'B' + A'B'C >= CC'B' + A'B'C'
=> A'B'C >= A'B'C'

So it is not possible to have, A'B'C' > A'B'C i.e. for case 3B, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

Case 4:
C' < C and B' = B

Here,
ABC >= ABC'
=> ABC + CC'A >= ABC' + CC'A
=> ABC + CC'A >= ABC'C
=> ABC + CC'A >= ABC + CC'B
=> CC'A >= CC'B
=> CC'A' >= CC'B
=> CC'A' + A'BC' >= CC'B + A'BC'
=> A'BC'C >= CC'B + A'BC'
=> A'BC + CC'B >= CC'B + A'BC'
=> A'BC >= A'BC'
=> A'BC >= A'B'C' (as B' = B)

So it is not possible to have, A'B'C' > A'BC i.e. for case 3B, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

Case 5:
B' > B and C' < C


Now,
A'BC >= A'BC' (From case 4)
=> A'BC + CC'B >= A'BC' + CC'B
=> A'BC'C >= A'BC' + CC'B
=> A'BC' + CC'A' >= A'BC' + CC'B
=> CC'A' >= CC'B
=> CC'A' >= CC'B'
=> CC'A' + A'B'C' >= CC'B' + A'B'C'
=> CC'B'A' >= CC'B' + A'B'C'
=> CC'B' + A'B'C >= CC'B' + A'B'C'
=> A'B'C >= A'B'C'

So it is not possible to have, A'B'C' > A'B'C i.e. for case 5, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C).

So, from all of the 5 possible cases, we found that, s(Pi + 1, Pi + 2, Pi + 3) = s(Pi + 1, B, C). So our algorithm will still work, if we start our iteration for first point i = i + 1, without changing the position 2nd and 3rd point. We will move first point from 0 to n - 1, and 2nd and 3rd point will also not move more than n times. So our algorithm in worst case is O (3 * n) = O(n)

My simple implementation in C++.


#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<vector>
#include<string>
#include<algorithm>
#include<set>
#include<map>
#include<cassert>
#include<sstream>
#include<queue>
#include<stack>
#include<iostream>
using namespace std;

struct Point {
int x;
int y;
Point(int _x, int _y) : x(_x), y(_y) {}
Point() : x(0), y(0) {}
bool operator <(const Point& p) const {
return x < p.x || (x == p.x && y < p.y);
}
};

bool left_turn(const Point& p1, const Point& p2, const Point& p3) {
return (p2.y - p1.y) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.y - p1.y) > 0;
}

// Returns list of points of convex hull in counter clockwise
// The last point and first point will be equal
vector<Point> convex_hull(vector<Point> p) {
vector<Point> ret;
sort(p.begin(), p.end());
// build lower hull
for (int i = 0; i < p.size(); ++i) {
while (ret.size() >= 2 &&
left_turn(ret[ret.size() - 2], ret[ret.size() - 1], p[i])) {
ret.pop_back();
}
ret.push_back(p[i]);
}
int lower_hull_size = ret.size();
// build upper hull
for (int i = p.size() - 2; i >= 0; --i) {
while (ret.size() - lower_hull_size >= 1 &&
left_turn(ret[ret.size() - 2], ret[ret.size() - 1], p[i])) {
ret.pop_back();
}
ret.push_back(p[i]);
}
return ret;
}

double triangle_area (const Point& p1, const Point& p2, const Point& p3) {
return abs((double) p1.x * p2.y +
(double) p2.x * p3.y +
(double) p3.x * p1.y -
(double) p1.y * p2.x -
(double) p2.y * p3.x -
(double) p3.y * p1.x) / 2.0;
}

int main (void)
{
int n;
while (scanf("%d", &n) && n != -1) {
vector<Point> p;
for (int i = 0; i < n; ++i) {
int x, y;
scanf("%d%d", &x, &y);
p.push_back(Point(x, y));
}
p = convex_hull(p); p.pop_back();
int i = 0;
int j = i + 1;
int k = j + 1;
double res = 0.;
while (true) {
double area = triangle_area(p[i], p[j], p[k]);
while (true) {
while (true) {
int nk = (k + 1) % n;
double narea = triangle_area(p[i], p[j], p[nk]);
if (narea >= area) {
k = nk;
area = narea;
} else {
break;
}
}
int nj = (j + 1) % n;
double narea = triangle_area(p[i], p[nj], p[k]);
if (narea >= area) {
j = nj;
area = narea;
} else {
break;
}
}
if (area > res) res = area;
i++;
if (i == j) j = (j + 1) % n;
if (j == k) k = (k + 1) % n;
if (i == n) break;
}
printf("%.2lf\n", res);
}
return 0;
}