Lanczos（ランチョス法）【ついでにSpline36】

SSDを交換して容量不足から解消され、久々にVisualStudioを入れることができたので再開しま。
今回はランチョス（Lanczos）法と、おまけでSpline法をご紹介します。

Lanczos法の概要

バイリニアは2×2、バイキュービックでは4×4の範囲のピクセルをサンプルし、重み付け関数によりそれぞれに重みを付けを計算しました。
今回紹介するランチョス（Lanczos）では更に高度な計算で重みを付けます。
なお、本稿ではバイキュービックよりサンプル範囲を広くした6×6ピクセルを対象としたLanczos3を採用します。

ちなみに、Lanczos3はGIMPの画像拡大縮小でも選択可能なアルゴリズムとなっており、バイキュービックよりも高画質であるとされています。

重み付けの関数は以下のようになります。

lanczos_weight(d) = sinc(d) * sinc(d/n)  (|d| <= n)
                  = 0                    (|d| > n)
sinc(x) = sin(PI * x) / (PI * x)

下記サイトを参考にさせていただきました。

Lanczos関数による画像の拡大縮小

sinc関数は窓関数という働きをします。信号処理の分野だと結構当たり前に利用されるものですが説明が面倒なので省略！
詳しく知りたい方はWikipediaへどうぞ。

窓関数 - Wikipedia

ja.wikipedia.org

そして、ここで出てくる*nですが、これはサンプルの範囲を示すものです。
Lanczos3の『3』というのはこのnのことでつまりは、カレントのピクセルから-3～+3の距離にあるピクセルをサンプルして重み付けを行うということです。
なのでnは2や4などでもできますが、クオリティとコストのバランスが良いのが3らしいです。

プログラム実装

では、バイキュービック編2回目で用いたソースを一部抜粋して、多少見やすくしたものに対してアルゴリズムの修正をしたのが以下のコードになります。

static double sinc(double x)
{
    return Math.Sin(x * Math.PI) / (x * Math.PI); 
}

static double lanczosWeight(double d, int n = 3)
{
    return d == 0 ? 1 : (Math.Abs(d) < n ? sinc(d) * sinc(d / n) : 0);
}

static Bitmap resizeLanczos(Bitmap bmp, int width, int height, int n = 3)
{
    int sw = bmp.Width, sh = bmp.Height;
    double wf = (double)sw / width, hf = (double)sh / height;

    var dst = new Bitmap(width, height, PixelFormat.Format32bppRgb);
    var bmpData = bmp.LockBits(new Rectangle(Point.Empty, bmp.Size), ImageLockMode.ReadOnly, PixelFormat.Format32bppRgb);
    var dstData = dst.LockBits(new Rectangle(Point.Empty, dst.Size), ImageLockMode.WriteOnly, PixelFormat.Format32bppRgb);

    unsafe
    {
        var bmpScan0 = (byte*)bmpData.Scan0;
        var dstScan0 = (byte*)dstData.Scan0;

        Parallel.For(0, height, iy =>
        {
            for (var ix = 0; ix < width; ix++)
            {
                var wfx = wf * ix;
                var wfy = hf * iy;
                var x = (int)wfx;
                var y = (int)wfy;
                double r = .0, g = .0, b = .0;

                for (int jy = y - n + 1; jy <= y + n; jy++)
                {
                    for (int jx = x - n + 1; jx <= x + n; jx++)
                    {
                        var w = lanczosWeight(wfx - jx, n) * lanczosWeight(wfy - jy, n);
                        if (w == 0) continue;
                        var sx = (jx < 0 || jx >= sw) ? x : jx;
                        var sy = (jy < 0 || jy >= sh) ? y : jy;
                        var sPos = 4 * sx + bmpData.Stride * sy;
                        b += bmpScan0[sPos + 0] * w;
                        g += bmpScan0[sPos + 1] * w;
                        r += bmpScan0[sPos + 2] * w;
                    }
                }

                var dPos = 4 * ix + dstData.Stride * iy;
                dstScan0[dPos + 0] = trimByte(b);
                dstScan0[dPos + 1] = trimByte(g);
                dstScan0[dPos + 2] = trimByte(r);
            }
        });
    }

    bmp.UnlockBits(bmpData);
    dst.UnlockBits(dstData);

    return dst;
}

サンプル範囲がかなり広くなったので処理時間も結構長くなったと思います。
それでは恒例の水蓮の画像を拡大した画像を見てみましょう。
（800x600pxから1280x960pxへ拡大）

…もはや、この画像じゃ役不足だ！（誤用）
良くなったのかはっきりしませんが、気持ち元から乗ってるノイズがより鮮明に見えるようになった気がします（逆効果）

まあ、今度もっとまともな画像データ使って比較しますよ！（汗）

Spline36への換装

AviUtlなどにお世話になっている方は『Spline36』というリサイズフィルタをよく利用したり、目にすると思います。
Lanczosより速く、アニメ映像などと相性が良いと聞きます。
そして実は今回のソースにおいてLanczos3の重み付け関数はそのSpline36のものにスッと差し替えることが可能です。

こちらのサイトを参考にさせていただきました。

Splineリサイズの係数式とその導出法

これによると重み付け関数はこんな感じになります。（dは絶対値）

spline36_weight(d) = 13/11*d^3 - 453/209*d^2 - 3/209*d + 1          (0 <= d <= 1)
                     -6/11*d^3 + 612/209*d^2 - 1038/209*d + 540/209 (1 < d <= 2)
                     1/11*d^3 - 159/209*d^2 + 434/209*d - 384/209   (2 < d <= 3)

という具合に係数が予め決まっています。
この係数を導き出すにはnの値が必要です。36という数字はLanczos3の3にあたるnを2倍（直径）の2乗（面積）したものだそうで、
n = 3ならば以下のように36になるわけです。

(3*2)^2 = 36

ちなみに参考サイトにあるとおり、係数はnを元に連立方程式を解く必要があり、簡単には計算式を実装することはできないのでnの値ごとに予め算出済みの定数として埋め込んでおくのがセオリーっぽいです。

プログラム実装

計算回数を減らすため、予め通分と展開を行った式を用意しました。

static double spline36Weight(double d)
{
    d = Math.Abs(d);
    if (d < 1.0)
    {
        return (((247.0 * d - 453.0) * d - 3.0) * d + 209.0) / 209.0;
    }
    if (d < 2.0)
    {
        return (((-114.0 * d + 612.0) * d - 1038.0) * d + 540.0) / 209.0;
    }
    if (d < 3.0)
    {
        return (((19.0 * d - 159.0) * d + 434.0) * d - 384.0) / 209.0;
    }
    return 0.0;
}

とりあえずこのコードを追加して、
lanczosWeight(d, n)を呼び出しているところをspline36Weight(d)に置き換えるだけです。

それじゃあ実行。
（800x600pxから1280x960pxへ拡大）

あい。。。わからん／(^o^)＼
とにかくLanczosよりかは速い！

おわりに

正直言ってここまでくると原理を完全に把握できないので、重み付けアルゴリズムそのものについてはご質問されても答えられないのでご容赦ください(^_^;)

しかしながら、いままで拡大ばっかりやってきたのですが、このコードで縮小してみるとどうなるのか…
実際にやってみればわかると思いますが、実写だとあまりわからなかったりしますが、イラスト等ではイマイチな結果になると思います。
バイキュービックやバイリニアもどっこいなものになると思います。

実は拡大と縮小ではサンプルの仕方を変えないといけないようです。
次回はその説明という名の他サイトへの丸投げと面積平均法についてご紹介したいと思います。