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

　手振れ補正のPythonコードをOpenCVの力を借りて書いてみます。上が手振れ、下が補正。こんなサンプルです。

手振れ補正
画像処理内容
ちょっとしたコツ
結果
まとめ
全コード

手振れ補正

　画像のふらふらを解消する手振れ補正。それを画像処理で行う、いわゆる電子式手振れ補正というやつをコーディングしてみようと思います。

　昨今のスマホでは標準装備ですね。これをコーディングしてみます。動画ファイルから1フレームずつ取り出してもよいのですが、扱いやすいので、今回ソースとするのは連番JPEG画像とします。

画像処理内容

　nフレーム目の画像とn+1フレーム目の画像とを位置合わせして、その位置合わせ結果からn+1フレーム目を幾何変換（今回はアフィン変換）を行うだけです。これだけでそれっぽくなります。似たようなコードもいろんなところに落ちてます。

　位置合わせには、nフレーム目の画像とn+1フレーム目の画像の特徴点をそれぞれ取り出し、その両者の突き合わせをして、対応点同士のマッチングを取り、最後アフィン変換します。全体の流れは、

特徴点取り出し（AKAZE）
特徴点同士突き合わせ（Brute-Force）
変換行列推定（estimateAffine）
変換（warpAffine）

　この繰り返しです。

特徴点取り出し

　ここでは、2枚の画像からそれぞれAKAZEという回転、拡大縮小耐性の強い、特徴点抽出手法を使って特徴点を取り出すことにします。（他にもいろいろ手段はありますが、権利的な部分でこれを採用。OpenCVを入れるだけで使えるのでお手軽。）

　コードにしてみると以下

##
# @brief AKAZEによる画像特徴量取得
# @param img 特徴量を取得したい画像（BGR順想定）
# @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_BGR2GRAY)
    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)

　指定した任意の矩形内に限定して、特徴点を取るようにしています。これにより画像の中心部をターゲットにした位置合わせができるようになります。この探索にはグレースケールの画像を用いています。

　この画像（2560×1440 pixel）をソースにすると、得られる特徴点は

　赤点で打った箇所がすべて特徴点です。このようにして特徴点列が取得できます。16571点も取れています。

特徴点同士突き合わせ

　続いて得られた特徴点同士のマッチングを取ります。

　コードにしてみると以下

##
# @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):
    bf = cv2.BFMatcher()
    matches = bf.knnMatch(des1, des2, k=2)
 
    good = []
    for m, n in matches:
        if m.distance < 0.7 * n.distance:
            good.append(m)
 
    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

　Keypoint同士の突き合わせを総当たり（Brute-Force, KNN）で行い、いい感じに一致するものを残します。queryIdxがkp1に対する座標。trainIdxがkp2に対する座標です。二つの座標対として結果を出力します。

　この2枚から得られた対応点は

　だいぶ減りましたが、アフィン行列を推定するには十分です。

変換行列推定

　この2つの対応点列から変換ルールを求めます。今回はアフィン変換とします。他にもいろいろ手段はあると思いますが、あまりに非線形な変換しまうと、フレームをまたぐにつれて画像に歪みが出てきてしまうので、アフィン変換程度が適当だと思われます。そのため変換ルールはいわゆるアフィン行列になります。逆にシーンが固定だったら、射影変換使ってびったびたに合わせるのもありだと思います。

　コードは

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

　これだけです。この関数の詳細はこちら。いい感じに外れ値を除いた行列を返してくれます。

　ただしこのアフィン変換行列は拡大、縮小成分も含んでいます。シーンが推移する動画において、あるフレームに写る画像に向かってサイズが固定されるのは望ましくありません。なので補正するのは回転成分とシフト成分だけにします。この推定した行列から、回転成分とシフト成分を取り出します。

　具体的には回転中心座標と回転角度が求まれば十分です。ココで取り出し方の詳細があります。それを持ってきます。基本的には連立方程式を解いているだけです。一応倍率も返すようにしています。

##
# @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

　このようにして得られた回転中心、回転角度の情報を用いて、

1	`mtx` `=` `cv2.getRotationMatrix2D((center[0], center[1]), angle,` `1)`

　から新たな回転行列を算出します。

変換

　こうして求めた行列をn+1フレーム目の画像に対して適用します。

cv2.warpAffine(img2, mtx, (w, h), borderMode=cv2.BORDER_TRANSPARENT, dst=img1)

　アフィン変換時のborderModeに前フレーム(nフレーム目の画像)を用います。これにより回転、シフトしたときに出てくる画像外の値をもっともらしくすることができます。borderModeに関してはコチラ。

　上がはみ出し部が黒になってしまった例。下が前フレームの画像にオーバーレイさせたもの。隣の手が邪魔ですが、それ以外は黒よりはマシになっていると思います。

　基本はこの繰り返しです。

ちょっとしたコツ

　馬鹿正直に毎フレーム前後の特徴点を取って、突き合わせるのでは処理時間がかかっていけません。なので一度取った特徴点は残しておき、次のフレームのマッチングで使いまわします。

　n+1フレーム目の変換を行う際に、nフレーム目の変換時に作ったkeypoint[n]を使いまわす。という流れです。

