球と線分の交点

計算
コード

計算

直線の式は
P1 = P0 + V*t

それぞれの成分に分解すると
X = P.x + V.x * t
Y = P.y + V.y * t
Z = P.z + V.z * t

球の式は
X^2 + Y^2 + Z^2 = R^2

それぞれの成分を球の式に代入してtで整理する

(P.x + V.x * t)^2 + (P.y + V.y * t)^2 + (P.z + V.z * t)^2 = R^2

(V.x^2 + V.y^2 + V.z^2) t^2 + 2(P.x * V.x + P.y * V.y + P.z * V.z) t + (P.x^2 + P.y^2 + P.z^2 – R^2)

解の公式を使う

ax^2 + bx + c = 0 のとき
x = (-b ± sqrt(b^2 – 4ac)) / 2a
※b^2 – 4acがマイナスのときは解なし

a = (V.x^2 + V.y^2 + V.z^2)
b = 2(P.x * V.x + P.y * V.y + P.z * V.z)
c = (P.x^2 + P.y^2 + P.z^2 – R^2)

コード

// 球の中心座標
vector offset = set(0, 0, 0);

// 直線の座標を入力
vector p0 = point(1,"P", 0) - offset;
vector p1 = point(1,"P", 1) - offset;

vector v = normalize(p0 - p1);

// 半径
float r = 1;

float a = v.x * v.x + v.y * v.y + v.z * v.z;
float b = 2 * (p0.x * v.x + p0.y * v.y + p0.z * v.z);
float c = p0.x * p0.x + p0.y * p0.y + p0.z * p0.z - r * r;

// 解があるかチェック
if((b * b - 4 * a * c) >= 0)
{
    float t = (-b + sqrt(b * b - 4 * a * c)) / 2 * a;
    vector cross0 = p0 + v * t;
    
    // 直線上にあるかチェック
    if(dot(p0 - cross0, p1 - cross0) < 0)
    {
        cross0 += offset;
        int pt = addpoint(0, cross0);
        setpointattrib(0, "N", pt, cross0);
        setpointattrib(0, "Cd", pt, set(1, 0.5, 0));
    }
    
    
    t = (-b - sqrt(b * b - 4 * a * c)) / 2 * a;
    vector cross1 = p0 + v * t;
    
    // 直線上にあるかチェック
    if(dot(p0 - cross1, p1 - cross1) < 0)
    {
        cross1 += offset;
        int pt = addpoint(0, cross1);
        setpointattrib(0, "N", pt, cross1);
        setpointattrib(0, "Cd", pt, set(1, 0.5, 0));
    }
}