AR

Unity5.4でComputeShader(HLSL)のコンパイルに失敗する

Boltzmannです.
ラベリングの並列実装を試みていたところ, あるアルゴリズムを実装したところでUnityに異変が生じました.
具体的には, エディタ(VisualStudio2015を使用)でComputeShader(HLSL)のスクリプトを編集した後, Unityに移行するとUnityがフリーズし, その後PCもフリーズするという状態です.
これが二回ほど起きて, さすがにこれはまずいと思い原因を調べ始めました.

結論から言えば, もっとも考えられる原因はこれでした.
[COMPUTE] HLSLCC TIMES OUT COMPILING COMPLEX SHADERS (BAD PERF ON COMPLEX SHADERS) -Unity Issue Tracker

もともと5.3で作業していて, このページによれば5.4.0で解消されたということで5.4.2にアップグレードして再度試しましたがやはりだめでした.

具体的な症状は, 「Unity Shader Compilerが大量にメモリを食ってからタイムアウトで終了する」というものです.
問題が発生した環境は以下の通りです.
Unity5.3.?(ちゃんとメモしなかった), 5.4.2
CPU:Intel Corei5-5200U
メモリ:8GB
(GPU:Intel HD Graphics 5500)

次にShaderComplierのタスクを以下に載せます. だいぶメモリを食っていることがわかるかと思います.
shadercompilerprofile3

これが5~10分ほど続いて, 最終的にタイムアウトで終了します. タイムアウトで終了するので, シェーダーはコンパイルできていません. この状態でPlayしても以下のようにカーネル関数を見つけることができないというエラーが出ます. (もちろんスクリプト上は存在するカーネル関数です. )
shadercompilertimeout

最後に, この問題が発生した際のスクリプトを記載します. スクリプト全体ではありませんが, この部分をコメントアウトすれば正常にコンパイルされる, という部分を記載します.
なお, 以下のスクリプトの参考文献を次に示します.

柴田直樹, 山本眞也. SumiTag : あまり目立たないARマーカーとGPGPUを利用した読み取り方法. 研究報告マルチメディア通信と分散処理(DPS). 2011, 2011-DPS-149巻, 7号

<br />
    bin = lerp(0, 1, step((float) grayscale, 0.4));<br />
    Input.GetDimensions(width, height);<br />
    index = width * thid.y + thid.x;</p>
