多边形交叉区域计算面积_多边形与圆相交面积计算

POJ 3675 裸的题了。直接上模板就行了。注意的是，求出的是有向面积，有可能是负数。

#include

#include

#include

#include

#include

#include

#include

#include

#include

#include

#include

#include

#define CL(arr, val) memset(arr, val, sizeof(arr))

#define REP(i, n)for((i) = 0; (i) < (n); ++(i))

#define FOR(i, l, h)for((i) = (l); (i) <= (h); ++(i))

#define FORD(i, h, l)for((i) = (h); (i) >= (l); --(i))

#define L(x)(x) << 1

#define R(x)(x) << 1 | 1

#define MID(l, r)(l + r) >> 1

#define Min(x, y)x < y ? x : y

#define Max(x, y)x < y ? y : x

#define E(x)(1 << (x))

#define iabs(x) (x) < 0 ? -(x) : (x)

#define OUT(x)printf("%I64d\n", x)

#define lowbit(x)(x)&(-x)

#define Read()freopen("data.in", "r", stdin)

#define Write()freopen("data.out", "w", stdout);

const double eps = 1e-8;

typedef long long LL;

const int inf = ~0u>>2;

using namespace std;

inline double max (double a, double b) { if (a > b) return a; else return b; }

inline double min (double a, double b) { if (a < b) return a; else return b; }

inline int fi (double a){

if (a > eps) return 1;

else if (a >= -eps) return 0;

else return -1;

}

class Vector{

public:

double x, y;

Vector (void) {}

Vector (double x0, double y0) : x(x0), y(y0) {}

double operator * (const Vector& a) const { return x * a.y - y * a.x; }

double operator % (const Vector& a) const { return x * a.x + y * a.y; }

Vector verti (void) const { return Vector(-y, x); }

double length (void) const { return sqrt(x * x + y * y); }

Vector adjust (double len) {

double ol = len / length();

return Vector(x * ol, y * ol);

}

Vector oppose (void) { return Vector(-x, -y); }

};

class point{

public:

double x, y;

point (void) {}

point (double x0, double y0) : x(x0), y(y0) {}

Vector operator - (const point& a) const { return Vector(x - a.x, y - a.y); }

point operator + (const Vector& a) const { return point(x + a.x, y + a.y); }

};

class segment{

public:

point a, b;

segment (void) {}

segment (point a0, point b0) : a(a0), b(b0) {}

point intersect (const segment& s) const {

Vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b;

double s1 = v1 * v2, s2 = v3 * v4;

double se = s1 + s2;

s1 /= se, s2 /= se;

return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1);

}

point pverti (const point& p) const {

Vector t = (b - a).verti();

segment uv(p, p + t);

return intersect(uv);

}

bool on_segment (const point& p) const {

if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 &&

fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true;

else return false;

}

};

double radius;

point polygon[110];

double kuras_area (point a, point b, double cx, double cy){

point ori(cx, cy);

int sgn = fi((b - ori) * (a - ori));

double da = (a - ori).length(), db = (b - ori).length();

int ra = fi(da - radius), rb = fi(db - radius);

double angle = acos(((b - ori) % (a - ori)) / (da * db));

segment t(a, b); point h, u; Vector seg;

double ans, dlt, mov, tangle;

if (fi(da) == 0 || fi(db) == 0) return 0;

else if (sgn == 0) return 0;

else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn;

else if (ra >= 0 && rb >= 0){

h = t.pverti(ori);

dlt = (h - ori).length();

if (!t.on_segment(h) || fi(dlt - radius) >= 0)

return radius * radius * (angle / 2) * sgn;

else{

ans = radius * radius * (angle / 2);

tangle = acos(dlt / radius);

ans -= radius * radius * tangle;

ans += radius * sin(tangle) * dlt;

return ans * sgn;

}

}

else{

h = t.pverti(ori);

dlt = (h - ori).length();

seg = b - a;

mov = sqrt(radius * radius - dlt * dlt);

seg = seg.adjust(mov);

if (t.on_segment(h + seg)) u = h + seg;

else u = h + seg.oppose();

if (ra == 1) swap(a, b);

ans = fabs((a - ori) * (u - ori)) / 2;

tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length()));

ans += radius * radius * (tangle / 2);

return ans * sgn;

}

}

const double pi = acos(-1.0);

int main (){

int n;

while(scanf("%lf",&radius)!=EOF){

scanf("%d",&n);

for(int i=0;i

scanf("%lf%lf",&polygon[i].x,&polygon[i].y);

}

double ans=0;

for(int i=0;i

ans+=kuras_area(polygon[i],polygon[(i+1)%n],0,0);

}

printf("%.2lf\n",fabs(ans));

}

}

HDU 3982。很明显，求出切割的半平面交之后就可以得出一个多边形。然后求多边形与圆面积相交部分即可。

#include

#include

#include

#include

#include

using namespace std;

const double eps = 1e-6;

typedef long long LL;

const int inf = ~0u>>2;

const int MAXN=4222;

inline double max (double a, double b) { if (a > b) return a; else return b; }

inline double min (double a, double b) { if (a < b) return a; else return b; }

inline int fi (double a){

if (a > eps) return 1;

else if (a >= -eps) return 0;

else return -1;

}

struct Point {

double x,y;

Point(double xx,double yy){x=xx; y=yy;}

Point(){}

}convex[MAXN],cherry;

class Vector{

public:

double x, y;

Vector (void) {}

Vector (double x0, double y0) : x(x0), y(y0) {}

double operator * (const Vector& a) const { return x * a.y - y * a.x; }//向量叉乘

double operator % (const Vector& a) const { return x * a.x + y * a.y; }//向量点乘

Vector verti (void) const { return Vector(-y, x); }

double length (void) const { return sqrt(x * x + y * y); }

Vector adjust (double len) {

double ol = len / length();

return Vector(x * ol, y * ol);

}

Vector oppose (void) { return Vector(-x, -y); }

};

class point{

public:

double x, y;

point (void) {}

point (double x0, double y0) : x(x0), y(y0) {}

Vector operator - (const point& a) const { return Vector(x - a.x, y - a.y); }

point operator + (const Vector& a) const { return point(x + a.x, y + a.y); }

void _copy(Point t){

x=t.x; y=t.y;

}

};

class segment{

public:

point a, b;

segment (void) {}

segment (point a0, point b0) : a(a0), b(b0) {}

point intersect (const segment& s) const {

Vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b;

double s1 = v1 * v2, s2 = v3 * v4;

double se = s1 + s2;

s1 /= se, s2 /= se;

return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1);

}

point pverti (const point& p) const {

Vector t = (b - a).verti();

segment uv(p, p + t);

return intersect(uv);

}

bool on_segment (const point& p) const {

if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 &&

fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true;

else return false;

}

};

double radius;

point polygon[MAXN];//多边形

double kuras_area (point a, point b, double cx, double cy){

point ori(cx, cy);

int sgn = fi((b - ori) * (a - ori));

double da = (a - ori).length(), db = (b - ori).length();

int ra = fi(da - radius), rb = fi(db - radius);

double angle = acos(((b - ori) % (a - ori)) / (da * db));

segment t(a, b); point h, u; Vector seg;

double ans, dlt, mov, tangle;

if (fi(da) == 0 || fi(db) == 0) return 0;

else if (sgn == 0) return 0;

else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn;

else if (ra >= 0 && rb >= 0){

h = t.pverti(ori);

dlt = (h - ori).length();

if (!t.on_segment(h) || fi(dlt - radius) >= 0)

return radius * radius * (angle / 2) * sgn;

else{

ans = radius * radius * (angle / 2);

tangle = acos(dlt / radius);

ans -= radius * radius * tangle;

ans += radius * sin(tangle) * dlt;

return ans * sgn;

}

}

else{

h = t.pverti(ori);

dlt = (h - ori).length();

seg = b - a;

mov = sqrt(radius * radius - dlt * dlt);

seg = seg.adjust(mov);

if (t.on_segment(h + seg)) u = h + seg;

else u = h + seg.oppose();

if (ra == 1) swap(a, b);

ans = fabs((a - ori) * (u - ori)) / 2;

tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length()));

ans += radius * radius * (tangle / 2);

return ans * sgn;

}

}

const double pi = acos(-1.0);

struct Line {

Point a,b;

double ang;

}line[MAXN/2],st[MAXN/2];

int n,cnt;

double operator *(const Point x,const Point y){

return x.x*y.y-x.y*y.x;

}

Point operator -(Point x,const Point y){

x.x-=y.x; x.y-=y.y;

return x;

}

Point operator * (const Line &x,const Line &y){//交点

double a1=(y.b-x.a)*(y.a-x.a),a2=(y.a-x.b)*(y.b-x.b);

Point r;

r.x=(x.a.x*a2+x.b.x*a1)/(a2+a1);

r.y=(x.a.y*a2+x.b.y*a1)/(a2+a1);

return r;

}

bool operator ==(const Point &a,const Point &b){

return fabs(a.x-b.x)

}

bool operator

if(fabs(x.ang-y.ang)

return (y.b-x.a)*(x.b-y.a)>eps;

return x.ang

}

bool JudgeOut(const Line &x,const Point &p){

return (p-x.a)*(x.b-x.a)>eps;

}

bool Parellel(const Line &x,const Line &y){

return fabs((x.b-x.a)*(y.b-y.a))

}

void HplaneI(){

sort(line,line+n);

int top=1,bot=0,tmp=1;

for(int i=1;i

if(line[i].ang-line[i-1].ang>eps) line[tmp++]=line[i];

}

n=tmp;

st[0]=line[0],st[1]=line[1];

for(int i=2;i

if(Parellel(st[top],st[top-1])||Parellel(st[bot],st[bot+1]))return ;

while(bot

while(bot

st[++top]=line[i];

}

while(bot

while(bot

if(top<=bot+1) return ;

st[++top]=st[bot]; cnt=0;

for(int i=bot;i

polygon[cnt++]._copy(st[i]*st[i+1]);

}

}

int main(){

int T,icase=0;

scanf("%d",&T);

while(T--){

scanf("%lf%d",&radius,&n);

for(int i=0;i

scanf("%lf%lf%lf%lf",&line[i].a.x,&line[i].a.y,&line[i].b.x,&line[i].b.y);

}

scanf("%lf%lf",&cherry.x,&cherry.y);

for(int i=0;i

if((line[i].a-cherry)*(line[i].b-cherry)<0){

swap(line[i].a,line[i].b);

}

line[i].ang=atan2(line[i].b.y-line[i].a.y,line[i].b.x-line[i].a.x);

}

line[n].a=Point(-radius,-radius); line[n].b=Point(radius,-radius); line[n++].ang=0;

line[n].a=Point(radius,-radius); line[n].b=Point(radius,radius); line[n++].ang=pi/2;

line[n].a=Point(radius,radius); line[n].b=Point(-radius,radius); line[n++].ang=pi;

line[n].a=Point(-radius,radius); line[n].b=Point(-radius,-radius); line[n++].ang=-pi/2;

cnt=0;

HplaneI();

double ans=0;

for(int i=0;i

//printf("%lf %lf \n",polygon[i].x,polygon[i].y);

ans+=kuras_area(polygon[i],polygon[(i+1)%cnt],0,0);

}

printf("Case %d: ",++icase);

printf("%.5lf",fabs(fabs(ans)-pi*radius*radius)

puts("%");

}

return 0;

}