球体与轴对齐立方体重叠检测方案优化问询

我需要用八叉树表示球体，核心工作是判断轴对齐立方体与球体的重叠关系，具体分为三类：

• 立方体完全在球内，标记为full；
• 立方体完全在球外，标记为empty；
• 立方体部分在球内，需对其细分后递归处理。

以下是我基于C23标准的当前实现，现寻求更低复杂度、更少运算量的简化方案：

#include <OpenGL/OpenGL.h>
#include <math.h>

/* A 3D point with type-punning to access coordinates by name or index */
typedef union point_3d {
    struct {
        GLdouble x, y, z;
    } coord;

    GLdouble t[ 3 ];
} point_3d;

/* An axis-aligned cube */
typedef struct cube {
    /* corner with smallest coordinates */
    point_3d p_min;

    /* side length */
    GLdouble len;
} cube_t;

typedef struct sphere {
    unsigned int radius;
    point_3d center;
} sphere_t;

/* Clearly defined values for our result */
typedef enum overlap {
    no_overlap = -1,
    partial_overlap = 0,
    total_overlap = 1,
} overlap_t;

static inline GLdouble sq( GLdouble val ) {
    return val * val;
}

static bool point_in_sphere( point_3d p, struct sphere s ) {
    GLdouble dx = p.t[ 0 ] - s.center.t[ 0 ];
    GLdouble dy = p.t[ 1 ] - s.center.t[ 1 ];
    GLdouble dz = p.t[ 2 ] - s.center.t[ 2 ];
    return sq( s.radius ) >= sq( dx ) + sq( dy ) + sq( dz );
}

overlap_t cube_overlaps_sphere( cube_t c, sphere_t s ) {
    point_3d p_min = c.p_min;
    point_3d p_max = { .t= {
        p_min.t[0] + c.taille,
        p_min.t[1] + c.taille,
        p_min.t[2] + c.taille,
    } };
    point_3d c_centre = { .t = { 
        p_min.t[0] + c.taille / 2,
        p_min.t[1] + c.taille / 2, 
        p_min.t[2] + c.taille / 2,
    } };
    point_3d v = { .t = {
        c_centre.t[ 0 ] - s.centre.t[ 0 ],
        c_centre.t[ 1 ] - s.centre.t[ 1 ],
        c_centre.t[ 2 ] - s.centre.t[ 2 ],}
    };
    point_3d far_corner = { .t = { 
        v.t[ 0 ] < 0 ? p_min.t[ 0 ] : p_max.t[ 0 ],
        v.t[ 1 ] < 0 ? p_min.t[ 1 ] : p_max.t[ 1 ],
        v.t[ 2 ] < 0 ? p_min.t[ 2 ] : p_max.t[ 2 ],
    } };

    if ( point_in_sphere( far_corner, s ) ) {
        return total_overlap;
    }

    /* Algorithm by Jim Arvo */
    GLdouble dmin = 0;
    for ( unsigned int i = 0; i < 3; i++ ) {
        if ( s.center.t[ i ] < p_min.t[ i ] ) {
            dmin += pow( s.center.t[ i ] - p_min.t[ i ], 2 );
        } else if ( s.center.t[ i ] > p_max.t[ i ] ) {
            dmin += pow( s.center.t[ i ] - p_max.t[ i ], 2 );
        }
    }
    if ( dmin <= pow( s.radius, 2 ) ) {
        return partial_overlap;
    } else {
        return no_overlap;
    }
}

注：当前方案先通过判断立方体最远角是否在球内识别完全重叠，否则使用Jim Arvo算法检测是否存在重叠。

优化后的实现方案 #

核心优化点 #

