視差画像算出方法とそれに伴うPythonコード

後日追記

「平行化」の章名を「位置合わせ」に変更し，合わせて各所の表現を修正しました．
また今回対象とする環境・データについての前提条件に関する記載を追加しました．

0. はじめに

この記事では，2枚の画像から1枚の視差画像を算出する方法の一例と，その場合に使用するPythonコードを簡単にまとめる．

用語の説明

Python

プログラミング言語．
文法が直感的かつ簡単で初心者でも扱いやすい．

Google Colaboratory

ブラウザ上でPythonが実行可能なGoogleのサービス．
機械学習の教育や研究のために無償で提供されている．
各種環境構築の必要がなく，手軽にPythonを実行可能．

OpenCV

画像処理用のライブラリ．
PythonやC++などで利用可能．

Numpy

配列演算用のライブラリ．
複雑な配列計算の高速処理が可能．

特徴量

画素や画像領域を特定するための指標のこと．
画像内の輝度の分布や勾配などをもとに算出されることが多い．
その算出方法は様々であり，それぞれのアルゴリズムに対応する検出器のことを特徴量記述子という．

特徴点 (key-point)

画像内で特に特徴量が大きいと判断された画素，あるいは領域．
この特徴点同士をマッチングして画像間の対応点関係を構築する．

視差 (dispairty)

視差とは，ある物体を2地点から撮影した際の像のズレのこと
下図の$d_p$に相当する．

人間の目も両目からの画像を脳が適切に重ね合わせ，視差を認識することで遠近感を把握する．
例えば遠くの山を見る時，手前にある建物などは2つにずれて見える．
これは脳が両目の画像中の山を重ね合わせた結果，山以外のものがずれて写るから．
簡単に言うとこのズレが視差．

この視差はカメラに近い物体ほど大きくなる．
この性質を利用して測量分野で使用されることが多い．

視差画像 (disparity image)

文字通り，視差を可視化した画像．
2枚の画像のうち片方の画像の各画素に写っている物体に対して視差を定義し画像化する．
詳しくは以下のMATLABのドキュメントを参照．

キャリブレーション (Calibration)

カメラレンズの僅かな歪みの影響を受けて，カメラに写される像も若干歪む．
歪んだままの画像を使用していると正しく画像処理ができない．
キャリブレーションとは，簡単に言うと画像の歪みを除去するためのパラメータを推定すること．
（厳密には歪み係数以外にも，画像主点位置などのカメラ内部パラメータなども推定される．）
実際にキャリブレーションを行う場合には，チェッカーボードなどの既知の模様を使用するカメラで異なる角度から数枚撮影し，それらの画像を解析することで歪み係数を推定する．

詳細な説明は以下を参照．

基本的な流れ

2枚の画像から視差画像を算出する基本的な流れは以下．

1. Google Colaboratoryの準備
2. カメラのキャリブレーション
3. 2枚の画像の撮影
4. 歪み補正
5. 特徴点抽出
6. 対応点マッチング
7. 位置合わせ
8. 視差画像算出

各ステップについて説明する．

1. Google Colaboratoryの準備

今回はPythonの実行環境としてGoogle Colaboratory (Colab)を想定．
Google Colaboratoryでは環境構築の必要がなく，手軽に利用しやすいため．
自分のPCで既にPythonの実行環境を構築できているのであればそちらを使用しても全く問題ない．

以下のリンクをブラウザで開くことでColabを利用可能．

以下の画面が開く．
「+ ノートブックを新規作成」をクリック．

今回はノートブック名を「Disparity_calculation」とする．

Google Colabは，PCのローカルに保存されている画像に直接アクセスできない．
そのため今回はGoogle Drive経由で画像を読み込む．

まずは画面左のファイルタブから「ドライブをマウント」を選択．

適当にアカウントなどを選ぶと，以下のようにMy driveが追加される．
これで自分のドライブ上の画像にアクセスできるようになった．

試しに1枚適当に画像をアップロードして確かめてみる．
画像のアップロードはアップロードしたいフォルダの右端のメニューからアップロードを選択して行う．
「削除されますよ」的な警告は無視．

アップロードしたファイルの右端のメニューから「パスをコピー」を選択．
（パスとは，ファイルの住所のようなもの．システムがファイルにアクセスする際に使用される．）

次に「+ コード」を選択して，新たなコードブロックを生成．

コードブロック内に以下のスクリプトを記入．
ただし，「file_path」の部分には先ほどコピーしたパスをペーストする．

.py

import cv2 as cv # OpenCV のインポート．
import matplotlib.pyplot as plt # Matplotlib のインポート．

img = cv.imread("file_path") # 画像の読込．
img = cv.cvtColor(img, cv.COLOR_BGR2RGB) # 画像を BGR から RGB に変換．
plt.show(plt.imshow(img)) # 画像の描画．

コードブロックの左上の実行ボタンを押すと，コードが実行され，以下のように画像が描画される．

2. カメラのキャリブレーション

今回はカメラとしてiPhoneのカメラを使用．
このカメラをキャリブレーションする．

キャリブレーションの基礎理論については以下を一読すると良い（再掲）．

チェッカーボードの撮影

まずキャリブレーションに使用するチェッカーボードを印刷．
OpenCVでは白黒のチェッカーボードや複数の円形模様などを用いたキャリブレーションが実装されているが，今回は以下のチェッカーボードを使用．

