対極幾何学の基礎

理論的背景

ピンホールカメラモデルでは、3次元空間から2次元画像への変換時に奥行き情報が失われる。このため、単眼画像からは各ピクセルがカメラからどれだけ離れているかを特定できない。このような深度情報を回復するには複数の視点が必要であり、両目視覚のように二つのカメラを用いる方法が立体視として知られている。

複数視点幾何学の基本概念を理解するために、まず対極幾何学について考察する。二つのカメラで同一シーンを撮影した場合の基本的な構成を考える。

片方のカメラのみでは、画像上の点に対応する3次元座標を特定できない。なぜなら、視線方向の直線上にあるすべての点が同じ2次元位置に投影されるからである。しかし、もう一方のカメラの情報を考慮すると、異なる3次元点は右側の画像平面上で異なる位置(x')に投影される。この原理により、二つの画像を使って正しい3次元点を三角測量で求めることができる。

対極制約

OX直線上の異なる点群が右側の画像平面上で直線(linel')を形成し、これを対極線(エピポーラルライン)と呼ぶ。この意味は、右側の画像上で対応点を探す際、この対極線上を探索すればよいというものである。これにより、画像全体を検索する必要がなくなり、より効率的かつ正確なマッチングが可能になる。すべての点は他方の画像に対応する対極線を持ち、その集合が対極面を構成する。

OとO'はそれぞれ左右のカメラ中心である。左側の画像上には右側カメラO'の投影が見られ、これが主点(エピポール)となる。主点はカメラ中心を通る直線と画像平面との交点であり、逆に右側のカメラに対しても同様に主点が存在する。場合によっては画像外に主点が存在することもある。

すべての対極線は主点を通るため、多数の対極線を求めそれらの交点を見つけることで主点の位置を決定できる。

基本行列と基礎行列

対極線や主点を求めるには、**基礎行列(F)と基本行列(E)**という二つの重要な要素が必要になる。基礎行列は二つのカメラ間の並進と回転情報を含み、グローバル座標系での相対配置を記述する。

基本行列は基礎行列と同じ情報を保持しつつ、二つのカメラの内部パラメータも含むため、二つのカメラの画素座標を関連付けることができる（較正済み画像を使用し、焦点距離で標準化された場合はF=Eとなる）。基本行列Fは一つの画像中の点を他方の画像中の直線にマッピングする。これは二つの画像の対応点から計算され、少なくとも8つの対応点が必要である（8点アルゴリズム使用時）。より多くの点を選択しRANSACを使用することで信頼性の高い結果を得られる。

対極線描画の実装

まず、基礎行列を求めるために二つの画像間の可能な限り多くの対応点を発見する必要がある。ここではSIFT特徴量とFLANNベースのマッチャーおよび比率テストを組み合わせて使用する。

以下は対極線を描画するPythonコードの実装例である：

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt

# 画像の読み込み
image_left = cv.imread('./opencv_data/left_image.jpg', 0)
image_right = cv.imread('./opencv_data/right_image.jpg', 0)

# SIFT特徴抽出器の初期化
feature_detector = cv.SIFT_create()

# キーポイントと記述子の計算
keypoints_left, descriptors_left = feature_detector.detectAndCompute(image_left, None)
keypoints_right, descriptors_right = feature_detector.detectAndCompute(image_right, None)

# FLANNマッチングパラメータ設定
flann_index_type = 1
index_params = dict(algorithm=flann_index_type, trees=5)
search_params = dict(checks=50)
matcher = cv.FlannBasedMatcher(index_params, search_params)

# 特徴量マッチング
all_matches = matcher.knnMatch(descriptors_left, descriptors_right, k=2)

# データ構造の初期化
valid_matches = []
left_points = []
right_points = []

# Loweの比率テストによるフィルタリング
for match_pair in all_matches:
    if len(match_pair) >= 2:
        first_match, second_match = match_pair
        if first_match.distance < 0.8 * second_match.distance:
            valid_matches.append(first_match)
            right_points.append(keypoints_right[first_match.trainIdx].pt)
            left_points.append(keypoints_left[first_match.queryIdx].pt)

# データ型変換
left_points = np.int32(left_points)
right_points = np.int32(right_points)

# 基礎行列の計算
fundamental_matrix, inlier_mask = cv.findFundamentalMat(left_points, right_points, cv.RANSAC, 4, 0.999)

# 内点のみ選択
left_inliers = left_points[inlier_mask.ravel() == 1]
right_inliers = right_points[inlier_mask.ravel() == 1]

def render_epipolar_lines(img_a, img_b, epipolar_lines, points_a, points_b):
    height, width = img_a.shape
    img_a_colored = cv.cvtColor(img_a, cv.COLOR_GRAY2BGR)
    img_b_colored = cv.cvtColor(img_b, cv.COLOR_GRAY2BGR)
    
    for line, point_a, point_b in zip(epipolar_lines, points_a, points_b):
        random_color = tuple(np.random.randint(0, 255, 3).tolist())
        y_intercept = int(-line[2] / line[1]) if line[1] != 0 else 0
        slope_factor = -(line[2] + line[0] * width) / line[1] if line[1] != 0 else 0
        
        img_a_colored = cv.line(img_a_colored, (0, y_intercept), (width, int(slope_factor)), random_color, 1)
        img_a_colored = cv.circle(img_a_colored, tuple(point_a), 5, random_color, -1)
        img_b_colored = cv.circle(img_b_colored, tuple(point_b), 5, random_color, -1)
        
    return img_a_colored, img_b_colored

# 右画像の対応点から対極線を左画像に描画
epipolar_lines_left = cv.computeCorrespondEpilines(right_inliers.reshape(-1, 1, 2), 2, fundamental_matrix)
epipolar_lines_left = epipolar_lines_left.reshape(-1, 3)
result_img_left, result_img_right = render_epipolar_lines(image_left, image_right, epipolar_lines_left, left_inliers, right_inliers)

# 左画像の対応点から対極線を右画像に描画
epipolar_lines_right = cv.computeCorrespondEpilines(left_inliers.reshape(-1, 1, 2), 1, fundamental_matrix)
epipolar_lines_right = epipolar_lines_right.reshape(-1, 3)
result_img_right_rev, result_img_left_rev = render_epipolar_lines(image_right, image_left, epipolar_lines_right, right_inliers, left_inliers)

plt.figure(figsize=(12, 6))
plt.subplot(121), plt.imshow(result_img_left), plt.title('左画像')
plt.subplot(122), plt.imshow(result_img_right_rev), plt.title('右画像')
plt.tight_layout()
plt.show()

左側の画像ではすべての対極線が右側の画像外の一点で収束していることが確認できる。これは主点の位置を示している。より良い結果を得るには、高解像度で多くの非平面ポイントを含む画像を使用すべきである。

ステレオ画像からの深度マップ生成

理論的基礎

対極制約などの基本概念を確認した。同一シーンの二つの画像があれば、直感的に深度情報を得ることが可能であることを示した。以下の簡単な数学的関係式でこの考え方が証明される：

視差 = x - x' = B×f / Z

ここでxとx'はカメラ中心からシーン中の3次元点に対応する画像平面上の点までの距離である。Bは二つのカメラ間の距離（既知）、fはカメラの焦点距離（既知）である。つまり、シーン中の点の深度は、対応する画像点とカメラ中心との距離差に反比例する。この情報を利用して、画像中のすべてのピクセルの深度を導き出すことができる。

OpenCVによる視差マップ生成

StereoBMアルゴリズムには調整が必要なパラメータがいくつかある：

texture_threshold: 信頼性の高いマッチングに不十分なテクスチャを持つ領域を除外 speckle_sizeとspeckle_range: ブロックベースのマッチングで発生するアーチファクトを除去するための後処理パラメータ num_disparities: スライディングウィンドウのピクセル数。値が大きいほど深度範囲が広くなるが計算負荷も増加 min_disparity: 左ピクセルのx位置からの探索開始オフセット uniqueness_ratio: 最適なマッチング視差が探索範囲内の他の視差に対して十分に優れているかを判定する後処理

以下は視差マップを生成するコード例：

import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt

# 画像読み込み
left_image = cv.imread('./opencv_data/stereo_left.png', 0)
right_image = cv.imread('./opencv_data/stereo_right.png', 0)

# StereoBMインスタンス作成
stereo_matcher = cv.StereoBM_create(numDisparities=32, blockSize=15)
computed_disparity = stereo_matcher.compute(left_image, right_image)

# 結果表示
display_images = [left_image, right_image, computed_disparity]
labels = ['左画像', '右画像', '視差']

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for idx, (img, label) in enumerate(zip(display_images, labels)):
    axes[idx].imshow(img, cmap='gray')
    axes[idx].set_title(label)
    axes[idx].axis('off')
plt.tight_layout()
plt.show()

視差マップから3次元再構成

以下は視差マップを用いて3次元点群を生成する完全な実装：

from __future__ import print_function
import cv2 as cv
import numpy as np

# PLYファイル形式のヘッダ定義
point_cloud_format = '''ply
format ascii 1.0
element vertex %(vertex_count)d
property float x
property float y
property float z
property uchar red
property uchar green
property uchar blue
end_header
'''

def save_point_cloud(filename, vertices, rgb_colors):
    """PLY形式で点群データを保存"""
    reshaped_vertices = vertices.reshape(-1, 3)
    reshaped_colors = rgb_colors.reshape(-1, 3)
    combined_data = np.hstack([reshaped_vertices, reshaped_colors])
    
    with open(filename, 'wb') as output_file:
        header_content = point_cloud_format % dict(vertex_count=len(combined_data))
        output_file.write(header_content.encode('utf-8'))
        np.savetxt(output_file, combined_data, fmt='%f %f %f %d %d %d')

def execute_stereo_reconstruction():
    """ステレオ再構成の実行"""
    print('画像を読み込んでいます...')
    
    # 処理速度向上のために画像サイズを縮小
    left_input = cv.pyrDown(cv.imread(cv.samples.findFile('./opencv_data/left_sample.jpg')))
    right_input = cv.pyrDown(cv.imread(cv.samples.findFile('./opencv_data/right_sample.jpg')))
    
    # SGBMパラメータ設定
    window_param = 7
    min_disparity_value = 16
    total_disparities = 128 - min_disparity_value
    
    stereo_algorithm = cv.StereoSGBM_create(
        minDisparity=min_disparity_value,
        numDisparities=total_disparities,
        blockSize=window_param,
        P1=8 * 3 * window_param ** 2,
        P2=32 * 3 * window_param ** 2,
        disp12MaxDiff=1,
        uniquenessRatio=15,
        speckleWindowSize=100,
        speckleRange=32
    )
    
    print('視差計算中...')
    disparity_result = stereo_algorithm.compute(left_input, right_input).astype(np.float32) / 16.0
    
    print('3次元点群生成中...')
    img_height, img_width = left_input.shape[:2]
    focal_length = 0.8 * img_width  # 焦点距離の推定
    
    # 再投影行列Qの定義
    reprojection_matrix = np.float32([
        [1, 0, 0, -0.5 * img_width],
        [0, -1, 0, 0.5 * img_height], 
        [0, 0, 0, -focal_length],
        [0, 0, 1, 0]
    ])
    
    # 3次元座標への再投影
    point_cloud_3d = cv.reprojectImageTo3D(disparity_result, reprojection_matrix)
    color_values = cv.cvtColor(left_input, cv.COLOR_BGR2RGB)
    
    # 有効な視差値を持つマスク作成
    valid_mask = disparity_result > disparity_result.min()
    extracted_points = point_cloud_3d[valid_mask]
    extracted_colors = color_values[valid_mask]
    
    output_filename = './output/pointcloud_output.ply'
    save_point_cloud(output_filename, extracted_points, extracted_colors)
    print(f'{output_filename} が保存されました')
    
    # 結果表示
    cv.imshow('左画像', left_input)
    normalized_disparity = (disparity_result - min_disparity_value) / total_disparities
    cv.imshow('視差マップ', normalized_disparity)
    cv.waitKey()
    print('処理完了')

if __name__ == '__main__':
    print(__doc__)
    execute_stereo_reconstruction()
    cv.destroyAllWindows()