simplestarの技術ブログ

目的を書いて、思想と試行、結果と考察、そして具体的な手段を記録します。

Math:Unityで対称行列の対角化と同時に固有値・固有ベクトルを求める方法

ゲーム業界でも固有値固有ベクトルという単語はよく出てきますが、対称行列を与えられた状態から具体的な解き方までご存知の方はそう多くいません。
答えを先に言うと、対称行列を対角化すると、その時点で対角行列の各要素は固有値、左右の回転行列はその固有値に対する固有ベクトルで構成されることになります。
数学カテゴリ第二回目の今回は、Unity での対称行列の対角化と、固有値固有ベクトルの求め方を示したいと思います。

今回示すプログラムを走らせると、次の動作を確認できます。
f:id:simplestar_tech:20170412084653g:plain
元の行列の形をまったく変えずに、左下の行列が対角行列に相似変換されていく様子と、それを実現する左右の回転行列(相似変換行列)が更新されていく様子が確認できます。

さて、固有値固有ベクトルの求め方だけでネットで調べると、固有方程式を立てて、高次方程式の解として固有値を得るというものが多く見つかります。
私が固有方程式を解く方法を嫌う理由としては、与えられた行列が解に近い形をしているとまず解けませんし、立てた高次方程式が結局解けないという問題や、重解となったときに対応する複数の固有ベクトル固有値から求められないという問題があるからです。(固有方程式は筆算でよく使われますが、数値計算プログラムとしては欠陥だらけですね)

もう少し資料を探してみると、ギブンス回転による相似変換(ヤコビ法)、ハウスホルダー変換、三重対角行列変換とQR法(ウィルキンソンの移動法)などが見つかります。
速度の速い順に手法を紹介しましたが、対照的に精度が上がっていきます。
ゲーム業界で使われるだろう次元数は最大で4くらいですので、この中で最速のギブンス回転による相似変換(ヤコビ法)が、ゲーム業界向けに最適な手法と言えます。

そこで、今回はギブンス回転による相似変換(ヤコビ法)によって、対称行列の対角化と固有値固有ベクトルを求めるプログラムを示します。

