動画のための手振れ補正２

　複数フレームの画像を元に手振れ補正のPythonコードをOpenCV使って書いてみる続きです。

　前回は前フレームにピッタリ位置合わせするパターンでコードを書いてみましたが、今回は現フレームの画像をある程度尊重し、シーンチェンジに堅牢な手振れ補正にしてみようと思います。

フィードバック率

　いろいろな意味がありそうですが、今回ここでは現フレームの角度、中心位置をどの程度変化させるか。という意味で使います。

　つまり前フレームにピッタリ位置を合わせると10度回す必要があるが、フィードバック率を80%とするならば8度の回転にとどめる。とかそんな意味に使います。中心位置が10画素ずれるのであれば同様に8画素にとどめる。とかそんな使い方です。特定のフレームの歪みに引きずられるのを回避するのが目的です。

　例えば、このようにスタートフレームをわざと斜めにしておくと、この動画をソースに手振れ補正のコードを書くと最後まで斜めになりっぱなしです。

　このような不良フレームがシーンの冒頭にあると

　このように最後まで引きずられてしまいます。今回はこれをある程度抑制させる仕組みをコーディングしてみます。

アルゴリズム

　そんなに難しく考えません。位置合わせで得られたアフィン行列から、移動、回転量が抑制された新たな行列を作ります。流れは

• 拡大率、回転角を求める
• 中心座標の移動先を求める
• 回転、拡大が抑制された新たに行列を作る
• 新たな行列での中心座標移動先を求める
• 抑制された中心位置を求める
• 行列のシフト量を変更する

　こんな感じです。

拡大率、回転角を求める

　入力の行列から拡大率、回転角度、を求めます。これは以前作った

Affine行列から回転角度を取り出す。

　アフィン変換の変換行列からその回転成分（回転角度、回転中心座標）だけ取り出してみます。ちょっとだけAffine変換の動きを復習したうえで試します。OpenCVとPythonを使います。Affine変換　アフィン変換そのも...

campkougaku.com

2020.07.14

　このアルゴリズムを使います。これにより拡大率、回転角度を求めます。詳細は上記に譲りますが、連立方程式を解いているだけです。

中心座標の移動先を求める

　これも簡単で、アフィン変換して座標を求めているだけです。入力された行列がmtxfitだとして、

# 中心座標の移動先を計算
mx = mtxfit[0, 0] * cx + mtxfit[0, 1] * cy + mtxfit[0, 2]
my = mtxfit[1, 0] * cx + mtxfit[1, 1] * cy + mtxfit[1, 2]

　cx,cyがオリジナルの中心。mx,myが中心の移動先。行列演算を手書きしてるだけです。

　こんなイメージ

回転、拡大が抑制された新たに行列を作る

　得られた回転角度の逆向きの回転方向で、フィードバック率を乗じた角度で回す行列を作ります。例えばフィードバック率が0.8で、入力が10度回す行列であれば、-2度回する行列を作るイメージです。

　拡大率も逆数を取ります。これは必ず倍率を1倍にするための対応で、90%縮小するのであれば、10/9倍拡大することになります。

# 倍率、回転量の抑制 (逆回し)
scale = 1/scale
angle = (feedback-1) * angle

　このようにして求めた逆回転、逆数倍拡大する行列（mtxback）を新たに作ります。cv2.getRotationMatrix2Dを使います。ついでに行列の積を取るために3×3に拡張しています。

# 抑制された行列 mtx の生成 (3x3)
mtxback = cv2.getRotationMatrix2D((mx, my), angle, scale)
mtxback = np.concatenate((mtxback, np.array([[0.0, 0.0, 1.0]])))

　そのうえでその行列を入力行列に対して左から乗じます。

mtx = np.dot(mtxback, mtxfit)

　これにより少しだけ元に戻った、拡大率1倍の新たな行列が得られます。

新たな行列での中心座標移動先を求める

　これも新たな行列でアフィン変換するだけ。中心位置はここで求めた座標に移動します。

# 中心座標の移動先を計算
mx = mtx[0, 0] * cx + mtx[0, 1] * cy + mtx[0, 2]
my = mtx[1, 0] * cx + mtx[1, 1] * cy + mtx[1, 2]

行列のシフト量を変更する

　新たな中心移動先を元に、フィードバック率を使って元の中心位置に近づけます。アフィン行列のシフト成分を逆算します。行列の

$$ Sh=\left( \begin{matrix} 1 & 0 & Sx \\ 0 & 1 & Sy \\ 0 & 0 & 1 \end{matrix} \right) $$

　このSx, Sy の部分を再計算してあげると、シフト成分だけ補正されます。

