無題

プログラミングとか研究とか

2017

特に脈絡のない記録.

2017年

1月 ~ 3月

本職のゲーマーとしてはOverwatchの国内最強決定戦という大会に出場した. Twitterで集めたメンバー参加したので,全体練習がほとんどしていない状態で本番に挑んだが, 予選決勝まで進むことが出来て中々いい結果だった.

RNG5に殲滅されたMotley Ordnancesでした.
“PlayStation 4 オーバーウォッチ国内最強決定戦”代表決定戦の模様を振り返る【動画あり】 - ファミ通エクストリームエッジ|ファミ通.com

  

副業の大学では12月に卒論のテーマが変わり, 4月から違う研究室に移ることもあってやる気が最低値で卒業諦めて就職しようかと何度も考えたが, 何とか1ヶ月で実験から卒論をほぼ完成まで書き上げた. 今ちょうど卒論やってる時期だと思うので折れそうになったときは僕で何とかなったので大丈夫だと自信を持って欲しい.

何年かぶりに家族全員で旅行に行ってエモだった.

2月

卒論の提出や,口頭試問があったはず. 疲弊しすぎて記録も記憶もない.

3月

この月で一番面白かったイベントは約2年半ぐらい付き合ってた彼女と RADWINPSのライブに行って,甘い歌を一緒に聞いた後に別れたことですね.
新しいラボに移動するための後処理などなど. Life is Tech!のHAL大阪キャンプでチームを持ってメンターしたりも. 新年度からどういう立ち回りをしようか模索していた. その残骸
yoshi-sa.hatenablog.com

4月

大学院への進学&新しい研究室
新しいテーマを決めるためにサーベイといくつか候補を練っていた.
それ以外では研究に絡めていけそうな強化学習やNN系の勉強もしていた.
Gの人が大学にインターンの募集にきていて,話をしたときに, 自分がエンジニアになりたいのか,このまま研究者になるのか,それ意外なのかいよいよ決定すべきことに気づいた. その後それが長い間精神疲労の原因になった.

5月

研究の方ではテーマがある程度固まったので実験環境の実装に入っていた.
他の時間では前項の理由により迷走してひたすら作ったり本読んだりしていた.
yoshi-sa.hatenablog.com

6月

迷走続き. yoshi-sa.hatenablog.com このとき書いた記事,DrawMeshInstancedIndirectって調べると公式リファレンスの次に出てくる様になっていて少し嬉しい.

7月

友達の紹介で京都のTomarubaというベンチャー企業iOSエンジニアとしてバイトを始めた. だいたい家とかラボで作業のリモートなのに時給2000円近くで今も楽しく働かせてもらっている. 社員さんたちもとてもよくしてくれて良い職場です.
この月は調子が戻ってきて,勉強していたことも色々繋がってきていい傾向が出てきた.
後教養をつけるためにCTFをやってみたりもした.

その他 yoshi-sa.hatenablog.com

8月

記憶障害, ひたすら何か作っていた気がする.

9月

ラボの色んな問題に衝突した. 一時は本気で大学院をやめようと思って就活していた.

10月

何かきっかけがあって研究をがっつりやりだした.10月の末には予備実験も終えた.

彼と同じ様な理由で競プロをちゃんと始めた.
shopetan.hatenablog.com でも実際やってるとそこを評価基準にするのは妥当だなぁと感じてきた.

11月

予備実験が終わってデータ解析用のプログラムを組んだり(一応OSSにしてある),
今のテーマのその先ともう1つのテーマに関して考える時間を多くとった.
この辺りから,自分の分野のジャーナルはアラートつけてアブストだけでも読むようになった.
変な焦りみたいなのがなくなって生産性がやっと戻った月だった.
もうすぐモンハンが発売されることもあって,PS4をPS4Proに買い換えたり,{i}の人たちと旅行にいけたりした.
内定もらっている会社だけじゃなくて他も見たほうがいいかーってことで逆求人とかに参加した.
2年前に作ったプロジェクションマッピングでなんかのコンペで最優秀賞とって20万円をいただいた.

12月

競プロ,10月から少しづつしかやってなかったので,やっとこの月でABCは全完できるようになってきた. 就活は内定2社でもう十分だと思ったし疲弊するのが嫌だったので,話は聞くけど余程じゃないとそれ以上はいかないというスタンスに切り替えた. 12月の後半は7日間で20人の実験をした(一人あたり2時間半ぐらいかかる).
色々サーベイしてて筋の良さそうな次のテーマを指導教員の方と話していてとても楽しかった.

2018年

研究室が変わって,色々あったけど,やっと論文が出せそうで,次もいい感じになりそうなので, 2018年も研究はいい精神状態で楽しんでいけそう.
エンジニアバイトは今のTomarubaから抜ける予定はないけれど他にもお誘いいただいてるので, 疲弊しない範囲でやっていきたい.

来年大事なのが進路選択で,就職するか博士に進むかですね. とりあえずやれるだけ両方やってみて選択します.

2018年もよろしくお願いいたします.

7月

GeneticTexture

遺伝的プログラミングをTextureに適用してみた.
報酬とかないのでよくわからないのしか生まれなかった.

gen from yoshi on Vimeo.

github.com

EMアルゴリズム, Variational Bayes, Gibbs samplingの実装

大学の課題でやった.
コード汚いので整理してあげる予定.
HaskellPythonで書いた.

Planet Understanding the amazon from space

Kaggleのコンペ.
時間なくてあんまりいいスコア出せなかった.
GPUマシンが導入されたので今後はもっと色々できる予定.

Planet: Understanding the Amazon from Space | Kaggle

CpawCTF

All Done

CpawCTF - Main page

久しぶりのiOSアプリ

Tomarubaという会社でiOSアプリ開発の仕事を始めた.

Rails tutorial

Rails怖い

Ruby on Rails チュートリアル:実例を使って Rails を学ぼう

読んだ本

Texture and Modeling
英語600ページ辛い.
ノイズの作り方やエイリアシング対策,フラクタル,自然現象のモデリング,ライティングモデル,レイトレーシングなどなど
サンプルコードも多くてProcedural approachの基本をさらうにはとてもいい本ではないかと思う.

https://www.amazon.co.jp/dp/B01E54DILE/ref=dp-kindle-redirect?_encoding=UTF8&btkr=1

Haskell

https://www.amazon.co.jp/gp/product/4798048062/ref=oh_aui_detailpage_o03_s00?ie=UTF8&psc=1

記号創発ロボティクス
立命の谷口先生の本

記号創発ロボティクス 知能のメカニズム入門 (講談社選書メチエ) | 谷口忠大 | 工学 | Kindleストア | Amazon

王様達のヴァイキング

https://www.amazon.co.jp/dp/B00F4TLKDI/ref=docs-os-doi_0

メイドインアビス

https://www.amazon.co.jp/gp/product/B00KYEH7GW/ref=oh_aui_d_detailpage_o00_?ie=UTF8&psc=1

レイン
昔めっちゃ読んでた小説のコミック版

https://www.amazon.co.jp/dp/B00LHBJ5O8/ref=docs-os-doi_0

他色々

研究進めながら他のことをするには1ヶ月の限界を感じてきた

6月

ComputeShader Study

これを作った.

demo4 from yoshi on Vimeo.

yoshi-sa.hatenablog.com

Sansanデータ解析チャレンジ

yoshi-sa.hatenablog.com

Unityからのパラレルポート制御

実験環境をUnityで作っていて脳波計測器との同期信号を出すためにDOIボードを取り付けてUnityから制御できるようにした.

youtu.be

読んだ本

森博嗣のXシリーズ6冊.シリーズ作前作コンプリートまでもうすぐ.

www.amazon.co.jp

石黒先生の「ロボットとは何か」

www.amazon.co.jp

ホリエモンの「多動力」あんまり面白くなかった.

どっかでオススメされたGoogleのやつ

www.amazon.co.jp

BLAME! 内容理解できなかった

https://www.amazon.co.jp/gp/product/B00ZHJU1N0/ref=oh_aui_d_detailpage_o01_

あげくの果てのカノン. 絵が好きで買ったけど微妙

www.amazon.co.jp

BIOMEGA BLAME!よりよかった

https://www.amazon.co.jp/gp/product/B00KNNORXU/ref=oh_aui_d_detailpage_o03_