以下の通り

        /// <summary>
        /// 対称行列の固有値・固有ベクトル分解
        /// </summary>
        /// <param name="A">[in]分解する対象のNxNの実対称行列</param>
        /// <param name="D">[out]実対称行列に対応する、対角成分を固有値とする行列</param>
        /// <param name="L">[out]その列の固有値に対する固有列ベクトルで構成される行列</param>
        /// <param name="e">収束判定用最小誤差</param>
        /// <param name="maxCount">最大ギブンス回転回数</param>
        /// <returns>最大回数以下で収束した時trueを返す</returns>
        public static bool EigenDecomposition(float[][] A, float[][] D, float[][] L, float e, int maxCount)
        {
            int sz = A.Length;
            for (int i = 0; i < sz; i++)
            {
                for (int j = 0; j < sz; j++)
                {
                    D[i][j] = A[i][j];
                    if (j != i)
                    {
                        L[i][j] = 0;
                    }
                }
                L[i][i] = 1;
            }

            int count = 0;
            while (count++ < maxCount)
            {
                int row = 0;
                int col = 0;
                AbsoluteMaxOffDiagonalElementPosition(D, out row, out col);
                if (e > Mathf.Abs(D[row][col]))
                    break;
                // pp = qq = (pp - qq) * sin2θ + 2 * pq * cos2θ = 0
                // αsin2θ + βcos2θ = 0
                float alpha = (D[row][row] - D[col][col]);
                float beta = - 2 * D[row][col];
                // (αcos2θ) * (αcos2θ) + (αsin2θ) * (αsin2θ) = α * α 
                // (α * α + β * β)cos2θ^2 = α * α
                // cos2θ = Abs(α) / Sqrt(α * α + β * β)
                float cos2 = Mathf.Abs(alpha) / Mathf.Sqrt(alpha * alpha + beta * beta);
                // cosθ
                float cos1 = Mathf.Sqrt((1 + cos2) / 2);
                // sinθ
                float sin1 = Mathf.Sqrt((1 - cos2) / 2);
                if (0 > alpha * beta)
                    sin1 *= -1.0f;

                GivensRotation(D, row, col, sin1, cos1);

                SimilarityTransformationLeft(L, row, col, sin1, cos1);
            }

            if (count >= maxCount)
                return false;

            return true;
        }

        /// <summary>
        /// ギブンス回転
        /// </summary>
        /// <param name="A">入出力の対角化する行列</param>
        /// <param name="row">非対角成分における絶対値最大要素の行</param>
        /// <param name="col">非対角成分における絶対値最大要素の列</param>
        /// <param name="sin1">その要素=0としたときの回転量sinθ</param>
        /// <param name="cos1">その要素=0としたときの回転量cosθ</param>
        private static void GivensRotation(float[][] A, int row, int col, float sin1, float cos1)
        {
            float sin1cos1 = sin1 * cos1; // sinθcosθ
            float sin1_2 = sin1 * sin1; // sinθ^2
            float cos1_2 = cos1 * cos1; // cosθ^2

            // pp * cosθ^2 + qq * sinθ^2 - 2 pq * sinθcosθ
            float tmp1 = A[row][row] * cos1_2 + A[col][col] * sin1_2 - 2 * A[row][col] * sin1cos1;
            // pp * sinθ^2 + qq * cosθ^2 + 2 pq * sinθcosθ
            float tmp2 = A[row][row] * sin1_2 + A[col][col] * cos1_2 + 2 * A[row][col] * sin1cos1;

            A[row][col] = A[col][row] = 0; // pq, qp = (pp - qq) / 2 * sin2θ + pq * cos2θ = 0
            A[row][row] = tmp1; // pp
            A[col][col] = tmp2; // qq

            int sz = A.Length;

            // A[*][p], A[*][q]
            for (int r = 0; r < sz; ++r)
            {
                if (r == row || r == col) continue;
                tmp1 = A[r][row] * cos1 - A[r][col] * sin1;
                tmp2 = A[r][row] * sin1 + A[r][col] * cos1;
                A[r][row] = tmp1;
                A[r][col] = tmp2;
            }
            // A[p][*], A[q][*]
            for (int c = 0; c < sz; ++c)
            {
                if (c == row || c == col) continue;
                tmp1 = A[row][c] * cos1 - A[col][c] * sin1;
                tmp2 = A[row][c] * sin1 + A[col][c] * cos1;
                A[row][c] = tmp1;
                A[col][c] = tmp2;
            }
        }

        /// <summary>
        /// 行列の相似変換(右側変換行列の積み上げ)
        /// </summary>
        /// <param name="L">入出力の相似変換行列</param>
        /// <param name="row">変換行</param>
        /// <param name="col">変換列</param>
        /// <param name="sin1">sinθ</param>
        /// <param name="cos1">cosθ</param>
        private static void SimilarityTransformationLeft(float[][] L, int row, int col, float sin1, float cos1)
        {
            int sz = L.Length;
            // for Right Matrix
            //for (int c = 0; c < sz; ++c)
            //{
            //    float d1 = cos1 * R[row][c] - sin1 * R[col][c];
            //    float d2 = sin1 * R[row][c] + cos1 * R[col][c];
            //    R[row][c] = d1;
            //    R[col][c] = d2;
            //}

            // for Left Matrix
            for (int r = 0; r < sz; ++r)
            {
                float d1 = L[r][row] * cos1 - L[r][col] * sin1;
                float d2 = L[r][row] * sin1 + L[r][col] * cos1;
                L[r][row] = d1;
                L[r][col] = d2;
            }
        }

        /// <summary>
        /// 対称行列の非対角成分における絶対値最大の要素の位置を返す
        /// </summary>
        /// <param name="A">対称行列</param>
        /// <param name="row">行の位置</param>
        /// <param name="col">列の位置</param>
        private static void AbsoluteMaxOffDiagonalElementPosition(float[][] A, out int row, out int col)
        {
            int sz = A.Length;
            row = 0;
            col = 1;
            float maxValue = Mathf.Abs(A[row][col]);

            for (int r = 0; r < sz; ++r)
            {
                for (int c = r + 1; c < sz; ++c)
                {
                    float val = Mathf.Abs(A[r][c]);
                    if (maxValue < val)
                    {
                        row = r;
                        col = c;
                        maxValue = val;
                    }
                }
            }
        }

この関数の使い方は次の通り

        int n = 6;
        SimpleMath.CreateMatrix(n, n, out _A);
        int row = 0;
        int col = 0;
        foreach (InputField child in panelA.GetComponentsInChildren<InputField>())
        {
            float random = Random.Range(-3, 3);
            child.text = random.ToString("0.00");
            _A[row][col++] = float.Parse(child.text);
            if (n == col)
            {
                col = 0;
                row++;
                if (n == row)
                {
                    break;
                }
            }
        }
        SimpleMath.CreateMatrix(n, n, out _AtA);
        SimpleMath.GetAutocorrelation(_A, _AtA);

        float[][] L;
        SimpleMath.CreateMatrix(n, n, out L);
        float[][] D;
        SimpleMath.CreateMatrix(n, n, out D);
        SimpleMath.EigenDecomposition(_AtA, D, L, 0.001f, n * n);

動画はギブンス回転の最大回数を n * n から 1 にして、一回転ずつの変換の様子を観察したものです。
だいたい 30~40 回転させると収束するようです。
これで n x n の対称行列の対角化、固有値固有ベクトルが得られるようになりました。

優秀な数値計算ライブラリに頼るのもアリですが、ちょこっとした固有値問題を解くケースで、上記コードを使ってみたいと思いまして作ってみました。
いつでも自分が参照できるように、ここに公開します。

以下、参考文献です。

ヤコビ法のロジックはここの資料を見て書きました。ギブンス回転により必ず対角化に向かうことの証明も書かれています。(大変助かりました、ありがとうございます!)
http://www.na.scitec.kobe-u.ac.jp/~yamamoto/lectures/appliedmathematics2004/chapter11.PDF