# 移動量を抑制 (画像中心に寄せる)
mtx[0, 2] = (mx * feedback + cx * (1-feedback)) - mtx[0, 0] * cx - mtx[0, 1] * cy
mtx[1, 2] = (my * feedback + cy * (1-feedback)) - mtx[1, 0] * cx - mtx[1, 1] * cy

　この部分が新たな中心座標です。線形補間しているだけです。

mx * feedback + cx * (1-feedback)

　このようにして新たな行列を作ってあげると、フィードバック率に応じて現フレームの向き、位置が尊重されます。新たな変換後の位置はある程度変換前の位置が保持されます。以下のようなイメージ。

　リミッタをここで儲けてもよいかもしれません。

結果

　このアルゴリズムで行くと、先ほどとは違い、徐々にフレームをまたぐにつれて向きが水平に戻っていきます。

　不良フレームが途中にある場合には悪い方向に働く場合もあると思います。シーン毎に最適なフィードバック率があると思います。

ちょっとしたコツ

　前回borderModeとして、前フレームの画像を後ろに透かして、はみ出した黒い領域を少しでもごまかそうとしました。ただ、今回のようにフィードバック率によっては前フレームと位置がずれて、はみ出した黒領域もギャップが出てしまいます。

　なので、前フレームの画像を現フレームの画像とピッタリ位置合わせするような変換行列を計算し、前フレームを変換し、

　これも行列の積を取っていくだけです。変換の順番としては、

・前フレームに対して、前回使った行列の逆行列を使って、元の向きにする
・位置合わせした行列の逆行列を使って、前フレームを変換前の現フレームとピッタリ合わせる
・今回作った行列で現フレームと同じ変換を行う。

　この順で行列を乗じることで、変換後の現フレームとピッタリ位置が合った前フレームの画像が作れます。

　画像そのものの変換は1回だけです。行列を3回乗じて3回分の変換を合成できます。

　これにより隅っこもそれっぽく補正されます。

最終結果

　最初のフレームに引きずられることなく、手振れ補正できています。四隅もまぁまぁ自然です。今回のようなシーンでは使えます。フィードバック率は適当に0.8としました。

まとめ

　2回かけて動画手振れ補正のコードを書いてみました。まぁ最近のカメラであれば内部でこれに近いことやってくれるので自分で書く必要はないですけどね。加えてジャイロセンサも使えば無敵でしょうし。今回はお勉強です。

　正直4辺は削った方がいいですね。カメラを動画モードにすると途端に画角が狭くなりますが、その程度は周りを削った方がよさそうです。

　シーン判定してこのフィードバック率を自動で求められるとさらに面白そうです。アルゴリズムが思いつきませんけど。

全コード

　上記のフィードバックの制御に加え、特徴量を取る範囲を中心位置に狭め処理速度を上げ、保存時に4辺を削ったコードにしています。

import cv2
import numpy as np
import math
import glob

##
# @brief AKAZEによる画像特徴量取得
# @param img 特徴量を取得したい画像（RGB順想定）
# @param pt1 特徴量を求める開始座標 tuple (default 原点)
# @param pt2 特徴量を求める終了座標 tuple (default None=画像の終わり位置)
# @return key points
def get_keypoints(img, pt1=(0, 0), pt2=None):
    if pt2 is None:
        pt2 = (img.shape[1], img.shape[0])

    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
    mask = cv2.rectangle(np.zeros_like(gray), pt1, pt2, color=1, thickness=-1)
    sift = cv2.AKAZE_create()

    # find the key points and descriptors with AKAZE
    return sift.detectAndCompute(gray, mask=mask)


##
# @brief 特徴記述子kp2/des2にマッチするような pointを求める
# @param kp1 合わせたい画像のkeypoint
# @param des1 合わせたい画像の特徴記述
# @param kp2 ベースとなる画像のkeypoint
# @param des2 ベースとなる画像の特徴記述
# @return apt1 kp1の座標　apt2 それに対応するkp2
def get_matcher(kp1, des1, kp2, des2):
    if len(kp1) == 0 or len(kp2) == 0:
        return None

    # Brute-Force Matcher生成
    bf = cv2.BFMatcher()
    matches = bf.knnMatch(des1, des2, k=2)

    # store all the good matches as per Lowe's ratio test.
    good = []
    for m, n in matches:
        if m.distance < 0.7 * n.distance:
            good.append(m)

    if len(good) == 0:
        return None

    target_position = []
    base_position = []
    # x,y座標の取得
    for g in good:
        target_position.append([kp1[g.queryIdx].pt[0], kp1[g.queryIdx].pt[1]])
        base_position.append([kp2[g.trainIdx].pt[0], kp2[g.trainIdx].pt[1]])

    apt1 = np.array(target_position)
    apt2 = np.array(base_position)
    return apt1, apt2


