计算几何及其应用——立体几何

 

  承接《计算几何及其应用——解析几何》一文中对平面几何的探索，这一片文章将开始对计算几何中的空间几何进行探索。

  我们先来看一个简单的小题目。(Problem source:hdu1140)

  题目大意：这里给出k个进攻卫星，和地球上m个点，进攻卫星可以进攻从该点出发做直线能够到达的任意点，问你这k个卫星公能够打到多少个点。
 
  数理分析：这里显然将地球抽象成球体，卫星是球外的任意一点。我们过卫星向地球做无数的切线，我们连接此点和球心，然后我们沿着垂直于这条直线的方向(这里可以旋转360°)看地球，无论如何，直线和球体相交都变成了直线和圆相交，而这个圆是球体中截面最大的那个圆。
  分析到这里，我们就不难发现，卫星能够触及的最远点就是平面中的切点，即：满足dis(卫星到目标) <= dis(卫星到球心) - r*r(地球半径)，卫星就可以打到该目标。

  编程实现上：考虑到精度问题，我们可以用反三角函数acos得到PI的精确值，随后再设置循环穷举所有情况。

  参考代码如下。
  

#include<math.h>
#include<stdio.h>

#define PI acos(-1.0)
using namespace std;
struct point
{
    double x , y ,z;
}p[110];
double dis(point a , point b)
{
     return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z));
}

int main()
{
    int m , k ;
    int i , j;
    double r = 40000.0/2/PI;
    while(scanf("%d%d" , &k , &m) && (k || m))
    {

          for(i = 0;i < k;i++)
             scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
          point target,o;
          o.x= 0 , o.y= 0 , o.z=0;

          int sum = 0;
          for(i = 0;i < m;i++)
          {
              scanf("%lf%lf%lf",&target.x,&target.y,&target.z);


                for(j = 0;j < k;j++)
                {
                   double dis1 = dis(o,p[j]);
                   double dis2 = dis(target,p[j]);
                     if(dis2*dis2 <= dis1*dis1 - r*r)
                         break;
                }
                if(j < k)
                     sum++;
          }
      printf("%d\n",sum);
    }
}

 

  在立体几何当中，我们常常讨论的一个几何体就是四面体，下面，我们来探讨一下如何通过四面体的六个边来得到该四面体的体积。由于该方法是欧拉开发，所以又叫做欧拉四面体公式。

  我们先从四面体的基本公式，V = 1/3sh开始进行转化，并引入向量来进行表示。

  我们得到了四面体体积用向量计算式，此时我们再用行列式的计算方法，对向量计算式转换成和点坐标有关的行列式。


    随后我们再将行列式平方，根据行列式转置(第i行变成第i列)后值不变的性质，对其中一个行列式进行转置，这是为了使点坐标在后面的计算中和棱长联系起来。
  随后在再需根据余弦定理和三阶行列式相乘的法则进行数学推导即可，最终将得到以下答案。

  值得注意的是，这里方程的推导是严格的边对应的，因为是在坐标系下给定点坐标的前提下进行推导，所以如果边的对应顺序不对，会引起错误。

  基于以上对欧拉四面体公式的推导，我们来对一个具体的问题进行实践。(Problem source : hdu 1411)

  很明显的需要用欧拉四面体公式的问题，唯一需要注意的一点就是输入边的对应。
  在编程实现时，既可以选择直接套用公式，也可以采取分步间接的方式来计算。这里体现了计算几何的一个特点，就是需要从一个不同角度看同一个问题，数学方法使问题用数学语言表示得间接，而编程实现就需要用计算机语言使问题变得间接。这里时间缘故，暂时不给出编程实现中后一种方法的代码分析。
  参考代码如下(给出的是后一种方法，前一种方法直接套用公式即可，但还是要注意边的对应)。
 

 #include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