相似変換は、ランク、トレース、行列式固有値、固有方程式を保存するという説明向けに紹介したいページです。
そう、ギブンス回転を何回行っても固有値は変わりません。(だから、今回の方法で固有値を求めることができます)
線形代数II/相似変換に対するトレース、行列式、固有値の保存 - 武内@筑波大

ギブンス回転の実装の参考にしたコードが書かれているページです。
式の説明に落とし穴がなく、すばらしい資料になっています。
ただ、ギブンス回転のコードには誤りがあるようで、そのまま利用するとなかなか対角行列に収束しませんでした。
fussy.web.fc2.com

ヤコビ法やハウスホルダー変換が簡潔に説明されている資料です。
上の資料との重複もあり、理解の確認用に使わせてもらいました。
http://www.aics.riken.jp/aicssite/wp-content/uploads/2014/03/sps14_enshu2-1.pdf

Math:Unityでの連立方程式の解き方

本ブログの数学カテゴリを作りまして、第一回目は連立方程式の解き方を示したいと思います。

みなさんご存じの通り、中学校で学ぶ連立方程式のことですが、これがなかなか簡単ではありません。
例えば、次の連立方程式がテストで出された時、変数 x, y, z, w, u, v を時間内にみなさんは求められましたでしょうか?

f:id:simplestar_tech:20170402193356j:plain

すみません、中学校ではここまで難易度の高い連立方程式は出題されませんでしたね。
ちなみに、上の式の答えは次の通り

x = 7.0
y = 5.0
z =-5.0
w = 3.0
u = 6.0
v = 4.0

さて、理系の大学を出て、そっち方面のお仕事に就くと、実際にこうした数値を計算する場面に遭遇します。
これをマイクロ秒オーダーで解くのはなかなかテクニックが必要でして、今日はその計算テクニックを動作確認済みのプログラムとして紹介します。

        /// <summary>
        /// 連立方程式を解く
        /// </summary>
        /// <param name="J">係数行列</param>
        /// <param name="b">定数項ベクトル</param>
        /// <param name="x">解ベクトル</param>
        /// <returns>true:成功,false:失敗</returns>
        public static bool SolveEq(float[][] J, float[] b, float[] x)
        {
            int jCount = J.Length;
            int n = x.Length;

            // 係数行列を正定値対称行列化
            float[][] JtJ;
            SimpleMath.CreateMatrix(n, n, out JtJ);
            for (int row = 0; row < n; row++)
            {
                for (int col = 0; col < n; col++)
                {
                    JtJ[row][col] = 0;
                    for (int i = 0; i < jCount; i++)
                    {
                        JtJ[row][col] += J[i][row] * J[i][col];
                    }
                }
            }

            // 定数項を対称行列化にともない修正
            float[] Jtb = new float[n];
            for (int row = 0; row < n; row++)
            {
                Jtb[row] = 0;
                for (int col = 0; col < jCount; col++)
                {
                    Jtb[row] += J[col][row] * b[col];
                }
            }

            // 修正 Cholesky 分解
            float[][] L;
            SimpleMath.CreateMatrix(n, n, out L);
            float[] d = new float[n];
            SimpleMath.ModifiedCholeskyDecomposition(JtJ, L, d);

            // 以下、連立方程式を解く

            // 前進代入
            float[] y = new float[n];
            for (int row = 0; row < n; row++)
            {
                y[row] = Jtb[row];
                for (int col = 0; col < row; col++)
                {
                    y[row] -= L[row][col] * y[col];
                }
                y[row] /= L[row][row];
            }
            // 後退代入
            for (int row = n - 1; row >= 0; row--)
            {
                x[row] = y[row];
                for (int col = n - 1; col > row; col--)
                {
                    x[row] -= (L[col][row] * d[row]) * x[col];
                }
            }
            return true;
        }

        /// <summary>
        /// 修正コレスキー分解(modified Cholesky decomposition)
        /// </summary>
        /// <remarks>
        /// 対称行列A(n×n)を下三角行列(L:Lower triangular matrix)と対角行列の積(LDL^T)に分解する
        /// l_ii * d_i = 1とした場合
        /// L: i > jの要素が非ゼロで対角成分は1
        /// </remarks>
        /// <param name="A">A n×nの対称行列</param>
        /// <param name="L">対角成分が1の下三角行列</param>
        /// <param name="d">対角行列(対角成分のみ)</param>
        /// <returns>true:成功,false:失敗</returns>
        public static bool ModifiedCholeskyDecomposition(float[][] A, float[][] L, float[] d)
        {
            int n = d.Length;
            if (n <= 0)
                return false;

            L[0][0] = A[0][0];
            d[0] = 1.0f / L[0][0];

            for (int row = 1; row < n; ++row)
            {
                for (int col = 0; col < n; ++col)
                {
                    L[row][col] = 0;
                    if (col <= row)
                    {
                        float A_SumLdL = A[row][col];
                        for (int nonZeroIdx = 0; nonZeroIdx < col; ++nonZeroIdx)
                        {
                            A_SumLdL -= L[row][nonZeroIdx] * d[nonZeroIdx] * L[col][nonZeroIdx];
                        }
                        L[row][col] = A_SumLdL;
                    }
                }
                d[row] = 1.0f / L[row][row];
            }

            return true;
        }