青の祓魔師

https://www.amazon.co.jp/gp/product/B06X3ZQR64/ref=oh_aui_d_detailpage_o02_

リーンカーネーションの花弁 圧倒的にオススメ.過去の偉人や罪人の能力使って戦うやつ.

www.amazon.co.jp

素足のメテオライト リーンカーネーションの花弁の作者が書いてるから買ったけど微妙

https://www.amazon.co.jp/gp/product/B009WNSA3K/ref=oh_aui_d_detailpage_o01_

僕らはみんな河合荘 良さしかない.ひたすらりっちゃんが可愛すぎる話

https://www.amazon.co.jp/gp/product/B00FZFBN38/ref=oh_aui_d_detailpage_o01_

みたもの

アニメ色々

インターンシップ Googleインターンにおっさん二人が参加する話.話はユーモアあって面白いけどこんな人たちはGoolgeのインターンに受からないのでただのコメディ

https://www.amazon.co.jp/gp/product/B00UMBCAKK/ref=oh_aui_d_detailpage_o00_?ie=UTF8&psc=1

Cross-Entropy MethodによるSVMのパラメータ最適化

パターン認識の講義の課題で以下のコンペに参加したのでその記録に. http://universityofbigdata.net/competition/5723788444434432

Hog特徴量

今回のテーマが与えられた名刺画像の記載項目の予測だったので,画像を横長に調整してからhog特徴量に変換をするようにした.

www26.atwiki.jp

DIR_IMAGES = 'sansan-001/images'
IMG_SIZE = 100

def load_hog(df):
    X = []
    for i, row in df.iterrows():
        img = Image.open(os.path.join(DIR_IMAGES, row.filename))
        img = img.crop((row.left, row.top, row.right, row.bottom))
        if img.size[0] < img.size[1]:
            img = img.transpose(Image.ROTATE_90)
        img = img.convert('L')
        img = np.array(img.resize((216, 72)))/255.0
        img = hog(img, orientations=6, pixels_per_cell=(12,12), cells_per_block=(1,1))
        x = np.asarray(img, dtype=np.float32)
        x = x.flatten()
        X.append(x)

    X = np.array(X)
    return X

Pipelineの生成とGridSearch

GridSearchにかけるパイプラインを生成.
今回は色々試したのとトレーニングデータの数が少なかったのでPCAかけてSVMを用いて精度を上げていくことにした.
9項目あったのでそれぞれに対してOne vs Restの判別器を作ってGridSearchで良さそうなパラメータの学習.
マシンスペックがよくないのでパターン数があまり増えないようにした.

    columns = ['company_name', 'full_name', 'position_name',
               'address', 'phone_number', 'fax',
               'mobile', 'email', 'url']

    '''''''''''''''
    Data Loading
    '''''''''''''''
    df_train = pd.read_csv('sansan-001/train.csv')
    X_train = load_hog(df_train)
    Y_train = df_train[columns].values

    '''''''''''''''
    GridSearch
    '''''''''''''''
    pipe_lr = Pipeline([('sc', StandardScaler()), ('pca', PCA()), ('svc', SVC(probability=True))])
    parameter =  {
        'pca__n_components': [100, 150, 200],
        'svc__C': [1, 10, 100, 1000],
        'svc__kernel': ['rbf'],
        'svc__gamma':[0.0001, 0.001, 0.01],
        'svc__probability':[True]
        }

    clf = GridSearchCV(pipe_lr, parameter, cv=5, n_jobs = -1)
    score_average = 0
    for n in range(Y_dev.shape[1]):
        clf.fit(X_train, Y_train[:, n])
        save_model(clf.best_estimator_, filename='best_svm-' + str(n) + '.pkl')
   

Cross-Entropy Methodによるパラメータ最適化

上記のGridSearchで得られたモデルを使ってテストデータに対して適応して提出してみたらこの時点では過学習もみられずAUCスコアが0.98ちょっとだったので, 更にパラメータを最適化して精度を上げるためにCross-Entropy Methodを実装してみた.

Cross-Entropy Method

Cross-Entropy Methodは分布推定アルゴリズムに類似したパラメータ学習法.
まずパラメータを持ったサンプルを p*N 個作る.
そのサンプルの分布を  U とし,Cross-Entropy  D_{CE}が最小となる分布gを選択する.
 D_{CE}(g||U) = \int g(x)\log \frac{U(x)}{g(x)}dx
この分布gを元に,サンプルをN個生成する.このサンプルを持ったSVMでテストデータに対して交差検証を行い,精度の良いサンプル p*Nをエリートサンプルとして取り出し,それを元に新たな分布Uとする.
これを繰り返してパラメータの学習を行う.

実装

パラメータの分布に関しては正規分布を前提として上位20%のパラメータの平均値を新しい分布の平均値とした.
分散に関しては初期値が最適解に近いと想定して分散一定にし,平均値だけを変動させてスコアの変動を確認した.

結果

色々あって実装するのが遅くなってコンペの開催期間中に結果出なかった辛い.
具体的な数字は出してないけれど,どの項目においてもテストデータの分類精度の向上はみられた.

終わりに

今回はじめて機械学習のコンペに参加して楽しかったので次はKaggleのコンペに参加してみている. https://www.kaggle.com/c/planet-understanding-the-amazon-from-space

機械学習詳しい人からのフィードバック欲しい…
まさかりでもいいので投げて欲しい…
機械学習できるマンになりたい…

ComputeShaderによるアニメーション中の位置計算とDrawMeshInstancedIndirectによるメッシュの描画

リポジトリ https://github.com/SakuragiYoshimasa/ShaderLab 以前ブログに書いたDrawMeshInstancedIndirectのサンプルの拡張.
ComputeShaderを書く練習をしたかったので,
任意のアニメーション中のSkinnedMeshに対してComputeShaderを使って
毎フレームメッシュの描画位置を再計算するのを実装してみた.

前バージョン.固定位置でしかない.

https://yoshidev.tumblr.com/post/161242251494/drawmeshinstancedindirect-study-2566000-tris

準備

最初に描画のためのメッシュの構築と,
アニメーション中のMeshの頂点位置のベイクを行う機構を作ります.

描画するメッシュの構築

今回用いたのはHairと名付けたメッシュ.
描画するときに円錐形を作りたいのでここではメッシュの各頂点の座標の代わりに
(x,y,z) = (位相, 0, 高さ) とした.

public class Hair : ScriptableObject {
    [SerializeField] public int divisions = 4;
    [SerializeField] public int segments = 32;
    [SerializeField] public Mesh mesh;

    public void RebuildMesh(){
        List<Vector3> vertices = new List<Vector3>();

        for(int i = 0; i < segments + 1; i++){
            for(int j = 0; j < divisions + 1; j++){
                float phai = Mathf.PI * 2.0f * (float)j / (float)divisions;
                vertices.Add(new Vector3(phai, 0, i));
            }
        }

        List<int> indices = new List<int>();
        int refi = 0;
        //インデックスの設定
        for (var i = 0; i < segments; i++){
           for (var j = 0; j < divisions; j++){
               indices.Add(refi);
               indices.Add(refi + 1);
               indices.Add(refi + 1 + divisions);

               indices.Add(refi + 1);
               indices.Add(refi + 2 + divisions);
               indices.Add(refi + 1 + divisions);

               refi++;
            }
            refi++;
       }

       mesh.SetVertices(vertices);
       mesh.SetIndices(indices.ToArray(), MeshTopology.Triangles, 0);
       mesh.UploadMeshData(true);
    }

    void OnEnable(){
         if (mesh == null){
           mesh = new Mesh();
           mesh.name = "Hair";
        }
    }
}

これをアセットとして保存

public class HairEditor : Editor {

    [MenuItem("Assets/Create/Hair")]
    static void CreteHair(){
            
        var path = AssetDatabase.GetAssetPath(Selection.activeObject);
        if (string.IsNullOrEmpty(path))
            path = "Assets";
        else if (Path.GetExtension(path) != "")
            path = path.Replace(Path.GetFileName(path), "");
        var assetPathName = AssetDatabase.GenerateUniqueAssetPath(path + "/Hair.asset");

            
        var asset = ScriptableObject.CreateInstance<Hair>();
        AssetDatabase.CreateAsset(asset, assetPathName);
        AssetDatabase.AddObjectToAsset(asset.mesh, asset);

        asset.RebuildMesh();

        AssetDatabase.SaveAssets();

        EditorUtility.FocusProjectWindow();
        Selection.activeObject = asset;
    }
}

