simplestarの技術ブログ

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

英語の勉強

今日はこの記事を訳してみます

www.theatlantic.com

Title:問題を起こす人
同一の特徴、それはドナルドトランプJr氏を彼の父親の政治的な相続人とする、は、彼を最新のスキャンダルの真っただ中に置いてしまうのかもしれない。

もう何十年も前、トランプ氏の長男がまだ少年だったころ、バーバラ記者(有名な女性記者)がドナルドトランプ氏を彼の家族の同席のもと取材したことがあった。
その時彼女は不動産王の彼に対し、どの子が家族の中で問題を起こす子だと思っていますか?と尋ねました。

トランプ氏はためらいなくすぐに"ドン"であると受け答えたと、当のドン自身の話によって明かされる。
ドンことドナルドトランプJr.氏(39歳のビジネスマン)は嬉々として語ってくれた。

ドン Jr. 氏は私が昨年の人物紹介で取材に訪れたときに、ある物語を歯を見せるような笑顔で語ってくれた。
"私は三人の中で最もヤンチャだったんですよ"と自身の兄妹について語った。
"いつも成績は良く、いい行いもしていましたが、その分たくさん遊んでいましたね"

元気よく、強い意志を持って、危機負担もある、そんな彼の特性はトランプの子供たちの中でドンJr.氏を最も目立つ存在としていたわけです、トランプ政権中は。
しかし、今週の暴露―2016年6月のメールの交換が火曜日に公開され、そのメールにはドンJr.氏が協力的なロシア政府からの政治的な援助運動の可能性が示された内容に肯定的な返事"そいつはいい!"をしていたという―これがこれまでの彼の特性を今までとは異なる場所へと投げ込んだことになった。
もう一度、ドンJr.氏は彼の父親にとって問題を起こす子だったわけだが、今回は遊びやゲームの度をはるかに越えた問題である。



…途中でやめられなくなって、30分以上作業してしまいました
ちょっと、政治的な単語知らな過ぎて読めない…
あと、Quality って今回は人の性格とか、特性って感じでとらえてよかったのかな?


うーむ、ひとまずの目標は30分で記事の最後までスラスラ訳せるようになることかな
明日も頑張ります。

追記:
って、あれ!元記事を読み直したら
以下の表現が変化してた
by Barbara Walters with his family present.

by Barbara Walters, along with his family.
確かに、下の方が誤解なく、"家族と一緒に"が伝わりますね。
修正前は、"家族の存在とともに"というより、場合によっては"家族への贈り物とともに"とも受け取れちゃうかも

英語の勉強

とりあえず一日30分は英語の記事を読んでみようと思います。

次の記事を読んでみた(30分だけ)
www.theatlantic.com

ちゃんとした訳ではなく、私が英語を読んで感じたことなので、信じないように…

タイトル:それらがすべて本当だとしたら何だというのか?
ここ数か月、トランプ大統領とその側近は主張を続けてきた、内容としてはいかなる提案の談合、相手がロシアで、そしてそれが偽物であると
しかし、もしトランプ政策が意欲的に情報をロシア政府から受け取りたいというスタンスだったなら、ほかに何か起こりえたのだろうか?

ここ数か月、噂、暗示、そして出どころの怪しい主張、それらは談合に関するもので、相手はトランプ政策とトランプ陣営とロシア政府というものがワシントンの周りを渦巻いていたわけで、時々大きな洪水が吹き出すように、そして常にゆっくり流れる小川のように
何回もドナルドトランプ氏とその味方は否定してきた
彼ら曰、そこに接触はなかった、選挙の前には



うーん、こりゃひどい
全然進まない…
もっとスラスラ読めるようになりたいですね
明日も頑張ります

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 の方でも、リアルタイムカメラ画像処理できることを確認しました。