Welzl のアルゴリズム

点集合の最小抱合円(点集合を周を含む内部に持つ最小の円)を O(N) の時間で求める Welzl のアルゴリズムというのがある。
論文が http://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf で読める。
論文通りに実装すると、点の数が多いときにスタックがあふれるので、再帰が深くならないようにしてみた。
言語は C。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define DIM 2
#define N 10000 /* テストデータの点の数。main でのみ参照 */
/* 集合 P のための領域を確保するとき余分に確保しておくサイズ */
/* 100万点でも300程度あれば十分(小さくても動くが、realloc が何度も呼ばれてしまう) */
#define DPSIZE 300
#define sq(x) ((x) * (x))

typedef struct {
    double r; /* 半径(計算の途中では使用しない) */
    double r2; /* 半径の自乗 */
    double center[DIM];
} Disk;

/* 集合 P のための構造体 */
typedef struct {
    const double **a; /* 点の座標へのポインタの配列を指す */
    int size; /* 確保したサイズ */
    int last; /* 配列の最も後ろの点のインデックス */
} p_type;

/* 0〜1 の一様乱数 */
double uniform_random(void)
{
    return rand() / (RAND_MAX + 1.0);
}

/* point が disk に含まれているか */
int inside(const double *point, const Disk *disk)
{
    return sq(point[0] - disk->center[0]) + sq(point[1] - disk->center[1]) <= disk->r2;
}

const Disk *make_disk(double r2, double x, double y)
{
    static Disk disk;
    disk.r2 = r2;
    disk.center[0] = x;
    disk.center[1] = y;
    return &disk;
}

/* 2点の最小抱合円 */
const Disk *calc_disk2(const double *const* r)
{
    return make_disk((sq(r[0][0] - r[1][0]) + sq(r[0][1] - r[1][1])) / 4,
    (r[0][0] + r[1][0]) / 2, (r[0][1] + r[1][1]) / 2);
}

/* 3点を周上に持つ円 */
const Disk *calc_disk3(const double *const* r)
{
    double ux, uy, vx, vy, zx, zy, c;
    ux = r[1][0] - r[0][0];
    uy = r[1][1] - r[0][1];
    vx = r[2][0] - r[0][0];
    vy = r[2][1] - r[0][1];
    c = (ux * vx + uy * vy) / (ux * vy - uy * vx);
    zx = (ux + vx + c * (uy - vy)) / 2;
    zy = (uy + vy + c * (vx - ux)) / 2;
    return make_disk(sq(zx) + sq(zy), r[0][0] + zx, r[0][1] + zy);
}

/* 集合 R から円を求める(rsize は R の要素数) */
const Disk *calc_disk(const double *const* r, int rsize)
{
    switch (rsize) {
        case 0:
            return make_disk(0, 0, 0);
        case 1:
            return make_disk(0, r[0][0], r[0][1]);
        case 2:
            return calc_disk2(r);
        default:
            return calc_disk3(r);
    }
}

/* 点を配列の末尾に移動 */
void move(int i, p_type *p)
{
    /* 必要なら配列のサイズを拡張 */
    if (++p->last == p->size)
        p->a = (const double**) realloc(p->a, sizeof(double**) * (p->size += DPSIZE));
    p->a[p->last] = p->a[i];
    p->a[i] = NULL; /* 抜き跡は NULL にしておく */
}

void shuffle(const p_type *p)
{
    int i, j;
    const double *temp;
    for (i = p->last; i > 0; i--) {
        j = (i + 1) * uniform_random();
        temp = p->a[i]; p->a[i] = p->a[j]; p->a[j] = temp;
    }
}

/* {p->a[k+1], …, p->a[p->last]} を集合 P、{r[0], …, r[rsize-1]} を集合 R として */
/* P∪R の最小抱合円 D(P,R) を求める (P, R の元は点の座標へのポインタとして実装している) */
const Disk *b_md(p_type *p, int k, const double **r, int rsize)
{
    const Disk *disk;
    int i;
    /* R の最小抱合円を disk とする */
    disk = calc_disk(r, rsize);
    if (rsize > DIM) return disk;
    /* 下のループカウンタ i の処理が終了した時点で disk は */
    /* R ∪ {p->a[i], …, p->a[p->last]} の最小抱合円になる */
    for (i = p->last; i > k; i--) {
        /* 下の !p->a[i] は点の抜き跡をスキップするため */
        if (!p->a[i] || inside(p->a[i], disk)) continue;
        r[rsize] = p->a[i];
        disk = b_md(p, i, r, rsize + 1);
        move(i, p); /* move-to-front heuristic */
    }
    return disk;
}

