Darkside（https対応しました）

　アクションカメラっぽい使い方での撮影を試してみた。手ぶれが極端に激しい場合に、ソフトウェア手ぶれ補正がどこまで機能するかを確認するためである。
　あれこれ問題が噴出しているのに撮影を進めなくても良さそうだが、撮るものは撮れるときに撮っておかねばならない。特に天候は、ままならない。
　計算時間が掛かり過ぎるのは問題だが、短時間のクリップになら現状でも適用可能である。ならば、適用して改良するのもまた意味はある。

　手ぶれが激しい映像に適用すると、すぐに問題が発覚した。まっとうに収束せず繰り返しループ内で無茶苦茶な推測値が続出したのだ。

　その余りに無茶苦茶な数字を眺めているうちに、気づいた。そうだ、異常値の判定は簡単じゃないか？
　ズレ量を推測し、その推測した値だけフレームを逆変換する。推測値が正確なら、逆変換後のフレームとのズレ量はゼロになる。ゼロにならなかったら、差を前の推測値に加算することで推測値をより正確にできる。これが基本アイデアだ。だが、推測はあくまで推測であり、不適切な推測値が出ることは珍しくない。どのような条件で異常値だと判定するかは、悩ましい問題だった。
　だが、推測値をフィードバックしてのクローズドループで収束を狙う場合、簡単だ。収束している場合、通常はループを回すごとにズレ量は小さくなる。１回前のループよりもズレ量が大きくなったら、不適切な推測が行われたと判定すれば良い。

　回転や拡大に関しては、異常な値を拝んだことがない。ほぼ問題はＸ軸Y軸の推測値だけで発生している。
　そこで、XズレとYズレの２乗加算を取り、それが１回前より大きくなっていたらループを打ち切ることにした。そして、１回前までの加算結果を採用する。
　ベンチマークに使っているスカイツリーのクリップと、今回アクションカメラ的に激しく手ぶれしているクリップと、その両方で試してみる。
　相変わらず膨大な時間を要するので、今日中に結果は出ない。その重い処理をしていて、GPU 使用率が０％のままというのが悲しい。

　だが、あれこれ手を広げ過ぎて手に負えなくなりつつあるのも確かだ。ここは、片付けられるものを１つずつ片付けよう。
　GPU 関係でいま最も解決に近いのが、ステッチである。位相限定相関法も極座標変換も登場せず、既に GPU により一部が高速化できている。numpy はなぜか遅いが、OpenCV は GPU 化で速くなっている。ステッチで残っているのは、倍率色収差低減部分の GPU 化である。RGB 配列の違いがあるので、うっかりすると別の色を補正してしまう。また、プログラムの記述自体も少し特殊らしい。
　完成一歩手前まで来ているが、慎重に作業せねばならない。

　静止画のステッチを使って動作確認し、これは何とか成功した。重要部分だけ↓

g_src = cv2.cuda_GpuMat()
g_dst = cv2.cuda_GpuMat()
im_bgr = cv2.cuda_GpuMat()
img_red = cv2.cuda_GpuMat((4096,5464), cv2.CV_8UC1)
im_dst = cv2.cuda_GpuMat((4096,5464), cv2.CV_8UC3)
　・・・
g_src.upload(img_right)
im_bgr = cv2.cuda.split(g_src)
M = cv2.getRotationMatrix2D((CRX, CRY), 0, MAG_RED)
img_red = cv2.cuda.warpAffine(im_bgr[1], M, (int(width/2), height), flags=WARP_FLAGS)
cv2.cuda.merge([im_bgr[0], img_red, im_bgr[2]], im_dst)
g_dst = cv2.cuda.remap(im_dst, cuxmapR, cuymapR, interpolation=WARP_FLAGS)
right = g_dst.download()

　ややこしいのはレイヤーごとの分割と合体である。特に merge は格納先 im_dst をサイズと型指定で用意しておかないといけない。

　numpy を cupy に書き直したが、不穏な警告が出る。

UserWarning: cuFFT plan cache is disabled on CUDA 11.1 due to a known bug, so performance may be degraded. The bug is fixed on CUDA 11.2+.
cache = get_plan_cache()

　をいをい、GPU 有効化した OpenCV は import 出来ずに苦しみまくり、ようやく少し古いバージョンの決め打ちで使えるようになったってのに！
　こともあろうか、CUDA 11.1 がバグでアウトだと！？

　最小二乗法の GPU 対応化がうまく行かず、その直前までの処理で実行時間を測ってみる。何と、正味の位相限定相関法だけで１．７５秒も掛かっている。cupy を使うことで、numpy より遅くなっている。もちろん OpenCV 関数との差は広がるばかり。これでは、最小二乗法を GPU 対応できたとしても無意味。
　新しい CUDA で OpenCV をビルドし直したら、速くなるのだろうか？
　それ以前の問題として、新しい CUDA でビルドし直した OpenCV は動くのだろうか？

