アルゴリズムって難しい

最小包含球に関して

教えている専門学校生は、ゲーム制作を主として勉強しているのだけれども、前回講義が終わってから質問があった。当たり判定を行う上で、キャラクターのポリゴンを包むエリアを球で表したいけれども、その最小の球の中心と半径を求めたいとな…。
ふむふむ、当然、そういうケースも起こりうるわなぁと。

「で、今はどんな方法を？」
「とりあえず、頂点の平均値から中心を求めているが、当然最小ではなくて、他の方法を探しています」
「平面であれば外接円を求めることになるよね？」
「それは分かりますが、空間で行う場合のアルゴリズムが分からないので…」

調べましたがな…

即答するのには、こちらも手持ちのネタがない。
「1週間待ってくれ。ヒントぐらいはあげられる筈だから…」（本当か？）と思いつつ、調査することに。

いろいろな方法が有るんだけれども、ゲームということは実行時間が短く、メモリをあまり食わない方法が最適。この条件で良さそうなものを探すことに。見つけたのが、
「点の集合を包含する球」という文献。
ここにも、様々な方法が紹介されているものの、最後のものが良さそう。

ということで、実際にコーディングして確かめることに。
条件は変なデータはないという前提でエラーチェックも最低限にして…。

#include <math.h>
#include <stdio.h>
#define INPUT_FILE "stars-s.txt"
#define DATANUM 30
typedef struct {
double x, y, z; /* coordinates */
} pos;
int n=DATANUM;
pos point[DATANUM];
pos center;
double distance(pos p,pos pp)
{
    double rr;
    return((pp.x-p.x)*(pp.x-p.x)+(pp.y-p.y)*(pp.y-p.y)+(pp.z-p.z)*(pp.z-p.z));
}
double solve_by_movement(void)
{
    int k,i,t;
    pos p;
    double move,max;
 
    p.x=50.0;
    p.y=50.0;
    p.z=50.0;
    move=0.5;
 
    while(move>1.0e-8){
        for(t=0;t<100;t++){
            max=0.0;
            for(i=0;i<n;i++)
                if(distance(p,point[i])>max){
                    max=distance(p,point[i]);
                    center = p;
                    k=i;
                }
            p.x += (point[k].x-p.x)*move;
            p.y += (point[k].y-p.y)*move;
            p.z += (point[k].z-p.z)*move;
        }
        move /=2.0;
    }
    return (max);
}
int main(void)
{
    int i;
double rr;
const char *path = INPUT_FILE;
if (freopen(path, "r", stdin) == NULL) {
perror(path);
return 1;
}
for (i = 0; i < n; i ++) {
if (scanf("%lf%lf%lf", &point[i].x,
 &point[i].y, &point[i].z) != 3) {
return 1;
}
}
rr=solve_by_movement();
printf("r=%.5f , (x,y,z)=(%.5f,%.5f,%.5f)\n", sqrt(rr),center.x,center.y,center.z);
return 0;
}

動作確認したが、取り敢えず正解っぽい。
まあ距離を2乗で求めるようにして、表示の時のみルートを取るということで簡略化。
あとは紹介されているアルゴリズム（コード）を流用してサクッと完了。

どの程度の処理か確認
動けばいいって問題ではないので、実際のコードの動きをチェック。

$ gcc smallest.c -lm -pg

で、プロファイルを取ることに。で、以下がその内容。（$ gprof --brief）

Flat profile:

Each sample counts as 0.01 seconds.
 no time accumulated

  %   cumulative   self              self     total        
 time   seconds   seconds    calls  Ts/call  Ts/call  name  
  0.00      0.00     0.00    87316     0.00     0.00  distance
  0.00      0.00     0.00        1     0.00     0.00  solve_by_movement

Call graph

granularity: each sample hit covers 4 byte(s) no time propagated

index % time    self  children    called     name
                     0.00    0.00   87316/87316       solve_by_movement [2]
[1]      0.0    0.00    0.00   87316                    distance [1]
-----------------------------------------------
                     0.00    0.00       1/1                    main [6]
[2]      0.0    0.00    0.00       1                         solve_by_movement [2]
                     0.00    0.00   87316/87316       distance [1]
-----------------------------------------------

Index by function name

   [1] distance                [2] solve_by_movement

距離を求める関数を87000回も呼び出してる…。
初期値を、空間の中心からスタートしているから、試行回数をもっと減らしても良いはず。ということで、t<100をt<20まで減らしても正確な答えを求められた。
すると、呼び出し回数は17456回に収まった。
あとは精度の問題でもっと荒くても良いのであれば、随分早くなるんだろうな。

PCが速すぎる

今回条件を変えたりしても、実質0.0秒以下で答えが求まってしまう。マシンのスペックが原因ですね。昔のような遅いマシンだったら、アルゴリズムや条件によって、答えが出るまでに随分差が有るんだろうけど。今時は多少無駄な方法を使っても、判断できないほど高速なので逆に困ってしまう。速くなったという実感すらわかない。
こういった、アルゴリズムを比較するためには、学習用の遅いマシンか、遅いエミュレータを用意する必要があるかも。

案外、教育現場なら需要があるのかな〜なんて思ったりして。誰か作りません？（笑）