チェッカーボードが歪んでいては元も子もないので，ホワイトボードやアクリル板など平面を維持できるものに貼り付けて撮影する．モニターに写したものでもOK（モニターの周波数次第では画像上に干渉縞が発生する場合もあるため注意）．

様々な角度から撮影した10~20枚程度の画像を用意．
できるだけ画像内に大きくチェッカーボードが写るようにすると良い．
今回は以下のような15枚の画像を使用．

拡張子やファイル名を変更したい場合には以下のコードを使用．
好みの問題で別に変更する必要はない．

image_rename.py

# 画像の拡張子とファイル名を変更するコード．
import cv2 as cv # OpenCV をインポート．

# 親ディレクトリのパスを指定．
folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc/20231211_calibration_image"

counter = 0
for i in range(279, 294) : # 画像ファイル名に従い繰り返し処理．
    img = cv.imread(f"{folder_path}/IMG_0{i}.JPG") # 画像の読込．
    cv.imwrite(f"{folder_path}/{counter}.png", img) # 画像の出力．
    counter += 1 # 1 を加算．

テスト

複数枚のキャリブレーションを行う前に，まずは1枚のチェッカーボード画像から各頂点を検出し，各種パラメータを推定するテストコードを動かしてみる．
ここで推定されるパラメータの精度は低い．
使用するコードが以下．

calibration_one_img.py

.py

import cv2 as cv
import numpy as np
import glob
from matplotlib import pyplot as plt

# 画像ファイルの親ディレクトリのパスを指定．
folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc/20231211_calibration_image"

CHECKERBOARD = (7, 7) # チェッカーボードの頂点数を(縦, 横)で指定．

