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