第2章：3D幾何学と交差判定

はじめに

レイトレーシングの核心は、レイと3Dオブジェクトの交差判定です。本章では、各プリミティブ（球、平面、円柱）との交差判定アルゴリズムを学びます。

---

1. 交差判定の基礎

1.1 問題の定式化

レイと物体の交差点を求める問題：

レイ: P(t) = O + t * D
     O = 始点
     D = 方向（正規化）
     t = パラメータ（t > 0 で前方）

求めるもの:
- 交差するか？
- 交差点 P（最も近い t）
- 法線 N（照明計算に必要）

1.2 交差情報の構造体

typedef struct s_hit {
    int     valid;      // 交差したか
    double  t;          // パラメータ値
    t_vec3  point;      // 交差点
    t_vec3  normal;     // 法線（外向き）
    t_color color;      // 物体の色
    // ボーナス用
    double  u;          // テクスチャ座標
    double  v;
} t_hit;

// 初期化
t_hit hit_miss(void) {
    return (t_hit){.valid = 0, .t = INFINITY};
}

---

2. 球との交差判定

2.1 数学的導出

球の方程式：

|P - C|² = r²

P = 球上の点
C = 球の中心
r = 半径

レイの方程式を代入：

|O + t*D - C|² = r²

ベクトル OC = O - C として:
|OC + t*D|² = r²

展開:
(OC + t*D) · (OC + t*D) = r²
|OC|² + 2t(OC · D) + t²|D|² = r²

|D| = 1（正規化済み）なので:
t² + 2t(OC · D) + |OC|² - r² = 0

これは t の二次方程式:
at² + bt + c = 0

a = 1
b = 2(OC · D)
c = |OC|² - r²

判別式: Δ = b² - 4ac
       = 4(OC · D)² - 4(|OC|² - r²)
       = 4[(OC · D)² - |OC|² + r²]

2.2 実装

typedef struct s_sphere {
    t_vec3  center;
    double  radius;
} t_sphere;

t_hit sphere_intersect(t_sphere *sphere, t_ray *ray) {
    t_vec3 oc = vec3_sub(ray->origin, sphere->center);

    // 二次方程式の係数
    double a = vec3_dot(ray->direction, ray->direction);  // = 1
    double half_b = vec3_dot(oc, ray->direction);
    double c = vec3_dot(oc, oc) - sphere->radius * sphere->radius;

    // 判別式
    double discriminant = half_b * half_b - a * c;

    if (discriminant < 0)
        return hit_miss();

    // 最も近い交差点を選択
    double sqrt_d = sqrt(discriminant);
    double t = (-half_b - sqrt_d) / a;

    // t が負なら遠い方を試す
    if (t < EPSILON) {
        t = (-half_b + sqrt_d) / a;
        if (t < EPSILON)
            return hit_miss();
    }

    // 交差情報を構築
    t_hit hit;
    hit.valid = 1;
    hit.t = t;
    hit.point = ray_at(ray, t);
    hit.normal = vec3_normalize(
        vec3_sub(hit.point, sphere->center)
    );

    return hit;
}

2.3 法線の方向

// 法線が常に外向きになるよう調整
void fix_normal(t_hit *hit, t_ray *ray) {
    if (vec3_dot(ray->direction, hit->normal) > 0) {
        // レイが内側から当たっている
        hit->normal = vec3_scale(hit->normal, -1);
    }
}

外側からの入射:        内側からの入射:
       ↘                    ↗
        \   N              /
         \  ↑             /
          \ |            /
           \|           /|
            ●          ● |
           球          球 N（反転後）

---

3. 平面との交差判定

3.1 数学的導出

平面の方程式：

(P - P₀) · N = 0

P  = 平面上の点
P₀ = 平面上の既知の点
N  = 平面の法線

レイの方程式を代入：

(O + t*D - P₀) · N = 0
(O - P₀) · N + t(D · N) = 0

t = [(P₀ - O) · N] / (D · N)

3.2 実装

typedef struct s_plane {
    t_vec3  point;   // 平面上の点
    t_vec3  normal;  // 法線（正規化済み）
} t_plane;

t_hit plane_intersect(t_plane *plane, t_ray *ray) {
    double denom = vec3_dot(ray->direction, plane->normal);

    // レイが平面と平行（または後ろ向き）
    if (fabs(denom) < EPSILON)
        return hit_miss();

    t_vec3 p0_o = vec3_sub(plane->point, ray->origin);
    double t = vec3_dot(p0_o, plane->normal) / denom;

    // 後ろの交差は無視
    if (t < EPSILON)
        return hit_miss();

    t_hit hit;
    hit.valid = 1;
    hit.t = t;
    hit.point = ray_at(ray, t);
    hit.normal = plane->normal;

    // 法線の向きを調整
    if (denom > 0)
        hit.normal = vec3_scale(hit.normal, -1);

    return hit;
}

