Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@cycloawaodorin

AviUtlの拡張編集用のLanczos3リサイズスクリプトを作ってみた

Last updated at Posted at 2019-12-12

#はじめに
Lanczos 関数による補間を用いた画像の拡大縮小アルゴリズムの存在については，聞いたことがある人も多いだろう．実際，ウェブには非常に多くの解説記事がある．以下はそのほんの一部の例．

AviUtl では，この方法で拡大・縮小を行うフィルタプラグインである Lanczos 3-lobed 拡大縮小 プラグインが非常に有名である．SIMD 命令を使用しているようで，比較的高速に処理することができる．

他方，これはフィルタプラグインであるため，拡張編集とともに用いるには，アルファチャンネルを持たないシーン全体に対して (フィルタオブジェクトの適用) しか使えない．シーン数が増大するとその管理が面倒で，中間ファイルを利用する場合，可逆圧縮形式の中間ファイルを生成するストレージ負荷が大きい．また，中間ファイルを生成するたびに不可逆圧縮を繰り返すのも避けたいところである．そこで，Lanczos 関数を用いた，拡張編集向けのリサイジングスクリプトを作ってみた，というのが本記事の内容である．

ちなみに拡張編集本家でもリサイズエフェクトや拡大率により拡縮は可能であるが，Lanczos 関数を使った方法ではない．以前拡大はベル関数，縮小は面積平均法を使っているとどこかで目にした記憶があるが，見た感じ拡張編集でも同様と考えられる (結果の比較等の検証はしていない)．

#使い方等
本記事は，スクリプトをどう作ったかの説明がメインであるが，使いたいだけという人のために，使い方も簡単にまとめておく．ただし，拡張編集そのものの使い方は説明しない．

##インストール
以下は例．わかっている人はお好きにどうぞ．

  1. AviUtl と拡張編集プラグインをインストールする (詳細は各自で)．
  2. KSAのリリース一覧 から最新版をダウンロード，展開し，好きな場所に配置する (またはgit cloneする)．
  3. [拡張編集のルート]/script以下に[KSAのルート]/scripts以下のファイルをコピーするか，[KSAのルート]/scriptsへの シンボリックリンクを生成する．

##使い方
拡張編集上でリサイズしたいオブジェクトにアニメーション効果をつけ，Lanczos3リサイズ@ksaを選択する．パラメータの意味は以下の通り．

ドット数でサイズ指定が未チェック
横を $(拡大率/100)\times (\mathrm{X}/100)$ 倍，縦を $(拡大率/100)\times (\mathrm{Y}/100)$ 倍する．
ドット数でサイズ指定がチェック済み
横を X，縦を Y のサイズにリサイズする．
設定 / スレッド数
並列処理に使うスレッド数を指定する．0 以下の値が指定された場合，$\max(\mbox{CPUスレッド数}+指定値,1)$ 個のスレッドを使う．
設定 / 上，下，左，右
リサイズ前に，指定した数値分クリッピングを行う．単位は常にピクセル．リサイズが比率指定の場合，クリッピング後のサイズに対して比率を適用する．

##使い所
本スクリプト は以下をすべて満たす場合に使うとよい．

  • アルファチャンネルを持つオブジェクトを拡縮したい．
  • アルファチャンネルが不要ならば Lanczos 3-lobed 拡大縮小 を勧める．
  • 本家の拡大率・リサイズの品質が不満である．
  • 拡大率・リサイズによる拡縮に比べて遅くなる，不満がなければ拡大率・リサイズをつかうとよい．
  • 動画を拡縮したい．
  • 静止画の拡縮は，事前に行うほうがよい．毎フレーム同じ拡縮を繰り返すのは無駄である．
  • 拡縮率を時間変化させたり，アニメーション後に拡縮したい場合は，本スクリプトは適している．

なお，クリッピングだけ行いたい場合は，本家のクリッピングエフェクトを使うこと．本スクリプトでは無駄に補間計算処理が入ってしまう．