係数行列をfloatの2次元配列として渡して、定数項ベクトルもfloat配列で渡すと、解ベクトルがfloat配列で返ってきます。
正しく動作するか確認するプログラムは次の通りです。

    public Transform panelA;
    public Transform panelx;
    public Transform panelb;

    float[][] _A;
    float[] _b;
    float[] _x;

    public void CreateProblem()
    {
        int n = 6;
        SimpleMath.CreateMatrix(n, n, out _A);
        int row = 0;
        int col = 0;
        foreach (InputField child in panelA.GetComponentsInChildren<InputField>())
        {
            float random = Random.Range(-9.0f, 9.0f);
            child.text = random.ToString("0.00");
            _A[row][col++] = float.Parse(child.text);
            if (n == col)
            {
                col = 0;
                row++;
                if (n == row)
                {
                    break;
                }
            }
        }

        _x = new float[n];
        
        int i = 0;
        foreach (InputField child in panelx.GetComponentsInChildren<InputField>())
        {
            float random = Random.Range(-9, 9);
            child.text = random.ToString("0.00");
            _x[i++] = float.Parse(child.text);
            if (n == i)
            {
                break;
            }
        }

        float[][] b;
        float[][] x;
        SimpleMath.CreateMatrix(n, 1, out b);
        SimpleMath.CreateMatrix(n, 1, out x);

        for (int k = 0; k < n; k++)
        {
            x[k][0] = _x[k];
        }

        if (SimpleMath.MatrixMultiply(_A, x, b))
        {
            _b = new float[n];
            for (int k = 0; k < n; k++)
            {
                _b[k] = b[k][0];
            }
            i = 0;
            foreach (InputField child in panelb.GetComponentsInChildren<InputField>())
            {
                child.text = _b[i++].ToString("0.00");
                if (n == i)
                {
                    break;
                }
            }
        }
    }

    public void SolveProblem()
    {
        if(SimpleMath.SolveEq(_A, _b, _x))
        {
            for (int i = 0; i < _x.Length; i++)
            {
                Debug.Log("x[" + i + "] = " + _x[i].ToString("0.000"));
            }
        }
    }

この連立方程式を解く関数はちょっとその辺の LU 分解タイプとは異なり、式の数が変数の数よりも多くても、もっともらしい解を返します。
疑似逆行列とか最小二乗法で調べてもらえると、原理の説明にたどり着けるでしょう。
そこで出てくる正定値対称行列のLU分解はコレスキー分解と呼ばれていて、そのコレスキー分解もまた高速化された修正版が存在します。
今回はその修正コレスキー分解を用いて、連立方程式を解きました。

変数の数が 6 より大きい場合も利用できますので、もし同じ問題でお困りの方がいましたら、使ってみてください。

以下、参考にしたサイトです。

こちらの記事で、コレスキー分解について学びました。
コレスキー分解 - PukiWiki for PBCG Lab

Unity:4画面テクスチャビューアの配置スクリプト

画像処理していると、複数枚のテクスチャを見比べたりしていきたいところなのですが
Android 環境にもっていくと、配置がずれていたりして困ったことになったので
簡単な配置スクリプトを作りました。

自分が再利用する目的で記事書きます。

f:id:simplestar_tech:20170402153656j:plain

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class UIQuadViewBehaviour : MonoBehaviour {

    public enum QuadLocation
    {
        LeftTop = 0,
        LeftBottom,
        RightTop,
        RightBottom
    }

    public QuadLocation MyLocation;

	// Use this for initialization
	void Start () {
        ResizeLocate();
    }
	
	// Update is called once per frame
	void Update () {
    }

    private void ResizeLocate()
    {
        int width = Screen.width / 2;
        int height = Screen.height / 2;

        Vector3 location = Vector3.zero;

        switch (MyLocation)
        {
            case QuadLocation.LeftTop:
                location.x = -width / 2;
                location.y = height / 2;
                break;
            case QuadLocation.LeftBottom:
                location.x = -width / 2;
                location.y = -height / 2;
                break;
            case QuadLocation.RightTop:
                location.x = width / 2;
                location.y = height / 2;
                break;
            case QuadLocation.RightBottom:
                location.x = width / 2;
                location.y = -height / 2;
                break;
            default:
                break;
        }

        RectTransform rc = GetComponent<RectTransform>();

        rc.SetSizeWithCurrentAnchors(RectTransform.Axis.Horizontal, width);
        rc.SetSizeWithCurrentAnchors(RectTransform.Axis.Vertical, height);
        rc.localPosition = location;
    }
}

UI の RawImage にこのスクリプトコンポーネントを付けると、起動時に画面4分割の配置に移動します。
Anchor Presets は center と middle の組み合わせに限ります。

PC でテストして Android で実行するとビューが壊れているなんていう問題は、一応これで回避
テクスチャの解像度やアスペクト比まで考えて、賢く配置してくれと言いたい文句もありますが、そこはテクスチャをセットする側が調整してあげて…

Unity:深度画像の取得

オブジェクトのエッジを強調したり、モデルベースの画像処理を行いたいときなどは、見ているシーンの深度画像が必要になります。(これは本当)