アニメーション中のSkinnedMeshのTextureBaking

アニメーション中のSkinnedMeshRendererから直接頂点位置を引っ張る方法がわからなかったので
keijiroさんのSkinner(https://github.com/keijiro/Skinner )を参考に

RenderTexture currPositionBuffer;
RenderTexture currNormalBuffer;

この2つのバッファにベイクします.
実装としてはカメラにFragmentShaderの出力として頂点の位置と法線を出力するShaderを
ReplacementShaderとして設定し,対象のMeshのみから情報を抽出するようにします.

その前準備として対象のSkinnedMeshのuv座標をそのindexを元に書き換えます.

[SerializeField] public SkinnedMeshRenderer targetSMR;

void RecreateMesh(){
    Mesh mesh = new Mesh();
    Mesh origMesh = targetSMR.sharedMesh;

    var vertices = new List<Vector3>(origMesh.vertices);
    var normals = new List<Vector3>(origMesh.normals);
    var tangents = new List<Vector4>(origMesh.tangents);
    var boneWeights = new List<BoneWeight>(origMesh.boneWeights);
    int[] indices = new int[origMesh.vertexCount];
            
    List<Vector2> uv = new List<Vector2>();
    for(int i = 0; i < origMesh.vertexCount; i++){
        uv.Add(new Vector2(((float)i+0.5f) / (float) origMesh.vertexCount, 0));
        indices[i] = i;
    }
            
    mesh.subMeshCount = 1;
    mesh.SetVertices(vertices);
    mesh.SetNormals(normals);
    mesh.SetTangents(tangents);
    mesh.SetIndices(indices, MeshTopology.Points, 0);
    mesh.SetUVs(0, uv);
    mesh.bindposes = origMesh.bindposes;
    mesh.boneWeights = boneWeights.ToArray();

    mesh.UploadMeshData(true);

    targetSMR.sharedMesh = mesh;
}

これをReplacementShaderでBufferの位置として利用して位置と法線をベイクします.

struct v2f {
    float4 position : SV_POSITION;
    float3 texcoord: TEXCOORD0;
    float3 normal: NORMAL;
    float psize : PSIZE;
};
            
struct fragout {
    float4 position: SV_TARGET0;
    float4 normal: SV_TARGET1;
};

v2f vert (appdata v) {
    v2f o;
    o.position = float4(v.texcoord.x * 2 - 1, 0, 0, 1);
    o.normal = UnityObjectToWorldNormal(v.normal);
    o.texcoord = mul(unity_ObjectToWorld, v.vertex).xyz;
    o.psize = 1;
    return o;
}
            
fragout frag (v2f i) : SV_TARGET0 {
    fragout o;
    o.position = float4(i.texcoord, 1.0);
    o.normal = float4(i.normal, 1.0);
    return o;
}

詳しくはSkinnerもしくは上記のリポジトリを参考にしてください.
Skinnerから多少改変してあります.

メッシュの描画位置の計算と描画

CSスクリプト

全体のコントローラとしてSkinnedMeshHairRenderer.csを記述する.
まずはShaderに流すパラメータやバッファ等の定義.

//描画用のマテリアルとメッシュ
[SerializeField] Hair _hair;
[SerializeField] Material _material; 
        
//Shaderに流すパラメータ
[SerializeField] int _instanceCount;
[SerializeField] float _scale = 1.0f;
[SerializeField] float _zScale = 0.07f;
[SerializeField] float _noisePower = 0.1f;
[SerializeField] float _frequency = 1.0f;
[SerializeField] float _radiusAmp = 1.0f;
[SerializeField] Vector3 _gravity = new Vector3(0f, -8.0f, 4.0f);
[SerializeField] Color _color = new Vector4(0,0,0,0);
[SerializeField] Color _gradcolor = new Vector4(0,0,0,0);
        
//バッファ
uint[] _drawArgs = new uint[5]{0, 0, 0, 0, 0};
ComputeBuffer _drawArgsBuffer;
ComputeBuffer positionBuffer;
ComputeBuffer rotationBuffer;
ComputeBuffer indicesBuffer;
ComputeBuffer weightsBuffer;
ComputeBuffer prevPositionBuffer;
        
//ベイク用のコンポーネントとその他参照
AnimatedSkinnedMesh asm;
SkinnedMeshRenderer _skinnedMeshR;
Mesh _targetmesh;
Bounds _bounds = new Bounds(Vector3.zero, Vector3.one * 4 * 32);
MaterialPropertyBlock _props;

//計算用のComputeShader
[SerializeField] ComputeShader _recalcGlownPositionsShader;

初期化.

InitBuffers()でバッファの初期化を行います.
SkinnedMeshの各トライアングルに対して面積が小さいものは除外して
描画するメッシュがなるべく重複しないように三角形を選んでバッファに格納していきます.
ついでにメッシュの半径と長さとランダムに選択して位置と回転のwに格納しています.

void Start(){
    //まずComputeShaderのDispatchに使う引数をセット
    _drawArgsBuffer = new ComputeBuffer(1, 5 * sizeof(uint), ComputeBufferType.IndirectArguments);
    _drawArgs[0] = (uint)_hair.mesh.GetIndexCount(0);
    _drawArgs[1] = (uint)_instanceCount;
    _drawArgsBuffer.SetData(_drawArgs);
            
    //マテリアルに基本パラメータを流す
    _props = new MaterialPropertyBlock();
        _props.SetFloat("_UniqueID", Random.value);
        _material.SetInt("_ArraySize", _hair.segments + 1);
        _material.SetInt("_InstanceCount", _instanceCount);
        _material.SetInt("_SegmentCount", _hair.segments);

    //ベイク用のコンポーネントから参照を引っ張る
    asm = GetComponent<AnimatedSkinnedMesh>();
    _skinnedMeshR = asm.targetSMR;
    _targetmesh = new Mesh();
            
    InitBuffers();
}
        
void InitBuffers(){
    if(positionBuffer != null) positionBuffer.Release();
    if(rotationBuffer != null) rotationBuffer.Release();
    if(indicesBuffer != null) indicesBuffer.Release();
    if(weightsBuffer != null) weightsBuffer.Release();
    if(prevPositionBuffer != null) prevPositionBuffer.Release();

    positionBuffer = new ComputeBuffer(_instanceCount, sizeof(float) * 4);
    prevPositionBuffer = new ComputeBuffer(_instanceCount, sizeof(float) * 4);
    rotationBuffer = new ComputeBuffer(_instanceCount, sizeof(float) * 4);
    indicesBuffer = new ComputeBuffer(_instanceCount * 3, sizeof(int));
    weightsBuffer = new ComputeBuffer(_instanceCount, sizeof(float) * 3);

    Vector4[] positions = new Vector4[_instanceCount];
    Vector4[] rotations = new Vector4[_instanceCount];
    Vector3[] weights = new Vector3[_instanceCount];
    int[] indices = new int[_instanceCount * 3];

    _skinnedMeshR.BakeMesh(_targetmesh);
    int targetTriangleCount = _targetmesh.triangles.Length / 3;
            
    HashSet<int> nouse = new HashSet<int>();
    int index = (int)Random.Range(0f, targetTriangleCount + 0.5f);
    float threshold = 0.04f;

    for (int i=0; i < _instanceCount; i++) {
        index = (index + (int)Random.Range(0f, targetTriangleCount + 0.5f)) % targetTriangleCount;

        while(nouse.Contains(index)){
            index = (index + (int)Random.Range(0f, targetTriangleCount + 0.5f)) % targetTriangleCount;
        }

        Vector3 v1 = _targetmesh.vertices[_targetmesh.triangles[index * 3]];
        Vector3 v2 = _targetmesh.vertices[_targetmesh.triangles[index * 3 + 1]];
        Vector3 v3 = _targetmesh.vertices[_targetmesh.triangles[index * 3 + 2]];
        
        Vector3 n1 = _targetmesh.normals[_targetmesh.triangles[index * 3]];
        Vector3 n2 = _targetmesh.normals[_targetmesh.triangles[index * 3 + 1]];
        Vector3 n3 = _targetmesh.normals[_targetmesh.triangles[index * 3 + 2]];
                
        float mag = ((v1 - v2).magnitude + (v2 - v3).magnitude + (v1 - v3).magnitude)/3.0f; 
        if(mag < threshold) {
            nouse.Add(index);
            i--;
            continue;
        }

        float p1 = Random.Range(0, 1.0f);
        float p2 = Random.Range(0, 1.0f - p1);
        float p3 = 1.0f - p1 - p2;
        Vector3 p = v1 * p1 + v2 * p2 + v3 * p3;
        Vector3 n = n1 * p1 + n2 * p2 + n3 * p3 + new Vector3(0, 1.0f,0);

        float radius = Random.Range(0.015f, 0.05f);
                
        positions[i] = new Vector4(p.x, p.y, p.z, radius);

        Vector3 rotation = Quaternion.LookRotation(n, Vector3.up).eulerAngles;
        rotations[i] = new Vector4(rotation.x / 180.0f * Mathf.PI, rotation.y / 180.0f * Mathf.PI, rotation.z / 180.0f * Mathf.PI, mag);

        indices[i * 3] = _targetmesh.triangles[index * 3];
        indices[i * 3 + 1] = _targetmesh.triangles[index * 3 + 1];
        indices[i * 3 + 2] = _targetmesh.triangles[index * 3 + 2];
        weights[i] = new Vector3(p1, p2, p3);
    }

    positionBuffer.SetData(positions);
    rotationBuffer.SetData(rotations);
    indicesBuffer.SetData(indices);
    weightsBuffer.SetData(weights);

    _material.SetBuffer("_PositionBuffer", positionBuffer);
    _material.SetBuffer("_RotationBuffer", rotationBuffer);
}

ComputeShaderのDispatchとDrawMeshInstancedIndirectによる描画

ComputeShaderと描画用のマテリアルにバッファを渡してDispatch, DrawMeshInstancedIndirectを呼ぶだけ.
多分こんなに毎回バッファ送らなくていい.

 
void Update(){
    int kernel =_recalcGlownPositionsShader.FindKernel("RecalcPositions");
    prevPositionBuffer = positionBuffer;
    var data = new Vector4[_instanceCount];
    positionBuffer.GetData(data);
    prevPositionBuffer.SetData(data);
            
    _recalcGlownPositionsShader.SetBuffer(kernel, "PositionBuffer", positionBuffer); 
    _recalcGlownPositionsShader.SetBuffer(kernel, "RotationBuffer", rotationBuffer); 
    _recalcGlownPositionsShader.SetBuffer(kernel, "IndicesBuffer", indicesBuffer);  
    _recalcGlownPositionsShader.SetBuffer(kernel, "WeightsBuffer", weightsBuffer);
    _recalcGlownPositionsShader.SetTexture(kernel, "CurrPositionBuffer", asm.CurrPositionBuffer);
    _recalcGlownPositionsShader.SetTexture(kernel, "CurrNormalBuffer", asm.CurrNormalBuffer);
    _recalcGlownPositionsShader.Dispatch(kernel, 64, 1, 1);
}

void LateUpdate(){
    _material.SetFloat("_Scale",_scale);
    _material.SetFloat("_ZScale",_zScale);
    _material.SetFloat("_NoisePower",_noisePower);
    _material.SetVector("_Gravity", _gravity);
    _material.SetFloat("_Scale",_scale);
    _material.SetFloat("_Frequency", _frequency);
    _material.SetFloat("_RadiusAmp", _radiusAmp);
    _material.SetColor("_MainColor", _color);
    _material.SetColor("_GradColor", _gradcolor);
    _material.SetBuffer("_PositionBuffer", positionBuffer);
    _material.SetBuffer("_RotationBuffer", rotationBuffer);
    _material.SetBuffer("_PrevPositionBuffer", prevPositionBuffer);
    Graphics.DrawMeshInstancedIndirect(_hair.mesh, 0, _material, _bounds, _drawArgsBuffer, 0, _props);
    }
}

