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

最近、BREW でアフィン変換を掛けながら転送する処理を書きました。
2,3回ぐらいに分けてそれを書いていきたいと思います。

固定小数点数

基本的に計算は全部固定小数で行っています。固定小数点数クラス fixed は以前書きました。

固定小数点数クラス - melpon日記 - HaskellもC++もまともに扱えないへたれのページ

また fixed_point というクラスが出てきますが、これはメンバに fixed x, y という 2 つの値を持っている構造体です。

アフィン変換用の行列を考える

アフィン変換は回転と拡大縮小とせん断と平行移動を組み合わせたりできるわけですが、その変換を表現するための行列を用意します。
2次元では平行移動を含めると 3x3 の行列が必要になるのですが 3 行目はそれぞれ 0,0,1 で固定でいいので、2x3 の 6 つの値で表現できます。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}}$
別に 3x2 の行列でもよかった気がするのですが、DoJa のメソッドは 2x3 になっていたので、そちらに合わせています。

この行列にそれぞれ変換用の行列を掛けてやることで、複数の変換を同時に行うことができるようです。
ちなみに自分は OpenGL が好きな人なので変換行列を後ろに掛けていきます。この場合、各変換は逆順に適用されていくようです。逆順が嫌な場合は手前に掛けていくようにすればいいかと。
でも逆順のほうが (-1,-1)-(1,1) 空間なんかに入れてやりたいときに最後に掛けてやる行列を事前に用意しておけるので、効率が良かったりするからこうなってるんじゃないかなとか思ったりするのですがよく知らないです誰か教えてください。

平行移動

平行移動を表現する行列です。
$\Large{\begin{pmatrix}1 & 0 & x \\ 0 & 1 & y \\ 0 & 0 & 1\end{pmatrix}}$
掛けると以下のようになります。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}1 & 0 & x \\ 0 & 1 & y \\ 0 & 0 & 1\end{pmatrix}=\begin{pmatrix}a_{11} & a_{12} & a_{11}x + a_{12}y + a_{13} \\ a_{21} & a_{22} & a_{21}x + a_{22}y + a_{23} \\ 0 & 0 & 1\end{pmatrix}}$

拡大縮小

拡大縮小（スケーリング）を表現する行列です。
スケール値 $x$ , $y$ をマイナスの値にすれば、左右反転や上下反転もできます。
$\Large{\begin{pmatrix}x & 0 & 0 \\ 0 & y & 0 \\ 0 & 0 & 1\end{pmatrix}}$
掛けると以下のようになります。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}x & 0 & 0 \\ 0 & y & 0 \\ 0 & 0 & 1\end{pmatrix}=\begin{pmatrix}a_{11}x & a_{12}y & a_{13} \\ a_{21}x & a_{22}y & a_{23} \\ 0 & 0 & 1\end{pmatrix}}$

せん断

せん断を表現する行列です。せん断をうまく使えば影の表現にも使えると思います。
$\Large{\begin{pmatrix}1 & x & 0 \\ y & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}}$
掛けると以下のようになります。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}1 & x & 0 \\ y & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}=\begin{pmatrix}a_{11}+a_{12}y & a_{11}x+a_{12} & a_{13} \\ a_{21}+a_{22}y & a_{21}x+a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}}$

回転

回転を表現する行列です。
$\Large{\begin{pmatrix}\cos\theta & \sin\theta & 0 \\ -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 1\end{pmatrix}}$
掛けると以下のようになります。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}\cos\theta & \sin\theta & 0 \\ -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 1\end{pmatrix}=\begin{pmatrix}a_{11}\cos\theta-a_{12}\sin\theta & a_{11}\sin\theta+a_{12}\cos\theta & a_{13} \\ a_{21}\cos\theta - a_{22}\sin\theta & a_{21}\sin\theta + a_{22}\cos\theta & a_{23} \\ 0 & 0 & 1\end{pmatrix}}$

実際に変換後の座標を求める

$x$ , $y$ にある点から変換後の座標を求めるには以下のように掛けます。
$\Large{\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}\begin{pmatrix}x \\ y \\ 1\end{pmatrix}=\begin{pmatrix}a_{11}x+a_{12}y+a_{13} \\ a_{21}x+a_{22}y+a_{23} \\ 1\end{pmatrix}}$
様々なアフィン変換を加えた後、一度の計算で座標が求まるのが行列のすごいところですね。

逆行列