OpenCV に現在フレームの深度画像を渡したいときは float 配列としてカメラからのZ距離が得られると文句はないのです。
この記事ではその float 配列の取得方法を動作確認込みで示します。

f:id:simplestar_tech:20170402120750j:plain

Unity で深度画像を得る方法をいくつか存じておりますが、私が知る中で最速な方法を示したいと思います。

前回、前々回と触ってきた Compute Shader を使います。(以下の方法が具体的にイメージできない人は読み返してね)

Compute Shader を新規作成したら、次のコードを記述します。

// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel CSMain

float n_f;
float f;
Texture2D<float> _zBuffer;
RWTexture2D<float> _cameraZ;

[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
#if SHADER_API_GLES3
	// for Android
	_cameraZ[id.xy] = 1.0f / (n_f * (1.0f - _zBuffer[id.xy]) + f);
#else
	// for Windows
	_cameraZ[id.xy] = 1.0f / (n_f * _zBuffer[id.xy] + f);
#endif
}

単なるシェーダーではなく Compute Shader であるため UNITY_REVERSED_Z による判定が行えません。
Android では UNITY_REVERSED_Z ではなかったので泣く泣く SHADER_API_GLES3 でコードを変更しています。

この Compute Shader を使うスクリプトは次の通り

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

[RequireComponent(typeof(Camera))]
public class VirtualCameraBehaviour : MonoBehaviour {

    private int width;
    private int height;

    Camera _virtualCamera = null;

    RenderTexture _color = null;
    RenderTexture _zBuffer = null;

    public RenderTexture _cameraZ = null;

    public ComputeShader Depth;
    public GameObject uiResult0;
    public GameObject uiResult1;
    public GameObject uiResult2;

    void Start () {

        width = Screen.width / 2;
        height = Screen.height / 2;

        _virtualCamera = GetComponent<Camera>();

        _color = new RenderTexture(width, height, 0, RenderTextureFormat.ARGB32);
        _zBuffer = new RenderTexture(width, height, 32, RenderTextureFormat.Depth);

        _virtualCamera.SetTargetBuffers(_color.colorBuffer, _zBuffer.depthBuffer);

        _cameraZ = new RenderTexture(width, height, 0, RenderTextureFormat.RFloat);
        _cameraZ.enableRandomWrite = true;
        _cameraZ.Create();

        float n_inv = 1.0f / _virtualCamera.nearClipPlane;
        float f_inv = 1.0f / _virtualCamera.farClipPlane;
        Depth.SetFloat("n_f", n_inv - f_inv);
        Depth.SetFloat("f", f_inv);
        Depth.SetTexture(0, "_zBuffer", _zBuffer);
        Depth.SetTexture(0, "_cameraZ", _cameraZ);

        uiResult0.GetComponent<UnityEngine.UI.RawImage>().texture = _color;
        uiResult1.GetComponent<UnityEngine.UI.RawImage>().texture = _cameraZ;
        uiResult2.GetComponent<UnityEngine.UI.RawImage>().texture = _zBuffer;
    }
	
	void Update () {
		
	}

    private void OnPostRender()
    {
        Depth.Dispatch(0, _zBuffer.width / 8, _zBuffer.height / 8, 1);
    }
}

オフスクリーンレンダリングで、解像度を指定して VirtualCamera 画像を作り、そのカラー画像や深度画像を float 32bit 深度 1チャンネル画像として取得しています。

CPU 側で float 配列を手に入れたい場合は次の一工夫を加えて完了です。

// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel CSMain

float n_f;
float f;
int width;
Texture2D<float> _zBuffer;
RWTexture2D<float> _cameraZ;

// for float Array
RWStructuredBuffer<float> _zArray;

[numthreads(8,8,1)]
void CSMain (uint3 id : SV_DispatchThreadID)
{
#if SHADER_API_GLES3 // insted of UNITY_REVERSED_Z
	// for Android
	_cameraZ[id.xy] = 1.0f / (n_f * (1.0f - _zBuffer[id.xy]) + f);
#else
	// for Windows
	_cameraZ[id.xy] = 1.0f / (n_f * _zBuffer[id.xy] + f);
#endif

	_zArray[id.x + id.y * width] = _cameraZ[id.xy];
}

以下の対応で CPU 側の float 配列に深度画像の画素値が得られることを確認しました。
これで OpenCV にこのデータを使ってもらうことにより、さまざまなモデルベース処理が行えるという夢が広がります。

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

[RequireComponent(typeof(Camera))]
public class VirtualCameraBehaviour : MonoBehaviour {

    private int width;
    private int height;

    Camera _virtualCamera = null;

    RenderTexture _color = null;
    RenderTexture _zBuffer = null;

    public RenderTexture _cameraZ = null;

    public ComputeShader Depth;
    ComputeBuffer _zArray;
    float[] _zArrayData;

    public GameObject uiResult0;
    public GameObject uiResult1;
    public GameObject uiResult2;

