Note -「计算几何」模板
2024-10-19 20:31:41
尚未完整测试，务必留意模板 bug！
/* Clearink */
#include <cmath>
#include <queue>
#include <cstdio>
#include <vector>
#include <algorithm>
namespace PCG {
const double PI = acos ( -1. ), EPS = 1e-9, INF = 2e9;
/* treat x as 0 <=> -EPS < x < EPS */
inline double dabs ( const double x ) { return x < 0 ? -x : x; }
inline double dmin ( const double a, const double b ) { return a < b ? a : b; }
inline double dmax ( const double a, const double b ) { return b < a ? a : b; }
inline int dcmp ( const double x, const double y = 0 ) {
return dabs ( x - y ) < EPS ? 0 : x < y ? -1 : 1;
}
struct Point {
double x, y;
inline Point ( const double tx = 0, const double ty = 0 ):
x ( tx ), y ( ty ) {}
inline Point operator + ( const Point& p ) const {
return Point ( x + p.x, y + p.y );
}
inline Point operator - ( const Point& p ) const {
return Point ( x - p.x, y - p.y );
}
inline Point operator - () const {
return Point ( -x, -y );
}
inline Point operator * ( const double v ) const {
return Point ( x * v, y * v );
}
inline Point operator / ( const double v ) const {
return Point ( x / v, y / v );
}
inline double operator * ( const Point& p ) const { // dot.
return x * p.x + y * p.y;
}
inline double operator ^ ( const Point& p ) const { // cross.
return x * p.y - y * p.x;
}
inline bool operator == ( const Point& p ) const {
return !dcmp ( x, p.x ) && !dcmp ( y, p.y );
}
inline bool operator != ( const Point& p ) const {
return !( *this == p );
}
inline bool operator < ( const Point& p ) const { // as a pair (x,y).
return dcmp ( x, p.x ) ? x < p.x : y < p.y;
}
inline double length () const { return sqrt ( *this * *this ); }
inline Point unit () const { return *this / this->length (); }
inline Point normal () const { return Point ( y, -x ); }
inline double angle () const { // [0,2pi).
double t = atan2 ( y, x );
return t < 0 ? t + 2 * PI : t;
}
inline Point rotate ( const double alpha ) const {
double c = cos ( alpha ), s = sin ( alpha );
return Point ( x * c - y * s, y * c + x * s );
}
friend inline double angle ( const Point& p, const Point& q ) { // [0,pi).
double t = atan2 ( p ^ q, p * q );
t < 0 && ( t += 2 * PI, 0 );
return t < PI ? t : 2 * PI - t;
}
friend inline double dist ( const Point& p, const Point& q ) {
return ( p - q ).length ();
}
friend inline double slope ( const Point& p, const Point& q ) {
return dcmp ( p.x, q.x ) ?
( p.y - q.y ) / ( p.x - q.x ) : p.y < q.y ? INF : -INF;
}
inline void read () { scanf ( "%lf %lf", &x, &y ); }
inline void _show ( const char ch = '\n' ) const {
#ifdef RYBY
printf ( "(%f, %f)%c", x, y, ch );
#endif
}
};
typedef Point Vector;
struct Line {
Point p, v;
inline Line (): p (), v () {}
inline Line ( Point a, Point b, const bool type = false ):
p ( a ), v ( type ? b - a : b ) {}
inline Line ( const double a, const double b, const double c ):
p ( dcmp ( a ) ? Point ( -c / a, 0 ) : Point ( 0, -c / b ) ),
v ( -b, a ) {}
inline Point A () const { return p; }
inline Point B () const { return p + v; }
inline bool operator < ( const Line& l ) const {
return dcmp ( atan2 ( v.y, v.x ), atan2 ( l.v.y, l.v.x ) ) < 0;
}
inline bool onLeft ( const Point& q ) const {
return dcmp ( v ^ ( q - p ) ) > 0;
}
inline bool onLine ( const Point& q ) const {
return !dcmp ( ( q - p ) ^ v );
}
inline bool onRay ( const Point& q ) const {
return onLine ( q ) && dcmp ( ( q - p ) * v ) >= 0;
}
inline bool onSegment ( const Point& q ) const {
if ( !onLine ( q ) ) return false;
return dcmp ( ( A () - q ) * ( B () - q ) ) <= 0;
}
friend inline bool sameSide ( const Line& l,
const Point& p, const Point& q ) {
return dcmp ( ( ( p - l.p ) ^ ( p - l.B () ) )
* ( ( q - l.p ) ^ ( q - l.B () ) ) ) > 0;
}
friend inline bool interSegment ( const Line& l1, const Line& l2 ) {
return ( !sameSide ( l1, l2.p, l2.B () ) )
&& ( !sameSide ( l2, l1.p, l1.B () ) );
}
friend inline bool interRay ( const Line& l1, const Line& l2 ) {
return dcmp ( ( ( l2.p - l1.p ) ^ l2.v ) / ( l1.v ^ l2.v ) ) > 0
&& dcmp ( ( ( l1.p - l2.p ) ^ l1.v ) / ( l2.v ^ l1.v ) ) > 0;
}
friend inline Point lineInter ( const Line& l1, const Line& l2 ) {
return l1.p + l1.v * ( ( l2.p - l1.p ) ^ l2.v ) / ( l1.v ^ l2.v );
}
};
inline std::vector<Point> getConvex ( std::vector<Point> vec,
const bool allowCol ) {
static std::vector<Point> ret;
ret.resize ( vec.size () << 1 );
std::sort ( vec.begin (), vec.end () );
int n = ( int ) vec.size (), top = 0;
for ( int i = 0; i < n; ++i ) {
for ( int d; top > 1; --top ) {
d = dcmp ( ( ret[top - 1] - ret[top - 2] )
^ ( vec[i] - ret[top - 2] ) );
if ( !( ( !allowCol && d <= 0 ) || ( allowCol && d < 0 ) ) ) break;
}
ret[top++] = vec[i];
}
for ( int tmp = top, i = n - 2; ~i; --i ) {
for ( int d; top > tmp; --top ) {
d = dcmp ( ( ret[top - 1] - ret[top - 2] )
^ ( vec[i] - ret[top - 2] ) );
if ( !( ( !allowCol && d <= 0 ) || ( allowCol && d < 0 ) ) ) break;
}
ret[top++] = vec[i];
}
if ( n > 1 ) --top;
return ret.resize ( top ), ret;
}
inline bool poleCmp ( const Point& p, const Point& q ) {
static int t;
return ( t = dcmp ( p ^ q ) ) > 0
|| ( !t && dcmp ( p.length (), q.length () ) < 0 );
}
inline Point polyCentroid ( const std::vector<Point>& poly ) {
double area = 0; Point ret;
int n = ( int ) poly.size ();
for ( int i = 0; i ^ poly.size (); ++i ) {
double s = poly[i] ^ poly[( i + 1 ) % n];
area += s;
ret.x += ( poly[i].x + poly[( i + 1 ) % n].x ) * s;
ret.y += ( poly[i].y + poly[( i + 1 ) % n].y ) * s;
}
ret.x /= 3, ret.y /= 3; // triangle's centroid.
ret.x /= area, ret.y /= area; // average.
return ret;
}
inline double polyArea ( const std::vector<Point>& poly ) {
double ret = 0; int n = ( int ) poly.size ();
for ( int i = 0; i < n; ++i ) {
ret += poly[i] ^ poly[( i + 1 ) % n];
}
return dabs ( ret / 2 );
}
inline std::vector<Point> halfPlaneInter ( std::vector<Line> lvec ) {
static std::vector<std::pair<double, int> > ord;
static std::deque<Line> que; que.clear ();
static std::deque<Point> ret; ret.clear ();
lvec.push_back ( Line ( Point ( -INF, -INF ), Point ( 1, 0 ) ) );
lvec.push_back ( Line ( Point ( -INF, INF ), Point ( 0, -1 ) ) );
lvec.push_back ( Line ( Point ( INF, -INF ), Point ( 0, 1 ) ) );
lvec.push_back ( Line ( Point ( INF, INF ), Point ( -1, 0 ) ) );
int n = ( int ) lvec.size (); ord.resize ( n );
for ( int i = 0; i < n; ++i ) {
ord[i].first = atan2 ( lvec[i].v.y, lvec[i].v.x );
ord[i].second = i;
}
std::sort ( ord.begin (), ord.end () );
que.push_back ( lvec[ord[0].second] );
for ( int i = 1; i < n; ++i ) {
const Line& l ( lvec[ord[i].second] );
for ( ; que.size () > 1 && !l.onLeft ( ret.back () );
que.pop_back (), ret.pop_back () );
for ( ; que.size () > 1 && !l.onLeft ( ret[0] );
que.pop_front (), ret.pop_front () );
if ( dcmp ( l.v ^ que.back ().v ) ) {
que.push_back ( l );
if ( que.size () > 1 ) {
ret.push_back (
lineInter ( que[que.size () - 2], que.back () ) );
}
} else if ( que.back ().onLeft ( l.p ) ) {
que.back () = l;
if ( que.size () > 1 ) {
ret.back () = lineInter ( que[que.size () - 2], que.back () );
}
}
}
for ( ; que.size () > 1 && !que[0].onLeft ( ret.back () );
que.pop_back (), ret.pop_back () );
if ( que.size () <= 2 ) return {};
if ( que.size () > 1 ) ret.push_back ( lineInter ( que[0], que.back () ) );
return std::vector<Point> ( ret.begin (), ret.end () );
}
inline double convexDiameter ( const std::vector<Point>& conv ) {
int n = ( int ) conv.size ();
if ( n == 1 ) return 0;
if ( n == 2 ) return dist ( conv[0], conv[1] );
double ret = 0;
for ( int i = 0, j = 2; i < n; ++i ) {
for ( ; dabs ( ( conv[j] - conv[i] )
^ ( conv[( i + 1 ) % n] - conv[i] ) )
< dabs ( ( conv[( j + 1 ) % n] - conv[i] )
^ ( conv[( i + 1 ) % n] - conv[i] ) ); j = ( j + 1 ) % n );
ret = dmax ( ret, dmax ( dist ( conv[i], conv[j] ),
dist ( conv[( i + 1 ) % n], conv[j] ) ) );
}
return ret;
}
inline int findPole ( const std::vector<Point>& vec ) {
int ret = -1, n = ( int ) vec.size ();
for ( int i = 0; i < n; ++i ) {
if ( !~ret || dcmp ( vec[ret].y, vec[i].y ) > 0
|| ( !dcmp ( vec[ret].y, vec[i].y )
&& dcmp ( vec[ret].x, vec[i].x ) > 0 ) ) {
ret = i;
}
}
return ret;
}
inline void getPoleOrdered ( std::vector<Point>& conv ) {
int pid = findPole ( conv ), n = ( int ) conv.size ();
pid = n - pid - 1; // reversed.
std::reverse ( conv.begin (), conv.end () );
std::reverse ( conv.begin (), conv.begin () + pid + 1 );
std::reverse ( conv.begin () + pid + 1, conv.end () );
}
inline std::vector<Point> convexSum ( std::vector<Point> A,
std::vector<Point> B ) {
static std::vector<Point> ret; ret.clear ();
getPoleOrdered ( A ), getPoleOrdered ( B );
// if use `getConvexG`, there's no need to `getPoleOrdered`.
int n = ( int ) A.size (), m = ( int ) B.size ();
Point ap ( A[0] ), bp ( B[0] );
ret.push_back ( ap + bp );
for ( int i = 0; i < n - 1; ++i ) A[i] = A[i + 1] - A[i];
A[n - 1] = ap - A[n - 1];
for ( int i = 0; i < m - 1; ++i ) B[i] = B[i + 1] - B[i];
B[m - 1] = bp - B[m - 1];
int i = 0, j = 0;
while ( i < n && j < m ) {
ret.push_back ( ret.back ()
+ ( dcmp ( A[i] ^ B[j] ) >= 0 ? A[i++] : B[j++] ) );
}
for ( ; i < n; ret.push_back ( ret.back () + A[i++] ) );
for ( ; j < m; ret.push_back ( ret.back () + B[j++] ) );
return ret;
}
} using namespace PCG;
int n;
std::vector<Point> pos, cnv;
int main () {
/*
example: https://www.luogu.com.cn/problem/P2742
*/
scanf ( "%d", &n ), pos.resize ( n );
for ( int i = 0; i < n; ++i ) pos[i].readn ();
cnv = getConvex ( pos, true ); // also, it could be `getConvex ( pos, false )`.
#ifdef RYBY
for ( auto p: cnv ) p._show ( ' ' );
putchar ( '\n' );
#endif
int s = cnv.size ();
double ans = 0;
for ( int i = 0; i < s; ++i ) {
ans += dist ( cnv[i], cnv[( i + 1 ) % s] );
}
printf ( "%.2f\n", ans );
return 0;
}
/*
try this data:
5
0 0
0 3
1 2
2 1
3 0
*/
大概是壬寅年新版本。（
/*+Rainybunny+*/
// #include <bits/stdc++.h>
#include <cmath>
#include <queue>
#include <cstdio>
#include <vector>
#include <cassert>
#include <iostream>
#include <algorithm>
#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)
namespace ComputingGeometry {
const double EPS = 1e-9, PI = acos(-1.), DINF = 1e18;
template <typename Tp>
inline void chkmin(Tp& u, const Tp& v) { v < u && (u = v, 0); }
template <typename Tp>
inline void chkmax(Tp& u, const Tp& v) { u < v && (u = v, 0); }
template <typename Tp>
inline Tp imin(const Tp& u, const Tp& v) { return u < v ? u : v; }
template <typename Tp>
inline Tp imax(const Tp& u, const Tp& v) { return u < v ? v : u; }
template <typename Tp>
inline Tp iabs(const Tp& u) { return u < 0 ? -u : u; }
inline int sign(const double x) { return iabs(x) <= EPS ? 0 : x < 0 ? -1 : 1; }
struct Point {
double x, y;
Point(): x(0.), y(0.) {}
Point(const double u, const double v): x(u), y(v) {}
inline void read() { scanf("%lf %lf", &x, &y); }
inline Point operator + (const Point& p) const {
return { x + p.x, y + p.y };
}
inline Point operator - () const { return { -x, -y }; }
inline Point operator - (const Point& p) const {
return { x - p.x, y - p.y };
}
inline Point operator * (const double k) const {
return { k * x, k * y };
}
inline Point operator / (const double k) const {
return { x / k, y / k };
}
inline bool operator == (const Point& p) const {
return !sign(x - p.x) && !sign(y - p.y);
}
inline bool operator != (const Point& p) const {
return !(*this == p);
}
inline double operator * (const Point& p) const {
return x * p.x + y * p.y;
}
inline double operator ^ (const Point& p) const {
return x * p.y - y * p.x;
}
inline double leng() const { return sqrt(*this * *this); }
friend inline double dist(const Point& u, const Point& v) {
return (u - v).leng();
}
inline Point norm() const { return { -y, x }; }
inline double angle() const {
double ret = atan2(y, x);
if (ret < 0) ret += 2 * PI;
return ret;
}
friend inline double angle(const Point& u, const Point& v) {
double ret = v.angle() - u.angle();
if (ret > PI) ret -= 2 * PI;
if (ret < -PI) ret += 2 * PI;
return ret;
}
inline Point rotate(const double alp) const {
double ca = cos(alp), sa = sin(alp);
return { x * ca - y * sa, x * sa + y * ca };
}
inline bool operator < (const Point& p) const {
return !sign(x - p.x) ? y < p.y : x < p.x;
}
};
typedef Point Vector;
typedef std::vector<Point> Polygon;
typedef Polygon Convex;
struct Ray {
Point p, v;
Ray() {}
Ray(const Point& a, const Point& b, const bool type = true) {
if (type) p = a, v = b - a;
else p = a, v = b;
}
Ray(const double a, const double b, const double c):
p(sign(a) ? Point(-c / a, 0) : Point(0, -c / b)), v(-b, a) {}
inline Point st() const { return p; }
inline Point ed() const { return p + v; }
inline void readSeg() {
scanf("%lf %lf %lf %lf", &p.x, &p.y, &v.x, &v.y);
v = v - p;
}
friend inline bool isInterSeg(const Ray& a, const Ray& b) { // *
return imin(a.st().x, a.ed().x) <= imax(b.st().x, b.ed().x)
&& imax(a.st().x, a.ed().x) <= imin(b.st().x, b.ed().x)
&& imin(a.st().y, a.ed().y) <= imax(b.st().y, b.ed().y)
&& imax(a.st().y, a.ed().y) <= imin(b.st().y, b.ed().y)
&& sign(a.v ^ (b.st() - a.st())) * sign(a.v ^ (b.ed() - a.st())) <=0
&& sign(b.v ^ (a.st() - b.st())) * sign(b.v ^ (a.ed() - b.st())) <=0;
}
friend inline Point lineInter(const Ray& a, const Ray& b) {
return a.p + a.v * ((b.p - a.p) ^ b.v) / (a.v ^ b.v);
}
friend inline double pSegDist(const Point& p, const Ray& s) {
if (sign(s.v * (p - s.p)) < 0) return dist(p, s.p);
if (sign(-s.v * (p - s.ed())) < 0) return dist(p, s.ed());
return iabs(s.v ^ (p - s.p)) / s.v.leng();
}
friend inline double sSegDist(const Ray& s, const Ray& t) {
return imin(imin(pSegDist(s.st(), t), pSegDist(s.ed(), t)),
imin(pSegDist(t.st(), s), pSegDist(t.ed(), s)));
}
};
typedef Ray Line;
typedef Ray Segment;
typedef std::vector<Line> PlaneCut;
inline Convex getConvex(Polygon P, const bool allw = false) {
int n = int(P.size()), top = 0, tmp; Convex ret(n << 1);
std::sort(P.begin(), P.end());
for (Point& p: P) {
for (int s; top > 1; --top) {
s = sign((ret[top - 1] - ret[top - 2]) ^ (p - ret[top - 2]));
if (s - !allw >= 0) break;
}
ret[top++] = p;
}
std::reverse(P.begin(), P.end()), tmp = top;
for (Point& p: P) {
for (int s; top > tmp; --top) {
s = sign((ret[top - 1] - ret[top - 2]) ^ (p - ret[top - 2]));
if (s - !allw >= 0) break;
}
ret[top++] = p;
}
if (n > 1) --top;
return ret.resize(top), ret;
}
inline double getArea(const Polygon& P) {
double ret = 0.; int n = int(P.size());
rep (i, 0, n - 1) ret += P[i] ^ P[(i + 1) % n];
return iabs(ret) * 0.5;
}
inline std::pair<Point, Point> convexDiameter(const Convex& C) {
int n = int(C.size());
if (n == 1) return { C[0], C[0] };
if (n == 2) return { C[0], C[1] };
double dia = 0.; std::pair<Point, Point> ans;
for (int i = 0, j = 1; i < n; ++i) {
while (((C[(i + 1) % n] - C[i]) ^ (C[j] - C[i]))
< ((C[(i + 1) % n] - C[i]) ^ (C[(j + 1) % n] - C[i])))
j = (j + 1) % n;
double d1 = dist(C[i], C[j]), d2 = dist(C[(i + 1) % n], C[j]);
if (d1 > dia) dia = d1, ans = { C[i], C[j] };
if (d2 > dia) dia = d2, ans = { C[(i + 1) % n], C[j] };
}
return ans;
}
inline double convicesDist(const Convex& A, const Convex& B) {
int n = int(A.size()), m = int(B.size()), p = 0, q = 0;
rep (i, 1, n - 1) if (A[i].y < A[p].y) p = i;
rep (i, 1, m - 1) if (B[i].y > B[q].y) q = i;
double ret = 1e100;
rep (i, 0, n - 1) {
while (sign((A[(p + 1) % n] - A[p]) ^ (B[(q + 1) % m] - B[q])) > 0)
q = (q + 1) % m;
chkmin(ret, sSegDist(Ray(A[p], A[(p + 1) % n]),
Ray(B[q], B[(q + 1) % m])));
p = (p + 1) % n;
}
return ret;
}
inline Convex halfPlaneInter(PlaneCut& vec, const bool apd = false) {
if (apd) {
vec.push_back(Ray(Point(-DINF, -DINF), Point(1, 0), 0));
vec.push_back(Ray(Point(DINF, -DINF), Point(0, 1), 0));
vec.push_back(Ray(Point(DINF, DINF), Point(-1, 0), 0));
vec.push_back(Ray(Point(-DINF, DINF), Point(0, -1), 0));
}
int n = int(vec.size());
std::vector<double> pol(n); std::vector<int> ord(n);
rep (i, 0, n - 1) ord[i] = i, pol[i] = atan2(vec[i].v.y, vec[i].v.x);
std::sort(ord.begin(), ord.end(),
[&](const int u, const int v)->bool {
return sign(pol[u] - pol[v]) ? pol[u] < pol[v]
: (vec[v].v ^ (vec[u].p - vec[v].p)) > 0;
}
);
std::vector<int> tmp; tmp.push_back(ord[0]);
rep (i, 1, n - 1) {
if (sign(pol[ord[i]] - pol[ord[i - 1]])) {
tmp.push_back(ord[i]);
}
}
ord.swap(tmp), n = ord.size();
std::deque<Line> deq;
deq.push_back(vec[ord[0]]), deq.push_back(vec[ord[1]]);
std::deque<Point> pnt;
pnt.push_back(lineInter(deq.front(), deq.back()));
rep (i, 2, n - 1) {
Line& l(vec[ord[i]]);
while (deq.size() > 1 && sign(l.v ^ (pnt.back() - l.p)) <= 0)
deq.pop_back(), pnt.pop_back();
while (deq.size() > 1 && sign(l.v ^ (pnt.front() - l.p)) <= 0)
deq.pop_front(), pnt.pop_front();
pnt.push_back(lineInter(deq.back(), l));
deq.push_back(l);
}
while (deq.size() > 1
&& sign(deq.front().v ^ (pnt.back() - deq.front().p)) < 0)
deq.pop_back(), pnt.pop_back();
while (deq.size() > 1
&& sign(deq.back().v ^ (pnt.front() - deq.back().p)) < 0)
deq.pop_front(), pnt.pop_front();
if (apd) rep (i, 0, 3) vec.pop_back();
if (deq.size() <= 2) return Convex();
pnt.push_back(lineInter(deq.front(), deq.back()));
while (pnt.size() > 1 && pnt.front() == pnt.back()) pnt.pop_back();
std::vector<Point> ret; ret.push_back(pnt[0]);
rep (i, 1, int(pnt.size()) - 1) {
if (pnt[i] != pnt[i - 1]) {
ret.push_back(pnt[i]);
}
}
return ret;
}
} using namespace ComputingGeometry;
int main() {
int n; scanf("%d", &n);
PlaneCut L;
while (n--) {
int m; scanf("%d", &m); Polygon P(m);
for (auto& p: P) p.read();
rep (i, 0, m - 1) L.push_back(Ray(P[i], P[(i + 1) % m]));
}
printf("%.3f\n", getArea(halfPlaneInter(L)));
return 0;
}
最新文章
- 再看 AspriseOCR - OCR应用开发 -20151124
- calc常用
- [CS231n-CNN] Backpropagation(反向传播算法)
- MapReduce：详解Shuffle过程
- UI：UITableView 编辑、cell重用机制
- 【剑指Offer学习】【面试题40：数组中仅仅出现一次的数字】
- canonical 标签介绍
- 去掉eclipse js 错误提示
- 如何制作一个类似Tiny Wings的游戏(2) Cocos2d-x 2.1.4
- 在ubuntu中为navicat创建快捷方式
- Java排序算法之堆排序
- 原创：TSP问题解决方案-----禁忌搜索算法C实现
- Windows 10新功能
- MO_GLOBAL包中一些过程和函数的使用
- UOJ#7. 【NOI2014】购票 点分治 斜率优化 凸包 二分
- python语法_字符编码
- confluence与jira账号对接、查看到期时间及问题总结
- go语言学习-接口
- Springmvc的原理和业务处理
- this 锁与 static 锁