はじめに
正直内容なんて読まなくて良いので実験結果だけでも覗いていってほしい。
SG法（Savitzky-Golay法）というと１次元データのピーク形状にある程度追従しつつノイズを削減することができるフィルタとしての印象が強い。
また、微分係数が同時に得られる点も大きな利点である。
以前このSGフィルタを２次元に拡張して画像平滑化フィルタとして使用してみる試みを投稿した。
当時は２次元SGフィルタの活用に思い至らなかったが、平滑化と同時に偏微分係数が得られる特徴はヘッセ行列の各要素を自然な流れから得ることができる強力なツールとして機能する。
この利点を、ピーク座標のサブピクセル推定タスクに適用してみる。
どのような理屈で説明がつくのか私にも解らないが、このタスクの精度に関しては他のどのような手法も右に出ないかもしれない。
サブピクセルピーク座標推定
テンプレートマッチングや天体画像解析など、画像上のピークの位置をピクセルより細かいスケールで特定したいシーンは多岐にわたる。
こういった場合に一般的に適用される方法として以下のような選択肢があるだろう。
- 重心計算
- パラボラフィッティング
- 等角直線フィッティング
- 多項式やガウシアンでフィッティング
- スプラインでフィッティング
これらは一長一短で、そのピーク構造の偏りや広がり、ノイズに応じて精度が出るとき出ないときがある。
例えば以下は位相限定相関法（POC）の位相相関マップにおけるピーク周辺を切り出した画像であるが、ピーク周辺の構造やノイズの乗り方が変化したときに、安定して数十分の1以下の精度で頂点座標を推定することが難しそうなことは直感的にわかるだろう。
ちなみに正解はここである。
2DSGフィルタによるピーク座標推定
詳しいSGフィルタ・2DSGフィルタの解説は過去の記事に譲るとして、簡単にSGフィルタについて説明すると、
等間隔なデータを注目点の周りで多項式フィッティングしたときの、フィッティング結果の注目点の座標における値を計算するには、データに依存しない定数のカーネルをかければ良いというものであり、計算の際に同時に偏微分係数を得られるカーネルも出てくるという特徴がある。
移動平均を0次の多項式フィッティングの結果と捉えれば、それを高次に拡張したデータ平滑化手法と言える。
イメージは英語版Wikipediaに掲載の以下のgifが大変わかりやすい。
1次元データに適用されることの多いSGフィルタであるが、2次元以上の高次にも拡張可能であり、3次以上のSGフィルタであれば2階以上の偏微分係数を得ることができるカーネルが得られるため、極大値探索に必要なX方向勾配（Ix）、Y方向勾配（Iy）、X方向曲率（Ixx）、Y方向曲率（Iyy）、Ixyがノイズに頑健に計算できる。
これら係数が安定に得られれば、ニュートン法によって2次収束でピーク座標に辿り着くことができる。
導出は行わないが、ニュートン法では以下の式で値を更新していくのであった。
$$\mathbf{x}_{k+1} = \mathbf{x}_k - H(\mathbf{x}_k)^{-1} \nabla f(\mathbf{x}_k)$$
ここで$H(\mathbf{x}_k)$はヘッセ行列である。
画像（$I(x,y)$）では、
\begin{pmatrix}
x_{k+1}\\
y_{k+1}
\end{pmatrix}
=
\begin{pmatrix}
x_{k}\\
y_{k}
\end{pmatrix}
-
\begin{pmatrix}
I_{xx} & I_{xy} \\
I_{xy} & I_{yy}
\end{pmatrix}^{-1}
\begin{pmatrix}
I_x \\ I_y
\end{pmatrix}
となる。これを解くと更新量は
$$\Delta x
= \frac{-I_x I_{yy} + I_y I_{xy}}{I_{xx}I_{yy}-I_{xy}^2},
\qquad
\Delta y
= \frac{-I_y I_{xx} + I_x I_{xy}}{I_{xx}I_{yy}-I_{xy}^2}$$
$$(x_{k+1},y_{k+1})
= (x_k + \Delta x,; y_k + \Delta y)$$
として求まる。
この結論は2次近似の結果であり、実際にこのステップで更新すると極値を通り越して収束しないことが多いため、適当な定数$\le1$をかけて更新する。
この必要な各係数を2DSG法で求めるのが提案手法である。
実装
シンプルにピーク推定位置の周辺をサンプリングし、2DSGフィルタで得られたステップ更新するだけのもの。
サンプリング手法をバイリニアなど線形サンプリングに変更すると精度が半分くらいになる。
import numpy
import cv2
def SGFilter2D(m, N):
mask = (np.c_[:m+1] + np.r_[:m+1]).flatten() <= m
X = np.zeros(((2*N+1)**2, mask.sum()))
for i, j in np.indices((2*N+1, 2*N+1)).T.reshape(-1, 2):
X[i*(2*N+1)+j] = ((i-N)**np.c_[:m+1] * (j-N)**np.r_[:m+1]).flatten()[mask]
coords = np.mgrid[:m+1, :m+1].transpose(1, 2, 0).reshape(-1, 2)[mask]
order = np.array(sorted(range(coords.shape[0]), key=lambda i: (coords[i, 0] + coords[i, 1], coords[i].min())))
C = np.linalg.pinv(X[:, order]).reshape(-1, 2*N+1, 2*N+1)
return C
def SG_peakSubPix(image: np.ndarray, win_size: int=5, degrees: int=6, max_iter: int=15, interpolation: int=cv2.INTER_CUBIC, eps: float=1e-7, lr=0.7) -> np.ndarray:
SGs = SGFilter2D(degrees, win_size//2)[1:6] # 偏微分フィルタのみ使用
peak = np.asarray(np.unravel_index(image.argmax(), image.shape))
mg = np.mgrid[-(win_size//2): win_size//2: 1j*win_size, -(win_size//2): win_size//2: 1j*win_size]
for _ in range(max_iter):
AOI = cv2.remap(image, *(mg + peak[:, None, None]).astype(np.float32)[::-1], interpolation=interpolation)
Ix, Iy, Ixx, Iyy, Ixy = np.einsum('ijk, jk -> i', SGs, AOI)
detH = (Ixx+1e-15) * (Iyy+1e-15) - Ixy**2
dx = - (Iyy * Ix - Ixy * Iy) / detH
dy = - (Ixx * Iy - Ixy * Ix) / detH
peak = peak + np.clip(np.array([dy, dx]).T, -1, 1) * lr
if np.hypot(dx, dy) < eps: break
return peak
比較テスト
テスト１ テンプレートマッチング
Describable Textures Dataset (DTD)より、どこを切り出しても満遍なく特徴がありそうなcracked_0129をチョイス。
- 480x484pxのcracked_0129（グレースケール化・0~1に正規化）のランダムな非整数座標から
cv2.remapやscipy.ndimage.map_coordinatesを用いて64x64pxの画像を切り出す。 -
cv2.matchTemplate（TM_CCOEFF_NORMED）で類似度マップを計算する。 - 提案手法と、比較対象として以下の各手法でそれぞれサブピクセル座標を推定する。
- 以上を1000回繰り返す。（各手法には同じ画像ペアを使用）
（提案手法以外の各手法はGeminiに書かせたもの）
| 手法 | 表記 | 簡単な説明 |
|---|---|---|
| 2次元SGフィルタ | 2DGS | 提案手法 |
| 重心 | Centroid | ピークを中心とした周辺の画素値を重みとし、加重平均（重心）を計算してサブピクセル位置を推定する。実装が簡単で計算コストが低い。 |
| パラボラフィッティング | Parabola | ピーク点と隣接する2点（1次元的に）を用いて2次関数（放物線）を当てはめ、その頂点（極値）をサブピクセル位置とする。 |
| 2次多項式フィッティング | Biquadratic | ピークを中心とする$3 \times 3$パッチの9点*2変数2次多項式を当てはめ、その極値を連立方程式を解いて求める。Parabolaの2次元拡張版。 |
| 1次元ガウシアンフィッティング | Gaussian1D | ピーク点とその隣接点（1次元的に）に対数ガウシアン関数（ガウス関数の対数）を当てはめ、その極値を計算する。Parabolaより丸い応答形状に適する。 |
| 2次元ガウシアンフィッティング | Gaussian2D | $3 \times 3$などのパッチに2次元ガウス関数を非線形最小二乗法で当てはめ、その中心を求める。精度が高いが、計算負荷も高い。 |
| 2次元ガウシアンフィッティング（Log） | Log-Gaussian | $3 \times 3$パッチの対数値に2変数2次多項式（ガウス関数の対数）を当てはめ、その極値を計算する。Gaussian2Dより計算が速い。 |
| スプラインフィッティング | Spline | ピーク周辺の点群にRectBivariateSpline（双三次スプラインなど）で滑らかな関数を補間し、最適化手法（L-BFGS-Bなど）を用いて補間関数の極値を探索する。 |
各手法の実装
import numpy as np
from scipy.interpolate import RectBivariateSpline
from scipy.optimize import minimize, curve_fit
# --- 共通ヘルパー ---
def get_peak_integer(response_map):
peak_coord = np.unravel_index(np.argmax(response_map), response_map.shape)
iy, ix = peak_coord
h, w = response_map.shape
is_edge = (iy == 0) or (iy == h - 1) or (ix == 0) or (ix == w - 1)
return (iy, ix), is_edge
# --- 1. 重心法 ---
def find_peak_centroid(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
patch = response_map[iy-1:iy+2, ix-1:ix+2]
patch = np.maximum(patch, 0)
total_weight = np.sum(patch)
if total_weight == 0:
return float(iy), float(ix)
y_indices, x_indices = np.indices((3, 3)) - 1
dy = np.sum(y_indices * patch) / total_weight
dx = np.sum(x_indices * patch) / total_weight
return iy + dy, ix + dx
# --- 2. 放物線当てはめ ---
def find_peak_parabolic(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
z0 = response_map[iy, ix]
dy, dx = 0.0, 0.0
ym1 = response_map[iy-1, ix]
yp1 = response_map[iy+1, ix]
denom_y = 2 * z0 - ym1 - yp1
if denom_y > 0:
dy = (ym1 - yp1) / (2 * denom_y)
xm1 = response_map[iy, ix-1]
xp1 = response_map[iy, ix+1]
denom_x = 2 * z0 - xm1 - xp1
if denom_x > 0:
dx = (xm1 - xp1) / (2 * denom_x)
if abs(dy) > 1.0 or abs(dx) > 1.0:
return float(iy), float(ix)
return iy + dy, ix + dx
# --- 3. 2次多項式補間 ---
def find_peak_biquadratic(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
patch = response_map[iy-1:iy+2, ix-1:ix+2]
z = patch.flatten()
x = np.array([-1, 0, 1, -1, 0, 1, -1, 0, 1])
y = np.array([-1, -1, -1, 0, 0, 0, 1, 1, 1])
A = np.vstack([x**2, y**2, x*y, x, y, np.ones(9)]).T
try:
coeffs = np.linalg.lstsq(A, z, rcond=None)[0]
except np.linalg.LinAlgError:
return float(iy), float(ix)
a, b, c, d, e, _ = coeffs
M = np.array([[2*a, c], [c, 2*b]])
v = np.array([-d, -e])
try:
det = 4*a*b - c*c
trace = 2*a + 2*b
if det <= 0 or trace >= 0:
return float(iy), float(ix)
delta = np.linalg.solve(M, v)
dx, dy = delta[0], delta[1]
except np.linalg.LinAlgError:
return float(iy), float(ix)
if abs(dy) > 1.0 or abs(dx) > 1.0:
return float(iy), float(ix)
return iy + dy, ix + dx
# --- 4. 対数ガウス当てはめ ---
def find_peak_log_gaussian(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
patch = response_map[iy-1:iy+2, ix-1:ix+2]
if np.any(patch <= 1e-6): # ゼロまたは負の値チェック
return find_peak_parabolic(response_map)
log_patch = np.log(patch.flatten())
x = np.array([-1, 0, 1, -1, 0, 1, -1, 0, 1])
y = np.array([-1, -1, -1, 0, 0, 0, 1, 1, 1])
A = np.vstack([x**2, y**2, x, y, np.ones(9)]).T
try:
coeffs = np.linalg.lstsq(A, log_patch, rcond=None)[0]
except np.linalg.LinAlgError:
return float(iy), float(ix)
a, b, c, d, _ = coeffs
dx, dy = 0.0, 0.0
if a >= 0 or b >= 0:
return float(iy), float(ix)
try:
dx = -c / (2 * a)
dy = -d / (2 * b)
except ZeroDivisionError:
return float(iy), float(ix)
if abs(dy) > 1.0 or abs(dx) > 1.0:
return float(iy), float(ix)
return iy + dy, ix + dx
# --- 5. 1Dガウス当てはめ ---
def find_peak_gaussian_1d(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
def parabolic_1d(v_m1, v_0, v_p1):
denom = 2 * v_0 - v_m1 - v_p1
if denom <= 0: return 0.0
return (v_m1 - v_p1) / (2 * denom)
# Y方向
ym1 = response_map[iy-1, ix]
y0 = response_map[iy, ix]
yp1 = response_map[iy+1, ix]
dy = 0.0
if ym1 > 1e-6 and y0 > 1e-6 and yp1 > 1e-6:
log_ym1 = np.log(ym1)
log_y0 = np.log(y0)
log_yp1 = np.log(yp1)
denom_y = 2 * (2 * log_y0 - log_yp1 - log_ym1)
if denom_y > 0:
dy = (log_yp1 - log_ym1) / denom_y
else:
dy = parabolic_1d(ym1, y0, yp1)
else:
dy = parabolic_1d(ym1, y0, yp1)
# X方向
xm1 = response_map[iy, ix-1]
x0 = response_map[iy, ix]
xp1 = response_map[iy, ix+1]
dx = 0.0
if xm1 > 1e-6 and x0 > 1e-6 and xp1 > 1e-6:
log_xm1 = np.log(xm1)
log_x0 = np.log(x0)
log_xp1 = np.log(xp1)
denom_x = 2 * (2 * log_x0 - log_xp1 - log_xm1)
if denom_x > 0:
dx = (log_xp1 - log_xm1) / denom_x
else:
dx = parabolic_1d(xm1, x0, xp1)
else:
dx = parabolic_1d(xm1, x0, xp1)
if abs(dy) > 1.0 or abs(dx) > 1.0:
return float(iy), float(ix)
return iy + dy, ix + dx
# --- 6. 2Dガウス当てはめ ---
def _gaussian_2d(xy, amplitude, x0, y0, sigma_x, sigma_y, theta, offset):
x, y = xy
x0 = float(x0)
y0 = float(y0)
a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
g = offset + amplitude * np.exp( - (a*((x-x0)**2) + 2*b*(x-x0)*(y-y0) + c*((y-y0)**2)))
return g.ravel()
def find_peak_scipy_gaussian(response_map, window_size=3):
(iy, ix), _ = get_peak_integer(response_map)
h, w = response_map.shape
hw = window_size // 2
if iy < hw or iy >= h - hw or ix < hw or ix >= w - hw:
return float(iy), float(ix)
y_min, y_max = iy - hw, iy + hw + 1
x_min, x_max = ix - hw, ix + hw + 1
patch = response_map[y_min:y_max, x_min:x_max]
y, x = np.indices(patch.shape)
xy_data = (x.ravel(), y.ravel())
z_data = patch.ravel()
try:
initial_guess = (patch[hw, hw], hw, hw, 1.0, 1.0, 0, np.min(patch))
popt, _ = curve_fit(_gaussian_2d, xy_data, z_data, p0=initial_guess, maxfev=500)
dx_local = popt[1]
dy_local = popt[2]
final_y = y_min + dy_local
final_x = x_min + dx_local
if abs(final_y - iy) > 1.5 or abs(final_x - ix) > 1.5:
return float(iy), float(ix)
return final_y, final_x
except Exception:
return float(iy), float(ix)
# --- 7. スプライン補間 (Bicubic Spline) ---
def find_peak_spline(response_map):
(iy, ix), is_edge = get_peak_integer(response_map)
if is_edge:
return float(iy), float(ix)
y_range = np.arange(iy - 1, iy + 2)
x_range = np.arange(ix - 1, ix + 2)
patch = response_map[iy-1:iy+2, ix-1:ix+2]
try:
interp_func = RectBivariateSpline(y_range, x_range, patch, kx=3, ky=3)
bounds = ((iy - 0.5, iy + 0.5), (ix - 0.5, ix + 0.5))
result = minimize(
lambda p: -interp_func(p[0], p[1], grid=False),
x0=[iy, ix],
method='L-BFGS-B',
bounds=bounds
)
if result.success:
y_sub, x_sub = result.x
if (iy - 1.0 < y_sub < iy + 1.0) and (ix - 1.0 < x_sub < ix + 1.0):
return y_sub, x_sub
return float(iy), float(ix)
except Exception:
return float(iy), float(ix)
実験結果を以下に示す。
上段が推定誤差の散布図、下段がその絶対値のヒストグラムである。
2次元多項式フィッティングとガウシアンフィッティング系は比較的安定して高い精度が出ている。
問題の2DSGフィルタを使用した方法であるが、他の手法と比較して文字通り桁違いの精度を達成している。
その精度は今回の実験で平均0.007px、標準偏差0.004pxであった。
画像加工プログラム
from scipy.ndimage import map_coordinates
def crop_random_patch(img, patch_size=30):
h, w = img.shape
cy = np.random.uniform(patch_size, h - patch_size)
cx = np.random.uniform(patch_size, w - patch_size)
yy, xx = np.mgrid[0:patch_size, 0:patch_size]
yy = yy + cy - patch_size/2
xx = xx + cx - patch_size/2
coords = np.vstack([yy.ravel(), xx.ravel()])
patch = map_coordinates(img, coords, order=3, mode='reflect').reshape(patch_size, patch_size)
return patch, (cy, cx)
patch, center = crop_random_patch(img, patch_size=64)
response = cv2.matchTemplate(img, patch, cv2.TM_CCORR_NORMED)
テスト２ 位相限定相関法
- 512x512pxのBarbara（0~1に正規化）に縦横それぞれ±10px以下のランダムな非整数平行移動を加え、更に標準偏差0.1のノイズを加えた画像を用意し、元画像と加工画像の位相限定相関を計算する。
- 得られた相関マップを入力として以下の各手法でピーク座標をサブピクセル推定する。
- 以上を1000回繰り返す。（各手法には同じ画像ペアを使用）
- おまけで
cv2.findTransformECC（輝度ベースの変位量推定）による推定とcv2.phaseCorrelate（位相限定相関法のOpenCV実装）も比較対象として計算する。
先程は今ひとつだった重心ベースの方法が良い成績を示しているが、こちらも2DSGが圧倒的に高い精度を示している。
推定誤差は平均0.016px、標準偏差0.008pxであった。
また、輝度ベースのECC（Enhanced Correlation Coefficient12、これのみ相関マップでなく画像を渡しているためヒントも多い）やOpenCVのphaseCorrelate関数と比較しても大幅に良い成績である。
画像加工プログラム
def phase_only_correlation(img1: np.ndarray, img2: np.ndarray) -> np.ndarray:
F1 = np.fft.fft2(img1.astype(np.float64))
F2 = np.fft.fft2(img2.astype(np.float64))
eps = np.finfo(np.float64).eps
cross_power = F1 * np.conj(F2)
cross_power /= np.abs(cross_power) + eps
response = np.fft.ifft2(cross_power)
response = np.fft.fftshift(np.real(response))
response /= np.max(response) + eps
return response
shift_x, shift_y = np.random.uniform(-10, 10, 2)
M = np.float32([[1, 0, shift_x], [0, 1, shift_y]])
im_shifted = cv2.warpAffine(img, M, (img.shape[1], img.shape[0]), flags=cv2.INTER_CUBIC) + np.random.normal(0, .05, img.shape)
response = phase_only_correlation(img, img_shifted)
おわりに
2DSGフィルタによるニュートン法で画像ピークのサブピクセル座標を推定すると、どういうわけかとんでもない精度を達成できる衝撃は伝わっただろうか？
ノイズを平滑化しつつ、広域の形状を考慮しつつ、でも局所的な形状にも鈍感にならない特徴が上手くハマったのだと思うが、もし数学的に解明できる方がいたら教えていただきたい。
精度の高さ・分散の小ささは製造業の位置合わせタスク等において極めて重要な性質であり、そのような分野の開発の一助となると嬉しい3。
-
Georgios D Evangelidis and Emmanouil Z Psarakis. Parametric image alignment using enhanced correlation coefficient maximization. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 30(10):1858–1865, 2008. ↩
-
https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga1aa357007eaec11e9ed03500ecbcbe47 ↩
-
勾配に応じて反復計算を行っているため、データによっては稀に変なところに飛んでいくことがある。大きなズレが許されないタスクでは、その点は注意して対策を打って活用して欲しい。 ↩