メッシュの描画位置の計算

RecalcGlownPositionsShader.computeに実装を書きます.
結構単純.

kernelとバッファの宣言.
ローポリでも位置が被らないように
メッシュを描画する位置としてベイクされた頂点を直接には使わない.

#pragma kernel RecalcPositions

//ベイクされた座標と法線
RWTexture2D<float4> CurrPositionBuffer; 
RWTexture2D<float4> CurrNormalBuffer; 

//メッシュを描画する位置の計算のためのインデックスと重み係数,初期化以降変更なし
StructuredBuffer<int> IndicesBuffer;
StructuredBuffer<float3> WeightsBuffer;

//メッシュを実際に描画する位置と回転情報
RWStructuredBuffer<float4> PositionBuffer;
RWStructuredBuffer<float4> RotationBuffer;

実際の関数. 位置と法線は三点から重みを元に3頂点の位置の線形和をとる.
さらにこの辺を参考に,法線からオイラー回転量を計算してる.
けどそのままだと微妙だったので2倍して散らしてる(かなり雑).

http://answers.unity3d.com/questions/467614/what-is-the-source-code-of-quaternionlookrotation.html
https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles

[numthreads(64,1,1)]
void RecalcPositions (uint id : SV_DispatchThreadID)
{
    int id1 = IndicesBuffer[id * 3];
    int id2 = IndicesBuffer[id * 3 + 1];
    int id3 = IndicesBuffer[id * 3 + 2];

    float3 weights = WeightsBuffer[id].xyz;

    float3 p1 = CurrPositionBuffer[int2(id1, 0)].xyz;
    float3 p2 = CurrPositionBuffer[int2(id2, 0)].xyz;
    float3 p3 = CurrPositionBuffer[int2(id3, 0)].xyz;

    PositionBuffer[id].xyz = weights.x * p1 + weights.y * p2 + weights.z * p3;
    
    float3 n1 = CurrNormalBuffer[int2(id1, 0)].xyz;
    float3 n2 = CurrNormalBuffer[int2(id2, 0)].xyz;
    float3 n3 = CurrNormalBuffer[int2(id3, 0)].xyz;
    
    float3 forward = (weights.x * n1.xyz + weights.y * n2.xyz + weights.z * n3.xyz).xyz;

    float3 up = float3(0, 1.0, 0);
    float3 v1 = normalize(forward);
    float3 v2 = normalize(cross(up, v1));
    float3 v3 = cross(v1, v2);

    float m00 = v2.x;
    float m01 = v2.y;
    float m02 = v2.z;
    float m10 = v3.x;
    float m11 = v3.y;
    float m12 = v3.z;
    float m20 = v1.x;
    float m21 = v1.y;
    float m22 = v1.z;

    float num8 = (m00 + m11) + m22;
    float4 quaternion = float4(0,0,0,0);
    if (num8 > 0){
         float num = sqrt(num8 + 1.0);
         quaternion.w = num * 0.5;
         num = 0.5 / num;
         quaternion.x = (m12 - m21) * num;
         quaternion.y = (m20 - m02) * num;
         quaternion.z = (m01 - m10) * num;
         
    }else if ((m00 >= m11) && (m00 >= m22)){
         float num7 = sqrt(1.0 + m00 - m11 - m22);
         float num4 = 0.5 / num7;
         quaternion.x = 0.5 * num7;
         quaternion.y = (m01 + m10) * num4;
         quaternion.z = (m02 + m20) * num4;
         quaternion.w = (m12 - m21) * num4;
    }else if (m11 > m22){
         float num6 = sqrt(1.0 + m11 - m00 - m22);
         float num3 = 0.5 / num6;
         quaternion.x = (m10+ m01) * num3;
         quaternion.y = 0.5 * num6;
         quaternion.z = (m21 + m12) * num3;
         quaternion.w = (m20 - m02) * num3;
    }else{
        float num5 = sqrt(1.0 + m22 - m00 - m11);
        float num2 = 0.5 / num5;
        quaternion.x = (m20 + m02) * num2;
        quaternion.y = (m21 + m12) * num2;
        quaternion.z = 0.5 * num5;
        quaternion.w = (m01 - m10) * num2;
    }
    float PI = 3.14159265;
    float4 q = quaternion;
    q = q / length(q);
    float rotz = atan2(2.0 * (q.x * q.y + q.z * q.w), 1.0 - 2.0 * (q.y * q.y + q.z * q.z));
    float roty = asin(2 * (q.x * q.z - q.w * q.y));
    float rotx = atan2(2.0 * (q.x * q.w + q.y * q.z), 1.0 - 2.0 * (q.z * q.z + q.w * q.w));
    
    float3 rot = float3(rotx * 2.0, roty * 2.0, rotz * 2.0);
    RotationBuffer[id].xyz = rot;
}