##
# @brief Affine行列から回転角、回転中心、拡大率を求める
# @param mat アフィン行列
# @return da 回転角度(度) center　回転中心　s 拡大率
def getRotateShift(mat):
    da = -math.atan2(mat[1, 0], mat[0, 0])  # ラジアン
    s = mat[0, 0] / math.cos(da)  # 拡大率

    m = np.zeros([2, 2])
    m[0, 0] = 1 - mat[0, 0]
    m[0, 1] = -mat[0, 1]
    m[1, 0] = mat[0, 1]
    m[1, 1] = m[0, 0]

    mm = np.zeros([2, 1])
    mm[0, 0] = mat[0, 2]
    mm[1, 0] = mat[1, 2]

    center = np.dot(np.linalg.inv(m), mm).reshape([2])

    return math.degrees(da), center, s


##
# @brief アフィン行列mtxfitの変化量を抑制する
# @param matfit アフィン行列 (3x3)
# @param cx  中心のx座標
# @param cy  中心のy座標
# @param feedback  フィードバック量（1.0で全フィードバック）
# @return mtx  mtxfitからfeedbackだけ変化が抑制された行列 (3x3)
def get_suppressed_mtx(mtxfit, cx, cy, feedback):

    angle, center, scale = getRotateShift(mtxfit)

    # 倍率、回転量の抑制 (逆回し)
    scale = 1/scale
    angle = (feedback-1) * angle

    # 中心座標の移動先を計算
    mx = mtxfit[0, 0] * cx + mtxfit[0, 1] * cy + mtxfit[0, 2]
    my = mtxfit[1, 0] * cx + mtxfit[1, 1] * cy + mtxfit[1, 2]

    # 抑制された行列 mtx の生成 (3x3)
    mtxback = cv2.getRotationMatrix2D((mx, my), angle, scale)
    mtxback = np.concatenate((mtxback, np.array([[0.0, 0.0, 1.0]])))
    mtx = np.dot(mtxback, mtxfit)

    # 中心座標の移動先を計算
    mx = mtx[0, 0] * cx + mtx[0, 1] * cy + mtx[0, 2]
    my = mtx[1, 0] * cx + mtx[1, 1] * cy + mtx[1, 2]

    # 移動量を抑制 (画像中心に寄せる)
    mtx[0, 2] = (mx * feedback + cx * (1-feedback)) - mtx[0, 0] * cx - mtx[0, 1] * cy
    mtx[1, 2] = (my * feedback + cy * (1-feedback)) - mtx[1, 0] * cx - mtx[1, 1] * cy

    return mtx


####################################
## main
if __name__ == '__main__':

    flist = glob.glob("E:\\tmp\\blog\\20200729\\*.jpg")  # ソースの写真群がある場所
    save = "E:\\tmp\\blog\\20200729\\dst\\"
    img1 = cv2.imread(flist.pop(0))
    cv2.imwrite(save + "{0:03d}".format(0) + ".jpg", img1[180:-180, 320:-320, :])

    h, w = img1.shape[:2]

    # 探索領域を絞る
    sx = w//4
    ex = w - sx
    sy = h//4
    ey = h - sy

    mtx1 = np.eye(3)
    kp1, des1 = get_keypoints(img1, (sx, sy), (ex, ey))

    for index, fname in enumerate(flist):
        img2 = cv2.imread(fname)

        kp2, des2 = get_keypoints(img2, (sx, sy), (ex, ey))
        pt1, pt2 = get_matcher(kp2, des2, kp1, des1)

        # アフィン行列の推定
        mtx = cv2.estimateAffinePartial2D(pt1, pt2)[0]
        mtx = np.concatenate((mtx, np.array([[0.0, 0.0, 1.0]])))

        # 現在フレームへの行列
        mtx2 = np.dot(mtx1, mtx)
        mtx2 = get_suppressed_mtx(mtx2, w//2, h//2, 0.8)

        # 前フレームへの行列
        mtx4 = np.dot(mtx2, np.linalg.inv(mtx))
        mtx4 = np.dot(mtx4, np.linalg.inv(mtx1))

        img1 = cv2.warpAffine(img1, mtx4[:2, :], (w, h))
        img1 = cv2.warpAffine(img2, mtx2[:2, :], (w, h), borderMode=cv2.BORDER_TRANSPARENT, dst=img1)

        mtx1 = mtx2
        kp1 = kp2
        des1 = des2

        cv2.imwrite(save + "{0:03d}".format(index+1) + ".jpg", img1[180:-180, 320:-320, :])