　ひたすら前途多難である。なぜ CUDA に phaseCorrelate 関数が無いんだよ！
　現状は、phaseCorrelate 関数が余りに遅いので GPU を有効にしようと頑張ったけど、意味ありませんでした・・・である。いや、ここは冷静に考えてみよう。余りに酷い目に遭いまくってるせいでキレちゃったが、さすがにおかしい。位相限定相関法の正味の処理は、それほど長いコードではない。

img_f1 = np.fft.fft2(img_g1)
img_f2 = np.fft.fft2(img_g2).conjugate()
img_f = img_f1 * img_f2
img_n = img_f / np.abs(img_f)
img_p = np.fft.ifft2(img_n).real
img_p = np.fft.fftshift(img_p)

　この部分を cupy に置換したら実行時間が１倍半を超えるなんてのは、幾ら何でもおかし過ぎる。かなり大きな行列であり、GPU への転送オーバーヘッドを考慮しても変だ。ここは CUDA11.1 のバグをまずは疑うべきだ。
　CUDA11.7 を使用するよう（そこだけ）パラメーターを変えて OpenCV をビルドし直す。エラー無くインストールできた。cupy も 11.7 対応を入れ直す。結果は・・・１．４７秒に高速化された！
　駄目じゃん、まだ numpy より遅いじゃん。

　位相限定相関法は非常に有用であることが判明しているものの、とにかく計算時間が長い。
　GPU を使えば大きな時間短縮効果が期待できるものの、GPU に対応したプログラムは容易に作れない。それが現状である。OpenCV だと GPU 非対応なら位相限定相関法は関数１つ呼ぶだけで完了するものの、GPU 対応版は無い。そこで、GPU 非対応 OpenCV にも関数が存在しなかった昔のプログラムから、サブピクセルを求めるための最小二乗法を流用させて貰う。scipy をインストールしなければならない。　

import numpy as np
import cv2
from scipy.optimize import leastsq
import time
# ステッチ中心
LX = 2114
LY = 2023
# 判定範囲 w や h の 1/2
LD2 = 1440
fitting_shape = (9, 9)
def pocfunc_model(alpha, delta1, delta2, r, u):
    N1, N2 = r.shape
    V1, V2 = list(map(lambda x: 2 * x + 1, u))
    return lambda n1, n2: alpha / (N1 * N2) * np.sin((n1 + delta1) * V1 / N1 * np.pi) * np.sin((n2 + delta2) * V2 / N2 * np.pi)\
                                            / (np.sin((n1 + delta1) * np.pi / N1) * np.sin((n2 + delta2) * np.pi / N2))
img1 = cv2.imread('103093.png')
img2 = cv2.imread('103094.png')
img_p1 = img1[int(LY-LD2):int(LY+LD2),int(LX-LD2):int(LX+LD2)]
img_p2 = img2[int(LY-LD2):int(LY+LD2),int(LX-LD2):int(LX+LD2)]
img_g1 = np.asarray(cv2.cvtColor(img_p1, cv2.COLOR_BGR2GRAY), 'float')
img_g2 = np.asarray(cv2.cvtColor(img_p2, cv2.COLOR_BGR2GRAY), 'float')
tt1 = time.perf_counter()
(x, y), e = cv2.phaseCorrelate(img_g1, img_g2)
tt2 = time.perf_counter()
print(tt2-tt1)
print(e, 'Y=' + str(y), 'x=' + str(x))
tt1 = time.perf_counter()
img_f1 = np.fft.fft2(img_g1)
img_f2 = np.fft.fft2(img_g2).conjugate()
img_f = img_f1 * img_f2
img_n = img_f / np.abs(img_f)
img_p = np.fft.ifft2(img_n).real
img_p = np.fft.fftshift(img_p)
m = np.floor(list(map(lambda x: x / 2.0, img_g1.shape)))
u = list(map(lambda x: x / 2.0, m))
max_pos = np.argmax(img_p)
peak = (int(max_pos / img_g1.shape[1]), max_pos % img_g1.shape[1])
max_peak = img_p[peak[0], peak[1]]
mf = np.floor(list(map(lambda x: x / 2.0, fitting_shape)))
fitting_area = img_p[int(peak[0] - mf[0]) : int(peak[0] + mf[0] + 1), int(peak[1] - mf[1]) : int(peak[1] + mf[1] + 1)]
cv2.imwrite('out.png', img_p, [int(cv2.IMWRITE_PNG_COMPRESSION), 1])
p0 = [0.5, -(peak[0] - m[0]) - 0.02, -(peak[1] - m[1]) - 0.02]
y, x = np.mgrid[-mf[0]:mf[0] + 1, -mf[1]:mf[1] + 1]
y = y + peak[0] - m[0]
x = x + peak[1] - m[1]
errorfunction = lambda p: np.ravel(pocfunc_model(p[0], p[1], p[2], img_p, u)(y, x) - fitting_area)
plsq = leastsq(errorfunction, p0)
tt2 = time.perf_counter()
print(tt2-tt1)
print(plsq[0][0], 'Y=' + str(plsq[0][1]), 'x=' + str(plsq[0][2]))

　まだ GPU は使用しておらず、OpenCV の関数一発と計算時間を比較。
　10年前の旧機だと、

新しいパソコンだと、