これを毎回dispatchして位置と回転の更新を行う.

メッシュの描画用のShader

最後に大事な描画.
基本構成はこれ.
vertの中身は後述.

Shader "Instanced/SkinnedInstancedHairShader" {
    Properties {
        _MainTex ("Albedo (RGB)", 2D) = "white" {}
        _MainColor ("Color", Color) = (1.0, 1.0, 1.0, 1.0)
        _GradColor ("Color", Color) = (1.0, 1.0, 1.0, 1.0)
        _Glossiness ("Smoothness", Range(0,1)) = 0.5
        _Metallic ("Metallic", Range(0,1)) = 0.0
        _Scale("Scale", Range(1.0, 10.0)) = 1.0
        
    }
    SubShader {
        Tags { "RenderType"="Opaque" }
    
        CGPROGRAM
        #pragma surface surf Standard vertex:vert
        #pragma instancing_options procedural:setup
        #include "SimplexNoise3D.cginc"
        
        sampler2D _MainTex;

        struct Input {
            float2 uv_MainTex;
            int seg;
        };

        #ifdef UNITY_PROCEDURAL_INSTANCING_ENABLED
        StructuredBuffer<float4> _PositionBuffer;
        StructuredBuffer<float4> _RotationBuffer;
        StructuredBuffer<float4> _PrevPositionBuffer;
        uint _ArraySize;
        uint _InstanceCount;
        uint _SegmentCount;
        float3 _Gravity;
        float _Scale;
        float _ZScale;
        float _NoisePower;
        float _Frequency;
        float _RadiusAmp;
        #endif

        void vert(inout appdata_full v, out Input data){
            UNITY_INITIALIZE_OUTPUT(Input, data);
        }

        void setup(){
            #ifdef UNITY_PROCEDURAL_INSTANCING_ENABLED
            #endif
        }

        half _Glossiness;
        half _Metallic;
        fixed4 _MainColor;
        fixed4 _GradColor;

        void surf (Input IN, inout SurfaceOutputStandard o) {
            #ifdef UNITY_PROCEDURAL_INSTANCING_ENABLED
            fixed4 c = _MainColor * (float)IN.seg / _SegmentCount + (1.0 -  (float)IN.seg / _SegmentCount) * _GradColor;
            o.Albedo = c.rgb;
            o.Metallic = _Metallic;
            o.Smoothness = _Glossiness;
            o.Alpha = c.a;
            #endif
        }
        ENDCG
    }
    FallBack "Diffuse"
}

vertの記述

ComputeShaderで計算したデータを元にメッシュの描画を行います.
まずは最初に定義したように各頂点の高さと位相を抽出し,
ComputeShaderの計算結果をInstanceIDを用いて取得します.

float phi = v.vertex.x;
int seg = (int)v.vertex.z;
data.seg = seg;

float id = (float)unity_InstanceID;
float4 base =  _PositionBuffer[unity_InstanceID];
float4 rotation = _RotationBuffer[unity_InstanceID];

//ついでに格納してあった半径と長さの取得
float radius = base.w;
float length = rotation.w;
            

次に抽出した位相等を元にメッシュの形状を作ります.

float2 xy = float2(cos(phi), sin(phi)) * radius * _RadiusAmp * ((float)(_SegmentCount - seg) / (float)_SegmentCount);
float4 offset = float4(xy, v.vertex.z * length * _ZScale, 1.0)
float3 n_normal = float3(cos(phi), sin(phi), 0);
            

ここまでの結果を見るために

v.vertex.xyz = offset.xyz;
v.normal.xyz = n_normal;

と一時的にかいて見ると

これが設定したインスタンスの数だけ一箇所に描画されます.
次はComputeShaderの解散結果を用いてアニメーション中のメッシュの位置に配置していきます.

v.vertex.xyz = (base.xyz + offset.xyz) * _Scale;
v.normal.xyz = n_normal;

こうなる .
base.xyzがComputeShaderでの位置計算の結果で,
offset.xyzはそれぞれの頂点のメッシュ内での相対位置に値します.

次に計算されたオイラー回転量を元に回転行列を生成しoffset.xyzとの乗算をとりposとします.

float a = rotation.x;
float b = rotation.y;
float c = rotation.z;
            
float4 low1 = float4(cos(a) * cos(b) * cos(c) - sin(a) * sin(c), -cos(a) * cos(b) * sin(c) - sin(a) * cos(c), cos(a) * sin(b), 0);
float4 low2 = float4(sin(a) * cos(b) * cos(c) + cos(a) * sin(c), -sin(a) * cos(b) * sin(c) + cos(a) * cos(c), sin(a) * sin(b), 0);
float4 low3 = float4(-sin(b) * cos(c), sin(b) * sin(c), cos(b), 0); 
float4 low4 = float4(0, 0, 0, 1);
float4x4 rotateMat;
rotateMat._11_12_13_14 = low1;
rotateMat._21_22_23_24 = low2;
rotateMat._31_32_33_34 = low3;
rotateMat._41_42_43_44 = low4;

float3 pos = mul(rotateMat, offset).xyz;

この段階で最終的な出力とすると,
先ほどのメッシュが回転できたことがわかります.

v.vertex.xyz = (base.xyz + pos.xyz) * _Scale;
v.normal.xyz = n_normal;

最後に回転前にノイズと重力による摂動を加えて終わりです.

       
if(seg != 0){
  offset.x += snoise(float3(
                    (float)seg * 0.02 + sin((_Time.x + id * 0.1) * _Frequency),
                    radius + sin(_Frequency *(_Time.x +(float)seg * 0.03)) * cos(_Frequency * (_Time.y -  id * 0.1)),
                    length + cos(_Frequency * (_Time.y)))) * _NoisePower;
    
  offset.y += snoise(float3(
                    (float)seg * 0.02 + cos(_Time.y * _Frequency), 
                    radius - sin(_Time.y *  _Frequency) * cos(_Frequency * (_Time.x +(float)seg * 0.03 + id * 0.1)), 
                    length + cos(_Frequency * (_Time.x +(float)seg * 0.03 + id *0.2)))) * _NoisePower;
}
            
/////////////
//ここで回転
.
.
.
float3 pos = mul(rotateMat, offset).xyz;
/////////////          
            
if(seg!=0) pos += _Gravity * (abs(pos.x) + abs(pos.y - 1.5)) * 0.03;

v.vertex.xyz = (base.xyz + pos.xyz) * _Scale;
v.normal.xyz = mul(rotateMat, n_normal);
      

以上.

結果

結果はこんな感じ.

終わりに

回転の部分をかなり雑にやったので時間があったら直します.
4000個とかメッシュの描画しても耐えるからすごい.

2017.06.25

追記 2017. 09 LT会での資料を追加しました.

www.slideshare.net

5月

4回生はOverwatchやりすぎて成長が皆無だったので自戒のためと何かとすぐ忘れるので毎月やったことを残していきたい.

作ったもの

CEM study

CartPole問題に対して強化学習のタスクにCrossEntropyMethodを適応してエージェントを学習.

github.com

CameraController

ここを参考にUnityで使う簡単なカメラのコントローラーを作った. https://pixta.jp/channel/?p=15795

https://camo.githubusercontent.com/1bab067e3a6ff4532e623925de14cde2b2ba28d1/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f32353638613036342d323062622d306534312d643739662d3766383234633730643330642e676966

github.com

ProceduralStage

背景用の画像生成.

https://yoshidev.tumblr.com/post/161015703784/test-a-procedural-stage-tool-id-like-to
github.com

NoiseViewer

ノイズの実験をするための可視化.

https://camo.githubusercontent.com/caf2855d857336f495f0f147cfcba26e71a7c9c4/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f34353533653862312d653462382d643264352d353838652d6461363561313138386537312e676966