　ただここで問題になるのが、一度取った特徴点は、その後アフィン変換がかかっているため場所がずれています。下図でいうところの、一度取得したkeypoint[n]はkeypoint[n]’へ動いています。

　図中keypoint[n]とkeypoint[n+1]でマッチングを取って、その結果得られた行列でn+1フレーム目を変換したところで、それはあくまで変換前のnフレーム目と位置が合うだけで、変換がかかったnフレーム目（図中下）の画像とは位置が合いません。

　真面目にやると、変換後の特徴点keypoint[n]’を求めて、それとkeipoint[n+1]をマッチング取る必要が出てきます。これはこれでめんどくさいです。

　そこで、nフレーム目で行ったアフィン変換の行列も残しておいて、その行列とn+1フレーム目で求めた行列の積を取ることで、2段階の変換を1つの行列で表現させます。アフィン変換の組み合わせは行列の積で表現できるので、このようにシンプルな記述ができます。この詳細はコチラ。このページでは射影変換との結合を書いてますが、内容的にはアフィン変換行列同士も同じです。

$$ \left( \begin{matrix} x’ \\ y’ \\ 1 \end{matrix} \right) = M_{before}\cdot M_{after}\left( \begin{matrix} x \\ y \\ 1\end{matrix} \right) $$

　コードにすると、

m = np.array([[0.0, 0.0, 1.0]])
 
for index, fname in enumerate(flist):
    img2 = cv2.imread(fname)
    kp2, des2 = get_keypoints(img2)
    pt1, pt2 = get_matcher(kp2, des2, kp1, des1)
 
    # アフィン行列の推定
    mtx = cv2.estimateAffinePartial2D(pt1, pt2)[0]
    angle, center, scale = getRotateShift(mtx)
    mtx2 = cv2.getRotationMatrix2D((center[0], center[1]), angle, 1)
    mtx2 = np.concatenate((mtx2, m))
 
    mtx1 = np.dot(mtx1, mtx2)
    img1 = cv2.warpAffine(img2, mtx1[:2, :], (w, h), borderMode=cv2.BORDER_TRARENT, dst=img1)
 
    kp1 = kp2
    des1 = des2

　こんな感じのループを書いてみました。

　読み込んだ画像のポイントと前回のポイントからマッチングを取り、得られた対から行列を推定したのちに、

1	`mtx1` `=` `np.dot(mtx1, mtx2)`

　を使って変換行列の積を求めています。アフィン行列が2×3の行列なので、積を取るために3×3に拡張しています。

m = np.array([[0.0, 0.0, 1.0]])
mtx2 = np.concatenate((mtx2, m))

　最後にkeypointを次で使うために残しています。

結果

　何も補正しない場合。手持ちしたカメラで連写設定で撮った写真です。歩いてます。グラグラです。

　今回のコードの結果。なめらかになっていますが、やはり四隅は少し落とした方がよいですね。

まとめ

　まずはとりあえずで作ってみました。ただこのままだと視点が最初のフレームに縛られます。このサンプルでも徐々に視点が下がってきてしまいます。動画を想定するとピッタリと位置合わせしすぎない方がシーンチェンジに適応できるはずです。

　今度は位置合わせにその縛りを入れて、より動画っぽい手振れ補正にしてみます。

　その２に続く

全コード

　上記に対して少しバカ除けの記述を追加しています。

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　回転中心　scale 拡大率
def getRotateShift(mat):
    da = -math.atan2(mat[1, 0], mat[0, 0])  # ラジアン
    scale = 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, scale
 
 
####################################
## main
if __name__ == '__main__':
 
    flist = glob.glob("img\\*.jpg")  # ソースの写真群がある場所
    save = "dst\\"
    img1 = cv2.imread(flist.pop(0))
    cv2.imwrite(save + "{0:03d}".format(0) + ".jpg", img1)
 
    h, w = img1.shape[:2]
    kp1, des1 = get_keypoints(img1)
    mtx1 = np.eye(3)
    m = np.array([[0.0, 0.0, 1.0]])
 
#    debug code (特徴点列)
#    for pos in kp1:
#        cv2.circle(img1, (int(pos.pt[0]), int(pos.pt[1])), 5, (0, 0, 255), thickness=-1)
#    cv2.imwrite("point0.jpg", img1)
 
    for index, fname in enumerate(flist):
        img2 = cv2.imread(fname)
 
        kp2, des2 = get_keypoints(img2)
        pt1, pt2 = get_matcher(kp2, des2, kp1, des1)
 
        # アフィン行列の推定
        mtx = cv2.estimateAffinePartial2D(pt1, pt2)[0]
        angle, center, scale = getRotateShift(mtx)
        mtx2 = cv2.getRotationMatrix2D((center[0], center[1]), angle, 1)
        mtx2 = np.concatenate((mtx2, m))
 
        mtx1 = np.dot(mtx1, mtx2)
        img1 = cv2.warpAffine(img2, mtx1[:2, :], (w, h), borderMode=cv2.BORDER_TRANSPARENT, dst=img1)
 
        kp1 = kp2
        des1 = des2
 
        cv2.imwrite(save + "{0:03d}".format(index+1) + ".jpg", img1)