    void Start () {

        width = Screen.width / 2;
        height = Screen.height / 2;

        _virtualCamera = GetComponent<Camera>();

        _color = new RenderTexture(width, height, 0, RenderTextureFormat.ARGB32);
        _zBuffer = new RenderTexture(width, height, 32, RenderTextureFormat.Depth);

        _virtualCamera.SetTargetBuffers(_color.colorBuffer, _zBuffer.depthBuffer);

        _cameraZ = new RenderTexture(width, height, 0, RenderTextureFormat.RFloat);
        _cameraZ.enableRandomWrite = true;
        _cameraZ.Create();

        float n_inv = 1.0f / _virtualCamera.nearClipPlane;
        float f_inv = 1.0f / _virtualCamera.farClipPlane;
        Depth.SetFloat("n_f", n_inv - f_inv);
        Depth.SetFloat("f", f_inv);
        Depth.SetTexture(0, "_zBuffer", _zBuffer);
        Depth.SetTexture(0, "_cameraZ", _cameraZ);

        _zArrayData = new float[width * height];
        _zArray = new ComputeBuffer(_zArrayData.Length, sizeof(float));
        Depth.SetBuffer(0, "_zArray", _zArray);
        Depth.SetInt("width", width);

        uiResult0.GetComponent<UnityEngine.UI.RawImage>().texture = _color;
        uiResult1.GetComponent<UnityEngine.UI.RawImage>().texture = _cameraZ;
        uiResult2.GetComponent<UnityEngine.UI.RawImage>().texture = _zBuffer;
    }
	
	void Update () {
		
	}

    private void OnPostRender()
    {
        Depth.Dispatch(0, _zBuffer.width / 8, _zBuffer.height / 8, 1);
        _zArray.GetData(_zArrayData);

        float av = 0;
        for (int i = 0; i < _zArrayData.Length; i++)
        {
            av += _zArrayData[i];
        }
        av /= _zArrayData.Length;
        Debug.Log("average z = " + av.ToString("0.0000"));
    }
}

Unity:Compute Shaderでカメラ画像処理する時に最初に書くコード

ひとつ前の記事で Compute Shader を Android で実行するサンプルを載せましたので、今回はその Compute Shader を使ってカメラ画像をリアルタイムに処理します。

f:id:simplestar_tech:20170326232644j:plain

できました!
これからリアルタイムカメラ画像処理を Compute Shader で行う際は、以下のコードをコピペして使っていこうと思います。

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class DeviceCameraBehaviour : MonoBehaviour {

    // Histogram Gather
    public GameObject resultCube;
    public ComputeShader RGBA2GRAY;
    public Texture2D Source;
    RenderTexture _result;
    int _kernelIndex = 0;
    ComputeBuffer _HistogramBuffer = null;
    uint[] _emptyBuffer = new uint[256];
    uint[] _histogramBuffer = new uint[256];

    // Device Camera
    string _activeDeviceName = "";
    WebCamTexture _activeTexture = null;

    // Use this for initialization
    void Start () {

        _HistogramBuffer = new ComputeBuffer(256, sizeof(uint));
        for (int i = 0; i < _emptyBuffer.Length; i++)
        {
            _emptyBuffer[i] = 0;
        }
        _HistogramBuffer.SetData(_emptyBuffer);

        _kernelIndex = RGBA2GRAY.FindKernel("KHistogramGather");
        RGBA2GRAY.SetBuffer(_kernelIndex, "_Histogram", _HistogramBuffer);
    }
	
	// Update is called once per frame
	void Update () {

        if (null != _activeTexture && _activeTexture.isPlaying && _activeTexture.didUpdateThisFrame)
        {
            if (null == _result)
            {
                _result = new RenderTexture(_activeTexture.width, _activeTexture.height, 0, RenderTextureFormat.ARGB32);
                _result.enableRandomWrite = true;
                _result.Create();
                RGBA2GRAY.SetTexture(_kernelIndex, "_Source", _activeTexture);
                RGBA2GRAY.SetTexture(_kernelIndex, "_Result", _result);
                RGBA2GRAY.SetVector("_SourceSize", new Vector2(_activeTexture.width, _activeTexture.height));
                _HistogramBuffer.SetData(_emptyBuffer);
                resultCube.GetComponent<Renderer>().material.mainTexture = _result;

                resultCube.transform.localScale = new Vector3((_activeTexture.width/(float)_activeTexture.height), 1, 1);
            }
            RGBA2GRAY.Dispatch(_kernelIndex, Mathf.CeilToInt(_activeTexture.width / 32f), Mathf.CeilToInt(_activeTexture.height / 32f), 1);
        }

        if (Input.GetKeyDown(KeyCode.Space))
        {
            Debug.Log("Time = " + Time.time.ToString("000.0000"));
            PlayDeviceCamera();

        }
        for (int i = 0; i < Input.touchCount; ++i)
        {
            if (Input.GetTouch(i).phase == TouchPhase.Began)
            {

                // Construct a ray from the current touch coordinates
                Ray ray = Camera.main.ScreenPointToRay(Input.GetTouch(i).position);
                // Create a particle if hit
                if (Physics.Raycast(ray))
                {
                    Debug.Log("Time = " + Time.time.ToString("000.0000"));
                    PlayDeviceCamera();
                }
            }
        }
    }

    private void PlayDeviceCamera()
    {
        WebCamDevice[] devices = WebCamTexture.devices;
        for (int i = 0; i < devices.Length; i++)
        {
            if (0 != string.Compare(_activeDeviceName, devices[i].name))
            {
                if (null != _activeTexture)
                {
                    _activeTexture.Stop();
                }
                _activeDeviceName = devices[i].name;
                _activeTexture = new WebCamTexture(devices[i].name);
                _activeTexture.Play();
                RGBA2GRAY.SetTexture(_kernelIndex, "_Source", _activeTexture);
                _result = null;
                Debug.Log(devices[i].name);
                break;
            }
        }
    }
}

