请求讲解raycasting中torus相交检测函数的实现方法

你已经实现了球体的射线相交检测函数，现在需要类似的圆环检测函数，下面一步步拆解实现思路和代码：

圆环的数学定义 #

圆环可以看作是半径为r的小圆沿着半径为R的大圆的圆心轨迹旋转而成的曲面。我们默认圆环位于局部XY平面，中心在原点（如果需要平移或旋转，可先对射线做坐标变换，将圆环转换到局部坐标系后再计算）。

射线与圆环相交的核心推导 #

射线的参数方程为：P(t) = O + t*D，其中：

• O 是射线原点
• D 是射线的单位方向向量（建议单位化，能大幅简化计算）
• t > 0 是我们要求的参数，代表射线原点到交点的距离

圆环的隐式方程为：

(sqrt(x² + y²) - R)² + z² = r²

展开后得到：

(x² + y² + z² + R² - r²)² = 4R²(x² + y²)

将射线参数方程代入隐式方程，会得到一个关于t的四次方程：At⁴ + Bt³ + Ct² + Dt + E = 0。解这个方程就能得到交点对应的t值，取最小的正t就是射线与圆环的最近交点；如果没有正根，说明射线与圆环不相交。

可编码的系数计算步骤 #

先计算几个中间变量简化运算：

• dot_O_O = dot(O, O)：射线原点到圆环中心的距离平方
• dot_O_D = dot(O, D)：射线原点向量与方向向量的点积
• dot_D_D = dot(D, D)：射线方向向量的模长平方（单位化后为1）
• dot_Oxy_Oxy：射线原点在XY平面上的分量模长平方
• dot_Oxy_Dxy：射线原点XY分量与方向向量XY分量的点积
• dot_Dxy_Dxy：射线方向向量XY分量的模长平方（单位化后小于等于1）

基于这些变量，四次方程的系数为：

A = dot_D_D * dot_D_D
B = 2 * dot_D_D * q ，其中 q = 2 * dot_O_D
C = q*q + 2*dot_D_D*p - 4*R²*dot_Dxy_Dxy ，其中 p = dot_O_O + R² - r²
D_coeff = 2*q*p - 2*s ，其中 s = 2*R²*dot_Oxy_Dxy
E = p*p - t_coeff ，其中 t_coeff = 2*R²*dot_Oxy_Oxy

完整代码实现 #

辅助函数：解四次方程 #

#include <cmath>
#include <limits>

using namespace std;

// 解四次方程 At⁴ + Bt³ + Ct² + Dt + E = 0，返回最小正根，无正根返回-1
float solve_quartic(float A, float B, float C, float D, float E) {
    const float eps = 1e-6;
    if (fabs(A) < eps) return -1.0f;

    // 归一化系数，让A=1
    float invA = 1.0f / A;
    B *= invA;
    C *= invA;
    D *= invA;
    E *= invA;

    // Ferrari方法转化为三次方程求解
    float p = C - 3*B*B/8;
    float q = B*B*B/8 - B*C/2 + D;
    float r = E - B*D/4 + B*B*C/16 - 3*pow(B,4)/256;

    float a3 = 1.0f;
    float b3 = p / 2.0f;
    float c3 = (p*p)/16.0f - r/4.0f;
    float d3 = -q*q/64.0f;

    float delta0 = b3*b3 - 3*a3*c3;
    float delta1 = 2*pow(b3,3) - 9*a3*b3*c3 + 27*pow(a3,2)*d3;
    float delta = delta1*delta1 - 4*pow(delta0,3);

    float y;
    if (delta >= -eps) {
        float sqrt_delta = sqrt(max(delta, 0.0f));
        float temp = cbrt((delta1 + sqrt_delta)/2.0f) + cbrt((delta1 - sqrt_delta)/2.0f);
        y = temp - b3/(3*a3);
    } else {
        float theta = acos(delta1/(2*sqrt(pow(delta0, 3))));
        y = 2*sqrt(delta0/3)*cos(theta/3) - b3/(3*a3);
    }

    // 计算四个根
    float sqrt_y1 = sqrt(max(y*y - r, 0.0f));
    float sqrt_y2 = sqrt(max(2*y - p, 0.0f));

    float t1 = (-B/4) + (sqrt_y1 + sqrt_y2)/2.0f;
    float t2 = (-B/4) + (sqrt_y1 - sqrt_y2)/2.0f;
    float t3 = (-B/4) - (sqrt_y1 + sqrt_y2)/2.0f;
    float t4 = (-B/4) - (sqrt_y1 - sqrt_y2)/2.0f;

    // 筛选最小正根
    float min_t = numeric_limits<float>::infinity();
    if (t1 > eps && t1 < min_t) min_t = t1;
    if (t2 > eps && t2 < min_t) min_t = t2;
    if (t3 > eps && t3 < min_t) min_t = t3;
    if (t4 > eps && t4 < min_t) min_t = t4;

    return min_t == numeric_limits<float>::infinity() ? -1.0f : min_t;
}

圆环射线检测函数 #

// vec3 需包含x/y/z成员，以及dot、operator-等方法
struct vec3 {
    float x, y, z;
    vec3(float x=0, float y=0, float z=0) : x(x), y(y), z(z) {}
    vec3 operator-(const vec3& other) const { return vec3(x-other.x, y-other.y, z-other.z); }
};

float dot(const vec3& a, const vec3& b) {
    return a.x*b.x + a.y*b.y + a.z*b.z;
}

// 检测射线与圆环的相交：返回最近交点的t值，无交点返回-1
// 参数：center=圆环中心，R=大圆半径，r=小圆半径，ray_origin=射线原点，ray_dir=射线方向（建议单位化）
float hit_torus(vec3 center, float R, float r, vec3 ray_origin, vec3 ray_dir) {
    // 转换到圆环局部坐标系
    vec3 O = ray_origin - center;
    vec3 D = ray_dir;

    // 计算中间变量
    float dot_O_O = dot(O, O);
    float dot_O_D = dot(O, D);
    float dot_D_D = dot(D, D);

    vec3 O_xy(O.x, O.y, 0.0f);
    vec3 D_xy(D.x, D.y, 0.0f);
    float dot_Oxy_Oxy = dot(O_xy, O_xy);
    float dot_Oxy_Dxy = dot(O_xy, D_xy);
    float dot_Dxy_Dxy = dot(D_xy, D_xy);

    float p = dot_O_O + R*R - r*r;
    float q = 2.0f * dot_O_D;
    float s = 2.0f * R*R * dot_Oxy_Dxy;
    float t_coeff = 2.0f * R*R * dot_Oxy_Oxy;

    // 四次方程系数
    float A = dot_D_D * dot_D_D;
    float B = 2.0f * dot_D_D * q;
    float C = q*q + 2.0f * dot_D_D * p - 4.0f * R*R * dot_Dxy_Dxy;
    float D_coeff = 2.0f * q * p - 2.0f * s;
    float E = p*p - t_coeff;

    return solve_quartic(A, B, C, D_coeff, E);
}

注意事项 #

• 射线方向单位化：如果射线方向向量已经单位化，dot_D_D=1，能简化系数计算，提升运行效率。
• 精度控制：代码中的eps阈值可以根据你的渲染需求调整，避免因浮点精度问题导致错误检测。
• 旋转/平移适配：如果需要检测非XY平面的圆环，可以先对射线做旋转变换，将圆环转换到局部XY平面后再调用上述函数。

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