https://camo.githubusercontent.com/52ef8bee53695e6aa8e82b5868278934b472bb20/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f62313762353931302d636633362d353233322d373766632d3139376462353130373138632e676966

https://camo.githubusercontent.com/6338de25c3c21deb4cf729a52ff0fe9568d5cff8/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f37373261626634612d646361342d393039622d633633612d3961336166386131346264372e676966

github.com

DrawMeshInstancedIndirect study

Uniteに行って聴いたDrawMeshInstancedIndirectを試してみた.

https://camo.githubusercontent.com/fd31638471b670cad66f58e1276d2e4f9edf240a/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f30613335613238392d313236642d346635652d633161612d6633623866323936616462302e676966

https://camo.githubusercontent.com/782772a390d950816761b52ccea41eda7859ddd3/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f39393132396632652d393832362d626338322d366233332d6566373735626239363163362e676966

https://camo.githubusercontent.com/03813b5e361e3b19df9913e621ea52956019c69e/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f31626665613438372d646239662d376462652d326336382d6435616232646538313230382e676966

https://camo.githubusercontent.com/636776d9d8b75bd2ac1ca7a310c9bef5420333b5/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f32343266623063302d353136662d346265312d626339322d3136323164303266653866652e676966

https://camo.githubusercontent.com/8ee7ef86566a614df62ac35abed4c77a6080c9a1/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e616d617a6f6e6177732e636f6d2f302f36383537302f35653166653766322d323836362d393634312d373263312d6166373037626465313864642e676966

https://yoshidev.tumblr.com/post/161242251494/drawmeshinstancedindirect-study-2566000-tris

github.com

研究計画

修士での研究の方針をだいたい決めた.

読んだ本

Python 機械学習プログラミング
4月に買って結局読んでなかったやつ.基本理論をわかっている人が実際に機械学習やってみようっていうのにいい本.

眠っていたプログラミングコンテストチャレンジブック

プログラミングコンテストチャレンジブック [第2版] ?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?

プログラミングコンテストチャレンジブック [第2版] ?問題解決のアルゴリズム活用力とコーディングテクニックを鍛える?

Q.E.D 38まで

蜜の島   面白いし4巻だけでオススメできる

ハピネス
押見修造の漫画読み漁ってる

その他色々

その他

Unite2017

前々から気になっていて初めて参加.
見に行ったセッション
・基調講演
・Unity最適化講座 ~スペシャリストが教えるメモリとCPU使用率の負担最小化テクニック~
・パフォーマンス向上のためのスクリプトのベストプラクティス
・シェーダープログラミング入門!カスタムシェーダー、作るで!
・実演!Timeline+Cinemachineによるカットシーン作り
・スクリプタブル・レンダーパイプラインのカスタマイズと拡張 ・Unityライティング最新情報
・インスタンシングを用いた美麗なグラフィックの実現方法
・ゲームの見た目も盛ったら変わる!!!!ヤバい!!ポストプロセス!!入門!!!!!!!!!

基本的にグラフィック周りを中心に見に行った.いつもGithubでお世話になっているkeijiroさんと話せてよかった.

SICF

東京行ったついでに見に行った.

METROCK

久々にフェス.Nulbarichにハマった.

www.youtube.com

Life is Tech! Leaders9期

今回で3回目の講師を終えた.

Overwatch楽しいよ  

ウィドウメーカーやりたすぎて3つめのアカウントを作った. もう発売されて1年経つのに遊び続けてるの珍しい.今無料体験期間だからやってみて欲しい.

www.youtube.com

RenderTextureとShaderを利用した集団同期現象のシミュレーションと視覚化

はじめに

本記事ではUnityを用いて以下のシミュレーション映像を作るためのプロセスをメモついでに書いていきます.

https://yoshidev.tumblr.com/post/158231100497/oscillation-sphere-simulating-kuramoto-model

なお,記述に関して,不正確な部分が生じると思いますので指摘頂ければ幸いです.
また今回の実装にあたって,Keijiro TakahashiさんのOSSを参考している部分があります.

同期現象と蔵本モデルについて

まず,同期現象とは,例えば蛍の集団の点滅周期が森の中の蛍の集団で同期が見られたり,
人間の脳内の領野間でも脳活動の同期によって脳の可塑性を実現している説等,自然界に様々に見られる異なる振動子のリズムが揃っていく現象のことです.

そこで蔵本モデルとは

蔵本由紀によって提案された同期現象を記述する数学モデルである。特に、相互作用のある非線形振動子集団の振る舞いを記述するモデルである。
(https://ja.wikipedia.org/wiki/蔵本モデル)

本質的なことは蔵本先生の講義ノートとか著書を是非読んでみてください.
また同期現象に関してもSYNCという書籍が読み物としても面白いのでオススメです.

今回はこの蔵本モデルをUnityを用いて簡単なシミュレーションを行います.

実装

環境

・Unity 5.6.0b3
MacBook Pro (15-inch, Late 2016)
・CPU:2.9 GHz Intel Core i7/ メモリ:16 GB
Radeon Pro 460 4096 MB
macOS Sierra 10.12.3

使用するスクリプト

・ModelSimulation.cs
・Floor.cs
・CameraController.cs
・RandomBoxMuller.cs (参考)
・Kernel.shader
Surface.shader

モデルの更新処理の実装

まず最初に以下の数式で表されるモデルの更新処理の実装を行います.

文字の意味は
・文字の上のドットは時間微分
・R,Θ:集団の同期度を表す複素パラメタ
・K:それぞれの振動子の相互作用の結合係数
・N:振動子数
・i:振動子インデックス
・ω:振動子の自然振動数
・φ:振動子の位相

この式に従って実装していきます.

ModelSimulation.csとKernel.shaderの役割

今回の実装では多くの振動子集団に対して上記の式に従い更新処理を行うため計算時間の肥大化が考えられます.
そのため,上記の速度の式の変形で振動子間の相互作用は陽には現れない,
を用いることで個々の振動子の位相およびその速度の更新についてはShaderによって並列計算させることを考えます.
同期度を表す複素パラメタの算出については振動子集合の角位相が影響するため ModelSimulation.csにおいて計算を行います.

ModelSimulation.cs

Shaderでの処理のために,RenderTextureを位相,速度のバッファに,Texture2Dを自然振動数のバッファとして利用します.
振動子の数はRenderTextureのサイズ制限により4096までとしています.

フィールドおよびプロパティ

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

public class ModelSimulation : MonoBehaviour {

    [SerializeField, Range(1, 4096)]
    private int _pointNum; //振動子の数
    public int pointNum {
        get { return _pointNum;}
        private set { _pointNum = value;}
    }

    [SerializeField] 
    private float _baseFreq; //振動数分布の中央値
    public float baseFreq {
        get { return _baseFreq;}
        private set { _baseFreq = value;}
    }

    [SerializeField]
    private float _connectionCoefficient; //結合係数
    public float connectionCoefficient {
        get { return _connectionCoefficient;}
        private set { _connectionCoefficient = value;}
    }

    [SerializeField] Shader kernelShader; //位相,速度更新用のShader
    private Material kernelMat; 

    private Texture2D naturalFreqBuffer; //自然振動数のバッファ
    private RenderTexture phaseBuffer; //位相のバッファ
    private RenderTexture velocityBuffer; //位相の速度バッファ

}

マテリアル,バッファの初期化のための関数

ここは素直に.
フォーマットを浮動小数点にしておきます.

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
    Material CreateMaterial(Shader shader)
    {
        var material = new Material(shader);
        material.hideFlags = HideFlags.DontSave;
        return material;
    }

    RenderTexture createRTBuffer(int size){
        var format = RenderTextureFormat.ARGBFloat;
        var buffer = new RenderTexture(size, 1, 0, format);
        buffer.hideFlags = HideFlags.DontSave;
        buffer.filterMode = FilterMode.Point;
        buffer.wrapMode = TextureWrapMode.Clamp;
        return buffer;
    }

    Texture2D createT2Buffer(int size){

        var texture = new Texture2D(size, 1, TextureFormat.RGBAFloat, false);

        for(int i = 0; i < size; i++){
            texture.SetPixel (i, 1, new Color(0, 0, 0, 0));
        }
        texture.Apply ();
        return texture;
    }   
}