#条件の整理
さて，Lanczos 関数による補間処理は重いため，Lua で直接書くわけには行かず，C/C++ による DLL を介した実装を用いることになる．拡張編集向けの DLL の作り方に関しては，先人が 入門記事 を書いているので，概要については割愛する．

なお，今回は DLL では C++ 言語を使用し，コンパイラは MSYS2 上の g++ を使用したので，以下ではそれを前提に話を進めていく．また，以下ではソースコードから重要部分を抜粋して議論するが，全体は GitHubのページ にあるので，必要に応じて参照されたい (ただしリサイズ以外のスクリプトも同梱されているので注意)．

##インターフェース部
まず，Lua 側，C++ 側のインターフェースをそれぞれ簡単に紹介しておこう．

Lua側
local ksa = require("ksa") -- モジュールのロード
local src, sw, sh = obj.getpixeldata("alloc") -- 入力画像の取得後の操作で破壊されないよう，"alloc" で取得する
obj.setoption("drawtarget", "tempbuffer", dw, dh) -- 出力用の領域(仮想バッファ)を生成
obj.copybuffer("obj", "tmp") -- copybuffer でオブジェクトのサイズを目標サイズに合わせる
local dest, dw, dh = obj.getpixeldata() -- 出力用の領域を取得
ksa.clip_resize(src, sw, sh, dest, dw, dh, ct, cb, cl, cr, n_th) -- DLL の関数を呼び出す(引数は説明しない)
obj.putpixeldata(dest) -- 結果を反映させる
Luaモジュール(ksa.lua)
local ksa = require("ksa_ext")

ksa.xxx = function(arg)
	-- Lua 言語でのモジュール関数を追加する
	-- これは今回の記事とは関係ないが，こうすると
	-- DLL と Lua 側のモジュールがまとめられてスッキリするなというご紹介
end

return ksa
C++側メイン
// Lua との連携用のヘッダ，予めダウンロードしておく
#include "lua/lua.hpp"

// 関数本体は別ファイルにして #include している
#include "ksa_ext.cpp"

// モジュール登録用の関数リストは ksa_ext.cpp からスクリプトで自動抽出している
// これも本題とは関係ないが，同じことを二度書くのが嫌いな人なので自動化
static constexpr luaL_Reg ksa_ext[] = {
#include "functions.c"
	{ nullptr, nullptr }
};

// DLL がロードされたときに呼び出される関数
// ここでは __declspec(dllexport) ではなく，
// 別に用意した def ファイルで export している
extern "C" {
int
luaopen_ksa_ext(lua_State *L)
{
	luaL_register(L, "ksa_ext", ksa_ext);
	return 1;
}
}
本体(ksa_ext.cpp)
namespace KSA { // 全体を名前空間で囲って内部リンケージを生成する

// 画素値用の構造体
using PIXEL_BGRA = struct pixel_bgra {
	unsigned char b;
	unsigned char g;
	unsigned char r;
	unsigned char a;
};

class ClipResize {
	// 処理本体に関わるクラス定義 (メンバ変数とコンストラクタ以外は省略)
public:
	const PIXEL_BGRA *src; // 入力側のデータは書き換えない
	PIXEL_BGRA *dest;
	XY x, y; // 縦軸・横軸で同じ処理を繰り返すので，それを格納するインナークラスで表現
};

// モジュールに登録する関数
static int
ksa_clip_resize(lua_State *L)
{
	// 引数受け取り
	std::unique_ptr<ClipResize> p(new ClipResize());
	int i=0;
	p->src = static_cast<PIXEL_BGRA *>(lua_touserdata(L, ++i));
	p->x.src_size = lua_tointeger(L, ++i);
	p->y.src_size = lua_tointeger(L, ++i);
	p->dest = static_cast<PIXEL_BGRA *>(lua_touserdata(L, ++i));
	// 以下省略
	
	// 本処理 (一旦省略)
	// p->dest を補間済みの画素値で書き換える処理が入る
	
	// (Lua側に) 返り値は返さない (出力領域は引数で受け取る)
	return 0;
}

};