1. 替换低效运算：用内联平方计算替代pow函数（pow是通用幂函数，计算平方的开销远高于直接乘法）；
2. 减少冗余变量：去掉不必要的临时结构体，直接计算坐标值，降低内存操作开销；
3. 提前预计算：一次性计算球半径的平方、立方体最大坐标等，避免重复计算；
4. 修正类型与笔误：将球体半径改为GLdouble避免类型转换开销，修正原代码中taille的笔误（对应结构体的len）；
5. 展开循环逻辑：将X/Y/Z轴的计算分开，减少循环控制开销，更利于编译器优化。

优化后代码 #

#include <OpenGL/OpenGL.h>
#include <stdbool.h>
#include <math.h>

/* 3D点：支持按名称或索引访问坐标 */
typedef union point_3d {
    struct {
        GLdouble x, y, z;
    } coord;
    GLdouble t[3];
} point_3d;

/* 轴对齐立方体 */
typedef struct cube {
    point_3d p_min;  // 最小坐标顶点
    GLdouble len;    // 边长
} cube_t;

typedef struct sphere {
    GLdouble radius;  // 统一为GLdouble，避免类型转换
    point_3d center;
} sphere_t;

/* 重叠类型枚举 */
typedef enum overlap {
    no_overlap = -1,
    partial_overlap = 0,
    total_overlap = 1,
} overlap_t;

// 内联计算平方，比pow高效数倍
static inline GLdouble sq(GLdouble val) {
    return val * val;
}

overlap_t cube_overlaps_sphere(cube_t c, sphere_t s) {
    const GLdouble half_len = c.len * 0.5;
    // 预计算立方体最大坐标
    const GLdouble p_max_x = c.p_min.t[0] + c.len;
    const GLdouble p_max_y = c.p_min.t[1] + c.len;
    const GLdouble p_max_z = c.p_min.t[2] + c.len;
    // 预计算立方体中心坐标
    const GLdouble cube_center_x = c.p_min.t[0] + half_len;
    const GLdouble cube_center_y = c.p_min.t[1] + half_len;
    const GLdouble cube_center_z = c.p_min.t[2] + half_len;

    // 计算球心到立方体最远顶点的距离平方
    const GLdouble dx = (cube_center_x < s.center.t[0]) ? (p_max_x - s.center.t[0]) : (s.center.t[0] - c.p_min.t[0]);
    const GLdouble dy = (cube_center_y < s.center.t[1]) ? (p_max_y - s.center.t[1]) : (s.center.t[1] - c.p_min.t[1]);
    const GLdouble dz = (cube_center_z < s.center.t[2]) ? (p_max_z - s.center.t[2]) : (s.center.t[2] - c.p_min.t[2]);
    const GLdouble max_dist_sq = sq(dx) + sq(dy) + sq(dz);
    const GLdouble radius_sq = sq(s.radius);

    // 判断完全重叠：最远顶点在球内则整个立方体在球内
    if (max_dist_sq <= radius_sq) {
        return total_overlap;
    }

    // Jim Arvo算法：计算球心到立方体的最小距离平方
    GLdouble min_dist_sq = 0.0;
    // X轴判断
    if (s.center.t[0] < c.p_min.t[0]) {
        min_dist_sq += sq(s.center.t[0] - c.p_min.t[0]);
    } else if (s.center.t[0] > p_max_x) {
        min_dist_sq += sq(s.center.t[0] - p_max_x);
    }
    // Y轴判断
    if (s.center.t[1] < c.p_min.t[1]) {
        min_dist_sq += sq(s.center.t[1] - c.p_min.t[1]);
    } else if (s.center.t[1] > p_max_y) {
        min_dist_sq += sq(s.center.t[1] - p_max_y);
    }
    // Z轴判断
    if (s.center.t[2] < c.p_min.t[2]) {
        min_dist_sq += sq(s.center.t[2] - c.p_min.t[2]);
    } else if (s.center.t[2] > p_max_z) {
        min_dist_sq += sq(s.center.t[2] - p_max_z);
    }

    return (min_dist_sq <= radius_sq) ? partial_overlap : no_overlap;
}

内容的提问来源于stack exchange，提问作者Ter Maxima