這個題的題意是給定一個凸多邊形表示的海島,求海島離大海最遠的距離。可以轉化為一個凸多邊形內部最多能夠放入一個多大的圓。
顯然可以對圓的半徑進行二分,但是怎麼確定圓心了。確定是否存在圓心,可以把原來的凸多邊形往內部移動r(圓的半徑)的距離之後,
再對新的多邊形求半平面交,如果半平面交存在(是個點即可),那麼當前大小的圓能夠放入。
求半平面交的算法可以用上一篇中的N*N復雜度的基本算法。本題還涉及到一個知識,就是如何把一條直線往逆時針方向或者順時針方向
移動R的距離。其實,可以根據單位圓那種思路計算。因為相當於以原來直線上的一點為圓心,以r為半徑做圓,而且與原來的直線成90的夾
角,那麼後來點的坐標是((x0 + cos(PI / 2 +θ )),(y0 + sin(PI / 2 + θ))),轉化一下就是(x0 - sinθ,y0 + cosθ)。那麼直接可以
求出dx = / (vp[i].y - vp[(i + 1) % nN].y) * fR / fDis,dy = (vp[(i + 1) % nN].x - vp[i].x) * fR / fDis,fDis是線段的長度。
代碼如下:
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <vector>
using namespace std;
const double fPre = 1e-8;
struct Point
{
double x,y;
Point(){}
Point(const Point& p){x = p.x, y = p.y;}
Point(double fX, double fY):x(fX), y(fY){}
Point& operator+(const Point& p)
{
x += p.x;
y += p.y;
return *this;
}
Point& operator+=(const Point& p)
{
return *this = *this + p;
}
Point& operator-(const Point& p)
{
x -= p.x;
y -= p.y;
return *this;
}
Point& operator*(double fD)
{
x *= fD;
y *= fD;
return *this;
}
};
typedef vector<Point> Polygon;
int DblCmp(double fD)
{
return fabs(fD) < fPre ? 0 : (fD > 0 ? 1 : -1);
}
double Cross(Point a, Point b)
{
return a.x * b.y - a.y * b.x;
}
double Det(double fX1, double fY1, double fX2, double fY2)
{
return fX1 * fY2 - fX2 * fY1;
}
double Cross(Point a, Point b, Point c)
{
return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}
Point Intersection(Point a1, Point a2, Point b1, Point b2)
{
Point a = a2 - a1;
Point b = b2 - b1;
Point s = b1 - a1;
return a1 + a * (Cross(b, s) / Cross(b, a));
}
Polygon Cut(Polygon& pg, Point a, Point b)
{
Polygon pgRet;
int nN = pg.size();
for (int i = 0; i < nN; ++i)
{
double fC = Cross(a, b, pg[i]);
double fD = Cross(a, b, pg[(i + 1) % nN]);
if (DblCmp(fC) >= 0)
{
pgRet.push_back(pg[i]);
}
if (DblCmp(fC * fD) < 0)
{
pgRet.push_back(Intersection(a, b, pg[i], pg[(i + 1) % nN]));
}
}
//printf("pgRet number:%d\n", pgRet.size());
return pgRet;
}
double Dis(Point a, Point b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
//返回半平面的頂點個數
int HalfPlane(Polygon& vp, double fR)
{
Polygon pg;
pg.push_back(Point(-1e9, -1e9));
pg.push_back(Point(1e9, -1e9));
pg.push_back(Point(1e9, 1e9));
pg.push_back(Point(-1e9, 1e9));
int nN = vp.size();
for (int i = 0; i < nN; ++i)
{
double fDis = Dis(vp[i], vp[(i + 1) % nN]);
double dx = (vp[i].y - vp[(i + 1) % nN].y) * fR / fDis;
double dy = (vp[(i + 1) % nN].x - vp[i].x) * fR / fDis;
Point a = vp[i], b = vp[(i + 1) % nN], c(dx, dy);
a += c;
b += c;
//printf("%f %f %f %f\n", a.x, a.y, b.x, b.y);
pg = Cut(pg, a, b);
if (pg.size() == 0)
{
return 0;
}
}
return pg.size();
}
int main()
{
int nN;
vector<Point> vp;
while (scanf("%d", &nN), nN)
{
vp.clear();
Point p;
for (int i = 0; i < nN; ++i)
{
scanf("%lf%lf", &p.x, &p.y);
vp.push_back(p);
}
double fMin = 0.0, fMax = 10000.0;
while (DblCmp(fMin - fMax))
{
double fMid = (fMin + fMax) / 2;
int nRet = HalfPlane(vp, fMid);
//printf("fMid:%f, nRet:%d\n", fMid, nRet);
if (nRet == 0)
{
fMax = fMid;
} www.2cto.com
else
{
fMin = fMid;
}
}
printf("%.6f\n", fMax);
}
return 0;
}