<p>    LabelBuffer[index] = bin * index;</p>
<p>    for (i = 0; i &lt; 1000; i++)<br />
    {<br />
        AllMemoryBarrierWithGroupSync();</p>
<p>        label = LabelBuffer[index];<br />
        orig = label;</p>
<p>        if (bin != 0)<br />
        {<br />
            if ((thid.x &gt; 0) &amp;&amp; (label &lt; LabelBuffer[index &#8211; 1]))<br />
                label = LabelBuffer[index &#8211; 1];<br />
            if ((thid.x &lt; width &#8211; 1) &amp;&amp; (label &lt; LabelBuffer[index + 1]))<br />
                label = LabelBuffer[index + 1];<br />
            if ((thid.y &gt; 0) &amp;&amp; (label &lt; LabelBuffer[index &#8211; width]))<br />
                label = LabelBuffer[index &#8211; width];<br />
            if ((thid.y &lt; height &#8211; 1) &amp;&amp; (label &lt; LabelBuffer[index + width]))<br />
                label = LabelBuffer[index + width];</p>
<p>            label = LabelBuffer[LabelBuffer[LabelBuffer[LabelBuffer[label]]]];</p>
<p>            if (label != orig)<br />
            {<br />
                InterlockedMax(LabelBuffer[orig], label);<br />
                InterlockedMax(LabelBuffer[index], label);<br />
            }<br />
        }<br />
    }</p>
<p>    AllMemoryBarrierWithGroupSync();</p>
<p>    switch (LabelBuffer[index] % 20)<br />
    {<br />
        case 0:<br />
            color = float3(1, 0, 0);<br />
            break;<br />
        case 1:<br />
            color = float3(0, 1, 0);<br />
            break;<br />
        case 2:<br />
            color = float3(0, 0, 1);<br />
            break;<br />
        case 3:<br />
            color = float3(1, 1, 0);<br />
            break;<br />
        case 4:<br />
            color = float3(1, 0, 1);<br />
            break;<br />
        case 5:<br />
            color = float3(0, 1, 1);<br />
            break;<br />
        case 6:<br />
            color = float3(0.1, 0.2, 0.5);<br />
            break;<br />
        case 7:<br />
            color = float3(0.2, 0.5, 0.1);<br />
            break;<br />
        case 8:<br />
            color = float3(0.5, 0.1, 0.2);<br />
            break;<br />
        case 9:<br />
            color = float3(0.5, 0.2, 0.1);<br />
            break;<br />
        case 10:<br />
            color = float3(0.1, 0.1, 0.1);<br />
            break;<br />
        case 11:<br />
            color = float3(0.2, 0.2, 0.2);<br />
            break;<br />
        case 12:<br />
            color = float3(0.23, 0.56, 0);<br />
            break;<br />
        case 13:<br />
            color = float3(0.9, 0.2, 0.5);<br />
            break;<br />
        case 14:<br />
            color = float3(0.4, 0.8, 0.2);<br />
            break;<br />
        case 15:<br />
            color = float3(0.4, 0.3, 0.6);<br />
            break;<br />
        case 16:<br />
            color = float3(0.1, 0.3, 0.2);<br />
            break;<br />
        case 17:<br />
            color = float3(0.5, 0.5, 0.2);<br />
            break;<br />
        case 18:<br />
            color = float3(0.5, 0.1, 0.2);<br />
            break;<br />
        case 19:<br />
            color = float3(0.8, 0.8, 0.3);<br />
            break;<br />
        default:<br />
            color = float3(1, 1, 1);<br />
            break;<br />
    }</p>
<p>    Result[thid.xy] = float4(color, bin);<br />

ちなみにswitch文のブロックのみをコメントアウトしても問題は解消されません.
なんとなく, スクリプトを編集しているときの感覚として, InterlockedMaxをfor文の中に組み込んだのがとどめの一撃になった気がします.

とりあえず実装を見直しますが, 実装を誤るたびにUnityが大量にメモリを食ってフリーズするのは勘弁してほしいです.

<2017/2/22 追記>
原因不明として放置してあったのでとりあえずあの後の調査で考えられる原因を書いておこうと思います.
原因として有力なのはコンパイラによるfor文の自動展開です.
GPUデバイスコードは命令数が多いほど並列性が高くなったりするのでコンパイラが可能な範囲で積極的にfor文を展開することがあります.
今回の場合は展開した後のコードがあまりに長くなったためオーバーフローしてフリーズしたのではないかと考えられます.
調べていないので何とも言えませんが, こういった自動展開に関しては必ずそれを制御する#pragmaがあるはずなので, それを使うことで解決できるのではないかと思います.
他の解決策としては, 展開できない程度にコードを複雑にする, というものがあげられます.
上記のスクリプトに関しても, 終了条件をfor (i = 0; i < 1000; i++)のようにべた書きするのではなく, tmp=1000のようにしてからfor (i = 0; i < tmp; i++)のようにすることで回復したように思います.
<追記終わり>

UnityでRawImageを使ってARをする際の備忘録

・Canvasを使うと自動的にカメラの描画範囲にフィットする上、親子関係を作って相対位置を固定する必要がないのでPlaneにはっつけつより楽
・Kinect等と連携してDepthMappingなどをする際にはRenderingPipeLineに入らないといけない気がする

どこをいじったかいまいち覚えていないので全部載っけておく

CanvasのInspector
ARCanvas

RawImageのInspector
ARRawImage

Background.cs

今回の実装はAndroidようだったので、取得したWebcamTextureを90°回転させている
CanvasのRenderCameraに、Canvasを表示させたいカメラを指定するのを忘れないように
UIのRectTransformはTransformに相当する
WebCamTexture.Play()をしないとテクスチャに反映されないため, wtex.widthやwtex.heightが機能しないのでさっさとPlay()しておくこと
今回はゲーム開始時にARモードをいきなり使うわけではないので, 一通りの設定が終わった後にStop()させている.
本当はこのセットアップは画面遷移時にやるべき.

以上.

UnityのInput.gyro.attitudeとInput.accelerationを使ってARをやる際の備忘録

前提:ジャイロセンサによってスマートフォンとカメラを同期、その上で加速度の方向をUnity内部と現実世界で合わせる

<br />
    void Update()<br />
    {<br />
        Quaternion orient = Input.gyro.attitude;<br />
        orient.x *= -1;<br />
        orient.y *= -1;<br />
        Quaternion sub = Quaternion.Euler(90, 0, 0);<br />
        transform.rotation = sub * orient;</p>
<p>        Vector3 accelVecCam = new Vector3( Input.acceleration.x, Input.acceleration.y, -Input.acceleration.z);<br />
        transform.GetChild(0).localPosition = accelVecCam*20; //確認用. 適当なGameObjectをカメラの子にする<br />
    }<br />

accelVecCamはカメラ座標系での加速度ベクトルであることに注意.
World座標にしたいときはtransform.TransformVectorで変換すること.

Input.gyro.attitudeのほう
参考元:http://vr-cto.hateblo.jp/entry/2016/05/02/070000
参考元の参考元:http://qiita.com/fuqunaga/items/b1a3e38af71f062f0781

Input.accelerationのほう
参考元:https://developer.android.com/reference/android/hardware/SensorEvent.html

ラベリング処理の並列化について少し考察する

Boltzmanです.
以前記事に書いた通り, ラベリング処理の実装が難航しています.
正直全然進捗がありませんが, やや行き詰ってきたところもあるので, ここで一度風呂場でアヒルに向かって話すつもりで記事にしようかと思います.
ちなみに”Rubber duck debugging”という言葉はchokopan氏に教わりました.

続きを読む

二値化の実装が完了しました

Boltzmanです.
実に2か月ぶりの進展です.
正直学校の授業に追われて全くと言っていいほど手を付けていなかった今プロジェクトですが、
夏休みに入ると同時に少しの進展がありました.

とはいっても実は前回のDebugWindowを用いたデバッグは1週間くらいで終わっていました.
デバッグがすぐに終わったので, 閾値の平滑化まで終わらせてから記事を書こうと思っていたら2か月がたっていました.
v0.3.1

DebugWindowに関して変更は
・ヒストグラムに(平滑化前の)閾値が赤線で表示されるようになった
・平滑化時の, 中央のピクセルに対する周りの重みをスライダーで調整できるようになった
・選択した領域は赤と緑で二値化するようにした
の3点です.
DebugWindow

平滑化の効果は以下の通りです.
平滑化

局所的に閾値を計算した後そのまま用いると, 上の図のように領域内に同じ色(濃度の近い色)しか存在しない場合でも
無理やり分けようとしてノイズのようなものが入ってしまいます.
平滑化することによってそれを軽減することができます.

今回平滑化の実装でネックだったのが, HLSLのRWTexture2DがRead&Writeではなかったという点です.
RWTeture2Dを読み取るには, 第一にTexture2Dのように, 各ピクセルに割り当てられるデータ配列が一次元である必要があります.
これだけであれば, まだ問題ないのですが, 結局私は上のように宣言してもうまく読み込めなかったため, 諦めてカーネル関数を二つにて対処しました.
多分RandomWriteの略だったんでしょう.

次はラベリングの実装ですが, ラベリングは致命的に並列計算効率が悪いので, C#側で実装するかもしれません.
ただC#側で実装するにしても, ComputeShaderで処理した画像はRenderTextureに入っているため, そのままでは色を抽出できず厄介です.

何とか頑張ります.

DebugWindowを作成しました

Boltzmanです. 前回記事更新から大分経ったかと思いましたが大体1か月くらいですね.
色々と忙しくて今回もなかなか進捗が上がりませんでした.

この一か月ずっと何をしていたのかというと, 「Group(Block)ごとに閾値を決定するスクリプトのデバッグをするためのウィンドウのデバッグ」です.
ん?本末転倒かな?? このままさまよい続けて私はどこへ行くのでしょうか…
まあしかし(デバッグをするためのウィンドウの)デバッグが終わったので心置きなくデバッグができますね!!

DebugWindow2

DebugWindow_
(× GropuID -> ○ GroupID)

何写したらいいかわからなかったのでとりあえずJ.J.Sakurai映しておきました.
二値化画像上で表示されている赤い四角形の部分が現在選択されているGroup(Block)です. このGroupのヒストグラムや閾値を表示します.
Groupの選択はマウスクリックでできるようになっています.
ヒストグラムの表示領域は256px*256pxでとってあるので, 1pxあたりに一つの階級という形になってしまい, ヒストグラムが毛みたいになっています.
閾値の取得はこれからなのでまだ確認できません. まあすぐに実装して見せるさ…

作ってみて気づいたことは如何にDebug.Logが重いかということです. いままでてっきりCPU-GPU間のデータ移動コストのせいで(デバッグ中に)遅くなっているのかと思っていましたが, そうではなくDebug.Logが重かったようです.

スクリプトはEditorWindowを使っていますが, 依存関係の問題で(?)Monobehaviour側からEditorWindow側を参照できなかったため, EditorWindowのスクリプトの他にDebugUtilsという, デバッグに必要なデータを管理するクラスを書いて両者を仲介させました. (これに気づくのに少々時間がかかりました. )

デバッグ用のウィンドウを作った一番のモチベーションは, もう少し感覚的にわかりやすい形で結果を確認したい, ということです.
どうしても結果の画像やヒストグラムの数値だけだと, 「なんとなくうまくいっていそう」という以上に確信が持てず, あとになってからバグが見つかるということが多いので, それが少しでも解消すればと思い作りました.

EditorWindow書いてて楽しかったです. 君も自分だけのデバッグ環境を構築しよう!!!

次回はこれを使って心置きなくデバッグします.

大津の手法による閾値自動決定アルゴリズムを用いた二値化をGPUを用いて実装しました

Boltzmanです. 前回の記事から1か月以上が経過しました. 極めて難航しましたが, 進捗が上がりましたので記事にします. (遅すぎ)

前回, 大津の手法による自動閾値決定2値化アルゴリズムをCPUで実装してあまりにも重くなったため, GPUを用いて, 具体的にはComputeShaderを用いて実装しました.
実行環境は私のノートPC(Intel Core i5-5200U/RAM 4.00GB/Intel HD Graphics 5500)です. 決して良いとは言えないスペック.
実行結果は
CPUのみで実装   (v0.1.0):MonoBehaviourが平均330msくらい
GPUを取り入れて実装(v0.2.3):Gfx.WaitForPresentが平均16~20msくらい?(ガタつきが多い)

プロファイル比較

改めてこうしてみるとCPUのみで実装したほうが理不尽に重いので, どこかでコーディングミスをしている可能性が大いにあります.

また大津の手法による自動閾値決定により2値化画像がどのように変化するかが以下です.

大津の手法結果

上を見るとわかる通り, マーカーは濃淡の変化が顕著になるように作られているためか, あまり見栄えは変わりません.
しかしその他の背景画像は大津の手法を実装したときのほうが見やすくなっています.

続きを読む

GPUで自動二値化を実装中です

Boltzmanです.
まだ完成してはいませんが, 中間報告です.
現在, GPUを用いた大津の手法による自動閾値決定二値化アルゴリズムを実装中です.
GPUプログラミングは初めてなので, どのようなコードがGPUと相性がいいのか模索しながらの実装となっています.
とりあえず現時点で一応動きはしたものを以下に載せます.

まだまだ改良の余地が多く, 現在全体を改めて見直している最中です. 特に, 上のコードではnumやsumを愚直に計算してしまっているため, この辺りを何とかしたいです.
また本当はif文を使いたくなかったのですが, m1, m2を計算する際にゼロ割が生じるケースがあり, それによって分離度がNaNになってしまうと最大値を求める際に支障をきたすため, NaNを消すためにif文を導入しました. 他に良い方法はないだろうか…

最大値の算出はmax.computeのようになっていますが, こちらもやはり全体を書き換えている最中です.
上のコードでは, 2n番目と2n+1番目を比較して, 大きいほうをn番目に格納し, それを繰り返すことで最大値がn=0に格納されるようにしています.
saidaiti
このコードで初めてHLSLの同期を使いました. 同期って他のやつの進行を待つことになるので, 一番遅い奴に足並みがそろうことになります. なのであまり使わないほうがいいのかと思いましたが, コードを調べてみたら割とがっつり使っているものが多かったので気にしなくてもいいのでしょうか…

上のコードを見直すと同時に, ラベリング処理の並列化についても少し調べています.
論文を見てみるとCPUを並列させたほうが速いらしいです. しかしながら, 今までの処理をGPUでやっているため, ラベリングをCPUで計算するとなるとデータをCPUへ移さなければならず, そのコストを考えて総合的に評価するとGPUのほうが速くなる可能性があります.
一度GPUでやると抜け出せなくなりますね…

そういえばGPUは流体力学と相性がいいそうです. 私は詳しくないのですが, 流体に詳しい人から少しだけ話を聞きました. まあ紹介程度に.

引き続き実装頑張ります.

シェーダーを使って二値化を再実装しました

Boltzmanです. 魚眼でやりたいAR, 4回目です.
少しですが進捗が出たので進捗報告です.
以前二値化を実装した際にあまりに重かったので, サークル仲間の助言に従ってシェーダーを使って実装しなおすことにしました.

使ったシェーダーはCompute Shaderです. Standerd Surface Shader, Unlit Shader, Image Effect Shader, Compute Shaderのうちどのシェーダーを使うのか選ぶ時点でかなり難航しました. 様々なサイトをめぐって情報をあさって, 最終的には消去法で選択したような記憶があります.
選んだシェーダーがあっていたのかはわかりませんが, とりあえず実装はできました.
以下にコードを記します. まずはCompute Shaderのほうから

<br />
Texture2D&lt;float4&gt; Input;<br />
float1 thresh;<br />
RWTexture2D&lt;float4&gt; Result;</p>
<p>[numthreads(8,8,1)]<br />
void CSMain (uint3 id : SV_DispatchThreadID)<br />
{<br />
	//二値化<br />
	Result[id.xy] = lerp(float4(1, 1, 1, 0), float4(0, 0, 0, 0),<br />
                      ceil( clamp( -(Input[id.xy].r+Input[id.xy].g+Input[id.xy].b)/3 + thresh, 0, 1) ));<br />
}<br />

(注:ここからなぜか口調が変わりますが書いている人は同じです. 本来なら直すべきですが, どうしても気力がわかないということと, 内容に影響はないということからそのままにします. 申し訳ない. )
今見返すとlerpの1番目と2番目の引数がfloat4であるのに対して3番目の引数がfloat1になっているのでどうして動いているのか不思議である. (自動的に変換しているのだろうか).
(2016/3/5追記:こちらのページに”スカラー昇格”なる記述があり, それによるものかと考えられます)
まずuint3 idは三つの符号なし整数をパラメーターとして持つ変数であるが, このような変数はx,y,zやr,g,bなどを用いて要素を取り出すことができる.(詳細)
これから, Input[id.xy]ではid.xyによってテクスチャの座標を指定して, その座標の色を求めている.(たぶん)
Input[id.xy]の値はfloat4でrgbaらしいので, rgbの和を取って3で割って濃度値を算出し, threshからその値を引くと, threshより大きい値の時は結果が負, threshより小さな値の時は結果が正になる. これをclampすると, 負の値は0に, 正の値は0より大きく1以下の値になる. さらにceilによってこれを0と1に二分する. lerpで線形補間すると, 0のときは後ろの値, 1の時は前の値が採用され, 画像を二値化することができる.
と, 私は思っているのだが, 実際どうなってコードがうまくいっている(ように見える)のかはわからない.

次にC#スクリプト

<br />
using UnityEngine;<br />
using System.Collections;</p>
<p>public class Test : MonoBehaviour {<br />
    public ComputeShader cs;<br />
    public int width = 1080;<br />
    public int height = 720;<br />
    public int FPS = 30;<br />
    [Range(0, 1)]<br />
    public float thresh = 0.5f;</p>
<p>    private RenderTexture text;<br />
    private WebCamTexture wtex;</p>
<p>	// Use this for initialization<br />
	void Start () {<br />
        //WebCam準備<br />
        var devices = WebCamTexture.devices;<br />
        if (devices.Length &gt; 0)<br />
        {<br />
            wtex = new WebCamTexture(width, height, FPS);<br />
            wtex.Play();<br />
        }<br />
        else<br />
        {<br />
            Debug.LogError(&quot;Could not Find WebCam.&quot;);<br />
            return;<br />
        }<br />
        text = new RenderTexture(wtex.width, wtex.height, 0);<br />
        text.enableRandomWrite = true;<br />
        text.Create();<br />
        transform.localScale = new Vector3(width, 1, height);<br />
        transform.GetComponent&lt;Renderer&gt;().material.mainTexture = text;<br />
        cs.SetTexture(0, &quot;Input&quot;, wtex);<br />
        cs.SetFloat(&quot;thresh&quot;, thresh);<br />
        cs.SetTexture(0, &quot;Result&quot;, text);<br />
    }</p>
<p>	// Update is called once per frame<br />
	void Update () {<br />
        if (!SystemInfo.supportsComputeShaders)<br />
        {<br />
            Debug.LogError(&quot;Computer Shader is not Supported.&quot;);<br />
            return;<br />
        }<br />
        cs.Dispatch(0, width, height, 1);<br />
    }<br />
}<br />

cs.Set~で値を設定する. WebcamTextureを入力画像として渡し, 結果をRenderTextureで受け取っている.
cs.DispatchでCompute Shaderを実行する. (Compute Shaderの詳細)
ちゃんとRandomWriteをtrueにしておくことが重要らしいとどこかで読んだ.
(2016/02/29追記:ソースを発見したので載せます->詳細)
途中まで実装したのが今月の頭くらいで, 完成したのが2,3日前とずいぶん間が空いてしまったために, どこから引っ張ってきた情報なのかわからず大変忍びない.
やっていて感じたのは, シェーダー関係は情報を探すのが大変だということ.
シェーダーその物もそうだが, CgやHLSLの書き方も公式のリファレンスに見つけられただけでCやC++、C#と比べるとずっと調べにくかったように思う.
まあ私の見落としが多いだけなのだろうが…

前回のコードが大津の手法まで実装したものであったため, 正確な比較ではありませんが, 前回のコードのMonoBehaviourが300ms以上かかっていたのに対し, 今回のコードではCompute ShaderのDispatshは0.01ms以下となっています. 目に見えてわかるほど大きく改善しました.

次回は(必要があれば)大津の手法をシェーダーで実装することを目標にします.

「大津の手法」を実装しました

どうも, Boltzmanです. 魚眼らしいことがまだ全然できていないAR, 3回目です.

今回の目標は「大津の手法」(判別分析法)による閾値自動決定の実装です. 大津の手法を紹介しているサイトはたくさんありますが, ここでも一通り紹介したいと思います. 正直書くことがない.

大津の手法は「分離度」と呼ばれる値が最大になるような閾値を採用する方法で, 分離度は\[\frac{\sigma_b^2}{\sigma_w^2}\]で表されます. ここで\(\sigma_b\)はクラス間分散, \(\sigma_w\)はクラス内分散です.

 

分散は要素がどれくらい散らばっているかの指標です.

bunsan1

上の図では, 右側のほうが左側に比べて要素が散らばっているため, 右のほうが分散の値は大きくなります. 分散\(\sigma\)は, 要素\(\boldsymbol{x}_i\)とそれらの平均\(\boldsymbol{\mu}\)を用いて\begin{equation}\sigma^2=\sum_{i}(\boldsymbol{x}_i-\boldsymbol{\mu})^2\end{equation}と表されます. つまり, 分散は要素の平均と各要素との距離の和と考えることができます.

 

ここで要素を複数のクラスに分けます.

bunsan2

上のデータを見るとなんとなく二つに分けられそうな気がします.
bunsan3

こんな風に. このようにして括った二つのクラスについて, 分散を考えます.

まずクラス内分散は各クラス内の要素の分散です. つまり, \begin{equation}\sigma_w^2=\sum_{k}\sum_{\boldsymbol{x}_i\in\chi_k}\frac{\omega_k}{\omega}(\boldsymbol{x}_i-\boldsymbol{\mu}_k)^2\end{equation}こんな感じですかね. \(\boldsymbol{\mu}_k\)はクラス内での平均, \(\omega_k\)はクラス内の要素数, \(\omega\)は全体の要素数を表します.

クラス間分散は, 各クラスがどれくらい散らばっているかを示します. つまり, \begin{equation}\sigma_b=\sum_{k}\frac{\omega_k}{\omega}(\boldsymbol{\mu}-\boldsymbol{\mu}_k)^2\end{equation}こうですかね. \(\boldsymbol{\mu}\)は全体の平均を表します.

なんか, ちゃんと統計やらなきゃだめだなって思います. 今回は実装だけなので原理はかなりスルーしてます. 適当なこと書いてたらすいません.

で, 分離度はクラス内分散\(\omega_w\)が小さいほど, クラス間分散\(\omega_b\)が大きいほど大きくなります. つまり, 各クラスを見たときに要素が集まっているほど, クラス同士を見たときに各クラスが散らばっているほど大きくなります.

ここで, クラス内分散\(\sigma_w\)とクラス間分散\(\sigma_b\), それに全分散\(\sigma\)について, 以下の関係式を導入します. \[\sigma^2=\sigma_w^2+\sigma_b^2\]この証明ができなくて困っています. どのサイト見ても自明っぽく書いてあるので, たぶん私が何か重要な公式を知らないんじゃないだろうか. この関係式を認めるとすると, 分離度は以下のように書きなおせます. \[\frac{\sigma_b^2}{\sigma^2-\sigma_b^2}\]よって, クラス間分散が最大になる時, 分離度は最大になります. なぜクラス間分散のほうを残したのかというと, クラス内分散のほうは計算の時に各クラスでの分散を求める必要があって面倒だからです.

 

今回はこれを使って, 濃淡のヒストグラムを二つに分割します.

histo

閾値でヒストグラムを二つのクラスに分けた時の分離度を考えます. クラス間分散は上の定義から, \begin{equation}\sigma_b^2=\frac{\omega_1(\mu_1-\mu)^2+\omega_2(\mu_2-\mu)^2}{\omega}\end{equation}ここで, \[\mu=\frac{\omega_1\mu_1+\omega_2\mu_2}{\omega}\]より, \begin{eqnarray}\mu_1-\mu&=&\mu_1-\frac{\omega_1\mu_1+\omega_2\mu_2}{\omega}\\&=&\frac{\omega_1(\mu_1-\mu_2)}{\omega}\end{eqnarray}となるので, \begin{equation}\sigma_b^2=\frac{\omega_1\omega_2(\mu_1-\mu_2)^2}{\omega^2}\end{equation}と書き直せます. \(\omega\)はクラスの分け方によらず一定なので, \(\omega_1\omega_2(\mu_1-\mu_2)^2\)が最大となる時\(\sigma_b\)も最大となり, 分離度が最大になります.

上の図で, ヒストグラムの式を\[y=f(x)(y\in\mathbb{Z}:要素数, x\in\mathbb{Z}:グレースケールの値(0\leq x\leq 255))\]とすると, 平均は\begin{equation}\mu_k=\frac{1}{\omega_k}\sum_{x_i\in \chi_k}x_if(x_i)\end{equation}であたえられます. (平均の定義を勘違いしてたのは私だけでいい)
結局コードは以下のようになりました.

 

</p>
<p>public int getThresh()<br />
{<br />
    int th = 0; //ループ用の閾値<br />
    double num1, num2; //画素数<br />
    double sum1, sum2; //和<br />
    double m1, m2; //平均<br />
    double va_b; //クラス間分散の二乗<br />
    double max = 0; //va_bの最大値を保存しておく変数</p>
<p>    //ヒストグラムを生成していなければ弾く<br />
    if (hist == null) return -1;</p>
<p>    //(th=0)の処理<br />
    num1 = hist[0];<br />
    sum1 = hist[0] * 0;<br />
    m1 = num1 != 0 ? sum1 / num1 : 0;</p>
<p>    num2 = 0;<br />
    sum2 = 0;<br />
    m2 = 0;<br />
    for (int i = 1; i &amp;lt; 256; i++)<br />
    {<br />
        num2 += hist[i];<br />
        sum2 += hist[i] * i;<br />
    }<br />
    m2 = num2 != 0 ? sum2 / num2 : 0;</p>
<p>    va_b = (num1 * num2 * (m1 &#8211; m2) * (m1 &#8211; m2));<br />
    max = va_b;<br />
    //Debug.Log(&quot;0:&quot;+va_b);</p>
<p>    //(th=1)以降の処理<br />
    for (th = 1; th &amp;lt; 255; th++)<br />
    {<br />
        num1 += hist[th];<br />
        num2 -= hist[th];<br />
        sum1 += hist[th] * th;<br />
        sum2 -= hist[th] * th;<br />
        m1 = num1 != 0 ? sum1 / num1 : 0;<br />
        m2 = num2 != 0 ? sum2 / num2 : 0;<br />
        va_b = (num1 * num2 * (m1 &#8211; m2) * (m1 &#8211; m2));</p>
<p>        //Debug.Log(th+&quot;:&quot;+va_b);<br />
        if (max &lt; va_b)<br />
        {<br />
            thresh = th;<br />
            max = va_b;<br />
        }<br />
    }</p>
<p>    return thresh;<br />
}</p>
<p>

これで一応動いていますが, 原理が完全には腑に落ちていないせいか, どうにも達成感がありません.
今回はここまでです. 次回はシェーダーを勉強しないといけないと思うので, 遅れると思います. では.