任意の解像度を指定するとか、細かいことをこの後入れていく感じですね。
同じコードで Android の方でも、リアルタイムカメラ画像処理できることを確認しました。

Unity:AndroidでComputeShaderを使ってみました

画像処理では、並列計算による高速化が求められます。(常に…)
Compute Shader を使った並列計算による高速化は実装が簡単で効果抜群なのですが、DirectX 依存なので(と思い込んでいた私は) PC 上での動作に限られていると思っていました。

しかし、ふと Unity の Compute Shader って Android 用にビルドした場合も動作するのかな?
だとしたら、Unity すばらしいなって思いまして、さっそく試してみました。

左は入力のカラー画像、右は出力のグレースケール画像(…イラストは私が描きました)

f:id:simplestar_tech:20170326171134j:plain

できた!Unityすばらしい!

使った Compute Shader のコードはこちら↓

// file head
RWStructuredBuffer<uint> _Histogram;
RWTexture2D<float4> _Result;

// Gathering pass
Texture2D<float4> _Source;

uint2 _SourceSize;

// groupshared uint gs_histogram[256]; // group shared memory can be used in a Windows PC, but it could not be used in an Android platform. why?

#pragma kernel KHistogramGather
[numthreads(32, 32, 1)]
void KHistogramGather(uint3 id : SV_DispatchThreadID, uint3 _group_thread_id : SV_GroupThreadID)
{
	const uint thread_id = _group_thread_id.y * 32 + _group_thread_id.x;

	//if (thread_id == 0)
	//{
	//	for (int i = 0; i < 256; i++)
	//	{
	//		gs_histogram[i] = 0;
	//	}
	//}

	// GroupMemoryBarrierWithGroupSync();

	if (id.x < _SourceSize.x && id.y < _SourceSize.y)
	{
		float3 color = saturate(_Source[id.xy].xyz);

		// Convert color to grayscale
		float luminance = dot(color.rgb, float3(0.2125, 0.7154, 0.0721));
		uint idx_l = (uint)round(luminance * 255.0);

		// InterlockedAdd(gs_histogram[idx_l], 1);
		InterlockedAdd(_Histogram[idx_l], 1); // it is an alternative method for using a group shared memory
		_Result[id.xy] = float4(luminance, luminance, luminance, 1);
	}

	// GroupMemoryBarrierWithGroupSync();

	//if (thread_id == 0)
	//{
	//	for (int i = 0; i < 256; i++)
	//	{
	//		InterlockedAdd(_Histogram[i], gs_histogram[i]);
	//	}
	//}
}

単にグレースケール化するだけじゃなくて、輝度ヒストグラムまで構築するスグレモノなのですが
残念なことに groupshared メモリを使った高速化がAndroid上だけ、できないことが分かったのでコメントアウトしています。

シェーダーの利用はスクリプト側で行いますので、利用スクリプトの方も以下に載せます。

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class DeviceCameraBehaviour : MonoBehaviour {

    public GameObject resultCube;
    public ComputeShader RGBA2GRAY;
    public Texture2D Source;
    public RenderTexture Result;

    int _kernelIndex = 0;
    ComputeBuffer _HistogramBuffer = null;
    uint[] _emptyBuffer = new uint[256];
    uint[] _histogramBuffer = new uint[256];

    // Use this for initialization
    void Start () {

        Result = new RenderTexture(Source.width, Source.height, 0, RenderTextureFormat.ARGB32);
        Result.enableRandomWrite = true;
        Result.Create();

        _HistogramBuffer = new ComputeBuffer(256, sizeof(uint));
        for (int i = 0; i < _emptyBuffer.Length; i++)
        {
            _emptyBuffer[i] = 0;
        }
        _HistogramBuffer.SetData(_emptyBuffer);

        _kernelIndex = RGBA2GRAY.FindKernel("KHistogramGather");
        RGBA2GRAY.SetBuffer(_kernelIndex, "_Histogram", _HistogramBuffer);
    }
	
	// Update is called once per frame
	void Update () {
        if (Input.GetKeyDown(KeyCode.Space))
        {
            DoComputeShader();

            Debug.Log("Time = " + Time.time.ToString("000.0000"));
            WebCamDevice[] devices = WebCamTexture.devices;
            for (int i = 0; i < devices.Length; i++)
                Debug.Log(devices[i].name);
        }
        for (int i = 0; i < Input.touchCount; ++i)
        {
            if (Input.GetTouch(i).phase == TouchPhase.Began)
            {

                // Construct a ray from the current touch coordinates
                Ray ray = Camera.main.ScreenPointToRay(Input.GetTouch(i).position);
                // Create a particle if hit
                if (Physics.Raycast(ray))
                {
                    DoComputeShader();

                    Debug.Log("Time = " + Time.time.ToString("000.0000"));
                    WebCamDevice[] devices = WebCamTexture.devices;
                    for (int j = 0; j < devices.Length; j++)
                        Debug.Log(devices[j].name);
                }
            }
        }
    }

    private void DoComputeShader()
    {
        RGBA2GRAY.SetTexture(_kernelIndex, "_Source", Source);
        RGBA2GRAY.SetTexture(_kernelIndex, "_Result", Result);
        RGBA2GRAY.SetVector("_SourceSize", new Vector2(Source.width, Source.height));
        _HistogramBuffer.SetData(_emptyBuffer);
        RGBA2GRAY.Dispatch(_kernelIndex, Mathf.CeilToInt(Source.width / 32f), Mathf.CeilToInt(Source.height / 32f), 1);
        _HistogramBuffer.GetData(_histogramBuffer);
        uint total = 0;
        for (int i = 0; i < _histogramBuffer.Length; i++)
        {
            total += _histogramBuffer[i];
        }
        Debug.Log("histogram total = " + total);
        resultCube.GetComponent<Renderer>().material.mainTexture = Result;
    }
}