/* n 点の最小抱合円を求める */
/* (coord[i][0], coord[i][1], ..., coord[i][DIM-1]) が i 番目の点の座標 */
Disk smallest_enclosing_disk(const double coord[][DIM], int n)
{
    p_type p;
    const double *r[DIM + 1];
    Disk disk;
    int i;
    p.a = (const double**) malloc(sizeof(double**) * (p.size = n + DPSIZE));
    for (i = 0; i < n; i++)
        p.a[i] = coord[i];
    p.last = n - 1;
    shuffle(&p);
    disk = *b_md(&p, -1, r, 0);
    free(p.a);
    disk.r = sqrt(disk.r2);
    return disk;
}

/* 原点中心の単位円が最小抱合円になるようなテストデータを生成 */
/* n (≧3) は点の数 */
void generate_test_data(double coord[][2], int n)
{
    double t[3], temp;
    int i, j;

    /* 単位円周上に、鋭角三角形の頂点になるような3点を取る */
    t[0] = 2 * uniform_random();
    t[1] = uniform_random();
    t[2] = uniform_random();
    if (t[1] < t[2]) {temp = t[1]; t[1] = t[2]; t[2] = temp;}
    t[1] += t[0];
    t[2] += t[0] - 1;
    for (i = 0; i < 3; i++) {
        coord[i][0] = cos(M_PI * t[i]);
        coord[i][1] = sin(M_PI * t[i]);
    }

    /* 単位円盤上に n-3 点をばらまく */
    for (i = 3; i < n; i++)
        do {
            for (j = 0; j < DIM; j++)
                coord[i][j] = 2 * uniform_random() - 1;
        } while (sq(coord[i][0]) + sq(coord[i][1]) > 1);
}

int main(void)
{
    static double coord[N][DIM];
    Disk disk;
    srand(time(NULL));
    generate_test_data(coord, N);
    disk = smallest_enclosing_disk(coord, N);
    printf("center = (%f, %f)\n", disk.center[0], disk.center[1]);
    printf("radius = %f\n", disk.r);
    return 0;
}

プログラムはテストデータの生成もしているので、このままコンパイルして実行できる。
中心が原点、半径が1という結果が帰ってくれば、(そのテストデータについては)正常に動いているということになる。

DIM は次元数。DIM=2、つまり平面のときしか動かない書き方になってるけど、Welzl のアルゴリズムは次元数に依存しないので、calc_disk (とその下請け)と inside のところを次元数に合せて適当に書き換えれば一般の次元数でも動くようになっている(はず)。

b_md が Welzl のアルゴリズムの本体で、スタックがあふれないようにループで実装しているので論文との対応が分かりにくいが、論文の B_MTFDISK に対応している。
やっているのは以下のようなこと。

点集合 P と R が与えられていて、R が P∪R の最小抱合円の円周上にあることが分かっているとする(計算の最初では R は空集合)。
P の元は配列に並べておく。
P∪R の最小抱合円 D(P,R) を次のようにして求める。

最初 disk を R の最小抱合円として、最終的に disk が P∪R の最小抱合円になるように disk を更新していく。
P の元 q を配列の最後尾から見ていき、以下を繰り返す。

q が disk の内部または周上にあるなら何もしない。
q が disk の外にあるなら、P' を q より(配列上で)後ろの点の集合、R'=R∪{q} とすると、R' は P'∪R' の最小抱合円の円周上にあるので*1、disk = D(P',R') とする。

この処理が P の先頭まで達したら D(P,R) = disk となる。

b_md の move(i, p) のところが move-to-front heuristic。この行は削っても動くが、円盤上にランダムにばら撒いたデータだと Welzl のアルゴリズムの部分の実行時間が1.5倍程度になる。

if (rsize > DIM) から始まる1行も普通のデータでは削れて効率にもほとんど影響しないが(但し move-to-front heuristic と一緒に削ると Welzl のアルゴリズムの部分の実行時間が2倍程度になる)、最小抱合円の円周上に多くの点が乗っているようなデータでおかしくなるので削らないほうがいい。

smallest_enclosing_disk の shuffle(&p) の行も削っても動く。ここは配列上の点の並びが不都合なものだと(例えば最小抱合円の周に近い点が配列の最初に固まっている場合)Welzl のアルゴリズムが線形時間の効率を発揮できなくなるため、点を配列上でシャッフルしているのだが、手元の環境では大抵のデータでシャッフルは Welzl のアルゴリズム自体と同程度の時間がかかる。データの並びが特に偏ったものでない限り、この行は削ってしまっていいかもしれない。

*1:証明は論文にある。最初読んだときこの事実を使っていることが分からなくて苦労した。