先人の記事と概ね同様であるが，私見をいくつか述べておく．まず，引数受け取りのlua_toXXXの第 2 引数では，数値リテラルを直接書かず，一つの変数のインクリメントで処理している．引数の数や順番を変えたくなったとき，リテラルでは修正が面倒であるが，この方法なら行の追加・削除・入れ替えで対応できる．

C++ を (本格的に) 使うのが初めてだったこともあり，コンパイルオプションに少々戸惑ったので，参考までに以下に示す．C++ 的な要素を含む場合に，リンカをg++にすること，-static -lstdc++ -lgcc -lwinpthreadをつけることなどでハマった．

Makefile(一部)
CXX = g++
CXXFLAGS = -Wall -Wextra -Wdeprecated-declarations -Wpointer-arith -Wwrite-strings -Wmissing-noreturn -Wno-unused-parameter -Wsuggest-attribute=format -Wsuggest-attribute=noreturn -Wunused-variable -Winit-self -Wshadow -Wlogical-op -Wconversion -Wold-style-cast -c -O3
LD = g++
LDFLAGS = -shared -Wl,--dll,--enable-stdcall-fixup
LIBS = lua/lua51.dll -static -lstdc++ -lgcc -lwinpthread

フィルタプラグインと拡張編集スクリプト用 DLL の違い

いろいろと違いはあるが，最大の違いは色表現の方法だろう．拡張編集スクリプト用 DLL では，上述のアルファチャンネル付きの，画像ファイルとしておなじみの 32bit カラーである．

アルファチャンネル付き32bitカラー
using PIXEL_BGRA = struct pixel_bgra {
	unsigned char b;
	unsigned char g;
	unsigned char r;
	unsigned char a;
};

他方，フィルタプラグインでは基本的に AviUtl の内部表現に触れることになり，画素は下記のPIXEL_YC構造体である．詳細は AviUtlの内部形式について に詳しい．

YC構造体
typedef struct {
    short   y;          //  画素(輝度    )データ (     0 ～ 4096 )
    short   cb;         //  画素(色差(青))データ ( -2048 ～ 2048 )
    short   cr;         //  画素(色差(赤))データ ( -2048 ～ 2048 )
                        //  画素データは範囲外に出ていることがあります
                        //  また範囲内に収めなくてもかまいません
} PIXEL_YC;

RGB とはアフィン関係にあるため，重み付き線形和による補間では色空間の違いが本質的な違いにはならないと思うが，量子化深度が違う点に注意が必要である1．また，拡張編集ではアルファチャンネルが存在するが，本体側であるフィルタプラグインでは扱えない．拡張編集スクリプトではアルファチャンネル，つまり不透明度に注意する必要がある．

##不透明度付きの補間処理
補間時の不透明度の扱いだが，単純に追加の重みとして処理すればよい．ピクセル $i$ の各色を $b_i, g_i, r_i, a_i$，重みを $w_i$ として加重平均を求めるには，以下のようにすればよい．

\begin{align}
b &= \dfrac{\sum_iw_ia_ib_i}{\sum_iw_ia_i}, \\
g &= \dfrac{\sum_iw_ia_ig_i}{\sum_iw_ia_i}, \\
r &= \dfrac{\sum_iw_ia_ir_i}{\sum_iw_ia_i}, \\
a &= \dfrac{\sum_iw_ia_i}{\sum_iw_i}.
\end{align}

不透明度については，通常通り $w_i$ の重みを用いる．不透明度の分子が 3 色の分母になる点が面白い．なお，周辺が完全透明の場合，$a_i$ が常に $0$ になることから，$\sum_i w_ia_i$ が $0$ となり，ナイーブな実装だとゼロ除算が発生する．不透明度を適切に処理していない別のエフェクト等で変なことが起こらないようにするためには，$a=0$ であるピクセルの $b, g, r$ の値は $0$ であることが望ましい．