# 繰り返し計算終了の基準を設定．
criteria = (cv.TERM_CRITERIA_EPS + cv.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# 各頂点の三次元座標格納用の配列を作成．
objp = np.zeros((1, CHECKERBOARD[0]*CHECKERBOARD[1], 3), dtype = np.float32)
objp[0,:,:2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2)
prev_img_shape = None
# print(f"objp = \n{objp}") # 作成した配列を確認したい場合はコメントアウトを外す．

file_name = glob.glob(f"{folder_path}/*.png")[0] # ディレクトリ内のファイル名を一括取得し，その中の 0 番目を取得．
img = cv.imread(file_name) # 画像の読込．
gray = cv.cvtColor(img,cv.COLOR_BGR2GRAY) # 画像を BGR カラー画像から白黒画像に変換．
# img_for_disp = cv.cvtColor(img, cv.COLOR_BGR2RGB) # 画像を BGR から RGB に変換．
# plt.show(plt.imshow(img_for_disp)) # 読み込んだ画像を確認したい場合はコメントアウトを外す．

# チェッカーボードの頂点検出．
# ret には指定された数の頂点が検出された場合に True が格納される．
# corners には検出された頂点の座標が格納される．
ret, corners = cv.findChessboardCorners(gray, CHECKERBOARD, 
          cv.CALIB_CB_ADAPTIVE_THRESH+cv.CALIB_CB_FAST_CHECK+cv.CALIB_CB_NORMALIZE_IMAGE)

objpoints = [] # 三次元座標格納用空配列．
imgpoints = [] # 二次元画像座標格納用空配列．

if ret == True : # 頂点が適切に検出された場合．
    objpoints.append(objp) # objpoints に objp を追加．
    
    # 頂点画像座標を高精度化．
    corners2 = cv.cornerSubPix(gray, corners, (5,5), (-1,-1), criteria)
    
    imgpoints.append(corners2) # imgpoints に corners2 を追加．
    
    # 画像に高精度化された頂点を描画．
    img = cv.drawChessboardCorners(img, CHECKERBOARD, corners2, ret)
    img = cv.cvtColor(img, cv.COLOR_BGR2RGB) # 画像を BGR から RGB に変換．
    plt.show(plt.imshow(img)) # 読み込んだ画像の確認．
else : # 頂点が適切に検出されなかった場合．
    print("Corners were not detected.")

ret, mtx, dist, rvecs, tvecs = cv.calibrateCamera(objpoints, # 各頂点の三次元座標．
                                                  imgpoints, # 各頂点の二次元画像座標．
                                                  gray.shape[::-1],
                                                  None,
                                                  None
                                                  )

print(f"ret = \n{ret}")
print(f"mtr = \n{mtx}") # カメラ行列の出力．
print(f"dist = \n{dist}") # 歪み係数の出力．
print(f"rvecs = \n{rvecs[0]}") # 回転ベクトルの出力．
print(f"tvecs = \n{tvecs[0]}") # 並進ベクトルの出力．

img = cv.imread(file_name) # 歪みを除去した画像の作成用に再読込．
h, w, color = img.shape # 画像の縦，横の画素数を取得．
newcameramtx, roi = cv.getOptimalNewCameraMatrix(mtx, dist, (w,h), 1, (w,h)) # 最適なカメラ行列を取得．
dst = cv.undistort(img, mtx, dist, None, newcameramtx) # 画像の歪みを除去．
plt.imshow(cv.cvtColor(dst, cv.COLOR_BGR2RGB)) # OpenCV は色がGBR順なのでRGB順に並べ替える
plt.show()

mean_error = 0
for i in range(len(objpoints)) : # 頂点の数について繰り返し．
    imgpoints2, _ = cv.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv.norm(imgpoints[i],imgpoints2, cv.NORM_L2)/len(imgpoints2)
    mean_error += error

print("total error: ", mean_error/len(objpoints)) # 再投影誤差を出力．

なお，使用している関数の説明は以下などを参照．

結果が以下．
高解像度の画像（iPhoneのカメラはかなりの解像度）を使用したため，描画される色付きの線が相対的に細くなってしまったが，きちんと頂点が検出されている．

歪み補正を行なった画像が以下．

検出した頂点（外側の1列を除く6☓6の正方形）内では各辺が直線に近くなるように歪みが補正されている．

キャリブレーション

複数枚の画像を用いたキャリブレーションを行う．
使用するコードは以下．
先ほどのテスト処理を画像の枚数だけ繰り返す構造となっている．
また，注意点としては画像によってはかなり処理に時間がかかるものもある（画像1枚に10分程度）．
これは何らかの原因により繰り返し計算の結果がなかなか収束しない場合に起こると考えられる．
どうしてもこれを避けたい場合はfor文直後の部分で適宜画像番号を指定することで指定した画像の処理をスキップできる．
ちなみに筆者の場合は15枚のうち3枚はかなり時間がかかった．
また，Colabが固まってしまった場合はコードブロック左上のストップマークを何回か押せば強制停止できる．

calibration.py

.py

import cv2 as cv
import numpy as np
import glob
from matplotlib import pyplot as plt

# 画像ファイルの親ディレクトリのパスを指定．
folder_path = "/content/drive/MyDrive/研究資料/20231211_calibration_image"
image_num = 15 # 画像の枚数．
CHECKERBOARD = (7, 7) # チェッカーボードの頂点数を(縦, 横)で指定．

# 繰り返し計算終了の基準を設定．
criteria = (cv.TERM_CRITERIA_EPS + cv.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# 各頂点の三次元座標格納用の配列を作成．
# チェッカーボードの三次元座標は不変なので，この配列は使い回す．
objp = np.zeros((1, CHECKERBOARD[0]*CHECKERBOARD[1], 3), np.float32)
objp[0, :, :2] = np.mgrid[0:CHECKERBOARD[0], 0:CHECKERBOARD[1]].T.reshape(-1, 2)
prev_img_shape = None
# print(f"objp = \n{objp}") # 作成した配列を確認したい場合はコメントアウトを外す．
objpoints = [] # 三次元座標格納用空配列．
imgpoints = [] # 二次元画像座標格納用空配列．

for img_i in range(image_num) : # 画像の枚数だけ繰り返し．
#    if img_i in [2, 10, 13] : # 除外する画像を指定．
#        continue
    file_name = glob.glob(f"{folder_path}/*.png")[img_i] # ディレクトリ内のファイル名を一括取得し，その中の 0 番目を取得．
    img = cv.imread(file_name) # 画像の読込．
    gray = cv.cvtColor(img,cv.COLOR_BGR2GRAY) # 画像を BGR カラー画像から白黒画像に変換．
    # img_for_disp = cv.cvtColor(img, cv.COLOR_BGR2RGB) # 画像を BGR から RGB に変換．
    # plt.show(plt.imshow(img_for_disp)) # 読み込んだ画像を確認したい場合はコメントアウトを外す．

    # チェッカーボードの頂点検出．
    # ret には指定された数の頂点が検出された場合に True が格納される．
    # corners には検出された頂点の座標が格納される．
    ret, corners = cv.findChessboardCorners(gray, CHECKERBOARD,
              cv.CALIB_CB_ADAPTIVE_THRESH+cv.CALIB_CB_FAST_CHECK+cv.CALIB_CB_NORMALIZE_IMAGE)

    if ret == True : # 頂点が適切に検出された場合．
        objpoints.append(objp) # objpoints に objp を追加．

        # 頂点画像座標を高精度化．
        corners2 = cv.cornerSubPix(gray, corners, (5,5), (-1,-1), criteria)

        imgpoints.append(corners2) # imgpoints に corners2 を追加．

        # 画像に高精度化された頂点を描画．
        img = cv.drawChessboardCorners(img, CHECKERBOARD, corners2, ret)
        img = cv.cvtColor(img, cv.COLOR_BGR2RGB) # 画像を BGR から RGB に変換．
        plt.show(plt.imshow(img)) # 読み込んだ画像の確認．
        plt.close()
        print(f"In No.{img_i} image, corners were detected correctly.")
    else : # 頂点が適切に検出されなかった場合．
        print(f"--- In No.{img_i} image, corners were not detected. ---")
    # 繰り返し終了．

# 各パラメータの推定．
ret, mtx, dist, rvecs, tvecs = cv.calibrateCamera(objpoints, # 各画像・各頂点の三次元座標．
                                                  imgpoints, # 各画像・各頂点の二次元画像座標．
                                                  gray.shape[::-1],
                                                  None,
                                                  None
                                                  )

print(f"ret = \n{ret}")
print(f"mtr = \n{mtx}") # カメラ行列の出力．
print(f"dist = \n{dist}") # 歪み係数の出力．
print(f"rvecs = \n{rvecs[0]}") # 回転ベクトルの出力．
print(f"tvecs = \n{tvecs[0]}") # 並進ベクトルの出力．

# 各パラメータをドライブ上に保存．
np.savez(f"{folder_path}/camera_parameter.npz", ret, mtx, dist, rvecs, tvecs)

# 歪みを除去した画像を出力．
for img_i in range(image_num) : # 画像の枚数だけ繰り返し．
    img = cv.imread(glob.glob(f"{folder_path}/*.png")[img_i]) # 歪みを除去した画像の作成用に再読込．
    h, w, color = img.shape # 画像の縦，横の画素数を取得．
    newcameramtx, roi = cv.getOptimalNewCameraMatrix(mtx, dist, (w,h), 1, (w,h)) # 最適なカメラ行列を取得．
    dst = cv.undistort(img, mtx, dist, None, newcameramtx) # 画像の歪みを除去．
    plt.imshow(cv.cvtColor(dst, cv.COLOR_BGR2RGB)) # OpenCV は色がGBR順なのでRGB順に並べ替える
    plt.show()
    plt.close()

# 再投影誤差を計算．
mean_error = 0
for i in range(len(objpoints)) : # 頂点の数について繰り返し．
    imgpoints2, _ = cv.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv.norm(imgpoints[i],imgpoints2, cv.NORM_L2)/len(imgpoints2)
    mean_error += error

print("total error: ", mean_error/len(objpoints)) # 再投影誤差を出力．

なお，使用する関数については以下を参照．

推定された各種パラメータが以下．

歪み補正された画像の例が以下．

以上の処理により，カメラに関するパラメータを推定することができ，このカメラで撮影した画像の歪みを取り除くことができるようになった．

3. 2枚の画像の撮影

撮影時の必須条件は以下．

• ある程度の大きさの平面が画像内に写るようにする．
• カメラ方向はその平面の法線方向とする（平面を正面から写す）．この平面は後に2枚の画像を位置合わせする際の基準となる．
• 上記の平面に対し，移動経路が平行となるように（カメラが平面の方向を向けたまま）カメラを少しづつ移動させる．

撮影時の注意点は以下．

• できるだけカメラからの距離のバリエーションが多い方が分かりやすい結果が出る．
• できるだけカメラの高さを変えない．
• カメラを水平方向（左から右）に少しずつ動かしながら何枚か撮影．
• 手ブレ防止のためスマホスタンドやタイマー機能を使用するのがおすすめ．
• 2枚の画像，といいつつこの時点では10~20枚程度撮影しておくのが無難．

見た目ではかなりわずかな視差でもしっかり検知してくれるので，水平方向へのずらし幅は1cm程度で良い．
逆に大きいと視差画像の算出に失敗することが多い．

撮影時の様子が以下．
スマホスタンドを1cmずつ程度ずらしながら20枚ほど撮影した．
今回は位置合わせの基準と成る平面としてホワイトボードを写した．

これ以降は以下の2枚の画像を使用して説明する．

なお今回撮影した画像は以下のURLのフォルダ内からダウンロードしサンプル画像として使用可能。上記画像は13.pngおよび14.pngに対応する。
キャリブレーションデータも下記フォルダ内の20231211_calibration_imageに置いてある。このフォルダ内のチェッカーボード画像を使用し自分でパラメータを推測しても良いし、こちらで推定したパラメータもcamera_parameter.npzとして保存しているのでこれを使用しても良い。

4. 歪み補正

※これ以降の処理の説明では対応するコードを逐一示すが，全処理をまとめたコードも記事末尾に掲載する．

先述のキャリブレーションによって求めたカメラのパラメータ情報を用いて画像の歪みを除去する．
使用するコードは以下．

.py

import cv2 as cv
import numpy as np
import glob
from matplotlib import pyplot as plt

# 画像ファイル，およびカメラパラメータファイルの親ディレクトリのパスを指定．
folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc"
calb_folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc/20231211_calibration_image"

npz = np.load(f"{calb_folder_path}/camera_parameter.npz") # カメラパラメータの読込．
ret = npz["arr_0"]
mtx = npz["arr_1"] # カメラ行列．
dist = npz["arr_2"] # 歪み係数．
rvecs = npz["arr_3"] # 回転行列．
tvecs = npz["arr_4"] # 並進ベクトル．

def plt_output(material, title) : # 描画用の関数を定義．
    plt.figure(figsize = (9, 6), dpi = 600)
    plt.imshow(cv.cvtColor(material, cv.COLOR_BGR2RGB)) # OpenCV は色がGBR順なのでRGB順に並べ替える
    plt.title(title)
    plt.show()
    plt.close()

def undistortion(img, mtx, dist, h, w) : # 歪み除去の関数を定義．
    h, w = img.shape # 画像の縦，横の画素数を取得．
    # 最適なカメラ行列を取得．
    newcameramtx, roi = cv.getOptimalNewCameraMatrix(mtx, dist, # カメラ行列と歪み係数．
                                                     (w,h),
                                                     0, # この引数を 1 にすると元画像の全画素が維持される．
                                                     (w,h))
    dst = cv.undistort(img, mtx, dist, None, newcameramtx) # 画像の歪みを除去．
    plt_output(dst, "Undistorted Image")
    return dst

img1 = cv.imread(f"{folder_path}/13.png", cv.IMREAD_GRAYSCALE) # 左画像を白黒で読込．
img2 = cv.imread(f"{folder_path}/14.png", cv.IMREAD_GRAYSCALE) # 右画像を白黒で読込．
h, w = img1.shape

img1 = undistortion(img1, mtx, dist, h, w)
img2 = undistortion(img2, mtx, dist, h, w)

以下のように歪みが除去された画像が出力される．

5. 特徴点抽出

歪みを除去した左右画像から特徴点を抽出する．
今回は特徴量記述子としてSIFTとAKAZEを用いた．

一般的にSIFTは画像の拡大縮小・回転・輝度変化に頑健である．
またAKAZEは画像の拡大縮小・回転・輝度変化・Blurに頑健であり，処理速度も早い．
検出できる特徴点の数はSIFTの方が多い．

使用するコードが以下．

.py

# 特徴量記述子とそれに対応する距離関数を定義．
creater = cv.AKAZE_create() ; distance_func = cv.NORM_HAMMING # AKAZE を使用する場合．
# creater = cv.SIFT_create() ; distance_func = cv.NORM_L1 # SIFT を使用する場合．

kp1, des1 = creater.detectAndCompute(img1, None) # 特徴点を検出．
kp2, des2 = creater.detectAndCompute(img2, None)
img1_with_kp = cv.drawKeypoints(img1, kp1, None, flags = 4) # 特徴点を描画．
img2_with_kp = cv.drawKeypoints(img2, kp2, None, flags = 4)
plt_output(img1_with_kp, "Left Image with Key-Points")
print(f"The number of all key-points on the left image was {len(kp1)}")
plt_output(img2_with_kp, "Right Image with Key-Points")
print(f"The number of all key-points on the left image was {len(kp2)}")

AKAZEを使用した場合の結果が以下．

拡大図．

SIFTを使用した場合の結果が以下．

拡大図．

円の中心が特徴点，直線が特徴量の向き，半径が特徴量の大きさを示している．
SIFTの方が小さい円が多い．
今回は左右画像それぞれで，AKAZEでは約16000点ずつ，SIFTでは約30000点ずつの特徴点を抽出できた．

6. 対応点マッチング

抽出した左右画像それぞれの特徴点のうち，対応するものを対応点として検出する．

使用するコードが以下．

.py

matcher = cv.BFMatcher(distance_func, crossCheck = True) # マッチングアルゴリズムの定義．
if creater == cv.SIFT_create() : # 特徴量記述子が SIFT の場合．
    des1 = des1.astype(np.uint8) ; des2 = des2.astype(np.uint8) # 配列の型を変換．
matches = matcher.match(des1, des2) # 特徴点同士を総当りでマッチング．
matches = sorted(matches, key = lambda x:x.distance) # マッチングコスト順にソート．
src_pts = np.float32([ kp1[m.queryIdx].pt for m in matches ]).reshape(-1,1,2) # 特徴点からマッチングに成功した点を抽出
dst_pts = np.float32([ kp2[m.trainIdx].pt for m in matches ]).reshape(-1,1,2)
good_maches = matches[:round(len(matches) * 0.05)] # マッチングコストの高い上位 3 割を抽出．
img_with_matches = cv.drawMatches(img1, kp1, img2, kp2, good_maches, None, flags = 2)
plt_output(img_with_matches, "Two Images with Matched-Points")

使用した関数などについては以下を参照．

今回は上記コードにもある通り全対応点のうちマッチングコスト上位5%を描画した．
数値を変えれば更に多くの対応点を描画できるが，見やすさのために数を絞った．

AKAZEを使用した際の結果が以下．

SIFTを使用した際の結果が以下．

対応点が直線で結ばれている．
AKAZEを使用した場合は約12000組，SIFTを使用した場合は約14000組の対応点が検出された．

7. 位置合わせ

視差画像を算出するためには2枚の画像を適切に位置合わせする必要がある．
例えば，人間の目も，遠くの山を見るときには両目の画像内の山の像を重ね合わせて遠近感を把握する．
これと同じことを今回は，射影変換を用いて行う．
今回は，ある基準平面が含まれることを前提としているため，先ほど検出した対応点情報を用いて右画像内の基準平面を左画像内の基準平面に位置合わせするための射影変換行列を推定し，右画像に射影変換行列を作用させて左画像に位置合わせする．
射影変換を表現する数式が以下．

s
\begin{pmatrix}
x' \\
y' \\
1  \\
\end{pmatrix}
=
\begin{pmatrix}
a & b & c \\
d & e & f \\
g & h & 1 \\
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
1 \\
\end{pmatrix}
=
H\
\begin{pmatrix}
x \\
y \\
1 \\
\end{pmatrix}

ただし，$(x, y, 1)^T$：変換前の右画像，$(x', y', 1)^T$：変換後の右画像，$s$：スケール項，$H$：射影変換行列である．

また，OpenCVにおいて射影変換行列を推定するcv.findHomography()という関数ではオプションでRANSACの有無を指定できる．
RANSACとはデータの中から外れ値を検出・除去するためのアルゴリズムである．
今回の場合は，対応点の中から外れ値となる対応点を除去する．

射影変換は平面からある平面への変換といえる．言い方を変えれば，右画像のある平面上の像（点）を左画像上のある平面上に投影することはできるが，逆に右画像上で同一平面上にない点を左画像上のある平面上に投影することはできない．

そのため，位置合わせをするためのもっともらしい射影変換行列を計算するためには対応点の中から同一平面上にある点を選ぶ必要がある．

cv.findHomography()関数内のRANSACはある平面上にある対応点以外の対応点を外れ値として除外することで，もっともらしい射影変換行列の推定を可能にしている．
これを使用しない場合にはcv.findHomography()の3つ目の引数を0とすれば良い．

今回はある基準平面（今回であればホワイトボード）を含むような画像を前提としているため，この平面以外の対応点が除去される．

使用するコードが以下．
なお，今回は位置合わせの結果を分かりやすく描画するために2トーンの画像に変換したうえでオーバーラップ画像を出力した．

.py

# mask は extracted_matches と同じ要素数かつ要素が 0 or 1 の配列．
# 行列計算に使用された点は 1 で，使用されていない点は 0 で表現される．
H, mask = cv.findHomography(dst_pts, src_pts, cv.RANSAC, 1.0) # 射影変換行列の推定．
Re_img2 = cv.warpPerspective(img2, H, (w, h)) # 右画像の位置合わせ．
plt_output(Re_img2, "Rectified Right Image")

used_matches = [] # 行列算出に使用された対応点情報格納用の空配列を準備．
for i in range(mask.shape[0]) : # mask の要素について繰り返し．
    if mask[i][0] == 1 : # H の算出に使用した点である場合．
        used_matches.append(matches[i])
print(f"The number of used matched-points was {len(used_matches)}")
img_with_used_matches = cv.drawMatches(img1, kp1, img2, kp2, used_matches, None, flags = 2)
plt_output(img_with_used_matches, "Two Images with Used Matched-Points")

w_img1 = cv.applyColorMap(img1, cv.COLORMAP_JET)
w_Re_img2 = cv.applyColorMap(Re_img2, cv.COLORMAP_HSV)
overlap_img = cv.addWeighted(w_img1,0.5,w_Re_img2,0.5,0) # addWeighted(画像1,重み,画像2,重み,γ値) で，画像を重ね合わせる．
plt_output(overlap_img, "Overlapped Two Images")

位置合わせされた右画像が以下．
AKAZE，SIFTどちらを使用した場合でも位置合わせ結果はほとんど変わらなかった．
位置合わせによってやや右側にずれ，左端には空白が生じている．

位置合わせのための行列計算に使用された対応点が以下．

AKAZEを使用した場合．

SIFTを使用した場合．

今回はホワイトボードを含む平面上の対応点が行列計算に使用されたと考えられる．
そのため，左右画像のホワイトボードの像がぴったりと重なるように位置合わせされている．
それ以外の平面の対応点はRANSACで除去されている．
AKAZEを使用した場合には約2500点，SIFTを使用した場合には約3000点の対応点が使用された．

オーバーラップ画像が以下．
ホワイトボードを含む平面ではずれが全くないものの，それ以外の部分では若干のズレがある．これが視差（厳密にはホワイトボードの視差との差（視差差））である．

8. 視差画像算出

位置合わせした2枚の画像から視差画像を算出する．

視差画像を算出する前にまずはいちあわせした2枚の画像の重複領域のマスク画像を作成する．
これにより左右画像の非重複領域（片方の画像にしか写っておらず視差が計算不可である領域）を除去できる．

使用コードが以下．

.py

    ret, thresh1 = cv.threshold(img1,    1, 255, cv.THRESH_BINARY) # img1 の画像領域の形をしたマスクを作成．
    ret, thresh2 = cv.threshold(Re_img2, 1, 255, cv.THRESH_BINARY) # img2 の画像領域の形をしたマスクを作成．
    thresh1 = thresh1 * 255
    thresh2 = thresh2 * 255
    trimmed_img1 = img1 * thresh2 # 右画像の存在範囲領域の左画像を抽出．
    trimmed_img2 = Re_img2 * thresh1 # 左画像の存在範囲領域の右画像を抽出．

続いて視差画像を算出する．
視差画像を算出するアルゴリズムとしてはSemi-Grobal matchingというアルゴリズムを使用している．
この方法では，水平方向の画素マッチング結果に加えて，非現実的な視差の不連続を軽減するための仕組みが組み込まれている．そのため単純な画素ごとのマッチングのみを用いて算出される視差画像よりもなめらかで現実的な視差画像が得られやすい．
OpenCVではSGMの改良版であるSemi-Grobal Block Matchingがcv.StereoSGBM_create()関数で実装されているため，今回はこれを使用する．
ちなみに単純な画素マッチングのみを使用する関数はcv.StereoBM_create()で実装されており，下記コードでもオプションとして含めている．

また，またcv.StereoSGBM_create()関数では視差の計算に失敗した画素にminDisparityで指定した値を代入する仕組みになっているが，この視差値には本質的に意味が無いのでNaN(Not a Number)に置換している．
描画する際にはNaNは空白（白色）で表現される．

使用コードは以下．

.py

# SGBM 関数のパラメータを定義．
minDisparity  = -50 ; numDisparities = 100 ; blockSize       = 11 ; cost1 =   1 ; cost2        = 4
disp12MaxDiff =   0 ; preFilterCap   =   0 ; uniquenessRatio =  0 ; sWS   = 600 ; speckleRange = 2
P1 = cost1 * 3 * ( blockSize ** 2) ; P2 = cost2 * 3 * ( blockSize ** 2)
stereo = cv.StereoSGBM_create(minDisparity, numDisparities, blockSize, P1, P2,
                disp12MaxDiff, preFilterCap, uniquenessRatio, sWS, speckleRange, mode = cv.STEREO_SGBM_MODE_SGBM_3WAY)
# stereo = cv.StereoBM_create(numDisparities = 96, blockSize = 15)
### DI 算出．
DI = stereo.compute(trimmed_img1, trimmed_img2).astype(np.float32) / 16.0 # DI を算出．
DI = np.where((thresh1 == 1) & (thresh2 == 1), DI, np.nan) # 左右画像の重複領域以外の領域の視差を NaN に置換．
DI[np.where(DI <= minDisparity)] = np.nan # 無意味な最低視差値を NaN に置換．

plt.figure(figsize = (9, 6), dpi = 300)
# plt.rcParams["font.family"] = "Times New Roman"
plt.imshow(DI, cmap = plt.cm.jet, vmin = minDisparity, vmax = minDisparity + numDisparities)
plt.tick_params(labelsize = 20)
cbar = plt.colorbar(aspect = 40, pad = 0.02, shrink = 0.7, orientation = "vertical", extend = "both")
cbar.ax.tick_params(axis = "y", labelsize = 20)
cbar.set_label("Disparity", labelpad = 15, size = 30)
plt.xlabel("Width", fontsize = 30)
plt.ylabel("Height", fontsize = 30)
plt.show()
plt.close()

使用した関数については以下を参照．

AKAZEを使用した場合の視差画像が以下．

SIFTを使用した場合の視差画像が以下．

どちらの視差画像を通しても，カメラに近い物体の視差は大きく，カメラから遠い物体の視差は小さくなっている．
今回はホワイトボードを基準として位置合わせしたため，ホワイトボードより奥の物体の視差は負に，ホワイトボードより手前の物体の視差は正となっている．
視差を計算するための画素同士のマッチングは輝度の大小関係を元に行われるため，白い壁や黒い布部分など，テクスチャが乏しい部分ではマッチングがうまくいかず視差が計算されていない（空白になっている）部分が多い．

また今回は特徴量記述子による結果の違いはなかったが，用いる画像によってはかなり結果が異なる場合もある．

うまくいかない場合の対処法

• チェッカーボードをしっかりとした平面に貼り付ける．
• できるだけチェッカーボードが大きく写るように撮影してキャリブレーションする．
• 写す平面の大きさを大きくしてみる．
• cv.StereoSGBM_create()のパラメータの値を調整する．特にminDisparityやnumDisparitesの値を調整すると，視差が計算されなかった画素でも計算されるようになることがある．
• AKAZE，SIFT以外にもORBなど他の特徴量記述子もあるので試してみる．
• 画像の撮影時には，左右のズレ幅を小さくする．今回は1cmとしたが，これで充分．
• 他の画像で試してみる．

付録（上記をまとめたコード）

歪み補正から視差画像算出までをまとめたコードが以下．
実際には複数組の画像ペアから複数枚の視差画像を算出して結果を比較したいケースが多いので，画像番号でfor文を回せるようにしている．for文の中身は上記の処理そのまま．

DI_calc.py

DI_calc.py

import cv2 as cv
import numpy as np
import glob
from matplotlib import pyplot as plt

# 画像ファイル，およびカメラパラメータファイルの親ディレクトリのパスを指定．
folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc"
calb_folder_path = "/content/drive/MyDrive/研究資料/20231211_disparity_calc/20231211_calibration_image"

# 特徴量記述子とそれに対応する距離関数を定義．
creater = cv.AKAZE_create() ; distance_func = cv.NORM_HAMMING # AKAZE を使用する場合．
# creater = cv.SIFT_create() ; distance_func = cv.NORM_L1 # SIFT を使用する場合．
matcher = cv.BFMatcher(distance_func, crossCheck = True) # マッチングアルゴリズムの定義．

npz = np.load(f"{calb_folder_path}/camera_parameter.npz") # カメラパラメータの読込．
ret = npz["arr_0"]
mtx = npz["arr_1"] # カメラ行列．
dist = npz["arr_2"] # 歪み係数．
rvecs = npz["arr_3"] # 回転行列．
tvecs = npz["arr_4"] # 並進ベクトル．

# SGBM 関数のパラメータを定義．
minDisparity  = -50 ; numDisparities = 100 ; blockSize       = 11 ; cost1 =   1 ; cost2        = 4
disp12MaxDiff =   0 ; preFilterCap   =   0 ; uniquenessRatio =  0 ; sWS   = 600 ; speckleRange = 2
P1 = cost1 * 3 * ( blockSize ** 2) ; P2 = cost2 * 3 * ( blockSize ** 2)
stereo = cv.StereoSGBM_create(minDisparity, numDisparities, blockSize, P1, P2,
                disp12MaxDiff, preFilterCap, uniquenessRatio, sWS, speckleRange, mode = cv.STEREO_SGBM_MODE_SGBM_3WAY)
# SGBM ではなく BM を使用する際には以下．
# stereo = cv.StereoBM_create(numDisparities = 96, blockSize = 15)

def plt_output(material, title) : # 描画用の関数を定義．
    plt.figure(figsize = (9, 6), dpi = 300) # 図のサイズと解像度を定義．
    plt.imshow(cv.cvtColor(material, cv.COLOR_BGR2RGB)) # OpenCV は色がGBR順なのでRGB順に並べ替える
    plt.title(title) # タイトルを描画．
    plt.show()
    plt.close()

def undistortion(img, mtx, dist, h, w) : # 歪み除去の関数を定義．
    h, w = img.shape # 画像の縦，横の画素数を取得．
    # 最適なカメラ行列を取得．
    newcameramtx, roi = cv.getOptimalNewCameraMatrix(mtx, dist, # カメラ行列と歪み係数．
                                                     (w,h),
                                                     0, # この引数を 1 にすると元画像の全画素が維持される．
                                                     (w,h))
    dst = cv.undistort(img, mtx, dist, None, newcameramtx) # 画像の歪みを除去．
    plt_output(dst, "Undistorted Image") # 歪みを除去した画像を描画．
    return dst

for img_i in [13,14,15] : # 左画像として使用したい画像番号ごとに繰り返し．
    ### 画像の読込．
    img1 = cv.imread(f"{folder_path}/{img_i}.png", cv.IMREAD_GRAYSCALE) # 左画像を白黒で読込．
    img2 = cv.imread(f"{folder_path}/{img_i + 1}.png", cv.IMREAD_GRAYSCALE) # 右画像を白黒で読込．
    h, w = img1.shape # 画像の高さ，幅を取得．

    ### 歪み除去．
    img1 = undistortion(img1, mtx, dist, h, w) # 画像の歪みを除去．
    img2 = undistortion(img2, mtx, dist, h, w)

    ### 特徴点抽出．
    kp1, des1 = creater.detectAndCompute(img1, None) # 特徴点を検出．
    kp2, des2 = creater.detectAndCompute(img2, None)
    img1_with_kp = cv.drawKeypoints(img1, kp1, None, flags = 4) # 特徴点を描画．
    img2_with_kp = cv.drawKeypoints(img2, kp2, None, flags = 4)
    plt_output(img1_with_kp, "Left Image with Key-Points")
    plt_output(img2_with_kp, "Right Image with Key-Points")
    print(f"The number of all key-points on the left image was {len(kp1)}")
    print(f"The number of all key-points on the left image was {len(kp2)}")

    ### 特徴点マッチング．
    if creater == cv.SIFT_create() : # 特徴量記述子が SIFT の場合．
        des1 = des1.astype(np.uint8) ; des2 = des2.astype(np.uint8) # 配列の型を変換．
    matches = matcher.match(des1, des2) # 特徴点同士を総当りマッチング．
    matches = sorted(matches, key = lambda x:x.distance) # マッチングコスト順にソート．
    print(f"The number of all matched-points was {len(matches)}")
    src_pts = np.float32([ kp1[m.queryIdx].pt for m in matches ]).reshape(-1,1,2) # 特徴点からマッチングに成功した点を抽出
    dst_pts = np.float32([ kp2[m.trainIdx].pt for m in matches ]).reshape(-1,1,2)
    good_maches = matches[:round(len(matches) * 0.05)] # マッチングコストの高い上位 5% を描画用に抽出．
    img_with_matches = cv.drawMatches(img1, kp1, img2, kp2, good_maches, None, flags = 2)
    plt_output(img_with_matches, "Two Images with Matched-Points")

    #### 位置合わせ．
    # mask は extracted_matches と同じ要素数かつ要素が 0 or 1 の配列．
    # 行列計算に使用された点は 1 で，使用されていない点は 0 で表現される．
    H, mask = cv.findHomography(dst_pts, src_pts, cv.RANSAC, 1.0) # 射影変換行列の推定．
    Re_img2 = cv.warpPerspective(img2, H, (w, h)) # 右画像の位置合わせ．
    plt_output(Re_img2, "Rectified Right Image")

    used_matches = [] # 行列算出に使用された対応点情報格納用の空配列を準備．
    for i in range(mask.shape[0]) : # mask の要素について繰り返し．
        if mask[i][0] == 1 : # H の算出に使用した点である場合．
            used_matches.append(matches[i]) # used_matches に matches の i 番目の要素を追加．
    print(f"The number of used matched-points was {len(used_matches)}")
    img_with_used_matches = cv.drawMatches(img1, kp1, img2, kp2, used_matches, None, flags = 2)
    plt_output(img_with_used_matches, "Two Images with Used Matched-Points")

    w_img1 = cv.applyColorMap(img1, cv.COLORMAP_JET) # カラーマップを調整した左画像を作成．
    w_Re_img2 = cv.applyColorMap(Re_img2, cv.COLORMAP_HSV) # カラーマップを調整した右画像を作成．
    overlap_img = cv.addWeighted(w_img1,0.5,w_Re_img2,0.5,0) # 画像を重ね合わせ．
    plt_output(overlap_img, "Overlapped Rectified Two Images")

    ### マスク画像作成．
    ret, thresh1 = cv.threshold(img1,    1, 255, cv.THRESH_BINARY) # img1 の画像領域の形をしたマスクを作成．
    ret, thresh2 = cv.threshold(Re_img2, 1, 255, cv.THRESH_BINARY) # img2 の画像領域の形をしたマスクを作成．
    thresh1 = thresh1 * 255
    thresh2 = thresh2 * 255
    trimmed_img1 = img1 * thresh2 # 右画像の存在範囲領域の左画像を抽出．
    trimmed_img2 = Re_img2 * thresh1 # 左画像の存在範囲領域の右画像を抽出．

    ### DI 算出．
    DI = stereo.compute(trimmed_img1, trimmed_img2).astype(np.float32) / 16.0 # DI を算出．
    DI = np.where((thresh1 == 1) & (thresh2 == 1), DI, np.nan) # 左右画像の重複領域以外の領域の視差を NaN に置換．
    DI[np.where(DI <= minDisparity)] = np.nan # 無意味な最低視差値を NaN に置換．

    plt.figure(figsize = (9, 6), dpi = 300)
    plt.imshow(DI, cmap = plt.cm.jet, vmin = minDisparity, vmax = minDisparity + numDisparities)
    plt.tick_params(labelsize = 20)
    cbar = plt.colorbar(aspect = 40, pad = 0.02, shrink = 0.7, orientation = "vertical", extend = "both")
    cbar.ax.tick_params(axis = "y", labelsize = 20)
    cbar.set_label("Disparity", labelpad = 15, size = 30)
    plt.xlabel("Column", fontsize = 30)
    plt.ylabel("Row", fontsize = 30)
    plt.show()
    plt.imsave(f"{folder_path}/{img_i}_{img_i + 1}_DI.png", DI)
    plt.close()