初期化

各バッファ,マテリアルの初期化を行います.
マテリアルに流し込むデータに関してはKernel.shaderの部分で記述します.

Graphics.Blit(Texture source, RenderTexture dest, Material mat, int pass = -1);

を用いることで,Shaderのパスを指定して演算結果をdestに格納します.
詳しくはhttps://docs.unity3d.com/jp/540/ScriptReference/Graphics.Blit.htmlを参考にしてください.

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
            
    void Start () {
        Initialize ();
    }   
        
    void Initialize(){

        naturalFreqBuffer = createT2Buffer (pointNum);
        phaseBuffer = createRTBuffer (pointNum);
        velocityBuffer = createRTBuffer (pointNum);

        InitializeMat (); //マテリアル生成とパラメタ初期化
        InitializePosition ();  //Kernel.shaderによる位相の初期化
        InitializeNaturalFreqs (); //μ=0のがウス分布から自然周波数の初期化
    }

    void InitializeMat(){

        kernelMat = CreateMaterial (kernelShader);
        kernelMat.SetInt ("_PointNum", pointNum);
        kernelMat.SetFloat ("_K", connectionCoefficient);
        kernelMat.SetFloat ("_BaseFreq", baseFreq);
        kernelMat.SetFloat ("_DeltaTime", 0);
    }

    void InitializePosition(){
        
        Graphics.Blit (null, phaseBuffer, kernelMat, 0);
        kernelMat.SetTexture ("_PhaseTex", phaseBuffer);
    }

    void InitializeNaturalFreqs(){

        RandomBoxMuller random = new RandomBoxMuller();

        for(int i = 0; i < pointNum; i++){
            naturalFreqBuffer.SetPixel (1, i, new Color(((float)random.next(0, 2.0, true) - 0.5f) * baseFreq , 0, 0, 0)); //ガウス分布に従う自然周波数
        }
        naturalFreqBuffer.Apply ();

        kernelMat.SetTexture ("_NaturalFreqTex", naturalFreqBuffer);
        surfaceMat.SetTexture ("_NaturalFreqTex", naturalFreqBuffer);
    }
}

更新処理

初期化と同じようにGraphics.Blitを用いてKernel.shaderによる更新処理を行います.

public class ModelSimulation : MonoBehaviour {
            ~~~~~略~~~~~
    
    void Update () {
        
        kernelMat.SetFloat ("_DeltaTime", Time.deltaTime);
        
        UpdateVelocity (); //速度バッファの更新
        UpdatePhase (); //位相バッファの更新
    }
                    
    void UpdateVelocity(){
    
        Graphics.Blit (null, velocityBuffer, kernelMat, 1);

        kernelMat.SetTexture ("_VelocityTex", velocityBuffer);
    }

    void UpdatePhase(){
        
        Graphics.Blit (null, phaseBuffer, kernelMat, 2);

        float[] param = calcParams ();

        kernelMat.SetTexture ("_PhaseTex", phaseBuffer);
        kernelMat.SetFloat ("_ParamTheta", param[0]);
        kernelMat.SetFloat ("_ParamR", param[1]);
    }

    //同期度を表す複素パラメタの算出
    float[] calcParams(){

        Texture2D tex = new Texture2D(phaseBuffer.width, phaseBuffer.height, TextureFormat.RGBAFloat, false);

        RenderTexture.active = phaseBuffer;
        tex.ReadPixels(new Rect(0, 0, phaseBuffer.width, phaseBuffer.height), 0, 0);
        tex.Apply();

        float real = 0;
        float imag = 0;

        for(int i = 0; i < pointNum; i++){
            float phai = tex.GetPixel (i, 0).r;
            real += Mathf.Cos (phai);
            imag += Mathf.Sin (phai);
        }

        real /= (float)pointNum;
        imag /= (float)pointNum;

        float paramTheta = Mathf.Atan (imag/real);
        float paramR = Mathf.Sqrt (real * real + imag * imag);

        radius = paramR * 0.7f + 0.5f;


        return new float[2]{paramTheta, paramR};
    }
}

特に難しいことはしていないですがこれでModelSimulation.csnのモデルの更新処理の実装は終わりです.
次は実際に位相の初期化,位相,速度の更新を行うKernel.shaderの実装をします.

Kernel.shader

プロパティ

Shaderに流し込むデータをここで定義します.

Shader "Hidden/Kernel"
{
    Properties
    {
       _PhaseTex ("-", 2D) = ""{} //位相バッファ
       _VelocityTex ("-", 2D) = ""{} //速度バッファ
       _NaturalFreqTex("-", 2D) = ""{} //自然振動数バッファ
       _K           ("-", Float) = 0 //相互作用の結合係数
       _ParamTheta  ("-", Float) = 0 //同期度を表す複素パラメタのΘ
       _ParamR      ("-", Float) = 0 //同期度を表す複素パラメタのR
       _BaseFreq    ("-", Float) = 0 //自然振動数の分布の中央値
       _DeltaTime   ("-", Float) = 0 //dt
       _PointNum    ("-", int) = 0 //振動子の数
    }

    CGINCLUDE

    #include "UnityCG.cginc"
   
    sampler2D _PhaseTex;
    sampler2D _VelocityTex;
    sampler2D _NaturalFreqTex;

    float _BaseFreq;
    float _ParamTheta;
    float _ParamR;
    float _K;
    float _DeltaTime;

    int _PointNum;
    
    ENDCG
}

パスの定義

SubShader{}

の中でパスの定義を行います.
上から順にGraphics.Blit()で指定するpassになります.
VertexShaderはUnityCG.cgingで定義されているvert_imgを用いて,FragmentShaderのみ自前で実装します.

Shader "Hidden/Kernel"
{
            ~~~~~略~~~~~
            
    CGINCLUDE
            ~~~~~略~~~~~
    ENDCG
                
    SubShader
    {
        Pass //位相の初期化 pass=0
        {
            CGPROGRAM
            #pragma target 3.0
            #pragma vertex vert_img
            #pragma fragment frag_init_phase
            ENDCG
        }


        Pass //速度の更新 pass=1
        {
            CGPROGRAM
           #pragma target 3.0
           #pragma vertex vert_img
           #pragma fragment frag_update_velocity
            ENDCG
  
        }

        Pass //位相の更新 pass=2
        {
            CGPROGRAM
            #pragma target 3.0
            #pragma vertex vert_img
            #pragma fragment frag_update_phase
            ENDCG
        }
    }
    
}

FragmentShaderの実装

振動子の状態の初期化,更新を上記の式の通りに実装するだけです.

Shader "Hidden/Kernel"
{
            ~~~~~略~~~~~
    CGINCLUDE
            ~~~~~略~~~~~

    float nrand(float2 uv, float salt)
    {
        uv += float2(salt, 0);
        return frac(sin(dot(uv, float2(12.9898, 78.233))) * 43758.5453);
    }

    float4 frag_init_phase(v2f_img i) : SV_Target {     

        float v = nrand(i.uv, 10) * _BaseFreq * 2.0;

        return float4(v,0,0,0);
    }


    float4 frag_update_velocity(v2f_img i) : SV_Target {

        //バッファのColor.rにが格納されている位相,速度を取り出して計算する.
        float p = tex2D(_PhaseTex, i.uv).x; 
        float nf = tex2D(_NaturalFreqTex, i.uv).x;
        float v = nf + _K * _ParamR * sin(_ParamTheta - p);

        return float4(v, 0, 0, 0);
    }

    float4 frag_update_phase(v2f_img i) : SV_Target {

        float p = tex2D(_PhaseTex, i.uv).x;
        float v = tex2D(_VelocityTex, i.uv).x;

        v = p + v * _DeltaTime;

        return float4(v, 0, 0, 0);
    }       
            
    ENDCG
                        
            ~~~~~略~~~~~           
}
 

モデルの更新処理は以上になります.
この状態でGameObjectに貼り付けて実行してみるとRが大きくなっていくのがわかると思います.
以降では可視化のための実装を行います.

可視化の実装

可視化にはModelSimulation.cs,Surface.shaderに加え,Blenderで簡単に頂点数と三角形をいい感じにした球体を用いました.
この球体は別に他のものでもいいのでどこかにあげることはしない予定です.

ModelSimulation.cs