int main()
{
    double a,b,c,d,e,f;
    double t1,t2,t3,t4,temp,ans;
    while(scanf("%lf%lf%lf%lf%lf%lf",&a,&b,&c,&d,&e,&f) != EOF)
    {
         t1 = acos((c*c+b*b-f*f)/(2.0*b*c));
         t2 = acos((a*a+c*c-e*e)/(2.0*a*c));
         t3 = acos((a*a+b*b-d*d)/(2.0*a*b));
         t4 = (t1+t2+t3)/2.0;
         temp = sqrt(sin(t4)*sin(t4-t1)*sin(t4-t2)*sin(t4-t3));
         ans = a*b*c*temp/3.0;
         printf("%.4lf\n",ans);
    }


}

    我们再来看一道利用欧拉四面体公式的问题。(Problem source : pku 2208)
 

    很典型的欧拉四面体的问题，关于边的对应也是和上面的题目一样，这里就不重复贴出代码。

 

    在平面集合当中，我们探讨了凸包、最小圆覆盖等最小的图形覆盖类问题，而拿到三维空间，同样有三维凸包‘最小球包含等问题一一对应，接下来我们就将一一进行探讨分析。
  我们先来看最小球包含。
                         
  这里我们介绍模拟退火算法，来解决这个问题。

  模拟退火算法，本质上就是模拟了一种具体存在的物理过程——固体退火，并将这个过程中各个物理量抽象成了其他量而实现的一种随机化搜索的算法，因此，对于了解固体退火过程中的物理量和抽象化的模拟退火算法中的量的对应关系，是理解并应用模拟退火算法的关键。
  对于固体退火的物理过程，摘抄一下来自百度百科的介绍。
  模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis准则，粒子在温度T时趋于平衡的概率为e(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。
    这是从物理的角度给出的关于固体退火过程的解释。
 
   而模拟退火算法的原理，也就是物理过程和算法变量之间的对应关系，则如下。用固体退火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件S。
  
  下面则是更抽象化的模拟退火算法的模型。
(1) 初始化：初始温度T(充分大)，初始解状态S(是算法迭代的起点)，每个T值的迭代次数L

(2) 对k=1，……，L做第(3)至第6步：

(3) 产生新解S′

(4) 计算增量Δt′=C(S′)-C(S)，其中C(S)为评价函数

(5) 若Δt′<0则接受S′作为新的当前解，否则以概率exp(-Δt′/T)接受S′作为新的当前解.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。

终止条件通常取为连续若干个新解都没有被接受时终止算法。

(7) T逐渐减少，且T->0，然后转第2步。
  这里我们最终的目标就是得到C(S) = 0时的S的值，也就是我们模拟退火算法最终的解，也是对应着固体退火过程的常温的平衡态。显然，在(4)中，如果△t’<0，说明当前生成的新解S‘比S更加接近最终解，因此此时保留S’,否则的话，以概率exp的形式保留。这里关于变成实现概率exp形式保留当前解的编程实现，需要一定的数理上的推导，这里暂且不给出，笔者将在以后的文章中给出较详尽的解释。
  
  参考代码如下。
  #include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
using namespace std;
const int MAXN = 50;
const double EPS = 1e-7;
struct P{double x ,y,z;}p[MAXN],s;
int n;
double delta,ans;

void input()
{
     for(int i = 0;i < n;i++)  scanf("%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z);
}
double dis(P a , P b)
{
     return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
void solve()
{
     s.x=s.y=s.z=0;        
     delta=100,ans=1e20;    //delta模拟温度，ans则是解空间
     while(delta > EPS)
     {
          int d  = 0;
           for(int i = 0;i < n;i++)   //目标函数
               if(dis(s,p[i]) > dis(s,p[d]))  d = i;
        double md  = dis(s,p[d]);
        ans = min(ans,md);
        s.x += (p[d].x-s.x)/md*delta;  //概率形式保留当前解
        s.y += (p[d].y-s.y)/md*delta;
        s.z += (p[d].z-s.z)/md*delta;
        delta *= 0.98;//降温
     }
     printf("%.5lf\n",ans);
}
int main()
{
    while(scanf("%d",&n),n)
    {
         input();
         solve();
    }
    return 0;
}