3.3 無限平面の描画

        N
        ↑
        |
━━━━━━━━●━━━━━━━━ 平面
       P₀

特徴:
- 無限に広がる
- チェッカーボードパターンでよく使われる

// ボーナス：チェッカーボードパターン
t_color checkerboard(t_hit *hit, double scale) {
    int x = (int)floor(hit->point.x * scale);
    int z = (int)floor(hit->point.z * scale);

    if ((x + z) % 2 == 0)
        return color_new(1, 1, 1);  // 白
    else
        return color_new(0, 0, 0);  // 黒
}

---

4. 円柱との交差判定

4.1 無限円柱

無限円柱の方程式（Z軸に沿った円柱）：

x² + y² = r²

一般的な向きの円柱：

|P - C - ((P - C) · A) * A|² = r²

C = 円柱の軸上の点
A = 軸の方向（正規化）
r = 半径

4.2 数学的導出

レイ: P(t) = O + t*D

V = O - C（始点から円柱中心へのベクトル）

円柱表面上の条件:
|V + t*D - ((V + t*D) · A) * A|² = r²

展開して整理すると二次方程式:
at² + bt + c = 0

a = |D - (D·A)A|²
b = 2(D - (D·A)A) · (V - (V·A)A)
c = |V - (V·A)A|² - r²

4.3 有限円柱の実装

typedef struct s_cylinder {
    t_vec3  center;     // 底面の中心
    t_vec3  axis;       // 軸の方向（正規化）
    double  radius;
    double  height;
} t_cylinder;

t_hit cylinder_intersect(t_cylinder *cyl, t_ray *ray) {
    // 軸方向の成分を除去
    t_vec3 v = vec3_sub(ray->origin, cyl->center);
    t_vec3 d_perp = vec3_sub(ray->direction,
        vec3_scale(cyl->axis, vec3_dot(ray->direction, cyl->axis)));
    t_vec3 v_perp = vec3_sub(v,
        vec3_scale(cyl->axis, vec3_dot(v, cyl->axis)));

    // 二次方程式の係数
    double a = vec3_dot(d_perp, d_perp);
    double b = 2.0 * vec3_dot(d_perp, v_perp);
    double c = vec3_dot(v_perp, v_perp) - cyl->radius * cyl->radius;

    double discriminant = b * b - 4 * a * c;
    if (discriminant < 0)
        return hit_miss();

    // 解を計算
    double sqrt_d = sqrt(discriminant);
    double t1 = (-b - sqrt_d) / (2 * a);
    double t2 = (-b + sqrt_d) / (2 * a);

    // 有限円柱の高さチェック
    t_hit hit = hit_miss();
    hit.t = INFINITY;

    // t1をチェック
    if (t1 > EPSILON) {
        t_vec3 p = ray_at(ray, t1);
        double h = vec3_dot(vec3_sub(p, cyl->center), cyl->axis);
        if (h >= 0 && h <= cyl->height) {
            hit.valid = 1;
            hit.t = t1;
        }
    }

    // t2をチェック（t1が無効な場合）
    if (!hit.valid && t2 > EPSILON) {
        t_vec3 p = ray_at(ray, t2);
        double h = vec3_dot(vec3_sub(p, cyl->center), cyl->axis);
        if (h >= 0 && h <= cyl->height) {
            hit.valid = 1;
            hit.t = t2;
        }
    }

    if (!hit.valid)
        return hit_miss();

    // 交差情報を構築
    hit.point = ray_at(ray, hit.t);

    // 法線を計算
    t_vec3 cp = vec3_sub(hit.point, cyl->center);
    double h = vec3_dot(cp, cyl->axis);
    t_vec3 axis_point = vec3_add(cyl->center,
                                  vec3_scale(cyl->axis, h));
    hit.normal = vec3_normalize(vec3_sub(hit.point, axis_point));

    return hit;
}

4.4 円柱のキャップ