　いずれにしろ関数１つで済む OpenCV が３倍以上も速い。numpy を GPU 対応にできたとしても、この差を逆転できないようでは意味がない。ズレ量はほぼ同じだがサブピクセルの値は違う。これは、そもそもサブピクセルの推定をどうやっているのか？という問題があり今は保留である。
　クローズドループによる収束を試したとき、どちらが収束し易いか？が勝負となる。もちろんそれ以前に、速さを逆転できなければ全く意味は無くなる。

　試験用画像として、拡大を考慮しないと収束しなかったフレームを使用。

　１フレーム差の２枚について、位相限定相関法を使ってみる。

　オリジナルは 4096×4096 で、中心付近の 2880×2880 を切り出して使う。

　かなりすったもんだの末に、何とか相関グラフ画像らしきものが得られた。見える絵にするため、値を1000倍している。

　左上隅を (0,0) とした場合に(3,0) にピークが来ていて、たぶん正しい。

　だがこの瞬間に、思い切り厄介な問題に気付く。
　まず、ズレ量をサブピクセルまで求めようとすれば、どのような計算が良いのか？
　次に、符号。ズレ量がマイナスの場合は、どのような画像になるのか？

　試しに２枚の画像を入れ替えたところ、右上隅にピークが出現した。なるほど、0〜127 が正で 128〜255 が負、みたいな感じか。

　しかし、サブピクセルの求め方は何とかしなければならない。

　プログラム自体は単純だが、このままでは実用にならない。

import numpy as np
import cv2
# ステッチ中心
LX = 2114
LY = 2023
# 判定範囲 w や h の 1/2
LD2 = 1440
img1 = cv2.imread('103093.png')
img2 = cv2.imread('103094.png')
img_p1 = img1[int(LY-LD2):int(LY+LD2),int(LX-LD2):int(LX+LD2)]
img_p2 = img2[int(LY-LD2):int(LY+LD2),int(LX-LD2):int(LX+LD2)]
img_g1 = np.asarray(cv2.cvtColor(img_p1, cv2.COLOR_BGR2GRAY), 'float')
img_g2 = np.asarray(cv2.cvtColor(img_p2, cv2.COLOR_BGR2GRAY), 'float')
img_f1 = np.fft.fft2(img_g1)
img_f2 = np.fft.fft2(img_g2).conjugate()
img_f = img_f1 * img_f2
img_n = img_f / np.abs(img_f)
img_p = np.fft.ifft2(img_n).real
print(img_p)
img_p = img_p * 1000
cv2.imwrite('out.png', img_p, [int(cv2.IMWRITE_PNG_COMP/RESSION), 1])

　撮影後のソフトウェア手振れ補正で、９割近くのフレームは収束する。収束してくれないのは１割強である。
　その１割強に対して最終２回の平均を取る方式で、動画全体を作り直してみた。予想された通り、ほぼ変化無し。区別を付けるのは困難だ。なら、もうこれでいいか・・・と思っていたが、重大なことを試し忘れていたことに気付いた。それは、拡大率である。

　回転不変位相限定相関法では、ＸＹ平面のズレと回転だけでなく、拡大率も推定できる。しかし動画の手ぶれ補正においては、拡大率は明らかに１固定である。マクロ撮影でカメラが前後に動く場合などは拡大率まで追随させても良いが、余りに特殊用途過ぎる。
　かくして、拡大率という概念が頭から完全に抜け落ちていたのだった。

　だが、収束させる過程においては拡大率も考慮すべきではないか？
　最終的に手ぶれ補正処理を行う際には拡大率１固定でも良いが、クローズドループ中では拡大率も管理するのがベターでは？
　そう考えて、試してみた。どうしても収束してくれないフレームが、しっかり収束した。正解だった！

　こうなると、GPU による高速化はどうしても試したい。１分のクリップを処理するのに１日掛かるようでは、試行錯誤もままならない。実用性も、失われる。
　手ぶれ量推測処理は画像出力がなく、ほぼ画像処理だけである。それだけ GPU の恩恵も大きいはずだ。
　CUDA への移植を進めようとしたが、いきなりフーリエ変換に直面。そうだ、numpy の関数でフーリエ変換しているが、これ CUDA に出来るのか？
　調べると、CuPy を使えばできるらしい。numpy を GPU(CUDA)対応にしたものが、CuPy だ。

　Python + OpenCV による画像処理はお手軽だが、GPU を使って高速化を図ると突然面倒になる。と言うか情報収集に苦労するようになる。そりゃあ GPU 対応の有無で OpenCV インストールの難易度が何次元も違うからなぁ・・・

　すぐに行き詰まった。numpy は良いとして、phaseCorrelate や warpPolar に対応する CUDA 関数が見つからないのだ。

　折角なので、phaseCorrelate を使わず CUDA で（回転不変ではない）位相限定相関法のプログラムを作ることろから研究してみよう。長い時間を要しそうなので、その間に GPU 使わない版（拡大率も利用して収束させる）を動かしておこう。
　動かし始めてすぐに、収束しないフレームが発生。これ最後まで走らせてみないと、拡大を考慮した効果は分からない。