レンズ歪みを扱う関数(distortion, undistortion)のChainer実装が欲しかったので探したが、見つけられなかったので自作した。なかなか大変だったのでメモを残す。
distortion
OpenCVに習って、歪み係数を $(k_1, k_2, p_1, p_2, k_3, k_4, k_5, k_6)$ としたときの歪み関数を以下の式で定める。
$$(x', y') = {\rm distort}(x, y)$$
$$ \begin{aligned}
r^2 &= x^2 + y^2 \\
x' &= x\frac{1+k_1r^2+k_2r^4+k_3r^6}{1+k_4r^2+k_5r^4+k_6r^6}+2p_1xy+p_2(r^2+2x^2)\\
y' &= y\frac{1+k_1r^2+k_2r^4+k_3r^6}{1+k_4r^2+k_5r^4+k_6r^6}+2p_2xy+p_1(r^2+2y^2)\\
\end{aligned}$$
distortionはこれをそのまま実装すれば良い。
undistortion
逆関数 $(x, y) = {\rm undistort}(x', y')$ を陽に求める事は出来ないので、数値解法によって求める事が必要。
初期近似解の導出
まず最初に出来るだけ良い初期値を求める。以下の文献によればradial distortionの部分(第一項目)は厳密に解く事ができる。と言っても無限級数としてだが。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4934233/
まず、
$$P(x) = 1+k_1x+k_2x^2+k_3x^3$$
と置く。本関数は$x\geq 0$において単調増加かつ連続である。これを用いるとradial distortionは
$$
\begin{pmatrix} x'\\y'\end{pmatrix} = P(r^2)\begin{pmatrix} x\\y\end{pmatrix} \\
$$
と表せる。同様に、逆方向のradial distortionを
$$\begin{aligned}
Q(x) &= 1+a_1x+a_2x^2+a_3x^3+\cdots\\
\begin{pmatrix} x\\y\end{pmatrix} &= Q(r'^2)\begin{pmatrix} x'\\y'\end{pmatrix}\\
\end{aligned}$$
と表すと、これらは逆関数なので
$$P(r^2)Q(r'^2) = 1$$
が成立する。また $r'^2 = r^2P(r^2)^2$ であるので
$$P(r^2)Q(r^2P(r^2)^2) = 1$$
が任意の $r\geq 0$ について成立する。$x=r^2 \geq 0$ と置換して
$$P(x)Q(xP(x)^2) = 1$$
となる。これを展開して係数を求めればよいが、面倒なのでsympyで計算する。
from sympy import Symbol, Poly, solve
DEGREE = 3
x = Symbol('x')
Ks = [Symbol('k%d'%i) for i in range(1, 4)]
As = [Symbol('a%d'%i) for i in range(1, DEGREE+1)]
P = sum(k * x**d for d, k in enumerate([1] + Ks))
Q = sum(a * x**d for d, a in enumerate([1] + As))
print('P(x) =', P)
print('Q(x) =', Q)
coeffs = Poly(P * Q.subs(x, x*P**2) - 1, x).coeffs()
answer = solve(coeffs[-len(As):], As)
for a in As:
print(a, '=', answer[a])
実行結果は以下の通り(とりあえず3次で打ち切った)。
P(x) = k1*x + k2*x**2 + k3*x**3 + 1
Q(x) = a1*x + a2*x**2 + a3*x**3 + 1
a1 = -k1
a2 = 3*k1**2 - k2
a3 = -12*k1**3 + 8*k1*k2 - k3
Newton-Raphson法での解改良
続いて、上で求めた初期解をNewton-Raphson法により改良する。
$$ f({\mathbf x}) = {\rm distort}({\mathbf x}) - {\mathbf x'}$$
と置く。これを微分すると
$$\nabla f({\mathbf x}) = \begin{pmatrix}
g(r^2) + 2x^2g'(r^2) + 2p_1y + 6p_2x & 2xyg'(r^2)+2p_1x+2p_2y\\
2xyg'(r^2) + 2p_1x + 2p_2y & g(r^2) + 2y^2g'(r^2) + 2p_2x + 6p_1y\\
\end{pmatrix}$$
となる。但し
$$g(x) = \frac{1+k_1x+k_2x^2+k_3x^3}{1+k_4x+k_5x^2+k_6x^3}$$
これを使って、前節で求めた値を初期値として用いて更新式
$${\mathbf x_{i+1}} = {\mathbf x_i}-\nabla f({\mathbf x_i})^{-1}({\rm distort}({\mathbf x_i})-{\mathbf x'})$$
を数回まわして解を改善する。
実装
こちら。
https://github.com/Idein/chainer-graphics
実験
下の画像の左がOpenCVの cv2.undistort レナ画像を変換したもの。右が chainer_graphics.camera.undistort_image で同じパラメータを用いて変換したもの。
正しく出来ていそう。
これを chainer_graphics.camera.distort_image で同じパラメータで逆変換すると元に戻る。こちらも正しく出来ていそう。