##端部処理
画像の端では，参照範囲の一部が元画像の外側に出てしまう．このような場合の処理として「外部は $0$ (黒，透明) として扱う」「鏡像があるものとして扱う」「ループしているものとして扱う」などがあるが，今回は「重みを含め，範囲外は参照しない」こととする．高度な並列処理を行う場合など，参照範囲が場合によって異なることが計算上面倒である場合もあるのだが，今回の実装では問題ではないため，この処理を用いることとした．なお，背景用などでループが前提となっている画像の場合はループ扱いで処理したほうがいいだろう．

#高速化
この章では，今回実装した 2 つの高速化について述べる．

##重みの事前計算
Lanczos 関数は sinc 関数を 2 回使うため，重みが必要になるたびに Lanczos 関数を呼び出すと，非常に計算が重くなる．しかし，画像の拡縮では，必要な重みの数はごく限られるため，予め必要な重みをすべて計算して配列に並べ，重みが必要なときはその配列を参照するだけにすれば，かなり計算量を削減できる．

では，実際に必要な重みというのはどれくらいになるのだろうか？まずは具体例として，8 ピクセルを 16 ピクセルに拡大する場合について考えてみる．
double.png
この図の通り，出力の 6 番の画素を計算するためには，入力の 0 番から 5 番を参照し，入力の 0 番の画素に掛ける重みは $\mathrm{Lanczos3(-2.75)}$ である．図に従い，距離 $-2.75$ から距離 $2.25$ までの 6 種類の重みが必要になる．同様に，出力の 7 番の画素を計算するためには $-2.25$ から $2.75$ までの 6 種類の重みが，8 番は $-2.75$ から $2.25$ までの重みが必要である．実質反転しているだけだからどの画素を計算するにも 6 種類の重みしか使わないのだが，反転を計算するのも面倒なので，簡単に 12 種類の重みが必要だと思ってすすめよう．お気づきだろうが，出力の偶数番目の画素を計算するには $-2.75$～$2.25$ のセットが，奇数番目には $-2.25$～$2.75$ のセットが必要になり，各 6 種類のセットが 2 セットで 12 種類の重みとなる．

この 2 セットの "2" とはどこから出てきた値だろうか． 2 倍の 2 だ，というのは無関係ではないが，完全ではない．1.5 倍のときに必要なのは，当然 1.5 セットではない．必要なセット数は $v=出力サイズ/\mathrm{gcd}(入力サイズ,,出力サイズ)$ である．ここで，$\mathrm{gcd}(a,,b)$ は $a$ と $b$ の最大公約数である．整数倍のときの最大公約数は入力サイズに等しくなるので，倍率がそのまま $v$ になる．少し考えればわかるので詳細は省くが，縮小の場合に必要なのも $v$ セットだ (参照範囲が増えるので，セットあたりの重みの種類数は拡大のときより多くなる)．

ところで，事前計算した重み配列から必要な重みを取り出すには，何番目の重みかを正確に計算できなければならない．拡大の場合は距離が $3$ 以内のものを，番号の小さい順に並べれば問題にならないが，縮小の場合，厳密の計算した場合に距離がちょうど $3$ になる画素が存在する場合があり，距離を浮動小数点数で計算して $3$ 以内かを判定すると，この距離がちょうど $3$ の画素を，距離が $3.0001$ に見えて使わないようにしようとしたり $2.9999$ に見えて使おうとしたりして，正しい重みを選択できない場合がある．実際最初の実装ではこのバグを入れてしまった2．

こうしたバグを防ぐためには，範囲内かどうかを整数演算で正確に計算する必要がある3．全部展開して直接の整数演算だけで済むようにアルゴリズムを整えてもよいが，せっかく C++ を使っているので，分数で表現する有理数クラスを作り，距離等を浮動小数点数から有理数に置き換えてバグを解消した．距離がちょうど $3$ の場合，重みは $0$ になるので，これは常に含まないようにした．実装はちょっとごちゃごちゃしているので ソースコード をあたっていただけると幸いである．

