カメラの較正（キャリブレーション）は意外と難しい．理由は，画像の撮影法で結果が毎回変わってしまうくらい不安定だから．既存のマーカのコーナー検出ライブラリも不安定で使いづらい．OpenCVのマニュアルだけでは分からないことがある．せめて分かっていることだけでもメモして残しておこうかと．

長いので目次つくりました．

カメラ較正のプログラムコード
較正のコツ
歪補正画像の作り方
- 方法1：getOptimalNewCameraMatrix関数を用いる
- 方法2：画像の隅に対応する半径を求める

カメラ較正のプログラムコード

import numpy as np
import cv2
import glob #フォルダ内の複数のファイル入力用
import matplotlib.pyplot as plt

マーカ関係の変数の準備

#マーカのコーナー数（内側の格子のみカウントする）
size = (10,7)

# コーナーの三次元座標の二次元配列：(0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
objp = np.zeros((size[0]*size[1],3), np.float32)
objp[:,:2] = np.mgrid[0:size[0],0:size[1]].T.reshape(-1,2)

# コーナーの三次元座標と画像上の座標の一次元配列
objpoints = [] # 三次元座標
imgpoints = [] # 画像座標

# 画像ファイル名リスト
dir = "org"
images = glob.glob(dir +"/l*.jpg")

コーナ検出

for fname in images:
    #ファイルの飲み込み
    img = cv2.imread(fname)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    # コーナー検出
    ret, corners = cv2.findChessboardCorners(img, (10,7),None)

    # 検出出来なければスルーする
    if ret == False:
        print("could not find corners in " + fname)
        continue
        
    # サブピクセル・コーナー検出
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
    corners2 = cv2.cornerSubPix(gray, corners,(3,3),(-1,-1), criteria)

    # 全コーナーの三次元位置を追加
    imgpoints.append(corners2)
    objpoints.append(objp)
    
    draw = img.copy()
    draw = cv2.drawChessboardCorners(draw, (10,7), corners,ret)
    dfname = fname.replace(dir, 'draw')
    cv2.imwrite(dfname, draw)

#　最後の画像だけ拡大
plt.figure(figsize=(10,10))
plt.imshow(draw)
plt.axis('off')

内部パラメータ推定

flags = (cv2.CALIB_FIX_ASPECT_RATIO |\
         cv2.CALIB_FIX_K4 | cv2.CALIB_FIX_K5 | cv2.CALIB_FIX_K6 |\
         cv2.CALIB_ZERO_TANGENT_DIST)
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

#レンズ歪補正画像生成してundistフォルダに保存
for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    dst = cv2.undistort(img, mtx, dist)
    dfname = fname.replace(dir, 'undist/')
    cv2.imwrite(dfname, dst)

#最後の画像を表示
img = cv2.imread("img/l-0035.jpg")
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
dst = cv2.undistort(img, mtx, dist)

plt.figure(figsize=(14,14))
plt.subplot(1, 2, 1)
plt.imshow(img)
plt.subplot(1, 2, 2)
plt.imshow(dst)
plt.axis('off')

CSVファイルに保存

内部パラメータ行列

np.savetxt('intrinsic_matrix.csv', mtx, fmt='%.5e', delimiter=',')
print(mtx)

[[1.22109326e+03 0.00000000e+00 8.15761285e+02]
 [0.00000000e+00 1.21712653e+03 6.31027660e+02]
 [0.00000000e+00 0.00000000e+00 1.00000000e+00]]

レンズ歪み係数

np.savetxt('intrinsic_dist_left.csv', dist, fmt='%.5e', delimiter=',')
print(dist)

[[-0.36750751  0.20992473  0.00206952 -0.00136743 -0.08158695]]

較正のコツ

カメラの較正は，理論やプログラミでは表現しずらい色々なコツがあるのでまとめておく． OpenCVのようなライブラリが無料で利用できる時代にあって，未だにカメラの較正は，昔と変わらず面倒で教科書を読んで理論とプログラムを理解しただけでは，コツを知ることはできない．

マーカは画角内いろんな位置・姿勢で

特に広角のレンズを用いる場合はレンズ歪みが大きい場合が多く，画像の隅々までマーカを移すことが重要です．カメラがその光学中心で回転しただけで得られるような画像だとだめ．それを防ぐには，カメラにマーカが正対するよう配置するだけでなくマーカを傾けて色んな姿勢で撮影しなければならない．

plt.figure(figsize=(14,7))
row = 6
col = 8
num = 0

for fname in images:
    #ファイルの飲み込み
    img = cv2.imread(fname)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
    
    # 表示（不要の場合はこれ以降削除）
    num += 1
    plt.subplot(row, col, num)
    plt.axis('off')
    plt.imshow(img)

