Transpose

transposeとは転置行列、つまり行列の行と列をひっくり返して新しい行列をつくることです。
画像をひとつの行列として見る場合、これをtransposeすると以下のようになります。

transeposeしたものに更にFlipHorizontal()をかければTurnRightになりますし、

FlipVerticalをかければTurnLeftになります。

VapourSynthにstd.Transpose()があるのに、TurnLeft/TurnRightがないのはこのためです。

さて、このtransposeは画像処理においては意外に重要な役割を果たします。
transposeすると、縦と横がひっくり返るため、フィルタ処理において縦と横でそれぞれ同じような処理をする場合(resizerとかconvolutionとか)、一方向の処理だけ書けば済むのです。
特にSSEとかCUDAとかでベクトル演算を利用するならば縦方向処理のほうが速いことが多いので、transposeの重要性は高まります。

というわけで、これは画像処理を嗜むものであれば自前でそこそこ速いやつを書けなければならんと思い、書いてみることにしました。

書いてみてわかったこと

1. SSE2を使う場合、8x8単位で処理するよりも16x16単位で処理するほうが速い
xmmレジスタの数的に32bit(8個)だと16x16はいろいろとペナルティが大きそうなのですが、それでも16x16のほうが速いようです。
もし64bitであれば、xmmレジスタの数は16個に増えるので、さらに有利になりそうな気がします。

2. VC++10は俺の書いたコードをろくにloop unrollしてくれない
このコードを書いた当初は配列とforループを使ったもうちょっとスッキリしたコードでした。でも、コンパイラがループアンロールしてくれないの...orz
手動アンロールしただけで数十fpsスピードが上がったのを目の当たりにしたときはもうなんというか力が抜けてしまいました。
こんな糞コンパイラが世にはびこっているせいでC99の普及が遅れたなんて、もうなんというかね...

Alpha Max Plus Beta Min Algorithm

画像処理フィルタを書いていると頻出する計算に次のようなものがあります。

X = sqrt(A * A + B * B)

要するに二つの数の二乗の和の平方根を求めるわけです。
普通にC等で書くならば特に問題ないわけですが、いざSIMD化しようとなるととたんに面倒になります。

まず整数演算の場合、掛け算は16bit同士からしかできません。
よって8bit同士の掛け算をしたい場合は、両者をまず16bitに変換しなければなりません。

__m128i xmm0 = _mm_setr_epi8(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);
__m128i zero = _mm_setzero_si128()
__m128i xmm1 = _mm_unpackhi_epi16(xmm0, zero);
xmm0 = _mm_unpacklo_epi16(xmm0, zero);

更に、16bit同士の掛け算は最大で32bitになる可能性があるため(65,535 x 65,535 = 4,294,836,225)、出力は32bitになります。
しかも素直に32bitで出力してくれるような命令は存在せず、上位16bitと下位16bitがわかれた状態で出力されるか、積和計算を利用するかしなければなりません。

__m128i xmm2 = _mm_mullo_epi16(xmm0, xmm0);
__m128i xmm3 = _mm_mulhi_epi16(xmm0, xmm0);
xmm0 = _mm_unpacklo_epi16(xmm2, xmm3);
xmm2 = _mm_unpackhi_epi16(xmm2m xmm3);

xmm3 = _mm_unpackhi_epi16(xmm1, zero);
xmm1 = _mm_unpacklo_epi16(xmm1, zero);
xmm1 = _mm_madd_epi16(xmm1, xmm1);
xmm3 = _mm_madd_epi16(xmm3, xmm3);

ただ二乗するだけで、最初は128bitのXMMレジスタ一本で済んでいたのがいつのまにやら4本まで膨れあがり、コードもすでに12行です。
この操作をもう一つの数に対しても行い全部で8本になったレジスタを足しあわせて再び4本のレジスタになったあと、今度はこれらの平方根を求めなければなりません。
しかし、平方根を求める命令は整数演算にはないので、32bit整数を単精度浮動小数点に変換してから行うことになります。

__m128 flt0 = _mm_cvtepi32_ps(xmm0)
flt0 = _mm_sqrt_ps(flt0)
xmm0 = _mm_cvtps_epi32(flt0)

これを4回行った後、最後に再び8bit整数に戻します。

xmm0 = _mm_packs_epi32(xmm0, xmm1);
xmm1 = _mm_packs_epi32(xmm2, xmm3);
xmm0 = _mm_packus_epi16(xmm0, xmm1);

なんでたかがsqrt(A*A+B*B)を行うためにこんな辛い目にあわねばならないのでしょう。
しかもこれで計算できるのは一度にたったの16個です(SSE/SSE2の場合)。
普通にCで書いて16回計算するのとスピードは変わらないか下手すると遅くなってしまいます。
例えばTEdgeMaskはtriticalプラグインでは珍しくクリティカルパートがSIMD化されていないのですが、おそらくはこの計算が出てくるためにやらなかったのでしょう。