// 上下のキャップ（円形の蓋）
t_hit cylinder_cap_intersect(t_cylinder *cyl, t_ray *ray, int is_top) {
    // キャップの平面
    t_vec3 cap_center = is_top ?
        vec3_add(cyl->center, vec3_scale(cyl->axis, cyl->height)) :
        cyl->center;
    t_vec3 cap_normal = is_top ? cyl->axis :
        vec3_scale(cyl->axis, -1);

    // 平面との交差
    double denom = vec3_dot(ray->direction, cap_normal);
    if (fabs(denom) < EPSILON)
        return hit_miss();

    t_vec3 p0_o = vec3_sub(cap_center, ray->origin);
    double t = vec3_dot(p0_o, cap_normal) / denom;

    if (t < EPSILON)
        return hit_miss();

    // 交差点が円内にあるか確認
    t_vec3 p = ray_at(ray, t);
    t_vec3 v = vec3_sub(p, cap_center);
    double dist2 = vec3_dot(v, v);

    if (dist2 > cyl->radius * cyl->radius)
        return hit_miss();

    t_hit hit;
    hit.valid = 1;
    hit.t = t;
    hit.point = p;
    hit.normal = cap_normal;

    return hit;
}

---

5. 交差判定の最適化

5.1 バウンディングボックス

複雑な物体の前に、単純な箱で判定：

typedef struct s_aabb {
    t_vec3 min;
    t_vec3 max;
} t_aabb;

int aabb_intersect(t_aabb *box, t_ray *ray) {
    double tmin = -INFINITY;
    double tmax = INFINITY;

    for (int i = 0; i < 3; i++) {
        double invD = 1.0 / (&ray->direction.x)[i];
        double t0 = ((&box->min.x)[i] - (&ray->origin.x)[i]) * invD;
        double t1 = ((&box->max.x)[i] - (&ray->origin.x)[i]) * invD;

        if (invD < 0) {
            double tmp = t0;
            t0 = t1;
            t1 = tmp;
        }

        tmin = t0 > tmin ? t0 : tmin;
        tmax = t1 < tmax ? t1 : tmax;

        if (tmax < tmin)
            return 0;
    }

    return tmax > 0;
}

5.2 バウンディング球

typedef struct s_bsphere {
    t_vec3 center;
    double radius;
} t_bsphere;

int bsphere_intersect(t_bsphere *bs, t_ray *ray) {
    t_vec3 oc = vec3_sub(ray->origin, bs->center);
    double b = vec3_dot(oc, ray->direction);
    double c = vec3_dot(oc, oc) - bs->radius * bs->radius;

    // b > 0 かつ c > 0 なら、レイは球の外で遠ざかっている
    if (b > 0 && c > 0)
        return 0;

    double discriminant = b * b - c;
    return discriminant >= 0;
}

5.3 早期終了

t_hit find_closest_hit(t_scene *scene, t_ray *ray) {
    t_hit closest = hit_miss();
    closest.t = INFINITY;

    for (int i = 0; i < scene->object_count; i++) {
        t_object *obj = &scene->objects[i];

        // バウンディング球で早期判定
        if (!bsphere_intersect(&obj->bsphere, ray))
            continue;

        t_hit hit;
        switch (obj->type) {
            case OBJ_SPHERE:
                hit = sphere_intersect(&obj->data.sphere, ray);
                break;
            case OBJ_PLANE:
                hit = plane_intersect(&obj->data.plane, ray);
                break;
            case OBJ_CYLINDER:
                hit = cylinder_intersect(&obj->data.cylinder, ray);
                break;
        }

        // より近い交差点を記録
        if (hit.valid && hit.t < closest.t) {
            closest = hit;
            closest.color = obj->color;
        }
    }

    return closest;
}

---

6. 数値安定性

6.1 イプシロン

浮動小数点誤差への対処：

#define EPSILON 1e-6

// 自己交差を防ぐ
t_ray create_shadow_ray(t_hit *hit, t_vec3 light_dir) {
    // 交差点をわずかに法線方向にオフセット
    t_vec3 origin = vec3_add(hit->point,
                             vec3_scale(hit->normal, EPSILON));

    return (t_ray){origin, light_dir};
}

6.2 t = 0 付近の問題

t_hit sphere_intersect_safe(t_sphere *sphere, t_ray *ray) {
    // ... 二次方程式を解く ...

    // 非常に小さいtは無視
    double t = (-half_b - sqrt_d) / a;
    if (t < EPSILON) {
        t = (-half_b + sqrt_d) / a;
        if (t < EPSILON)
            return hit_miss();
    }

    // ... 残りの処理 ...
}

---

まとめ

本章で学んだこと：

• 交差判定の基礎: レイと物体の交点、法線
• 球との交差: 二次方程式、判別式
• 平面との交差: 線形方程式
• 円柱との交差: 有限円柱、キャップ
• 最適化: バウンディングボリューム
• 数値安定性: イプシロン、自己交差

次章では、ライティングとシェーディングを学びます。