2次元アフィン変換の高速転送(3) - melpon日記 - HaskellもC++もまともに扱えないへたれのページ

穴あきの問題は解決したので、次は速度が遅いという問題について考えてみます。

前回のプログラム

前回のコードは

// 逆行列を求める
matrix inv(mat);
inv.invert();
for (int y = lt.y; y <= rb.y; y++)
{
    for (int x = lt.x; x <= rb.x; x++)
    {
        fixed_point p = inv.transform(fixed_point(x, y)); // x, y から対応する src の位置を求める
        dst->set_pixel(x, y, src->get_pixel(p.x, p.y));
    }
}

こうやって書いていました。
ループの一番深い場所を高速化しないといけないわけですが、以前書いたように set_pixel, get_pixel はポインタで行う操作（とクリッピング操作）を簡略化しただけなので気にしないことにします。
問題は inv.transform です。transform の実装は 2次元アフィン変換の高速転送(1) で書きました。コードを見てみると、

fixed_point transform(const fixed_point& p) const
{
    return fixed_point(
        m11 * p.x + m12 * p.y + m13,
        m21 * p.x + m22 * p.y + m23);
}

このように乗算を4回も行っています。
これを毎ループ行うのは BREW には少々辛いので、これを無くしてみましょう。

差分で計算

上記プログラムの内側のループ、つまり X 方向のループ中で 1 ループあたり x が 1 ずつ増えていきますが、それに対して src の座標（プログラム中の p）の値がどれぐらい変化するのかを考えてみます。
$\begin{pmatrix}x^\prime \\ y^\prime \\ 1\end{pmatrix}=\begin{pmatrix}a_{11}x + a_{12}y + a_{13} \\ a_{21}x + a_{22}y + a_{23} \\ 1\end{pmatrix}$
で、x が 1 増えると、
$\begin{pmatrix}x^{\prime\prime} \\ y^{\prime\prime} \\ 1\end{pmatrix}=\begin{pmatrix}a_{11}(x+1) + a_{12}y + a_{13} \\ a_{21}(x+1) + a_{22}y + a_{23} \\ 1\end{pmatrix}$
となります。この差をとってみると、
$\Large \begin{matrix} & \begin{pmatrix}x^{\prime\prime} \\ y^{\prime\prime} \\ 1\end{pmatrix} &-& \begin{pmatrix}x^\prime \\ y^\prime \\ 1\end{pmatrix} \\ =& \begin{pmatrix}a_{11}(x + 1) + a_{12}y + a_{13} \\ a_{21}(x + 1) + a_{22}y + a_{23} \\ 1\end{pmatrix} &-& \begin{pmatrix}a_{11}x + a_{12}y + a_{13} \\ a_{21}x + a_{22}y + a_{23} \\ 1\end{pmatrix} \\ =& \begin{pmatrix}a_{11}(x + 1) - a_{11}x \\ a_{21}(x + 1) - a_{21}x \\ 0 & \end{pmatrix} \\ =& \begin{pmatrix}a_{11} \\ a_{21} \\ 0\end{pmatrix} \end{matrix}$
となります。
つまり、x が 1 つ増えるたびに p.x は $a_{11}$ 、p.y は $a_{21}$ だけ増えていることになります。毎ループこれを足すだけで、乗算を必要とせずに X ループが行えそうです。

for (int y = lt.y; y <= rb.y; y++)
{
    fixed_point p = inv.transform(fixed_point(lt.x, y)); // X 方向スキャン用の初期位置
    for (int x = lt.x; x <= rb.x; x++)
    {
        dst->set_pixel(x, y, src->get_pixel(p.x, p.y));
        p.x += inv.m11;
        p.y += inv.m21;
    }
}

これで X ループ中では座標変換を行わないので乗算が消えました。

同様に、あまり効果はないですが Y ループ中でも座標変換を行わないようにしてみましょう。
上記プログラムの外側のループ、つまり Y 方向のループ中で 1 ループあたり y が 1 ずつ増えていきますが、それに対して src の座標（プログラム中の p）の値がどれぐらい変化するのかを考えてみます。
$\begin{pmatrix}x^\prime \\ y^\prime \\ 1\end{pmatrix}=\begin{pmatrix}a_{11}x + a_{12}y + a_{13} \\ a_{21}x + a_{22}y + a_{23} \\ 1\end{pmatrix}$
で、y が 1 増えると、
$\begin{pmatrix}x^{\prime\prime} \\ y^{\prime\prime} \\ 1\end{pmatrix}=\begin{pmatrix}a_{11}x + a_{12}(y+1) + a_{13} \\ a_{21}x + a_{22}(y+1) + a_{23} \\ 1\end{pmatrix}$
となります。この差をとってみると、
$\Large \begin{matrix} & \begin{pmatrix}x^{\prime\prime} \\ y^{\prime\prime} \\ 1\end{pmatrix} &-& \begin{pmatrix}x^\prime \\ y^\prime \\ 1\end{pmatrix} \\ =& \begin{pmatrix}a_{11}x + a_{12}(y+1) + a_{13} \\ a_{21}x + a_{22}(y+1) + a_{23} \\ 1\end{pmatrix} &-& \begin{pmatrix}a_{11}x + a_{12}y + a_{13} \\ a_{21}x + a_{22}y + a_{23} \\ 1\end{pmatrix} \\ =& \begin{pmatrix}a_{12}(y+1) - a_{12}y \\ a_{22}(y+1) - a_{22}y \\ 0 & \end{pmatrix} \\ =& \begin{pmatrix}a_{12} \\ a_{22} \\ 0\end{pmatrix} \end{matrix}$
となります。
つまり、y が 1 つ増えるたびに p.x は $a_{12}$ 、p.y は $a_{22}$ だけ増えていることになります。毎ループこれを足すだけで、乗算を必要とせずに Y ループが行えそうです。

