球と線分の交点

計算
コード
1. 関数

計算

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

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

球の式は
X² + Y² + Z² = R²

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

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

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

解の公式を使う

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

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

コード

// 球の中心座標
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));
    }
}

関数

//
// 円とポリラインの交差判定
//
int isCircleCross(int geometry; int primnum; vector position; float radius; export vector cross)
{
    int result = -1;
    
    int pts[] = primpoints(geometry, primnum);
    
    for(int i = 0; i < len(pts)-1; i++)
    {
        // 球の中心
        vector offset = position;
        
        vector p0 = point(geometry, "P", pts[i]) - offset;
        vector p1 = point(geometry, "P", pts[i+1]) - offset;
            
        vector v = normalize(p0 - p1);
        
        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 - radius * radius;
        
        // 解があるかチェック
        if((b * b - 4 * a * c) >= 0)
        {
            // 判別式±の+の場合
            float t = (-b + sqrt(b * b - 4 * a * c)) / 2 * a;
            cross = p0 + v * t;
            
            // 直線上にあるかチェック
            if(dot(p0 - cross, p1 - cross) < 0)
            {
                cross += offset;
                result = 1;
                break;
            }
            
            // 判別式±の-の場合
            t = (-b - sqrt(b * b - 4 * a * c)) / 2 * a;
            cross = p0 + v * t;
            
            // 直線上にあるかチェック
            if(dot(p0 - cross, p1 - cross) < 0)
            {
                cross += offset;
                result = 1;
                break;
            }
        }
    }
    
    return result;
}