OpenCVを使わず幾何変換する
OpenCVが使えない, 使いたくない, その他の場合1.
import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt
import cv2
参考までにここで紹介する OpenCV を使わない方法より OpenCV を素直に使った方が数十倍速い.
画像データ
画像をプログラム上で扱う時, 通常は2～3次元配列のデータとして扱うことになる.
インデクスの [i, j, k] は前の2つ i, j がそれぞれ $y$, $x$ 座標 (ピクセルの位置) を表し, 3つ目 k は RGB や CMYK 等の画素値ベクトルのインデクスを示す.
2次元配列の場合, モノトーン画像等の画素値が1変数となる画像データを表す.
C 系のプログラミング言語では for ループの入れ子構造とデータアクセスの順序の関係で, 第1インデクスの i が縦方向の位置を表すと考えるため, (i, j) - ($y$, $x$) という対応付けがアルファベット順と逆になることに注意.
幾何変換
画像の幾何変換とは写像 $f$ によって座標 $[y, x]$ を他の座標 $[w, z] = f([y, x])$ に移すことで, 新たな画像に変換することを指す.
変換後画像の画素 $\mathrm{Dst}$ と元画像の画素 $\mathrm{Src}$ は $\mathrm{Dst}[w, z] = \mathrm{Src}[y, x]$ で対応する.
例えばアファイン変換なら座標の配列$[y, x]$ を行列 $A$ とベクトル $b$ で $f([y,x]) = A[y, x] + b$ と変換して新たな座標が計算される2.
$[y, x]$ として入力する画像データの全画素の座標 (インデクス) を渡した時, 変換後の座標 $[w, z]$ は画像データを構成できるものになるとは限らない.
例えば小数座標や負の座標は配列のインデクスにはならないし, 配列を構成するのに必要な要素のインデクスが得られないこともある.
よって, 変換後の画像データを作成するためには, 変換後の全座標 (インデクス) $[w, z]$ に対応する変換前の位置座標 $[y, x]$ を逆算する必要がある.
この時, 逆算された $[y, x]$ もまた元の画像データを構成するインデクスになるとは限らないので, 元画像の画素値を使って $[y, x]$ における値を補間しなければいけない.
画像データのような2次元に配置されたデータの補間手法には最近傍補間, 双線型 (bilinear) 補間, 双三次 (bicubic) 補間等があり, 詳細は下のサイト等を参照.
cf. 画素の補間（Nearest neighbor,Bilinear,Bicubic） 画像処理ソリューション
また, $[y, x]$ が補間できない座標の場合 (双線型, 双三次補間での負の座標等), 外挿も行うことになる.
補間が「周囲のデータの中間値を取る」という方向性が決まっているのに対して外挿は「画像内に写ってない箇所の画素値を予測する」という自由度の高い問題なので, 定数で埋める, 端の画素値を延長する, 元画像の画素値を据え置く等, 補間と比較してバリエーションが多い.
OpenCV での幾何変換
OpenCV を使わない幾何変換の前に OpenCV を使う方法を魚眼変換を例に解説する.
OpenCV (-Python) では cv2.remap に幾何変換の実装があり,
cv2.remap(img, x, y, interpolation, borderMode, borderValue)
で元画像 img, 変換前座標 x, y (ここはアルファベット順に書かれる), 補間方法 interpolation, 外挿方法 borderMode, borderValue を指定する.
次のコードでは interpolation=cv2.INTER_LINEAR で双線型補間を行って borderMode=cv2.BORDER_CONSTANT, borderValue=0 として定数0で外挿している.
画像は SIDBA (標準画像データベース) より引用したマンドリル.
# 画像読み込み, 元画像と同じサイズの変換後画像のインデクス配列生成
# 画像は https://sipi.usc.edu/database/ から
img = cv2.imread('mandrill.tiff')[..., ::-1]
wz = np.indices(img.shape[:2]).astype(np.float32)
# 魚眼変換の逆算
# cf. http://www.persfreaks.jp/main/aov/fisheye/
center = np.max(wz, axis=(1, 2), keepdims=True)/2
radius = np.linalg.norm(wz-center, axis=0)
radius = radius/np.max(radius)
theta = 3/8*np.pi*radius
yx = 0.7*np.tan(theta)*(wz-center)/radius+center
# OpenCV での幾何変換
y, x = yx
img_cv = cv2.remap(img, x, y, interpolation=cv2.INTER_LINEAR,
borderMode=cv2.BORDER_CONSTANT, borderValue=0)
# オリジナルのマンドリルと魚眼マンドリル
plt.imshow(img)
plt.show()
plt.imshow(img_cv)
plt.show()
OpenCV を使わない幾何変換
OpenCV を使わない幾何変換では変換前座標 $[y, x]$ を算出するところまでは同じだが, 補間と外挿で新たな座標での画素値を計算する必要がある.
ここでは NumPy, SciPy を使った方法を紹介するが, 行列演算を for 文で展開すれば Python 以外の言語でも実装できる.
NumPy による双線型補間・0埋め外挿
双線型補間は周囲4ヶ所の格子点 (座標が整数となる点) の重み付き平均で補間点の画素値を決める補間となる.
重みは各格子点と補間点との近さで決められる.
詳細は上で挙げたサイト等を参照.
NumPy を使うと行列演算で全座標について一括して書くことができる.
ただし前述のように補間点と周囲の格子点が元画像のインデクスの範囲を外れる可能性もあるので, パディングして範囲を拡張しておく必要がある.
この時, 外挿方法として0埋めを行うのであれば0パディングと双線型補間だけで追加の処理が不要になる.
# 画像読み込み, 元画像と同じサイズの変換後画像のインデクス配列生成
img = cv2.imread('mandrill.tiff')[..., ::-1]
wz = np.indices(img.shape[:2])
# 魚眼変換の逆算
center = np.max(wz, axis=(1, 2), keepdims=True)/2
radius = np.linalg.norm(wz-center, axis=0)
radius = radius/np.max(radius)
theta = 3/8*np.pi*radius
yx = 0.7*np.tan(theta)*(wz-center)/radius+center
# 双線型補間, 0埋め外挿
# 0パディング
yx_fl = np.floor(yx).astype(int)
yx_lim = np.stack([np.minimum(np.min(yx_fl, axis=(1, 2)), 0),
np.maximum(np.max(yx_fl, axis=(1, 2))+1, np.array(img.shape[:2])-1)])
pad = np.array([[*(-yx_lim[0]), 0], [*(yx_lim[1]-np.array(img.shape[:2])+1), 0]]).T
img_pad = np.pad(img, pad, mode='constant', constant_values=0)
# 0パディングで拡張した画像データを使って双線型補間
# この場合は0埋め外挿のための追加処理は不要
weights = np.stack([yx_fl+1-yx, yx-yx_fl])
yx_sh = yx_fl-yx_lim[0, :, np.newaxis, np.newaxis]
delta = np.array([[[0, 1]], [[0, 0]]])+np.array([[[0], [0]], [[0], [1]]])
indices = yx_sh[:, np.newaxis, np.newaxis]+delta[..., np.newaxis, np.newaxis]
img_mat = img_pad[indices[0], indices[1]]
img_np = np.einsum('ihw, ijhwc, jhw -> hwc', weights[:, 0], img_mat , weights[:, 1]).astype(int)
# NumPY 版魚眼マンドリル
plt.imshow(img_np)
plt.show()
scipy.interpolate による補間・外挿
Python では SciPy の scipy.interpolate に補間と簡単な外挿の実装がある.
画像データは座標がグリッド状に並ぶのでscipy.interpolate.interpn が使える.
欠損値ありや非矩形領域の場合等, グリッドでないデータを扱うなら代わりにscipy.interpolate.griddata を使う3.
scipy.interpolate.interpn の引数は変換後座標 points, 補間の参照値 value (ここでは画素値), 変換前座標 xi, 補間方法 method, 外挿に関するオプション bounds_error, fill_value.
bounds_error が True の場合は外挿せずエラーを返し, False ならば fill_value の値で定数埋め, または fill_value=None で線型外挿を行う.
次のコードでは method='linear' で双線型補間を, bounds_error=Fasle, fill_value=0 で0埋め外挿を行う.
# 画像読み込み, 元画像と同じサイズの変換後画像のインデクス配列生成
img = cv2.imread('mandrill.tiff')[..., ::-1]
wz = np.indices(img.shape[:2])
# 魚眼変換の逆算
center = np.max(wz, axis=(1, 2), keepdims=True)/2
radius = np.linalg.norm(wz-center, axis=0)
radius = radius/np.max(radius)
theta = 3/8*np.pi*radius
yx = 0.7*np.tan(theta)*(wz-center)/radius+center
# 補間, 外挿
# OpenCV の場合と渡すデータ形式が少し変わる
warr = np.arange(img.shape[0])
zarr = np.arange(img.shape[1])
img_sp = scipy.interpolate.interpn((warr, zarr), img, np.transpose(yx, [1, 2, 0]),
method='linear', bounds_error=False, fill_value=0).astype(int)
# SciPy 版魚眼マンドリル
plt.imshow(img_sp)
plt.show()