後で説明しますが、転送を行うために逆行列を作る必要が出てきます。
3x3 で逆行列を計算するための公式は結構複雑なのですが、3行目が固定値なので結構簡略化できます。
簡略化した後は以下のようになります。
$\Large\begin{pmatrix}a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ 0 & 0 & 1\end{pmatrix}^{-1}=\frac{1}{detA}\begin{pmatrix}a_{22} & -a_{12} & a_{12}a_{23} - a_{13}a_{22} \\ -a_{21} & a_{11} & a_{13}a_{21} - a_{11}a_{23} \\ 0 & 0 & 1\end{pmatrix}$ $\large(detA\neq0)$
$\Large detA=a_{11}a_{22}-a_{12}a_{21}$

アフィン変換用の行列クラスを作る

上記の変換を実際にプログラムに直してみると以下のようになりました。

class matrix
{
public:
    fixed m11;
    fixed m12;
    fixed m13;
    fixed m21;
    fixed m22;
    fixed m23;

    matrix() : m11(1), m12(0), m13(0), m21(0), m22(1), m23(0) { }
    matrix(
        const fixed& v11, const fixed& v12, const fixed& v13,
        const fixed& v21, const fixed& v22, const fixed& v23)
        : m11(v11), m12(v12), m13(v13), m21(v21), m22(v22), m23(v23)
    {
    }
    matrix(const matrix& rhs)
        : m11(rhs.m11), m12(rhs.m12), m13(rhs.m13)
        , m21(rhs.m21), m22(rhs.m22), m23(rhs.m23)
    {
    }

    const matrix& operator=(const matrix& rhs)
    {
        m11 = rhs.m11;
        m12 = rhs.m12;
        m13 = rhs.m13;
        m21 = rhs.m21;
        m22 = rhs.m22;
        m23 = rhs.m23;
        return *this;
    }

    // 恒等行列に設定
    void reset()
    {
        m11 = 1; m12 = 0; m13 = 0;
        m21 = 0; m22 = 1; m23 = 0;
    }
    // 平行移動
    void translate(const fixed& x, const fixed& y)
    {
        m13 += m11 * x + m12 * y;
        m23 += m21 * x + m22 * y;
    }
    // 回転
    void rotate(const fixed& rotate)
    {
        fixed s = math::sin(rotate);
        fixed c = math::cos(rotate);
        fixed t = m11;
        m11 = t   * c + m12 * s;
        m12 = m12 * c - t   * s;
        t = m21;
        m21 = t   * c + m22 * s;
        m22 = m22 * c - t   * s;
    }
    // スケーリング
    void scale(const fixed& sx, const fixed& sy)
    {
        m11 *= sx;
        m12 *= sy;
        m21 *= sx;
        m22 *= sy;
    }
    // せん断
    void shear(const fixed& sx, const fixed& sy)
    {
        fixed t = m11;
        m11 += m12 * sy;
        m12 += t   * sx;
        t = m21;
        m21 += m22 * sy;
        m22 += t   * sx;
    }

    // 座標変換
    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);
    }

    // 逆行列
    void invert()
    {
        fixed det = m11 * m22 - m12 * m21;
        if (det == 0) throw invert_not_exist();
        fixed rdet = 1 / det;
        fixed t = m13;
        m13 = (m12 * m23 - t   * m22) * rdet;
        m23 = (t   * m21 - m11 * m23) * rdet;
        m12 = -m12 * rdet;
        m21 = -m21 * rdet;
        t = m11;
        m11 = m22 * rdet;
        m22 = t   * rdet;
    }
};

そのまま定義どおりです。

行列クラスを使ってみる

こんな感じで使います。

matrix mat;
matrix.translate(50, 50);   // 4
matrix.scale(1.5f, 0.5f);   // 3
matrix.rotate(45);          // 2
matrix.translate(-25, -25); // 1

fixed_point points[] = {
    fixed_point(0, 0),
    fixed_point(0, 50),
    fixed_point(50, 50),
    fixed_point(50, 0),
};

points[0] = mat.transform(points[0]);
points[1] = mat.transform(points[1]);
points[2] = mat.transform(points[2]);
points[3] = mat.transform(points[3]);

graphics->fill_polygon(points, color::blue());

各変換は 1, 2, 3, 4 の順番に適用されていきます。
なので、この場合だと

初期位置が (0,0)-(50,50) の位置にある矩形に 1 を適用して平行移動するので、

このようになり、次は 2 を適用して（原点を中心に）回転するので

このようになり、次は 3 を適用して（原点を中心に）拡大縮小するので

このようになり、最後に 4 を適用して平行移動を行い、最終的に表示されるのは、

このようになります。

これで行列によるアフィン変換はできるようになりました。
次はこれを使って実際に転送処理を行っていこうと思います。
→書きました。

2次元アフィン変換の高速転送(2)