##並列化
実装したもう 1 つの高速化手段は，並列化である．このリサイズ処理は，互いの計算が影響しあわない，並列性の高いアルゴリズムのため，並列化は容易である．並列化に関係する部分だけを抜粋すると，以下の通りである．

ksa_ext.cppの並列化周り
#include <memory>
#include <thread>
namespace KSA {
class ClipResize {
	// 省略
public:
	static void
	invoke_interpolate(ClipResize *p, int y_start, int y_end)
	{
		for (int dy=y_start; dy<y_end; dy++) {
			for (int dx=0; dx<(p->x.dest_size); dx++) {
				// 出力側の (dx, dy) 座標の画素値を内挿計算する
				p->interpolate(dx, dy);
			}
		}
	}
};
static int
ksa_clip_resize(lua_State *L)
{
	// 引数受け取りとか重み計算とかを省略
	
	// スレッドの配列を確保
	std::unique_ptr<std::unique_ptr<std::thread>[]> threads(
		new std::unique_ptr<std::thread>[n_th] );
	// スレッド番号ごとに引数をズラして，それぞれのスレッドをたてる
	for (int t=0; t<n_th; t++) {
		threads[t].reset(new std::thread(ClipResize::invoke_interpolate, p.get(),
			( t*(p->y.dest_size) )/n_th, ( (t+1)*(p->y.dest_size) )/n_th));
	}
	// 全部のスレッドの終了を待ち合わせる
	for (int t=0; t<n_th; t++) {
		threads[t]->join();
	}
	return 0;
}
};

出力側の画素値を1点ずつClipResize::interpolate()によって計算して置き換えるのだが，その呼び出しループにおいて，出力画像の縦方向に単純にn_th分割して並列に走らせるという実装である．この実装においてp.get()で各スレッドに渡しているClipResizeオブジェクトは同一であるため，書き込みに衝突がある場合はスレッドセーフではないが，画素値の書き込みしかなく，書き込み先の座標がスレッドごとに異なるため，スレッドセーフ性は特に考えなくても問題ないという楽ちん実装である．

##高速化のあとがき
Lanczos 3-lobed 拡大縮小 ばりに SIMD 命令を使うとか，画像がでかいなら GPU を使うとかできればもっと速くなるだろうが，そういうのはよく知らないのでとりあえずこれで．上の 2 つの高速化を入れただけで「動画出力時にはちょっと重いけど，編集中は気にならない」程度には速くなったのでとりあえず満足．

(三角関数を使わない) Spline36 のほうが Lanczos より速い などと言われることもあるが，重みの事前計算を行えば，入力サイズと出力サイズが互いに素などの不幸なケース4でもない限り，重み計算そのものがネックになることはあまりない．

#おまけ：なんで sinc 関数なのか
この人 のように，sinc 関数をコンパクト化した Lanczos 関数を使うと「理論的にいい」らしいけど，なんで sinc 関数なのかわからない，という人が少なくない．実際，Lanczos 関数を使ったリサイズ処理を実装したという記事は数あれど，そこまで踏み込んだ説明はあまり見ない．あまり詳しくは述べないが，sinc 関数を使った補間が何を意味するのか，簡単に紹介しよう．キーワードだけ並べる感じなので，気になった人はキーワードを元に教科書を探して読んでみてほしい．また，以下の記述は間違いがあるかも知れないので，その場合は指摘していただけるとありがたい．

まず，sinc 関数がどこから出てきたのかというと，これは矩形関数の (逆)フーリエ変換によって得られる関数である．距離に応じた重みを付けて加重平均を取ることをスライドさせて繰り返す，という計算は，数学で畳み込みとよばれる演算である．したがって，フーリエ変換における畳込み定理により，「sinc 関数を畳み込む」計算は「フーリエ変換して，矩形関数を掛けて，フーリエ逆変換する」操作に等しい．ここで，「フーリエ変換して，矩形関数を掛ける」という操作は「周波数領域で矩形のローパスフィルタを掛ける」ことを意味する．このローパスフィルタのカットオフ周波数は，sinc 関数の幅に反比例する．

このため「拡大」における sinc 関数の畳み込みは，「サンプリング定理に基づいて『完璧に再現』した連続関数から，新しいサンプリングレートで再サンプリングする」操作に等しく，「縮小」における sinc 関数の畳み込みは，「サンプリング定理に基づいて『完璧に再現』した連続関数を，新たなサンプリングレートにおいて再現できない高周波の信号を除去し，再現可能な形状に変形してから再サンプリングする」操作に等しい．

このサンプリング定理に基づく前提は，「距離」をどう測るのが「正しい」かも規定する．斜めの「距離」はユークリッド距離で測り，重みは $\mathrm{sinc}\left(\sqrt{x^2+y^2}\right)$ とするのが「正しい」のでは，と考えてしまう人が時々いる (その疑問は尤もである) が，縦横それぞれの軸で独立にサンプリングしている，という立場からは，重みは $\mathrm{sinc}\left(x\right)\mathrm{sinc}\left(y\right)$ を用いるのが「正しい」のである5．距離を，サンプリングレートが低い側，つまり拡大なら入力側，縮小なら出力側の座標で測るのが「正しい」こともこの立場から説明できる．

ところで，sinc 関数が「なんの意味で理想的」で，現実の画像が「その意味で理想的なのか？」という点に注意しなければならない．元の画像が，ナイキスト周波数をカットオフ周波数とする矩形ローパスフィルタで処理して，(デルタ関数で) サンプリングしたものなのかというと，当然そんなことはない．「縮小」において，矩形のローパスフィルタを掛けるのが，「画質」の意味で理想的かというと，必ずしもそうではない．エイリアシングを防いだ上で，情報をできるだけ落とさない，という意味では理想的ではあるが，それが「人の目で見て最も自然」かどうかは定かではない．Lanczos の範囲を拡大して sinc 関数に近づけたら，却って変な画像になった，というのはこのあたりの「そもそもの前提が間違っている」ことの現れである．

ちなみに，sinc 関数のコンパクト近似において Lanczos 関数が何らかの意味で理想的なのかどうかについては私は知らない．

  1. つまり，複数のフィルタプラグインを掛けるのと同様の処理をすべて拡張編集上のスクリプトで行う場合，中間表現が「粗く」なるため，結果の精度が低下する原因となる．他方，1 回だけなら最終的な出力の量子化深度が変わるわけではないので，問題にならない． ↩

  2. 距離がちょうど $3$ になる場合が存在する倍率の縮小でないと現れないバグであり，発見に少々手間取った． ↩

  3. EPSILON を加えて誤差の影響を軽減する方法もあり，一時的に実装もしてみたが，EPSILON の大きさがどれくらいなら正確に判定できるかを計算するのが面倒なうえ，バグを埋め込まない自信がなかったので，整数演算に置き換えたほうが安心だと判断した． ↩

  4. その場合でも，1 画素あたり拡大で 36 点の画素を参照するが，重みの計算が 72 回 ( 縦 36 回 + 横 36 回 ) から 12 回 ( 縦 6 回 + 横 6 回 ) に減らせるので，事前計算は有用である．縮小の場合は参照範囲が増加するため，更に効果が大きい． ↩

  5. 独立でないと考えるならば，斜め方向の「正しいサンプリングレート」に基づいて sinc 関数の幅を決めなければならない．通常のデジタル画像が等間隔であるために，ランダム間隔である人間の目やフィルムカメラのフィルムに比べて「悪い」のもこの辺が関係している． ↩

0
3
7

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?