1231 - Super Star

Aizu Online Judge

問題概要

空間上にn個の点がある。それらをすべて包含するような、最小の球の半径を求めよう。球面上に点があっても包含しているものとする。
ただし 4<=n<=30 とする。

コード

#include<iostream>
#include<vector>
#include<cmath>

using namespace std;

typedef double elem;

const elem eps = 1.0e-8;
const elem inf = 1.0e20;

bool eq(elem a, elem b){return abs(b-a)<eps;}
bool leq(elem a, elem b){return eq(a,b)||a<b;}

struct point{
  elem x,y,z;
  point(elem x, elem y, elem z):x(x),y(y),z(z){}
  point():x(0),y(0),z(0){}

  const point operator+(const point &t)const{
    return point(x+t.x,y+t.y,z+t.z);
  }
  const point operator-(const point &t)const{
    return point(x-t.x,y-t.y,z-t.z);
  }
  const point operator*(elem t)const{
    return point(t*x,t*y,t*z);
  }
  const point operator/(elem t)const{
    return point(x/t,y/t,z/t);
  }

  elem length()const{
    return sqrt(x*x+y*y+z*z);
  }
  elem dot(const point &t)const{
    return x*t.x+y*t.y+z*t.z;
  }
};

ostream &operator<<(ostream &os, const point &t){
  os << '(' << t.x << ',' << t.y << ',' << t.z << ')';
  return os;
}
const point far_point(inf,inf,inf);

typedef point vec;

struct sphere{
  point center;
  elem r;
  sphere ():center(0.0,0.0,0.0),r(0.0){}
  sphere (point c, elem r):center(c),r(r){}
  bool inside(const point &t)const{
    return leq( (t-center).length(), r );
  }
};

inline point centerPoint(point a, point b)
{
  return (a+b)/2;
}

inline point circumscribedCircleCenter(point a, point b, point c)
{
  vec va=b-a;
  vec vb=c-a;
  elem A=va.length();
  elem B=vb.length();
  elem C=va.dot(vb)/(A*B);
  elem ds=2*A*pow(C,2.0)-2*A;
  elem dt=2*B*pow(C,2.0)-2*B;
  if( eq(ds,0.0) || eq(dt,0.0) ) return far_point;
  elem s=(B*C-A)/(2*A*pow(C,2.0)-2*A);
  elem t=(A*C-B)/(2*B*pow(C,2.0)-2*B);
  return a+va*s+vb*t;
}

inline point circumscribedTetrahedronCenter(point a, point b, point c, point d)
{
  vec B=b-a;elem la=B.length();
  vec C=c-a;elem lb=C.length();
  vec D=d-a;elem lc=D.length();
  elem delta = B.x*(C.y*D.z-D.y*C.z)+C.x*(D.y*B.z-B.y*D.z)+D.x*(B.y*C.z-C.y*B.z);
  if( eq( delta, 0.0 ) ) return far_point;
  elem x=(B.y*C.z-C.y*B.z)*pow(lc,2.0)+(D.y*B.z-B.y*D.z)*pow(lb,2.0)+(C.y*D.z-D.y*C.z)*pow(la,2.0);
  elem y=(C.x*B.z-B.x*C.z)*pow(lc,2.0)+(B.x*D.z-D.x*B.z)*pow(lb,2.0)+(D.x*C.z-C.x*D.z)*pow(la,2.0);
  elem z=(B.x*C.y-C.x*B.y)*pow(lc,2.0)+(D.x*B.y-B.x*D.y)*pow(lb,2.0)+(C.x*D.y-D.x*C.y)*pow(la,2.0);
  x/=2;y/=2;z/=2;
  return a+point(x,y,z)/delta;
}

elem solve(const vector<point> &vp)
{
  elem answer = inf;
  // nC2
  for(int i = 0; i < vp.size(); ++i){
    for(int j = i + 1; j < vp.size(); ++j){
      bool isContainingBall = true;
      sphere s(centerPoint(vp[i],vp[j]),(vp[i]-vp[j]).length()/2);
      for(int k = 0; k < vp.size(); ++k){
	if( !s.inside( vp[k] ) ){
	  isContainingBall = false;
	  break;
	}
      }
      if( isContainingBall ){
	answer = leq( s.r, answer ) ? s.r : answer;
      }
    }
  }

  // nC3
  for(int i = 0; i < vp.size(); ++i){
    for(int j = i + 1; j < vp.size(); ++j){
      for(int k = j + 1; k < vp.size(); ++k){
	bool isContainingBall = true;
	point cent = circumscribedCircleCenter(vp[i],vp[j],vp[k]);
	sphere s(cent, (vp[i]-cent).length() );
	for(int l = 0; l < vp.size(); ++l){
	  if( !s.inside( vp[l] ) ){
	    isContainingBall = false;
	    break;
	  }
	}
	if( isContainingBall ){
	  answer = leq( s.r, answer ) ? s.r : answer;
	}
      }
    }
  }
  
  // nC4
  for(int i = 0; i < vp.size(); ++i){
    for(int j = i + 1; j < vp.size(); ++j){
      for(int k = j + 1; k < vp.size(); ++k){
	for(int l = k + 1; l < vp.size(); ++l ){
	  bool isContainingBall = true;
	  point cent = circumscribedTetrahedronCenter(vp[i],vp[j],vp[k],vp[l]);
	  sphere s(cent, (vp[i]-cent).length() );
	  for(int m = 0; m < vp.size(); ++m){
	    if( !s.inside( vp[m] ) ){
	      isContainingBall = false;
	      break;
	    }
	  }
	  if( isContainingBall ){
	    answer = leq( s.r, answer ) ? s.r : answer;
	  }
	}
      }
    }
  }
  
  return answer;
}

int main()
{
  while(true){
    int n;
    vector<point> vp;

    cin >> n;
    if( n == 0 ) break;

    for(int i = 0; i < n; ++i){
      elem x,y,z;
      cin >> x >> y >> z;
      vp.push_back( point(x,y,z) );
    }
    
    printf("%.7lf\n", solve(vp));

  }
  return 0;
}

解法

全探索する。
任意の異なる2点を直径とする球、3点の外心を大円とする球、4点が成す四面体の外接球を求める。
2点を直径とする球は自明。
3点の外心を大円とする球の中心の位置ベクトルpは、各頂点の位置ベクトルを0,a,bとすれば、
p=sa+tb, s = (BC-A)/(2AC2-2A), t = (AC-B)/(2BC2-2B)
A=nrm(a),B=nrm(b),C=a・b/AB, nrm(x) は x の大きさ
となる。
また、四面体の方は、上記同様、各頂点の位置ベクトルを0,a,b,cとすれば
正方行列 A=(a b c) として
nrm(p) = nrm(p-a) = nrm(p-b) = nrm(p-c)
より
tAp=0.5t(|a|2 |b|2 |c|2)
p=0.5 tA-1 t(|a|2 |b|2 |c|2)
より比較的簡単な以下の式が得られる。
http://twitpic.com/2ac056/full
また、4点のうちを1点を原点に置かず、代数的な方法で計算すると以下のような式になる。
http://twitpic.com/29xcbg/full
参考: http://homepage2.nifty.com/PAF00305/math/triangle/node7.html

その他の解法