Transpose

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

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

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

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

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

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

	#include <stdint.h>
	#include <malloc.h>
	#include <windows.h>
	#include <emmintrin.h>
	#include "avisynth.h"
	static const AVS_Linkage* AVS_linkage = 0;
	static __forceinline __m128i mm_movehl_si128(const __m128i& xmm0)
	{
	return _mm_castps_si128(_mm_movehl_ps(_mm_castsi128_ps(xmm0),
	_mm_castsi128_ps(xmm0)));
	}
	/* MSVC++10 does not loop unroll to my code...FUCK! */
	static void __stdcall
	proc_8x8(const uint8_t* srcp, uint8_t* dstp, int width, int height, int src_pitch,
	int dst_pitch)
	{
	int hmod8 = (height / 8) * 8;
	int sp = src_pitch / 16;
	int dp = dst_pitch / 16;
	for (int y = 0; y < hmod8; y += 8) {
	__m128i* d = (__m128i*)(dstp + y);
	for (int x = 0; x < width; x += 8) {
	const __m128i* s = (__m128i*)(srcp + x);
	__m128i xm0, xm1, xm2, xm3, xm4, xm5, xm6, xm7;
	xm0 = _mm_loadl_epi64(s);
	xm1 = _mm_loadl_epi64(s + sp);
	xm2 = _mm_loadl_epi64(s + sp * 2);
	xm3 = _mm_loadl_epi64(s + sp * 3);
	xm4 = _mm_loadl_epi64(s + sp * 4);
	xm5 = _mm_loadl_epi64(s + sp * 5);
	xm6 = _mm_loadl_epi64(s + sp * 6);
	xm7 = _mm_loadl_epi64(s + sp * 7);
	xm0 = _mm_unpacklo_epi8(xm0, xm1);
	xm1 = _mm_unpacklo_epi8(xm2, xm3);
	xm2 = _mm_unpacklo_epi8(xm4, xm5);
	xm3 = _mm_unpacklo_epi8(xm6, xm7);
	xm4 = _mm_unpacklo_epi16(xm0, xm1);
	xm5 = _mm_unpackhi_epi16(xm0, xm1);
	xm6 = _mm_unpacklo_epi16(xm2, xm3);
	xm7 = _mm_unpackhi_epi16(xm2, xm3);
	xm0 = _mm_unpacklo_epi32(xm4, xm6);
	xm1 = _mm_unpackhi_epi32(xm4, xm6);
	xm2 = _mm_unpacklo_epi32(xm5, xm7);
	xm3 = _mm_unpackhi_epi32(xm5, xm7);
	xm4 = mm_movehl_si128(xm0);
	xm5 = mm_movehl_si128(xm1);
	xm6 = mm_movehl_si128(xm2);
	xm7 = mm_movehl_si128(xm3);
	_mm_storel_epi64(d, xm0);
	d += dp;
	_mm_storel_epi64(d, xm4);
	d += dp;
	_mm_storel_epi64(d, xm1);
	d += dp;
	_mm_storel_epi64(d, xm5);
	d += dp;
	_mm_storel_epi64(d, xm2);
	d += dp;
	_mm_storel_epi64(d, xm6);
	d += dp;
	_mm_storel_epi64(d, xm3);
	d += dp;
	_mm_storel_epi64(d, xm7);
	d += dp;
	}
	srcp += src_pitch * 8;
	}
	for (int y = hmod8; y < height; y++) {
	uint8_t* d = dstp + y;
	for (int x = 0; x < width; x++) {
	*d = srcp[x];
	d += dst_pitch;
	}
	srcp += src_pitch;
	}
	}
	static __forceinline void
	MM_TRANSPOSE4_SI128(__m128i& x0, __m128i& x1, __m128i& x2, __m128i& x3)
	{
	__m128 f0 = _mm_castsi128_ps(x0);
	__m128 f1 = _mm_castsi128_ps(x1);
	__m128 f2 = _mm_castsi128_ps(x2);
	__m128 f3 = _mm_castsi128_ps(x3);
	_MM_TRANSPOSE4_PS(f0, f1, f2, f3);
	x0 = _mm_castps_si128(f0);
	x1 = _mm_castps_si128(f1);
	x2 = _mm_castps_si128(f2);
	x3 = _mm_castps_si128(f3);
	}
	static __forceinline void mm_uplh_epi16(__m128i& x0, __m128i& x1)
	{
	__m128i t = _mm_unpackhi_epi16(x0, x1);
	x0 = _mm_unpacklo_epi16(x0, x1);
	x1 = t;
	}
	static __forceinline void mm_uplh_epi8(__m128i& x0, __m128i& x1)
	{
	__m128i t = _mm_unpackhi_epi8(x0, x1);
	x0 = _mm_unpacklo_epi8(x0, x1);
	x1 = t;
	}
	static void __stdcall
	proc_16x16(const uint8_t* srcp, uint8_t* dstp, int width, int height,
	int src_pitch, int dst_pitch)
	{
	int hmod16 = height / 16 * 16;
	int dp = dst_pitch / 16;
	int sp = src_pitch / 16;
	for (int y = 0; y < hmod16; y += 16) {
	__m128i* d = (__m128i*)(dstp + y);
	for (int x = 0; x < width; x += 16) {
	const __m128i* s = (__m128i*)(srcp + x);
	__m128i xm0, xm1, xm2, xm3, xm4, xm5, xm6, xm7,
	xm8, xm9, xm10, xm11, xm12, xm13, xm14, xm15;
	xm0 = _mm_load_si128(s);
	xm1 = _mm_load_si128(s + sp);
	xm2 = _mm_load_si128(s + sp * 2);
	xm3 = _mm_load_si128(s + sp * 3);
	xm4 = _mm_load_si128(s + sp * 4);
	xm5 = _mm_load_si128(s + sp * 5);
	xm6 = _mm_load_si128(s + sp * 6);
	xm7 = _mm_load_si128(s + sp * 7);
	xm8 = _mm_load_si128(s + sp * 8);
	xm9 = _mm_load_si128(s + sp * 9);
	xm10 = _mm_load_si128(s + sp * 10);
	xm11 = _mm_load_si128(s + sp * 11);
	xm12 = _mm_load_si128(s + sp * 12);
	xm13 = _mm_load_si128(s + sp * 13);
	xm14 = _mm_load_si128(s + sp * 14);
	xm15 = _mm_load_si128(s + sp * 15);
	mm_uplh_epi8(xm0, xm1);
	mm_uplh_epi8(xm2, xm3);
	mm_uplh_epi8(xm4, xm5);
	mm_uplh_epi8(xm6, xm7);
	mm_uplh_epi8(xm8, xm9);
	mm_uplh_epi8(xm10, xm11);
	mm_uplh_epi8(xm12, xm13);
	mm_uplh_epi8(xm14, xm15);
	mm_uplh_epi16(xm0, xm2);
	mm_uplh_epi16(xm1, xm3);
	mm_uplh_epi16(xm4, xm6);
	mm_uplh_epi16(xm5, xm7);
	mm_uplh_epi16(xm8, xm10);
	mm_uplh_epi16(xm9, xm11);
	mm_uplh_epi16(xm12, xm14);
	mm_uplh_epi16(xm13, xm15);
	MM_TRANSPOSE4_SI128(xm0, xm4, xm8, xm12);
	MM_TRANSPOSE4_SI128(xm1, xm5, xm9, xm13);
	MM_TRANSPOSE4_SI128(xm2, xm6, xm10, xm14);
	MM_TRANSPOSE4_SI128(xm3, xm7, xm11, xm15);
	_mm_store_si128(d, xm0);
	d += dp;
	_mm_store_si128(d, xm4);
	d += dp;
	_mm_store_si128(d, xm8);
	d += dp;
	_mm_store_si128(d, xm12);
	d += dp;
	_mm_store_si128(d, xm2);
	d += dp;
	_mm_store_si128(d, xm6);
	d += dp;
	_mm_store_si128(d, xm10);
	d += dp;
	_mm_store_si128(d, xm14);
	d += dp;
	_mm_store_si128(d, xm1);
	d += dp;
	_mm_store_si128(d, xm5);
	d += dp;
	_mm_store_si128(d, xm9);
	d += dp;
	_mm_store_si128(d, xm13);
	d += dp;
	_mm_store_si128(d, xm3);
	d += dp;
	_mm_store_si128(d, xm7);
	d += dp;
	_mm_store_si128(d, xm11);
	d += dp;
	_mm_store_si128(d, xm15);
	d += dp;
	}
	srcp += src_pitch * 16;
	}
	for (int y = hmod16; y < height; y++) {
	uint8_t* d = dstp + y;
	for (int x = 0; x < width; x++) {
	*d = srcp[x];
	d += dst_pitch;
	}
	srcp += src_pitch;
	}
	}
	class Transpose: public GenericVideoFilter {
	int bwidth;
	uint8_t *buff;
	int num_planes;
	void (_stdcall *proc)(const uint8_t* srcp, uint8_t* dstp, int width,
	int height, int src_pitch, int dst_pitch);
	public:
	Transpose(PClip child, int pf, IScriptEnvironment* env);
	~Transpose();
	PVideoFrame __stdcall GetFrame(int n, IScriptEnvironment* env);
	};
	Transpose::Transpose(PClip c, int pf, IScriptEnvironment* env)
	: GenericVideoFilter(c), num_planes(3)
	{
	if (!vi.IsPlanar() || vi.IsYV16() || vi.IsYV411()) {
	env->ThrowError("Transpose: unsupported format.");
	}
	if (vi.IsY8()) {
	num_planes = 1;
	}
	int t = vi.width;
	vi.width = vi.height;
	vi.height = t;
	bwidth = ((vi.width + 15) / 16) * 16;
	int bheight = ((vi.height + 15) / 16) * 16;
	buff = (uint8_t *)_aligned_malloc(bwidth * bheight, 16);
	if (!buff) {
	env->ThrowError("Transpose: failed to allocate buffer.");
	}
	proc = pf == 0 ? proc_16x16 : proc_8x8;
	}
	Transpose::~Transpose()
	{
	_aligned_free(buff);
	}
	PVideoFrame __stdcall Transpose::GetFrame(int n, IScriptEnvironment* env)
	{
	const int planes[] = {PLANAR_Y, PLANAR_U, PLANAR_V};
	PVideoFrame src = child->GetFrame(n, env);
	PVideoFrame dst = env->NewVideoFrame(vi);
	for (int i = 0; i < num_planes; i++) {
	int p = planes[i];
	const uint8_t* srcp = src->GetReadPtr(p);
	if ((intptr_t)srcp & 15) {
	env->ThrowError("Transepose: invalid memory alignment.");
	}
	proc(srcp, buff, src->GetRowSize(p), src->GetHeight(p),
	src->GetPitch(p), bwidth);
	env->BitBlt(dst->GetWritePtr(p), dst->GetPitch(p), buff, bwidth,
	dst->GetRowSize(p), dst->GetHeight(p));
	}
	return dst;
	}
	static AVSValue __cdecl
	create_transpose(AVSValue args, void* user_data, IScriptEnvironment* env)
	{
	if (!(env->GetCPUFlags() & CPUF_SSE2)) {
	env->ThrowError("Transpose: requires SSE2 capable CPU.");
	}
	int proc = args[1].AsInt(0); // 0:16x16 1:8x8
	if (proc < 0 || proc > 1) {
	env->ThrowError("Transpose: proc must be set to 0 or 1.");
	}
	return new Transpose(args[0].AsClip(), proc, env);
	}
	extern "C" __declspec(dllexport) const char* __stdcall
	AvisynthPluginInit3(IScriptEnvironment* env, const AVS_Linkage* const vectors)
	{
	AVS_linkage = vectors;
	env->AddFunction("Transpose", "c[proc]i", create_transpose, 0);
	return "TRANSPOSE";
	}