簡単にコードを解説してみますと

RGBA2GRAY というのが、Compute Shader で、先ほど示したコードで作ったシェーダーを Unity のエディタから紐づけます。
あとは、画像データを Source テクスチャとして用意して渡したら Disptach にて 32 分の 1 の幅と高さのスレッドグループ数を指定し実行しています。
なぜ 32 分の 1 にしなければいけないかというと、Compute Shader 側で[numthreads(32, 32, 1)]と、32 x 32のタイル領域ごとに処理をすると宣言しているからです。
Dispach で縦横のタイル数を指定してあげることで全画素についての処理が完成するというイメージを持っていただけると、理解しやすいかなって思います。

最後に Andorid のビルド設定ですが、以下のようにして、端末は2017年1月に購入した Phab 2 Pro で動作することを確認しました。(自身の携帯である Xperia Z3 では Compute Shader は動かなかった…)

f:id:simplestar_tech:20170326174354j:plain

今回の記事は以下のページを参考に作成しました。

まず Compute Shader を一度も触ったことがない人は、こちらの記事を参考にするとサクッとPC上で動作確認までできるようになると思います。
[Unity] UnityでComputeShaderを使う解説をしているページを訳してみた その2 - Qiita

次にテクスチャを入力に、Compute Shader で何かしら画素値を処理してテクスチャを出力する場合については次の記事を参考にします。
DirectCompute tutorial for Unity 3: Textures | Cheney Shen

ヒストグラムの計算については、次のプロジェクトの、このプルリクエストを参考にしました。(ヒストグラムなのに並列計算で高速化!?できないと思っていたので驚きましたよ、考えた人グッジョブです!)
Unity-Technologies / cinematic-image-effects / Pull request #31: [tcg] Histogram compute shader optimizations (~4x) — Bitbucket

Unity:uGUIにDebug.Logの内容を表示する方法1

UnityでDebug.Logした内容を、uGUIのUI画面でゲーム実行中に確認したい場合があります。
今回はその要望に応える形で次のようにログが流れる仕組みを作ったので、今後自分が再利用するために公開します。

f:id:simplestar_tech:20170319195427j:plain

Debug.Log をハンドリングするコードは次の通り

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.UI;

public class UILogBehaviour : MonoBehaviour {

    public Text TextPrefab;
    private RectTransform _myPanel;

	// Use this for initialization
	void Start () {
        _myPanel = GetComponent<RectTransform>();
    }
	
	// Update is called once per frame
	void Update () {
		
	}

    void OnEnable()
    {
        Application.logMessageReceived += Log;
    }

    void OnDisable()
    {
        Application.logMessageReceived -= Log;
    }

    public void Log(string logString, string stackTrace, LogType type)
    {
        Text logLine = Instantiate(TextPrefab, Vector3.zero, Quaternion.identity, _myPanel);
        logLine.name = "LogLine";
        logLine.text = logString;

        if (20 < _myPanel.transform.childCount)
        {
            Transform child = _myPanel.GetChild(1);
            GameObject.Destroy(child.gameObject);
        }
    }
}

使い方:このスクリプトを VerticalLayoutGroup コンポーネントを付けた Panel に追加して、最初の子要素の Text を TextPrefab に渡します。

すると 20 件以上ログが流れたら、古いログのラインから順番に消えていきます。
テスト用に Space キーを押したらログが流れる仕組みを次のように書いて、テストしました。

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class DeviceCameraBehaviour : MonoBehaviour {

	// Use this for initialization
	void Start () {
		
	}
	
	// Update is called once per frame
	void Update () {
        if (Input.GetKeyDown(KeyCode.Space))
        {
            Debug.Log("Time = " + Time.time.ToString("000.0000"));
        }
	}
}

Android などの実機で、デバイスの設定がどうなっているかなどを簡易にに確認する方法として使えると思い、作ってみました。
以上です。