fixed_point base = inv.transform(fixed_point(lt.x, lt.y)); // XY 方向スキャン用の初期位置
for (int y = lt.y; y <= rb.y; y++)
{
    fixed_point p = base;
    for (int x = lt.x; x <= rb.x; x++)
    {
        dst->set_pixel(x, y, src->get_pixel(p.x, p.y));
        p.x += inv.m11;
        p.y += inv.m21;
    }
    base.x += inv.m12;
    base.y += inv.m22;
}

実際のコード

あとは

fixed や fixed_point を使っている部分を int で計算してやる（最適化レベルに影響してくるので）
set_pixel とか get_pixel とか書いてる部分をちゃんとポインタを使った実装にしてやる
正確な座標を取るのと誤差対策のために dst を 0.5pix ずらして考える
ちゃんとクリッピングを行ってやる

とかいろいろあるのですが、それらは本題とは違うので省略します。
実際に書いたコードはこんな感じです。

// 8ビットインデックスのイメージから RGB565 のイメージへアフィン変換を掛けながら転送する
// パレット 0 番のインデックスは描画しない。
void draw_affine(
    uint16* dst,             // 転送先イメージ
    uint32 dpitch,           // dst の 1line のバイト数
    const fixed_rect& dclip, // 転送先クリップ領域
    const byte* src,         // 転送元イメージ
    uint32 spitch,           // src の 1line のバイト数
    const fixed_rect& sclip, // 転送元クリップ領域
    const uint16* pal,       // パレットデータ
    const matrix& mat)       // 座標変換用行列
{
    // (0, 0, sclip.w, sclip.h) の矩形をアフィン変換する
    array<fixed_point, 4> ps(mat.transform(sclip.w, sclip.h));

    // アフィン変換した四角形から、その四角形を包含する最小の矩形を計算
    fixed_point lt(
        std::min(std::min(std::min(ps[0].x, ps[1].x), ps[2].x), ps[3].x),
        std::min(std::min(std::min(ps[0].y, ps[1].y), ps[2].y), ps[3].y));
    fixed_point rb(
        std::max(std::max(std::max(ps[0].x, ps[1].x), ps[2].x), ps[3].x),
        std::max(std::max(std::max(ps[0].y, ps[1].y), ps[2].y), ps[3].y));
    // lt - rb で幅と高さを取得できるようにしたいので 1 を足しておく
    rb += fixed_point(1, 1);

    // クリッピングして左右上下のはみ出した量を求める
    fixed cl = std::max(fixed(0), dclip.x - lt.x);
    fixed cr = std::max(fixed(0), rb.x - (dclip.x + dclip.w));
    fixed cu = std::max(fixed(0), dclip.y - lt.y);
    fixed cd = std::max(fixed(0), rb.y - (dclip.y + dclip.h));

    // 描画する場所がない
    if (rb.x - lt.x - cl - cr < 0 ||
        rb.y - lt.y - cu - cd < 0) return;

    // dst を描画開始位置に設定、src はクリップ領域の左上が (0,0) になるように設定
    dst = (uint16*)((byte*)dst + ((lt.y + cu).get() >> 16) * dpitch + ((lt.x + cl).get() >> 16));
    src = src + (sclip.y.get() >> 16) * spitch + (sclip.x.get() >> 16);

    // 逆行列を求めてスキャン開始位置を調べる
    matrix inv(mat);
    inv.invert();
    fixed_point base = inv.transform(lt);

    // 転送先 X座標 1pixel に対する、転送元の X,Y の移動量
    fixed_point deltaX(inv.m11, inv.m21);
    // 転送先 Y座標 1line に対する、転送元の X,Y の移動量
    fixed_point deltaY(inv.m12, inv.m22);

    // クリッピングによる調整と、誤差調整
    base += deltaX * cl + (deltaX >> 1);
    base += deltaY * cu + (deltaY >> 1);

    // 転送
    int dw = (rb.x - lt.x - cl - cr).get() >> 16;
    int dh = (rb.y - lt.y - cu - cd).get() >> 16;
    int sw = sclip.w.get() >> 16;
    int sh = sclip.h.get() >> 16;
    int bu = base.x.get();
    int bv = base.y.get();
    int dxu = deltaX.x.get();
    int dxv = deltaX.y.get();
    int dyu = deltaY.x.get();
    int dyv = deltaY.y.get();
    for (int j = 0; j < dh; j++)
    {
        int u = bu;
        int v = bv;
        for (int i = 0; i < dw; i++)
        {
            int x = u >> 16;
            int y = v >> 16;
            if ((uint32)x < (uint32)sw && (uint32)y < (uint32)sh)
            {
                int n = src[y * spitch + x];
                if (n != 0) dst[i] = pal[n];
            }
            u += dxu;
            v += dxv;
        }
        bu += dyu;
        bv += dyv;
        dst = (uint16*)((byte*)dst + dpitch);
    }
}

これで完成です。

もっと速くする

まず、場合分けで速くする方法が考えられます。
例えば平行移動や拡大縮小、90, 180, 270度回転の場合は矩形が維持されるため src の座標が外側にはみ出さないので、一番深いループ中の if 文が省略できます。
平行移動のみの場合はもっと高速に処理することができます。これぐらいの処理であれば 8 ループ展開も効いてくるでしょう。
spitch が 2 のべき乗サイズであれば、アドレスの計算で乗算を使う必要も無くなります。これはうまくやらないとメモリ使用量が増えてしまうのでよく考える必要がありそうです。

また、それぞれの転送ルーチンはインラインアセンブラで書くという手もあります。
それぞれネックになっている部分を調べて、必要な箇所だけインラインアセンブラで書いていけばいいと思います。