Skip to content

Instantly share code, notes, and snippets.

@tanzby
Last active September 10, 2020 15:11
Show Gist options
  • Save tanzby/a9d5aab6b4a2c8024983d958c7b2d630 to your computer and use it in GitHub Desktop.
Save tanzby/a9d5aab6b4a2c8024983d958c7b2d630 to your computer and use it in GitHub Desktop.
program for computing area of two overlap triangles
#include <iostream>
#include <Eigen/Eigen>
using namespace std;
class TriangleOverlap
{
public:
using Vector2 = Eigen::Vector2d;
struct Triangle
{
Vector2 s;
Vector2 alpha, beta;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
vector<Vector2> get_vertices() const
{
return { s, s+alpha, s+beta};
}
};
protected:
bool inside(double x1, double y1, double x2, double y2, double xk, double yk )
{
return (x1 == x2 || (min(x1, x2) <= xk && xk <= max(x1, x2)))
&& (y1 == y2 || (min(y1, y2) <= yk && yk <= max(y1, y2)));
}
void update(Vector2& ans, double xk, double yk)
{
if (xk < ans[0] || (xk == ans[0] && yk < ans[1])) {
ans = {xk, yk};
}
}
public:
bool get_line_segment_intersection(Vector2 start1, Vector2 end1, Vector2 start2, Vector2 end2, Vector2& out_intersect)
{
// reference from https://leetcode-cn.com/problems/intersection-lcci/solution/jiao-dian-by-leetcode-solution/
double x1 = start1[0], y1 = start1[1];
double x2 = end1[0], y2 = end1[1];
double x3 = start2[0], y3 = start2[1];
double x4 = end2[0], y4 = end2[1];
if ((y4 - y3) * (x2 - x1) == (y2 - y1) * (x4 - x3))
{
if ((y2 - y1) * (x3 - x1) == (y3 - y1) * (x2 - x1))
{
if (inside(x1, y1, x2, y2, x3, y3)) {
update(out_intersect, x3, y3);
}
if (inside(x1, y1, x2, y2, x4, y4)) {
update(out_intersect, x4, y4);
}
if (inside(x3, y3, x4, y4, x1, y1)) {
update(out_intersect, x1, y1);
}
if (inside(x3, y3, x4, y4, x2, y2)) {
update(out_intersect, x2, y2);
}
return true;
}
}
else
{
double t1 = (x3 * (y4 - y3) + y1 * (x4 - x3) - y3 * (x4 - x3) - x1 * (y4 - y3)) /
((x2 - x1) * (y4 - y3) - (x4 - x3) * (y2 - y1));
double t2 = (x1 * (y2 - y1) + y3 * (x2 - x1) - y1 * (x2 - x1) - x3 * (y2 - y1)) /
((x4 - x3) * (y2 - y1) - (x2 - x1) * (y4 - y3));
if (t1 >= 0.0 && t1 <= 1.0 && t2 >= 0.0 && t2 <= 1.0)
{
out_intersect = {x1 + t1 * (x2 - x1), y1 + t1 * (y2 - y1)};
return true;
}
}
return false;
}
vector<Vector2> get_intersection_points(const Triangle& t1, const Triangle& t2)
{
vector<Vector2> ret;
auto v1 = t1.get_vertices();
auto v2 = t2.get_vertices();
vector<pair<Vector2, Vector2>> t1_line_segments {
{v1[0],v1[1]}, {v1[1],v1[2]}, {v1[2],v1[0]}
};
vector<pair<Vector2, Vector2>> t2_line_segments {
{v2[0],v2[1]}, {v2[1],v2[2]}, {v2[2],v2[0]}
};
Vector2 out_intersect = Eigen::Vector2d::Zero();
for (auto&&[p1, q1]: t1_line_segments)
{
for (auto&&[p2, q2]: t2_line_segments)
{
if (get_line_segment_intersection(p1, q1, p2, q2, out_intersect))
{
ret.emplace_back(out_intersect);
}
}
}
return ret;
}
bool is_point_in_triangle(Vector2 p, Triangle t)
{
Eigen::Matrix2d tri_base;
tri_base.col(0) = t.alpha;
tri_base.col(1) = t.beta;
Vector2 uv = tri_base.inverse() * (p - t.s);
return (uv[0] >= 0 && uv[1] >= 0 && uv[0]+uv[1] <= 1.0);
}
vector<Vector2> get_cover_points(Triangle& t1, Triangle& t2)
{
vector<Vector2> ret;
if (is_point_in_triangle(t1.s, t2)) ret.emplace_back(t1.s);
if (is_point_in_triangle(t1.s + t1.alpha, t2)) ret.emplace_back(t1.s + t1.alpha);
if (is_point_in_triangle(t1.s + t1.beta, t2)) ret.emplace_back(t1.s + t1.beta);
return ret;
}
double get_convex_area(vector<Vector2> points)
{
assert(points.size() >= 3);
Vector2 center = Vector2::Zero();
for (auto& p: points)
{
center += p;
}
center /= points.size();
sort(points.begin(), points.end(), [&](const Vector2& a, const Vector2& b)
{
Vector2 u = a - center;
Vector2 v = b - center;
return atan2(u.y(), u.x()) > atan2(v.y(), v.x());
});
points.emplace_back(points.front());
double area = 0;
for (size_t i = 1; i < points.size(); ++i)
{
Vector2 u = points[i-1] - center;
Vector2 v = points[i ] - center;
area += abs(u.x() * v.y() - u.y() * v.x()) * 0.5;
}
return area;
}
TriangleOverlap() = default;
~TriangleOverlap() = default;
double get_overlap_area(Triangle t1, Triangle t2)
{
vector<Vector2> all_points;
vector<Vector2> intersect_points = get_intersection_points(t1, t2);
all_points.insert(all_points.end(), intersect_points.begin(), intersect_points.end());
vector<Vector2> cover_points = get_cover_points(t1, t2);
all_points.insert(all_points.end(), cover_points.begin(), cover_points.end());
cover_points = get_cover_points(t2, t1);
all_points.insert(all_points.end(), cover_points.begin(), cover_points.end());
double area = get_convex_area(all_points);
return area;
}
};
int main()
{
TriangleOverlap::Triangle t1;
t1.s = {0.0, 0.0};
t1.beta = {3.0, 4.0};
t1.alpha = {5.0, 1.0};
TriangleOverlap::Triangle t2;
t2.s = {2.0, 0.0};
t2.beta = {0.0, 6.0};
t2.alpha = {7.0, 0.0};
cout << "Vertices of Triangle 1: " << endl;
for (auto& v: t1.get_vertices()) { cout << "(" << v.transpose() << ") "; }
cout << endl;
cout << "Vertices of Triangle 2: " << endl;
for (auto& v: t2.get_vertices()) { cout << "(" << v.transpose() << ") "; }
cout << endl << endl;
TriangleOverlap app;
cout << "All intersection points of all line segments:" << endl;
for (auto& p: app.get_intersection_points(t1, t2)) { cout << "(" << p.transpose() << ") "; }
cout << endl << endl;
cout << "All Cover points of two triangles:" << endl;
for (auto& p: app.get_cover_points(t1, t2)) { cout << "(" << p.transpose() << ") "; }
cout << endl;
for (auto& p: app.get_cover_points(t2, t1)) { cout << "(" << p.transpose() << ") "; }
cout << endl;
cout << "intersection area: " << app.get_overlap_area(t1, t2) << endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment