博主这里曾经学过计算几何(下文简称jj),所以没有证明或者说明某些算法,不适合初学者食用
辛普森积分Simpson~
用一道例题及黄学长的代码来理解:
黄学长代码~
#include<map>
#include<cmath>
#include<ctime>
#include<queue>
#include<cstdio>
#include<vector>
#include<bitset>
#include<cstring>
#include<iostream>
#include<algorithm>
#define eps 1e-8
#define pa pair<int,int>
#define ll long long
#define inf 1000000000
using namespace std;
struct P{
double x, y;
P(){}
P(double _x, double _y):x(_x), y(_y){}
friend P operator-(P a, P b){
return P(a.x - b.x, a.y - b.y);
}
friend double operator*(P a, P b){
return a.x * b.y - a.y * b.x;
}
}p[105];
struct L{
P a, b;
}l[105];
vector<double> ve;
int n;
double R;
void add(double X, L l)
{
if(l.a.x < X - eps && l.b.x < X - eps)return;
if(l.a.x > X + eps && l.b.x > X + eps)return;
if(fabs(X - l.a.x) < eps && fabs(X - l.b.x) < eps)return;
double t;
if(fabs(l.b.x - X) < eps)t = l.b.y;
else t = (X - l.a.x) / (l.b.x - X);
ve.push_back(l.a.y + (l.b.y - l.a.y) * t / (t + 1));
}
double F(double X)
{
if(abs(X) > R)return 0;
ve.clear();
for(int i = 1; i <= n; i++)
add(X, l[i]);
double U = sqrt(R * R - X * X), ANS = 0;
sort(ve.begin(), ve.end());
bool flag = 0;
for(int i = 0; i < (int)ve.size() - 1; i++)
{
flag ^= 1;
if(ve[i] > U)break;
if(flag)
ANS += min(ve[i + 1], U) - ve[i];
}
return ANS;
}
double cal(double fl, double fm, double fr, double L)
{
return (fl + fm * 4 + fr) * L / 6;
}
double simpson(double l, double r, double fl, double fm, double fr, double S, int dep)
{
double mid = (l + r) / 2;
double m1 = (l + mid) / 2, m2 = (r + mid) / 2;
double f1 = F(m1), f2 = F(m2);
double g1 = cal(fl, f1, fm, mid - l), g2 = cal(fm, f2, fr, r - mid);
if(fabs(g1 + g2 - S) < 1e-12 && dep >= 10)return g1 + g2;
return simpson(l, mid, fl, f1, fm, g1, dep + 1) + simpson(mid, r, fm, f2, fr, g2, dep + 1);
}
double solve(double l, double r)
{
if(l > r)return 0;
double mid = (l + r) / 2;
return simpson(l, r, 0, F(mid), 0, (F(l) + F(r) + F(mid) * 4) * (r - l) / 6, 0);
}
int main()
{
scanf("%d%lf", &n, &R);
for(int i = 1; i <= n; i++)
scanf("%lf%lf", &p[i].x, &p[i].y);
p[n + 1] = p[1];
for(int i = 1; i <= n; i++)
l[i].a = p[i], l[i].b = p[i + 1];
double ANS = solve(-R, -2) + solve(-2, 1) + solve(1, R);
printf("%.10lf\n", ANS);
return 0;
}
也就是说这道题只需要积分求面积就好了,那么为什么这个F函数这么可怕(难以理解)呢?
因为题目说多边形 所以啊,他有可能长这个样:
我们再想想我们的积分的本质——把每个x对应的
f(x)
f
(
x
)
求和 即
∑nx=1f(x)
∑
x
=
1
n
f
(
x
)
当然我只是这么表示一下方便脑补。也就是说对应到图像上我们要求的是图像内每条竖线的长度的和,举个栗子吧:
假设我们要求
f(X)
f
(
X
)
即粉线处的函数值,那么我们只需要求
L红线+L绿线+L黄线
L
红
线
+
L
绿
线
+
L
黄
线
这也是黄学长代码的意图,我分步剖析一下:
先说ve数组,这里面存的是x=X与多边形的交点
然后说add:我们枚举每条线段,看是否有交点,及在何处
void add(double X, L l)
{
if(l.a.x < X - eps && l.b.x < X - eps)return;
if(l.a.x > X + eps && l.b.x > X + eps)return;
if(fabs(X - l.a.x) < eps && fabs(X - l.b.x) < eps)return;
double t;
if(fabs(l.b.x - X) < eps)t = l.b.y;
else t = (X - l.a.x) / (l.b.x - X);
ve.push_back(l.a.y + (l.b.y - l.a.y) * t / (t + 1));
}
上方的三个return 分别对应 :该线段在X左边,X右边,垂直x轴。这就比较好理解了,因为这三种状况不可以贡献面积(第三种可以用直线的面积为0解释)
然后就是这个点坐标,看图解释吧
假设我们要求R的纵坐标 即绿色长度,我们设为
LG
L
G
,代码中l.a.x,l.a.y,l.a.x,l.b.y与图中
x1,y1,x2,y2
x
1
,
y
1
,
x
2
,
y
2
一一对应;
然后我们做辅助线(米黄色),我们发现上方是两个相似三角形!所以我们利用x的比值来确定
LG−LR
L
G
−
L
R
,代码中pickline两侧的米黄线比值为t 那么
LG−LR
L
G
−
L
R
=
(y2−Y1)∗tt+1
(
y
2
−
Y
1
)
∗
t
t
+
1
这个不理解的话自己玩两组数据就好了。
扫描线
继续拿题举例子
BZOJ 1845 [Cqoi2005] 三角形面积并
题目大意:给定n个三角形,求面积并 n<=100
怎么做呢?
1.simpson积分,与第一个例题相似,我们可以求交点,然后求f(x)做(比上一个例题还好理解一点)。
2.扫描线
看图233
我们对每两条线段的交点打标记,然后过标记点做扫描线。
我们看 第一条扫描线与第二条之间只需要 求一个三角形,而第二三条之间是一个三角加一个梯形,那我们就可以用梯形公式求积啦~
我们只需要特判一下平行x轴的边,因为它有可能无贡献(留个读者考虑)
或者用中位线考虑(比较繁琐)
旋转卡壳
首先说一下读音(毕竟有
24
2
4
种读法),xuan zhuan qia ke (音调:2 3 3 2
http://blog.csdn.net/wang_heng199/article/details/74477738
这篇博客讲的真的好,不用我xbb了~
UPD:抱歉我没看到该文章是转载的,原文地址404 not found了