Ritterの境界球

3Dの当たり判定やcullingをする際に、境界球が必要となる場面は多く、この分野ではRitterの境界球という有名なアルゴリズムがあります。
原作はACMのGraphics Gemsのサイトから取得できるのですが、記法がとても古く、マルチスレッドも考慮されていないので、そのまま使うのはナンセンスです。
どこかに現代語訳(?)した物が無いかと探してみたのですが、無いようなので下記に作成しました。

/*
An Efficient Bounding Sphere
by Jack Ritter
from "Graphics Gems", Academic Press, 1990
*/

/* Routine to calculate tight bounding sphere over    */
/* a set of points in 3D */
/* This contains the routine find_bounding_sphere(), */
/* the struct definition, and the globals used for parameters. */
/* The abs() of all coordinates must be < BIGNUMBER */
/* Code written by Jack Ritter and Lyle Rains. */

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

const double BIGNUMBER = 100000000.0; // Please set very large coord that is out of the world.

struct Point3d {
    double x, y, z;
};

inline double sq(double x) { return x * x; }

void CalcBoundingSphere(const Point3d* pPoints, int nPoints,
        Point3d* pCenter, double* radius) {

    Point3d xmin = {0};
    Point3d xmax = {0};
    Point3d ymin = {0};
    Point3d ymax = {0};
    Point3d zmin = {0};
    Point3d zmax = {0};

    // FIRST PASS: find 6 minima/maxima pPoints
    xmin.x = ymin.y = zmin.z = BIGNUMBER; // initialize for min/max compare
    xmax.x = ymax.y = zmax.z = -BIGNUMBER;
    for (int i=0; i < nPoints; i++)
    {
        if (pPoints[i].x < xmin.x)
            xmin = pPoints[i]; // New xminimum point
        if (pPoints[i].x > xmax.x)
            xmax = pPoints[i];
        if (pPoints[i].y < ymin.y)
            ymin = pPoints[i];
        if (pPoints[i].y > ymax.y)
            ymax = pPoints[i];
        if (pPoints[i].z < zmin.z)
            zmin = pPoints[i];
        if (pPoints[i].z > zmax.z)
            zmax = pPoints[i];
    }

    // Set span = distance between the 2 pPoints *min & *max (squared)
    double xspan = sq(xmax.x - xmin.x) + sq(xmax.y - xmin.y) + sq(xmax.z - xmin.z);
    double yspan = sq(ymax.x - ymin.x) + sq(ymax.y - ymin.y) + sq(ymax.z - ymin.z);
    double zspan = sq(zmax.x - zmin.x) + sq(zmax.y - zmin.y) + sq(zmax.z - zmin.z);

    // Set pPoints dia1 & dia2 to the maximally separated pair
    Point3d dia1 = xmin;
    Point3d dia2 = xmax; // assume xspan biggest
    double maxspan = xspan;
    if (yspan > maxspan)
    {
        maxspan = yspan;
        dia1 = ymin; dia2 = ymax;
    }
    if (zspan > maxspan)
    {
        dia1 = zmin; dia2 = zmax;
    }


    // dia1,dia2 is a diameter of initial sphere
    // calc initial pCenter
    pCenter->x = (dia1.x+dia2.x)/2.0;
    pCenter->y = (dia1.y+dia2.y)/2.0;
    pCenter->z = (dia1.z+dia2.z)/2.0;

    // calculate initial radius**2 and radius

    double rad_sq = sq(dia2.x-pCenter->x) + sq(dia2.y-pCenter->y) + sq(dia2.z-pCenter->z);
    *radius = sqrt(rad_sq);

    /* SECOND PASS: increment current sphere */

    for (int i = 0; i < nPoints; i++)
    {
        double old_to_p_sq = sq(pPoints[i].x-pCenter->x) + sq(pPoints[i].y-pCenter->y) + sq(pPoints[i].z-pCenter->z);

        if (old_to_p_sq > rad_sq)   /* do r**2 test first */
        {   /* this point is outside of current sphere */
            double old_to_p = sqrt(old_to_p_sq);
            /* calc radius of new sphere */
            *radius = (*radius + old_to_p) / 2.0;
            rad_sq = sq(*radius);   /* for next r**2 compare */
            double old_to_new = old_to_p - *radius;
            /* calc pCenter of new sphere */
            pCenter->x = (*radius*pCenter->x + old_to_new * pPoints[i].x) / old_to_p;
            pCenter->y = (*radius*pCenter->y + old_to_new * pPoints[i].y) / old_to_p;
            pCenter->z = (*radius*pCenter->z + old_to_new * pPoints[i].z) / old_to_p;
        }
    }
}

// テスト用
void main(void) {
    Point3d p[] = { {0, 0, 0}, {1, 1, 1}, {1, 0, 1}, {1, -1, 0} };

    Point3d center;
    double r = 0;
    CalcBoundingSphere(p, 4, &center, &r);

    printf("Result (%lf,%lf,%lf) %lf\n", center.x, center.y, center.z, r);
}

上のコードのライセンスは元のGraphicsGemsのそれと同様です。
(自己責任という条件で、商用非商用問わず使用改変再配布自由。)

Copyright by MUTATE Systems 2018
管理者宛メール