え？findChessboardCorners関数は内部で自動的に画像を二値化するんじゃないの？と思うかもしれないが，これが中々うまく行かない．人の目で見る限りかなりきれいに撮影された画像でも検出に失敗することがある．画角が広いと周辺減光が生じ周辺が暗くなるのでグラデーションが生じて画像全体の2値化だとうまく行かないことがある．ここではGoogle Photoでコントラスト調整するだけで検出がうまく行った例を示す．

plt.figure(figsize=(14,14))

#元の画像
img= cv2.imread("org/r-0017.jpg")
ret, corners = cv2.findChessboardCorners(img, (10,7),None)
print("original image: ", ret)
draw = cv2.drawChessboardCorners(img.copy(), (10,7), corners, ret)
plt.subplot(1, 2, 1)
plt.imshow(draw)
plt.axis('off')

#Google Photoでコントラストを調整した画像
img = cv2.imread("right/r-0017.jpg")
ret, corners = cv2.findChessboardCorners(img, (10,7),None)
print("adjusted image: ", ret)
draw = cv2.drawChessboardCorners(img.copy(), (10,7), corners,ret)
plt.subplot(1, 2, 2)
plt.imshow(draw)
plt.axis('off')

original image:  False
adjusted image:  True

この例では，img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)がコメントアウトされており，BGR画像が表示されている．このコメントを外すとコーナーは検出されるようになる．RとBを入れ替えるだけで結果が変わるくらい入力画像に敏感だということ．

レンズ歪みモデルを使い分ける

OpenCVのレンズ歪みモデルは，一部のパラメータを0に固定することで，パラメータ数が少ないシンプルなモデルになるようなっている．パラメーター数が多いとオーバーフィッティングになる可能性があるので，パラメータ数が少ないモデルから試して上手く行かなければ，徐々にパラメータが多いモデルに変えるのがいい．

calibrateCamera()関数のCALIB_FIX_K2などCALIB_FIX_XXXというオプションを用いればパラメータを固定できる．特に広角でもない普通のカメラなら歪みモデルは 𝑘1 のみのシンプルなモデルで十分な場合が多い．

上手く出来たかどうかの確認方法

正しく歪み補正されていれば，三次元空間中で直線のものは画像上でも直線になるはず．そこで歪み補正前後でコーナーに直線をひいて確認する．この確認で使う画像は，キャリブレーションに使っていない画像でないと意味がない．

#最後の画像を表示
img = cv2.imread("img/l-0035.jpg")
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
dst = cv2.undistort(img, mtx, dist)

#コーナー検出
ret, corners = cv2.findChessboardCorners(gray, (10,7),None)
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
corners2 = cv2.cornerSubPix(gray,corners,(3,3),(-1,-1), criteria)

#コーナーの歪補正
corners2 = corners2[:,0,:]
#undstpoints = corners2.copy()
undstpoints = cv2.undistortPoints(corners2, mtx, dist, P=mtx)
undstpoints = undstpoints[:,0,:]

#直線の描画
img_draw = img.copy() #
dst_draw = dst.copy() #
cv2.line(img_draw, (corners2[0,0], corners2[0,1]), (corners2[9,0],corners2[9,1]), (255,0,0), 5)
cv2.line(img_draw, (corners2[0,0], corners2[0,1]), (corners2[10*6,0],corners2[10*6,1]), (255,0,0), 5)
cv2.line(dst_draw, (undstpoints[0,0], undstpoints[0,1]), (undstpoints[9,0], undstpoints[9,1]), (255,0,0), 5)
cv2.line(dst_draw, (undstpoints[0,0], undstpoints[0,1]), (undstpoints[10*6,0], undstpoints[10*6,1]), (255,0,0), 5)

#表示
plt.figure(figsize=(14,14))
plt.subplot(1, 2, 1)
plt.imshow(img_draw)
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(dst_draw)
plt.axis('off')

他に方法ないの？

あります．でもインストールの手間まで考えると現状上の方法が一番楽かなと思う．マーカをロバストに検出する手法はいくつも研究されているけど一度しかしなくていいキャリブレーションのために色々環境を構築するのは面倒なので結局OpenCVを使ってしまうことが多い．

歪補正画像の作り方

レンズ歪みが補正された画像を逆変換（補正後の画像の座標から補正前の座標へ変換）で生成する場合，下の画像の左上や左下ように元の画像の中央部分が繰り返し現れることがある．この部分を現れないように補正画像を作る方法を２つメモしておく．

mtx_new = mtx.copy()
mtx_new[0,0] = mtx_new[0,0]*0.5
mtx_new[1,1] = mtx_new[1,1]*0.5
dst = cv2.undistort(gray, mtx, dist, newCameraMatrix=mtx_new);
plt.figure(figsize=(14,14))
plt.axis('off')
plt.imshow(dst, cmap="gray")