view raw transpose.cpp hosted with ❤ by GitHub

書いてみてわかったこと

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++がよくわからない自分には理解できない)書いてみたわけですが…

	#include <windows.h>
	#include "avisynth.h"
	static const AVS_Linkage* AVS_linkage = 0;
	class Footprint {
	int index;
	int width;
	int height;
	BYTE* map;
	DWORD* xy;
	const BYTE* altp;
	BYTE* dstp
	int a_pitch;
	int d_pitch;
	public:
	Footprint(int w, int h, const BYTE* a, BYTE* d, int ap, int dp,
	IScriptEnvironment* env) : width(w), height(h), altp(a),
	dstp(d), a_pitch(ap, d_pitch(d),index(-1) {
	map = (BYTE*)calloc(width * height, 1);
	xy = (DWORD*)malloc(width * height * sizeof(DWORD));
	if (!map || !xy) {
	env->ThrowError("failed to allocate Footprint");
	}
	}
	~Footprint() {
	free(map);
	free(xy);
	map = 0;
	xy = 0;
	}
	int get_width() {
	return width;
	}
	int get_height() {
	return height;
	}
	void push(int x, int y) {
	dstp[x + y * a_pitch] = altp[x + y * d_pitch];
	map[x + y * width] = 0xFF;
	xy[++index] = ((WORD)x << 16) | (WORD)y;
	}
	void pop(int* x, int* y) {
	*x = xy[index] >> 16;
	*y = xy[index--] & 0x0000FFFF;
	}
	bool is_empty() {
	return index == -1;
	}
	int passed(int x, int y) {
	return (int)map[x + y * width];
	}
	};
	class Hysteresis : public GenericVideoFilter {
	PClip child2;
	bool chroma;
	int num_planes;
	public:
	Hysteresis(PClip c, PClip c2, bool chroma, IScriptEnvironment* env);
	~Hysteresis() {}
	PVideoFrame __stdcall GetFrame(int n, IScriptEnvironment* env);
	};
	Hysteresis::Hysteresis(PClip c, PClip c2, bool ch, IScriptEnvironment* env)
	: GenericVideoFilter(c), child2(c2), chroma(ch)
	{
	if (!vi.IsPlanar()) {
	env->ThrowError("not planar format");
	}
	const VideoInfo& vi2 = child2->GetVideoInfo();
	if (vi.IsSameColorspace(vi2) == false) {
	env->ThrowError("not the same csp");
	}
	if (vi.width != vi2.width || vi.height != vi2.height) {
	env->ThrowError("not the same resolution");
	}
	num_planes = vi.IsY8() ? 1 : 3;
	}
	static void proc(Footprint* fp, int x, int y, const BYTE* altp, int alt_pitch)
	{
	fp->push(x, y);
	int posx, posy;
	while (!fp->is_empty()) {
	fp->pop(&posx, &posy);
	int x_min = posx > 0 ? posx - 1 : 0;
	int x_max = posx < fp->get_width() - 1 ? posx + 1 : posx;
	int y_min = posy > 0 ? posy - 1 : 0;
	int y_max = posy < fp->get_height() - 1 ? posy + 1 : posy;
	for (int yy = y_min; yy <= y_max; yy++) {
	for (int xx = x_min; xx <= x_max; xx++) {
	if (altp[xx + yy * alt_pitch] && !fp->passed(xx, yy)) {
	fp->push(xx, yy);
	}
	}
	}
	}
	}
	PVideoFrame __stdcall Hysteresis::GetFrame(int n, IScriptEnvironment* env)
	{
	const int planes[] = {PLANAR_Y, PLANAR_U, PLANAR_V};
	PVideoFrame base = child->GetFrame(n, env);
	PVideoFrame alt = child2->GetFrame(n, env);
	PVideoFrame dst = env->NewVideoFrame(vi);
	for (int i = 0; i < num_planes; i++) {
	if (i != 0 && !chroma) {
	return dst;
	}
	int p = planes[i];
	int width = base->GetRowSize(p);
	int height = base->GetHeight(p);
	int base_pitch = base->GetPitch(p);
	int alt_pitch = alt->GetPitch(p);
	int dst_pitch = dst->GetPitch(p);
	const BYTE* basep = base->GetReadPtr(p);
	const BYTE* altp = alt->GetReadPtr(p);
	BYTE* dstp = dst->GetWritePtr(p);
	Footprint *fp = new Footprint(width, height, altp, dstp, alt_pitch,
	dst_pitch,env);
	memset(dstp, 0, dst_pitch * height);
	for (int y = 0; y < height; y++) {
	for (int x = 0; x < width; x++) {
	if (basep[x + y * base_pitch] && altp[x + y * alt_pitch] &&
	!fp->passed(x, y)) {
	proc(fp, x, y, altp, alt_pitch);
	}
	}
	}
	delete fp;
	}
	return dst;
	}
	AVSValue __cdecl
	create_hysteresis(AVSValue args, void* user_data, IScriptEnvironment* env)
	{
	return new Hysteresis(args[0].AsClip(), args[1].AsClip(),
	args[2].AsBool(false), env);
	}
	extern "C" __declspec(dllexport) const char* __stdcall
	AvisynthPluginInit3(IScriptEnvironment* env, const AVS_Linkage* const vectors)
	{
	AVS_linkage = vectors;
	env->AddFunction("Hysteresis", "cc[chroma]b", create_hysteresis, 0);
	return "hysteresis mask filter";
	}

view raw histeresis.cpp hosted with ❤ by GitHub

再帰を使うと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使うことにしよ。