可視化とシミュレーションを同じスクリプトに書くのはどうかと思ったのですが,
ファイル数を減らしたかったのでそのままModelSimulation.csに可視化の実装も書いちゃいました.

フィールドおよびプロパティ

Shaderに流し込む経過時間,色とかの設定のためのフィールドとプロパティをつらつらと.

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    [SerializeField] Shader surfaceShader;
    [SerializeField] private Mesh mesh; //入力メッシュ
    private Material surfaceMat;
            
    [SerializeField] private Color _color1;
    public Color color1 {
        get{ return _color1;}
        set{ _color1 = value;}
    }

    [SerializeField] private Color _color2;
    public Color color2 {
        get{ return _color2;}
        set{ _color2 = value;}
    }

    [SerializeField, Range(0f, 1.0f)] private float _metallic = 0.5f;
    public float metallic {
        get{ return _metallic;}
        set{ _metallic = value;}
    }

    [SerializeField, Range(0f, 1.0f)] private float _smoothness = 0.5f;
    public float smoothness {
        get{ return _smoothness;}
        set{ _smoothness = value;}
    }

    [SerializeField, Range(0.5f, 2.0f)] private float _radius = 1.0f;
    public float radius {
        get{ return _radius;}
        set{ _radius = value;}
    }
        
    private float elapsedTime;
    
            ~~~~~略~~~~~
}

初期化

例によってここでマテリアルの初期化を行います.
また,それぞれの頂点に対して振動子を割り当てるためにUVの初期化を行います.
頂点数と振動子の数を合わせる必要はありませんがなるべく近い方がいいかもしれない.

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    void Initialize(){

            ~~~~~略~~~~~
            
        elapsedTime = 0f;
        int vNum = mesh.vertexCount;

        Vector2[] uv = new Vector2[vNum];

        for(int i = 0; i < vNum; i++){
            uv [i] = new Vector2 (((float)i + 0.5f)/ (float)vNum, 0.5f);
        }

        mesh.uv = uv;
    }
    
    void InitializeMat(){

        surfaceMat = CreateMaterial (surfaceShader);
            ~~~~~略~~~~~   
    }
    
    
            ~~~~~略~~~~~
}

更新処理

そのままですね,毎フレームマテリアルにバッファ等を流し込んで

Graphics.DrawMesh(Mesh mesh, Matrix4x4 matrix, Material material, int layer)

で描画を行うだけです.

public class ModelSimulation : MonoBehaviour {

            ~~~~~略~~~~~
            
    void Update () {
        

            ~~~~~略~~~~~
        elapsedTime += Time.deltaTime;
        DrawMesh ();
    }   
    
            ~~~~~略~~~~~
            
    void DrawMesh(){
        
        surfaceMat.SetColor ("_Color1", color1);
        surfaceMat.SetColor ("_Color2", color2);
        surfaceMat.SetTexture ("_PhaseTex", phaseBuffer);
        surfaceMat.SetFloat ("_Metallic", metallic);
        surfaceMat.SetFloat ("_Smoothness", smoothness);
        surfaceMat.SetFloat ("_Radius", radius);
        surfaceMat.SetFloat ("_ElapsedTime", elapsedTime);
        surfaceMat.SetFloat ("_BaseFreq", baseFreq);

        Graphics.DrawMesh (mesh, transform.localToWorldMatrix, surfaceMat, gameObject.layer);
    }           
}

以上でModelSimulation.csは完成です.
次に描画のためのShaderを書いていきます.

Surface.shader

プロパティ

Shader "Custom/Surface" {
    Properties
    {
        _PhaseTex ("-", 2D) = ""{}
        _NaturalFreqTex("-", 2D) = ""{}

        _Color1 ("-", Color) = (1, 1, 1, 1)
        _Color2 ("-", Color) = (1, 1, 1, 1)

        _Metallic ("-", Float) = 0
        _Smoothness ("-", Float) = 0

        _Radius ("-", Float) = 0
        _PointNum ("-", Float) = 0
    }

    CGINCLUDE

    sampler2D _PhaseTex;
    sampler2D _NaturalFreqTex;
    float4 _PositionTex_TexelSize;
    float4 _NaturalFreqTex_TexelSize;

    half4 _Color1;
    half4 _Color2;

    half _Metallic;
    half _Smoothness;

    float _Radius;
    float _PointNum;

    float _ElapsedTime;
    float _BaseFreq;
    
    ENDCG

}

Vertex, Framgment Shaderの記述

振動子の位相によって頂点の座標を変換します.
特に他に書くことはなさそうなのでこれぐらいにしておきます.

Shader "Custom/Surface" {
            ~~~~~略~~~~~
    CGINCLUDE
            ~~~~~略~~~~~
            
    struct Input {
        half color;
    };

    float nrand(float2 uv, float salt)
    {
        uv += float2(salt, 0);
        return frac(sin(dot(uv, float2(12.9898, 78.233))) * 43758.5453);
    }

    void vert(inout appdata_base v, out Input data)
    {
        UNITY_INITIALIZE_OUTPUT(Input, data);

        float2 uv = v.texcoord;
        float p = tex2Dlod(_PhaseTex, float4(uv , 0, 0)).x;
        float theta = _BaseFreq * _ElapsedTime + p;

        v.vertex.xyz = v.vertex.xyz + v.normal * (sin(theta) + 0.3) * (length(v.vertex.xyz)) * _Radius;

        float ln = (sin(theta) + 1.0) / 2.0;
        data.color = ln;
    }
    ENDCG
    
     SubShader
    {
        Tags { "RenderType"="Opaque" }

        CGPROGRAM

        #pragma surface surf Standard vertex:vert nolightmap addshadow
        #pragma target 3.0

        void surf(Input IN, inout SurfaceOutputStandard o)
        {
            o.Albedo = lerp(_Color1, _Color2, IN.color);
            o.Metallic = _Metallic;
            o.Smoothness = _Smoothness;
        }

        ENDCG
    }                   
}

Surface.shaderの実装は以上です.
現時点でこんな感じにしてやれば振動子の振る舞いが可視化されると思います.
パラメータを変えて振る舞いを鑑賞します.笑

スクリーンショット 2017-03-11 18.43.53.png

次の項目はおまけです.

見た目を整える

背景

Quadを6枚内側に向けて以下の画像のように配置.(Walls)

スクリーンショット 2017-03-11 18.53.42.png

新しくマテリアルWallMat,Wall.csを作ってそれぞれのWallに以下のように貼り付ける.

Wall.csの中身は以下.

public class Wall : MonoBehaviour {

    Texture2D emmitionTex;
    int texSize = 1024;

    void Start () {
        emmitionTex = new Texture2D (texSize, texSize);

        for(int i = 0; i < texSize; i ++){
            for(int j = 0; j < texSize; j++){
                if (i % 64 == 0 || j % 64 == 0) {
                    emmitionTex.SetPixel (i, j, new Color (1.0f, 1.0f, 1.0f)); 
                } else {
                    emmitionTex.SetPixel (i, j, new Color (32f / 255f, 32f / 255f, 32f / 255f));    
                }
            }
        }

        emmitionTex.Apply ();

        Material fMat = gameObject.GetComponent<Renderer> ().material;
        fMat.SetTexture ("_EmissionMap", emmitionTex);
        fMat.SetColor ("_EmissionColor", new Color(1.0f, 1.0f, 1.0f, 1.0f));
    }
}

Cameraとライト

Cameraの動き

テキトーにぐるぐるさせます.

public class CameraController : MonoBehaviour {

    float radius = 3.5f;
    float time = 0;

    void Update () {
        time += Time.deltaTime;
        this.transform.position = new Vector3 (radius * Mathf.Cos (time / 2.0f), radius/ 1.5f * Mathf.Sin (time / 3.0f), (radius + 1.0f) * Mathf.Sin (time / 2.0f));
        this.transform.LookAt (Vector3.zero, Vector3.up);
    }
}

ImageEffectとライト

いい感じになるようにカメラにImageEffectをつける.
いい感じになるようにPointLightを配置する.

終わりに

長くなりましたが以上になります.
この手のモデルの可視化が好きなので試してほしいモデルなどありましたら紹介ください.

また,意見や質問,歓迎しますのでよろしくお願いします.
ありがとうございました.

2017.03.11