方法1：getOptimalNewCameraMatrix関数を用いる

mtx_new, roi = cv2.getOptimalNewCameraMatrix(mtx, dist, (img.shape[0],img.shape[1]), 1.0)

undst = cv2.undistort(gray, mtx, dist, newCameraMatrix=mtx_new);
plt.figure(figsize=(14,14))
plt.imshow(undst, cmap="gray")

確かに無駄な部分は消えましたがなんだか結果はいまいち．左上のびよ～～んと歪んだ部分が残っている．あと，こうやってみるとマーカをプリントした紙の辺が歪んていることが分かる．画像の隅々までは歪み補正できていないということ．OpenCVのコーナー検出がマーカ全体が写っていないと検出されない仕組みなのでこうなりやすい．

方法2：画像の隅に対応する半径を求める

getOptimalNewCameraMatrixでなんとなく無駄な部分を削除することができたけど，まだ余計な部分が残っている．ここではシェーダで描画する場合などのように歪み無し画像の座標でどこまでを描くかを主点からの半径rで指定する方法を紹介する．ここでいうrとはOpenCVのマニュアル「Camera Calibration and 3D Reconstruction」の$r$のこと．中央から半径$r$方向外側に見ていった時，対応する歪有りの座標が(0,0)に最も近づいたときに，打ち切れるようにそのような半径$r$を求める．

import scipy.optimize as opt
cmat = (mtx[0,0], mtx[1,1], mtx[0,2], mtx[1,2])
(fx, fy, cx, cy) = cmat
(k1, k2, p1, p2, k3) = tuple(dist.ravel())

まず歪み無し(xp, yp)から歪みあり(xpp, ypp)への変換関数を定義する．これは後のコードを短くするため．

def dist_func(xp, yp, cmat, dist):
    r2 = xp*xp + yp*yp
    r4 = r2 * r2
    r6 = r4 * r2
    (fx, fy, cx, cy) = cmat
    (k1, k2, p1, p2, k3) = tuple(dist.ravel())
    xpp = xp*(1 + k1*r2 + k2*r4 + k3*r6) + 2*p1*xp*yp + p2*(r2 + 2*xp*xp)
    ypp = yp*(1 + k1*r2 + k2*r4 + k3*r6) + p1*(r2 + 2*yp*yp) + 2*p2*xp*yp
    return (xpp, ypp)

非線形最適化で歪みあり座標が(xpp, ypp)=(x0,y0)となる歪み無し座標(xp, yp)を求めるための評価関数を定義する．

def err(r, target, shape, cmat, dist):
    (fx, fy, cx, cy) = cmat
    w = shape[1]
    h = shape[0]
    x0 = target[1]
    y0 = target[0]
    xp = - r/np.sqrt(w**2 + h**2) * w
    yp = - r/np.sqrt(w**2 + h**2) * h
    (xpp, ypp) = dist_func(xp, yp, cmat, dist)
    return (fx*xpp + cx - x0)*(fx*xpp + cx - x0) + (fy*ypp + cy - y0)*(fy*ypp + cy - y0)

評価関数がどんな関数なのか確認するために描いてみる．

target = (0, 0)
shape = img.shape
r = np.linspace(0, 2, 100)
e = err(r, target, shape, cmat, dist)

plt.plot(r,e)

上の評価関数を最小化することで(x0,y0)=(0,0)となる半径rを求める．

r_opt = opt.fmin(err, 0.0, args = (target, shape, cmat, dist),disp=False)
print(r_opt)

[1.17065625]

これでr2の範囲を「r2< (1.17065625)**2」のように指定すれば余計な箇所が描画されずに済む

udimg = np.zeros((2000,2000),np.uint8)

for j in range(udimg.shape[0]):
    for i in range(udimg.shape[1]):
        ip = i - int((udimg.shape[1]-gray.shape[1])/2)
        jp = j - int((udimg.shape[0]-gray.shape[0])/2)
        xp = (ip - cx)/fx
        yp = (jp - cy)/fy
        r = xp*xp + yp*yp
        if r_opt < r:
            continue
        (xpp, ypp) = dist_func(xp, yp, cmat, dist)
        ipp = int(fx*xpp + cx)
        jpp = int(fy*ypp + cy)
        if ipp < 0 or 1200 <= ipp or jpp < 0 or 1200 <= jpp:
            continue
        udimg[j,i] = gray[jpp,ipp]

plt.axis('off')
plt.imshow(udimg, cmap="gray")

とこのように，無駄な部分を半径だけで判定して描画しないようにすることができた．今回は左上としてtarget=(0,0)を基準にしたけれど，画像の4隅ならどこでもいい．

武相荘

2021年1月3日日曜日

カメラ較正のテクニック色々メモ