というわけで、もしどうしてもスピードがほしいなら普通はこんなに馬鹿正直な計算はしません。

sqrt(A*A+B*B) -> abs(A)+abs(B)

正直、手抜きにも程があるとは思います。
しかし、実用上これで問題が出るケースは意外と少ないので、これで済ませてしまったりするのですね。

とまあ、コストとパフォーマンスの兼ね合いという人類普遍のテーマについてつらつらと考えていたところ、prunedtree氏が"alpha max plus beta min algorithm"を教えてくれました。

sqrt(A*A+B*B) -> max(abs(A),abs(B))*alpha+min(abs(A),abs(B))*beta

alphaとbetaの選び方によっては最大誤差3.96%、平均誤差1.3%……これはイイ！
流石はgradfun2dbやFRfunシリーズを書いた人です。こういうテクニックをよく知ってらっしゃる。
抜け道というものは探せばあるものなんだなぁと思いました。


ちなみにprunedtree氏は、昔はmarcFDという名前でした。
なんでも若気の至りによって犯した過ち(DeenやaWarpSharpのことだと思われる)を悔いて改名したそうです。

Hysteresis mask

Doom9のほうで「GenericFiltersにmt_hysteresisの同等品を追加してくれないか？あれはhandyでいいものだよ」というリクエストがありました。

mt_histeresisねぇ…そりゃ、あんたにとってはhandyかもしれんが、書く方の身にもなれよ。
あれはIIRフィルタだからSIMD化も出来ないし、遅いものしか書けないんじゃないかなぁ。
だいたいあれってgenericなfilterじゃないだろ、specificなfilterじゃねえか。

てなことをブツブツ言いながらもとりあえず書いてみることにしたわけです。

ちなみにmt_hysteresisは、二値化する際の閾値を変えた二つのエッジマスクからノイズの少ない一つのエッジマスクを作り出すフィルタです。

まずはアルゴリズムやロジックを理解するため、avisynthで同等のものが書けるかどうかを試してみます。
vapoursynthプラグインは9/10/16bit対応とか可変解像度/フォーマットとかのことも考えないといけないので、avisynthよりも面倒なのです。
で、Masktools1の方を参考に(Masktools2のコードはmanao氏の趣味なのかテンプレートやSTL使いまくりのメタメタコードなのでC++がよくわからない自分には理解できない)書いてみたわけですが…

再帰を使うとstack overflowを起こすし、かといってSTLの使い方もよくわからんので自分でスタックを実装するはめになりました(Masktools2はSTLのlistとpairを使っている)。
結果として解像度に縦横ともに65535までの制限ができたりメモリの使用量が多くなったり(1920x1080で約8MB)しましたが、まあそこらへんは現時点ではそれほど問題ではないと目をつぶることにします。

さて、このプラグイン、出力はmt_hysteresisと変わらないわけですが最大の問題はスピードです。
さっそくベンチマークだ！

#benchmark.avs
MPEG2Source("1440x1080_6360frames.d2v").ConvertToY8()
base = last.mt_edge(thY1=30, thY2=30)
alt = last.mt_edge(thY1=5, thY2=5)
ret= mt_hysteresis(base, alt) # Masktools2a48
# ret = Hysteresis(base, alt) # 今回の自作フィルタ
return ret

$ for i in {1..3}; do avs2pipemod -benchmark benchmark.avs; done

      mt_hysteresis  Hysteresis
1st     24.255fps     48.490fps
2nd     24.242fps     48.523fps
3rd     24.310fps     48.583fps
--------------------------------
avg     24.269fps     48.525fps

まさかのダブルスコアです。

原因がMSVC++のSTL実装が糞すぎるのか(Masktools2の配布バイナリはVC++10で筆者と同じ)、それともこんなクリティカルな処理にSTLを使うのが間違いなのかはわかりませんが、C++大好き人間は御託を並べる前にもう少し基本に立ち返るべきなのではなかろーかとか思いました。
世はC++11サイコーとか騒がしいですが、肝腎な部分でラクなものに走れば大事なものを失うことになりかねません。
つーか、Masktools2遅すぎるだろ。

さて、次はvapoursynth版書かなきゃならんのか…。

追記：
TurboPascal7氏よりmt_hysteresisが遅い件について指摘を受けました。
彼曰く、「あれはmanaoがlistを使ってるから遅いんだよ、vector使えば2、3倍速くなるのは確認済みだよ」
なるほど、vectorですか。失礼しました。
俺もmallocで一度に最大分確保するのやめてrealloc使うことにしよ。