计算几何(空间几何 + 三分) - Dome of Circus - UVA 1473

计算几何(空间几何 + 三分) - Dome of Circus - UVA 1473

题意：

$空间中给出n个点，要求用一个底面在xOy平面，顶点在z轴正半轴的圆锥包含这n个点。$

$要求这圆锥的体积最小，计算此时圆锥的底面半径和高。$

输入：

$多组测试数据，$

$每组首行包括一个正整数n，表示点的个数，$

$接着n行每行包括3个浮点数x,y,z，表示点的坐标。$

输出：

$两个浮点数，分别表示圆锥的高和底面半径。$

$小数点后保留三位。$

Sample Input

1
1.00 0.00 1.00
2
1.00 0.00 1.00
0.00 1.50 0.50
3
1.00 0.00 1.00
0.00 1.50 0.50
-0.50 -0.50 1.00

Sample Output

3.000 1.500
2.000 2.000
2.000 2.000

数据范围：

$1\le n\le 10000$

---

分析：

$不难发现，半径越大、高越大，就越容易包含所有点。$

$圆锥的体积应当满足一个单峰的凹函数，故我们考虑用三分法来求这个最小值。$

$我们三分半径，对于每次三分得到的两个半径r_1和r_2，$

$我们通过比例关系计算出，每个半径对应的，能够覆盖所有点的圆锥的高h_1和h_2，如下图$

$对于每一个点P(x,y,z)，利用相似三角形的性质，有如下等式：$

$$\frac{\sqrt{x^2+y^2}}{r}=\frac{h-z}{h}$$

$解得：$

$$h=\frac{z\cdot r}{r-\sqrt{x^2+y^2}}$$

那 么 h = m a x i = 1 n { z i ⋅ r r − x i 2 + y i 2 } 那么h=max_{i=1}^n\{\frac{z_i·r}{r-\sqrt{x_i^2+y_i^2}}\} 那么h=maxi=1n​{r−xi2​+yi2​ ​zi​⋅r​}

$最后三分输出极值点处的底面半径和高即可。$

代码：

#include<cstring>
#include<iostream>
#include<cmath>
#include<cstdio>
#include<algorithm>

using namespace std;

const double eps=1e-10;
const double pi=cos(-1.0);

int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return x<0 ? -1 : 1;
}

struct Point3
{
    double x,y,z;
    Point3(double x=0,double y=0,double z=0) : x(x), y(y), z(z) {}
    bool operator < (const Point3 &t) const
    {
        return z<t.z;
    }
};

int n;
Point3 V[10010];

double check(double r)
{
    double h=0;
    for(int i=0;i<n;i++)
    {
        double d=sqrt(V[i].x*V[i].x+V[i].y*V[i].y);
        h=max(h,V[i].z/(r-d)*r);
    }
    return h;
}

int main()
{
    while(~scanf("%d",&n))
    {
        double L=0,R=1e20;
        for(int i=0;i<n;i++) scanf("%lf%lf%lf",&V[i].x,&V[i].y,&V[i].z);

        while(R-L>eps)
        {
            double lmid=(L+R)/2, rmid=(lmid+R)/2;
            double h1=check(lmid), h2=check(rmid);
            double v1=1.0/3.0*pi*lmid*lmid*h1, v2=1.0/3.0*pi*rmid*rmid*h2;
            if(v1>v2) L=lmid;
            else R=rmid;
        }
        
        printf("%.3lf %.3lf\n",check(L),L);
    }
    return 0;
}