なおしのこれまで、これから

学んだこと・感じたこと・やりたいこと

近傍探索の実装とまとめ

f:id:vxd-naoshi-19961205-maro:20210116180008g:plain

始めに

今回の内容は次の記事を参考にして実装した内容になります。解説する内容は大体同じになってしまうかもしれません。

qiita.com


Boidsアルゴリズムを改良したいと考えていた際にこの記事を見つけました。どんなものか調べるために自分で実装してみました。


また今回のプロジェクトをGitHubの方に上げました。 バージョンはUnity2019.4.1fです。

github.com


近傍探索の処理

まず、サンプルとして一定の範囲内を動き続けるパーティクルを作成します。

f:id:vxd-naoshi-19961205-maro:20210116153620g:plain



近傍探索で行う処理は以下の通りです。

  • 各パーティクルが属する格子を求める
  • 格子のインデックスをもとにパーティクルのソート
  • 各格子の始点と終点を求める
  • パーティクルを格子のインデックス順に並び替える
  • 近傍探索

各パーティクルが属する格子を求める

f:id:vxd-naoshi-19961205-maro:20210116161247p:plain

各パーティクルの座標からどの格子に属しているか調べてパーティクルIDと格子Indexをペアで保存します。

ここで各パーティクルがどの格子に属しているかはindex = (position - gridMinBound) / gridSizeで求めることが出来ます。

図は2次元ですが3次元でも同様にして求めます。



格子のインデックスをもとにパーティクルのソート

f:id:vxd-naoshi-19961205-maro:20210116161846p:plain

先ほどのパーティクルIDと格子indexを保存した配列を格子indexでソートします。

ここでのソートアルゴリズムはBitonicSortを使用します。 BitonicSortは過去記事でまとめたので興味があれば参考にしてください。

shitakami.hatenablog.com



各格子の始点と終点を求める

まず、格子の個数分の配列を用意します。初めは始点終点ともに無限で初期化します。

f:id:vxd-naoshi-19961205-maro:20210116170118p:plain
次に先ほどソートしたパーティクルIDと格子indexの配列をもとにして各格子の始点と終点を求めます。 注意として begin <= index < end として保存しています。

f:id:vxd-naoshi-19961205-maro:20210116170140p:plain



パーティクルを格子のインデックス順に並び替える

ComputeShaderではランダムアクセスが遅いという性質があります。(参考:仮眠プログラマーのつぶやき : UnityでGPGPUその2 並列計算と時間測定、コアレスアクセスなど)

なので、一度パーティクルの情報が入った配列もソートされた順番に並び替えます。

この処理がないと近傍探索が遅くなってしまいます。



近傍探索

f:id:vxd-naoshi-19961205-maro:20210116165619p:plain

あるパーティクルの周囲にあるものを調べる際は先ほどの始点と終点を保存した配列をもとに調べます。

例えば図のようにパーティクルIDが2の要素の周囲を調べる際は自身が属する格子とその周りの格子を調べることで周りにあるパーティクルを調べることができます。



プログラム

かなり大きなプログラムになっているので注意してください。

Particle.cs

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System.Runtime.InteropServices;

public class Particle : MonoBehaviour
{

    [SerializeField]
    private ParticleRenderer m_particleRenderer;
    [SerializeField]
    private Material m_instanceMaterial;


    [Space(20)]
    [SerializeField]
    private int m_instanceCount;
    [SerializeField]
    private ComputeShader m_particleCalculator;

    [Header("パーティクル設定")]
    [Space(20)]
    [SerializeField]
    private float m_particleRange;

    [SerializeField]
    private float m_scale;
    [SerializeField]
    private Color m_color;
    [SerializeField]
    private float m_particleVelocity;
    [SerializeField]
    private float m_searchRange;
    [SerializeField]
    private int m_searchIndex;

    private int m_gridRange;

    private int m_calcParticlePositionKernel;
    private int m_initializeGridInfoKernel;
    private int m_calcParticleGridKernel;
    private int m_buildGridIndexRangeKernel;
    private int m_rearrengeBufferKernel;
    private int m_copyBufferKernel;
    private int m_searchNeighborKernel;
    private int m_bitonicSortKernel;
    
    private ComputeBuffer m_particleBuffer;
    private ComputeBuffer m_particleFieldBuffer;
    private ComputeBuffer m_gridInfoBuffer;
    private ComputeBuffer m_sortedParticleBuffer;
    
    private Vector3Int m_gridGroupSize;
    private Vector3Int m_particleGroupSize;

    private readonly int m_NumThreads = 64;

    private readonly int m_DeltaTimePropId = Shader.PropertyToID("_DeltaTime");
    private readonly int m_swapBitPropId = Shader.PropertyToID("_SwapBit");
    private readonly int m_upperBitPropId = Shader.PropertyToID("_UpperBit");

    /// <summary>
    /// ComputeBufferをReleaseする用
    /// </summary>
    private System.Action m_releaseBufferAction;

    struct ParticleInfo {
        public Vector3 position;
        public Vector4 color;
        public float scale;
        public Vector3 velocity;
    }

    // Start is called before the first frame update
    void Start()
    {
        m_instanceCount = Mathf.ClosestPowerOfTwo(m_instanceCount);
        m_particleGroupSize = new Vector3Int(m_instanceCount/m_NumThreads, 1, 1);
        m_particleRenderer.InitializeParticleRenderer(m_instanceCount, m_instanceMaterial);

        // マス目の個数を計算
        // Mathf.CeilToInt:引数以上の整数を返す
        // 参考 : https://docs.unity3d.com/ja/current/ScriptReference/Mathf.CeilToInt.html
        m_gridRange = Mathf.CeilToInt((m_particleRange * 2) / m_searchRange);
        m_particleCalculator.SetInt("_GridTotalRange", m_gridRange*m_gridRange*m_gridRange);

        InitializeParticleBuffer();
        InitializeParticleFieldBuffer();
        InitializeGridInfoBuffer();
        InitializeSortedParticleBuffer();

        SetUpCalcParticlePosition();
        SetUpCalcParticleGrid();
        SetUpBitonicSort();
        SetUpInitializeGridInfo();
        SetUpBuildGridIndexRange();
        SetUpSearchNeighbor();
        SetUpRearrengeBuffer();
        SetUpCopyBuffer();
    }


    void Update() {

        // パーティクルのGridIndexを求める
        m_particleCalculator.Dispatch(m_calcParticleGridKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);

        // GridIndexでソート
        BitonicSort();

        // Grid情報を初期化
        m_particleCalculator.Dispatch(m_initializeGridInfoKernel,
                                     m_gridGroupSize.x,
                                     m_gridGroupSize.y,
                                     m_gridGroupSize.z);


        // Grid内にあるパーティクルのbeginとendを求める
        m_particleCalculator.Dispatch(m_buildGridIndexRangeKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);

        // _SortedParticleBufferにGridIndex順でパーティクル情報を移す
        m_particleCalculator.Dispatch(m_rearrengeBufferKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);

        // m_SortedParticleBufferをm_ParticleBufferに移す
        m_particleCalculator.Dispatch(m_copyBufferKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);


        // 近傍探索
        m_particleCalculator.SetInt("_SearchIndex", m_searchIndex);
        m_particleCalculator.Dispatch(m_searchNeighborKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);

        // パーティクルの移動
        m_particleCalculator.SetFloat(m_DeltaTimePropId, Time.deltaTime);
        m_particleCalculator.Dispatch(m_calcParticlePositionKernel,
                                     m_particleGroupSize.x,
                                     m_particleGroupSize.y,
                                     m_particleGroupSize.z);


    }


    private void InitializeParticleBuffer() {

        ParticleInfo[] particles = new ParticleInfo[m_instanceCount];

        for(int i = 0; i < m_instanceCount; ++i) {
            particles[i].position = RandomVector(-m_particleRange, m_particleRange);
            particles[i].color = m_color;
            particles[i].scale = m_scale;
            // Random.onUnitySphere:半径1の球面上のランダムな点を返す
            // つまり、大きさm_particleVelocityのランダムなベクトルを計算
            particles[i].velocity = Random.onUnitSphere * m_particleVelocity;
        }

        m_particleBuffer = new ComputeBuffer(m_instanceCount, Marshal.SizeOf(typeof(ParticleInfo)));
        m_releaseBufferAction += () => m_particleBuffer?.Release();
        m_particleBuffer.SetData(particles);

        m_instanceMaterial.SetBuffer("_ParticleBuffer", m_particleBuffer);
    }


    private void InitializeParticleFieldBuffer() {

        m_particleFieldBuffer = new ComputeBuffer(m_instanceCount, Marshal.SizeOf(typeof(uint))*2);
        m_releaseBufferAction += () => m_particleFieldBuffer?.Release();

    }


    private void InitializeGridInfoBuffer() {

        // Gridはマス目の個数で初期化する
        m_gridInfoBuffer = new ComputeBuffer(m_gridRange*m_gridRange*m_gridRange, Marshal.SizeOf(typeof(uint))*2);
        m_releaseBufferAction += () => m_gridInfoBuffer?.Release();

    }

    private void InitializeSortedParticleBuffer() {

        m_sortedParticleBuffer = new ComputeBuffer(m_instanceCount, Marshal.SizeOf(typeof(ParticleInfo)));
        m_releaseBufferAction += () => m_sortedParticleBuffer?.Release();

    }

    private void SetUpCalcParticlePosition() {

        m_calcParticlePositionKernel = m_particleCalculator.FindKernel("CalcParticlePosition");

        m_particleCalculator.GetKernelThreadGroupSizes(m_calcParticlePositionKernel,
                                                      out uint x,
                                                      out uint y,
                                                      out uint z);

        m_particleCalculator.SetFloat("_ParticleRange", m_particleRange);
        m_particleCalculator.SetBuffer(m_calcParticlePositionKernel, "_Particle", m_particleBuffer);


    }


    private void SetUpCalcParticleGrid() {


        m_particleCalculator.SetInt("_GridRange", m_gridRange);
        float minGrid = -m_particleRange;
        m_particleCalculator.SetFloat("_MinGrid", minGrid);
        m_particleCalculator.SetFloat("_GridSize", m_searchRange);

        m_calcParticleGridKernel = m_particleCalculator.FindKernel("CalcParticleGrid");

        m_particleCalculator.SetBuffer(m_calcParticleGridKernel, 
                                      "_Particle", 
                                      m_particleBuffer);

        m_particleCalculator.SetBuffer(m_calcParticleGridKernel, 
                                      "_ParticleField", 
                                      m_particleFieldBuffer);

    }


    private void SetUpBitonicSort() {

        m_bitonicSortKernel = m_particleCalculator.FindKernel("ParallelBitonic");
        m_particleCalculator.SetBuffer(m_bitonicSortKernel,
                                      "_ParticleField", 
                                      m_particleFieldBuffer);

    }


    private void SetUpInitializeGridInfo() {

        int totalGridRange = m_gridRange*m_gridRange*m_gridRange;
        m_initializeGridInfoKernel = m_particleCalculator.FindKernel("InitializeGridInfo");
        m_particleCalculator.GetKernelThreadGroupSizes(m_initializeGridInfoKernel, out uint x, out uint y, out uint z);
        m_gridGroupSize = new Vector3Int(Mathf.CeilToInt((float)totalGridRange/x), (int)y, (int)z);

        m_particleCalculator.SetBuffer(m_initializeGridInfoKernel, "_GridInfo", m_gridInfoBuffer);
        
    }


    private void SetUpBuildGridIndexRange() {

        m_buildGridIndexRangeKernel = m_particleCalculator.FindKernel("BuildGridIndexRange");
        m_particleCalculator.SetBuffer(m_buildGridIndexRangeKernel,
                                      "_ParticleField",
                                      m_particleFieldBuffer);
        m_particleCalculator.SetBuffer(m_buildGridIndexRangeKernel,
                                      "_GridInfo",
                                      m_gridInfoBuffer);
        m_particleCalculator.SetInt("_ParticleCount", m_instanceCount);
    }


    private void SetUpSearchNeighbor() {

        m_searchNeighborKernel = m_particleCalculator.FindKernel("SearchNeighbor");
        m_particleCalculator.SetBuffer(m_searchNeighborKernel,
                                      "_Particle",
                                      m_particleBuffer);
        m_particleCalculator.SetBuffer(m_searchNeighborKernel,
                                      "_GridInfo",
                                      m_gridInfoBuffer);
        m_particleCalculator.SetBuffer(m_searchNeighborKernel,
                                      "_ParticleField",
                                      m_particleFieldBuffer);
        m_particleCalculator.SetBuffer(m_searchNeighborKernel,
                                      "_SortedParticle",
                                      m_sortedParticleBuffer);
        m_particleCalculator.SetFloat("_SearchRange", m_searchRange);
        m_particleCalculator.SetInt("_SearchIndex", m_searchIndex);
    }

    private void SetUpRearrengeBuffer() {

        m_rearrengeBufferKernel = m_particleCalculator.FindKernel("RearrengeBuffer");
        m_particleCalculator.SetBuffer(m_rearrengeBufferKernel,
                                      "_ParticleField",
                                      m_particleFieldBuffer);
        m_particleCalculator.SetBuffer(m_rearrengeBufferKernel, 
                                      "_SortedParticle", 
                                      m_sortedParticleBuffer);
        m_particleCalculator.SetBuffer(m_rearrengeBufferKernel, 
                                      "_Particle", 
                                      m_particleBuffer);

    }


    private void SetUpCopyBuffer() {

        m_copyBufferKernel = m_particleCalculator.FindKernel("CopyBuffer");
        m_particleCalculator.SetBuffer(m_copyBufferKernel,
                                      "_Particle",
                                      m_particleBuffer);
        m_particleCalculator.SetBuffer(m_copyBufferKernel,
                                      "_SortedParticle",
                                      m_sortedParticleBuffer);

    }

    private void BitonicSort() {

        int logCount = (int)Mathf.Log(m_instanceCount, 2);

        for(int i = 0; i < logCount; ++i) {

            for(int j = 0; j <= i; j++) {

                int swapBit = 1 << (i - j);
                int upperBit = 2 << i;

                m_particleCalculator.SetInt(m_swapBitPropId, swapBit);
                m_particleCalculator.SetInt(m_upperBitPropId, upperBit);
                m_particleCalculator.Dispatch(m_bitonicSortKernel,
                                             m_particleGroupSize.x,
                                             m_particleGroupSize.y,
                                             m_particleGroupSize.z);
            }
        }
    }


    private Vector3 RandomVector(float min, float max) {

        return new Vector3(
            Random.Range(min, max),
            Random.Range(min, max),
            Random.Range(min, max)
            );

    }

    // 領域の解放
    private void OnDisable() {

        m_releaseBufferAction();

    }


}


ParticleRenderer.cs

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

public class ParticleRenderer : MonoBehaviour
{

    [Header("DrawMeshInstancedIndirectのパラメータ")]
    [SerializeField]
    private Mesh m_mesh;

    [SerializeField]
    private Bounds m_bounds;

    [SerializeField]
    private ShadowCastingMode m_shadowCastingMode;

    [SerializeField]
    private bool m_receiveShadows;

    private ComputeBuffer m_argsBuffer;

    private Material m_instanceMaterial;

    public void InitializeParticleRenderer(int instanceCount, Material instanceMaterial) {

        uint[] args = new uint[5] { 0, 0, 0, 0, 0 };

        uint numIndices = (m_mesh != null) ? (uint) m_mesh.GetIndexCount(0) : 0;

        args[0] = numIndices;
        args[1] = (uint)instanceCount;
        
        m_argsBuffer = new ComputeBuffer(1, args.Length * sizeof(uint), ComputeBufferType.IndirectArguments);

        m_argsBuffer.SetData(args);
        m_instanceMaterial = instanceMaterial;
    }


    // Update is called once per frame
    void LateUpdate()
    {

        Graphics.DrawMeshInstancedIndirect(
            m_mesh,
            0,
            m_instanceMaterial,
            m_bounds,
            m_argsBuffer,
            0,
            null,
            m_shadowCastingMode,
            m_receiveShadows
        );

    }


    private void OnDestroy() {

        m_argsBuffer?.Release();

    }

}


ParticleCalclator.compute

#pragma kernel CalcParticlePosition
#pragma kernel CalcParticleGrid
#pragma kernel InitializeGridInfo
#pragma kernel BuildGridIndexRange
#pragma kernel ParallelBitonic
#pragma kernel SearchNeighbor
#pragma kernel RearrengeBuffer
#pragma kernel CopyBuffer

struct Particle {
    float3 position;
    float4 color;
    float scale;
    float3 velocity;
};


RWStructuredBuffer<Particle> _Particle;
RWStructuredBuffer<int2> _ParticleField;   // x:particle id,  y:grid index
RWStructuredBuffer<int2> _GridInfo;        // x:begin,  y:end
RWStructuredBuffer<Particle> _SortedParticle;

float _ParticleRange;
float _DeltaTime;

uint _ParticleCount;
int _GridRange;
uint _GridTotalRange;
float _MinGrid;
float _GridSize;
float _SearchRange;

uint _SearchIndex;

int CalcGridIndex(float3 particlePosition) {

    int GridX = (particlePosition.x - _MinGrid) / _GridSize;
    int GridY = (particlePosition.y - _MinGrid) / _GridSize;
    int GridZ = (particlePosition.z - _MinGrid) / _GridSize;

    return GridX + GridY*_GridRange + GridZ*_GridRange*_GridRange;

}

int3 CalcGridIndex3(float3 particlePosition) {

    int x = (particlePosition.x - _MinGrid) / _GridSize;
    int y = (particlePosition.y - _MinGrid) / _GridSize;
    int z = (particlePosition.z - _MinGrid) / _GridSize;

    return int3(x, y, z);

}

int CalcIndex3ToIndex(int x, int y, int z) {

    return x + y*_GridRange + z*_GridRange*_GridRange;

}


// パーティクルの座標をフィールドのindexにする計算
[numthreads(64, 1, 1)]
void CalcParticleGrid(uint id : SV_DISPATCHTHREADID) {

    int index = CalcGridIndex(_Particle[id].position);
    _ParticleField[id] = int2(id, index);
    
}


// インデックス情報の初期化
[numthreads(64, 1, 1)]
void InitializeGridInfo(uint id : SV_DISPATCHTHREADID) {

    if(id >= _GridTotalRange)
        return;

    _GridInfo[id] = int2(0xffffff, 0xffffff);

}


[numthreads(64, 1, 1)]
void BuildGridIndexRange(uint id : SV_DISPATCHTHREADID) {

    uint prevId = (id == 0) ? _ParticleCount - 1 : id - 1;
    uint nextId = (id == _ParticleCount - 1) ? 0 : id + 1;

    int index = _ParticleField[id].y;
    int prevIndex = _ParticleField[prevId].y;
    int nextIndex = _ParticleField[nextId].y;

    // 前の格子indexと異なればidがbeginとなる
    if(index != prevIndex) {
        _GridInfo[index].x = id;
    }

    // 後の格子indexと異なればid + 1がendとなる
    if(index != nextIndex) {
        _GridInfo[index].y = id + 1;
    }

}


// パーティクルの移動
[numthreads(64, 1, 1)]
void CalcParticlePosition(uint id : SV_DISPATCHTHREADID) {

    float3 velocity = _Particle[id].velocity;
    float3 pos = _Particle[id].position + velocity * _DeltaTime;

    if(abs(pos.x) > _ParticleRange) {
        velocity.x *= -1;
        pos.x = _Particle[id].position.x + velocity.x * _DeltaTime;
    }

    if(abs(pos.y) > _ParticleRange) {
        velocity.y *= -1;
        pos.y = _Particle[id].position.y + velocity.y * _DeltaTime;
    }

    if(abs(pos.z) > _ParticleRange) {
        velocity.z *= -1;
        pos.z = _Particle[id].position.z + velocity.z * _DeltaTime;
    }

    _Particle[id].position = pos;
    _Particle[id].velocity = velocity;
}


int _SwapBit;
int _UpperBit;

// パーティクルを格子のインデックスでソート
[numthreads(64, 1, 1)]
void ParallelBitonic(uint id : SV_DISPATCHTHREADID) {

    int low = id & (_SwapBit - 1);
    int swapPosX = (id << 1) - low;
    int swapPosY = swapPosX | _SwapBit;

    bool isUpper = (swapPosX & _UpperBit) == 0;

    if((_ParticleField[swapPosX].y > _ParticleField[swapPosY].y) == isUpper) {
        int2 temp = _ParticleField[swapPosX];
        _ParticleField[swapPosX] = _ParticleField[swapPosY];
        _ParticleField[swapPosY] = temp;
    }

}


[numthreads(64, 1, 1)]
void RearrengeBuffer(uint id : SV_DISPATCHTHREADID) {

    uint sortedId = _ParticleField[id].x;
    _SortedParticle[id] = _Particle[sortedId];

}


[numthreads(64, 1, 1)]
void CopyBuffer(uint id : SV_DISPATCHTHREADID) {

    _Particle[id] = _SortedParticle[id];

}


[numthreads(64, 1, 1)]
void SearchNeighbor(uint id : SV_DISPATCHTHREADID) {

    float3 pos = _Particle[id].position;
    int3 gridIndex3 = CalcGridIndex3(pos);
    int pX = gridIndex3.x;
    int pY = gridIndex3.y;
    int pZ = gridIndex3.z;

    int gridIndex = CalcIndex3ToIndex(pX, pY, pZ);

    _Particle[id].color = float4(0, 0, 0, 1);

    if(_SearchIndex != gridIndex)
        return;

    // 自身の格子と周りの格子を調べる
    for(int z = max(pZ - 1, 0); z <= min(pZ + 1, _GridRange - 1); ++z) {
        for(int y = max(pY - 1, 0); y <= min(pY + 1, _GridRange - 1); ++y) {
            for(int x = max(pX - 1, 0); x <= min(pX + 1, _GridRange - 1); ++x) {

                int index = CalcIndex3ToIndex(x, y, z);
                int begin = _GridInfo[index].x;
                int end = _GridInfo[index].y;

                /* 格子で色をつける */
                if(_SearchIndex == index) {

                    for(int i = begin; i < end; ++i) {
                        int pI = _ParticleField[i].x;
                        _Particle[pI].color = float4(1, 0, 0, 1);
                    }

                }
                else {

                    for(int i = begin; i < end; ++i) {
                        int pI = _ParticleField[i].x;
                        _Particle[pI].color = float4(0, 1, 1, 1);
                    }

                }

            }
        }
    }

}


ParticleShader.shader

Shader "Unlit/Particle"
{
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 100

        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag

            #include "UnityCG.cginc"

            struct Particle
            {
                float3 position;
                float4 color;
                float scale;
                float3 velocity;
            };

            struct appdata
            {
                float4 vertex : POSITION;
            };

            struct v2f
            {
                float4 vertex : SV_POSITION;
                float4 color : COLOR;
            };

            StructuredBuffer<Particle> _ParticleBuffer;

            v2f vert (appdata v, uint instanceId : SV_InstanceID)
            {
                Particle p = _ParticleBuffer[instanceId];

                v2f o;

                float3 pos = (v.vertex.xyz * p.scale) + p.position;
                o.vertex = mul(UNITY_MATRIX_VP, float4(pos, 1.0));

                o.color = p.color;

                return o;
            }

            fixed4 frag (v2f i) : SV_Target
            {
                return i.color;
            }

            ENDCG
        }
    }
}



解説

各パーティクルが属する格子を求める

各パーティクルの座標をもとに格子のindexを求めてパーティクルIDとペアで保存します。 また、3次元のindexを1次元のindexに直して保存しています。(参考:多次元配列を1次元配列に変換する計算 - なおしのこれまで、これから)

int CalcGridIndex(float3 particlePosition) {

    int GridX = (particlePosition.x - _MinGrid) / _GridSize;
    int GridY = (particlePosition.y - _MinGrid) / _GridSize;
    int GridZ = (particlePosition.z - _MinGrid) / _GridSize;

    return GridX + GridY*_GridRange + GridZ*_GridRange*_GridRange;

}

// パーティクルの座標をフィールドのindexにする計算
[numthreads(64, 1, 1)]
void CalcParticleGrid(uint id : SV_DISPATCHTHREADID) {

    int index = CalcGridIndex(_Particle[id].position);
    _ParticleField[id] = int2(id, index);
    
} 



格子のインデックスをもとにパーティクルのソート

BitonicSortを使ってソートを行います。実装は次の記事を参考にさせて頂きました。(仮眠プログラマーのつぶやき : UnityでGPGPU応用編 バイトニックソートを高速化)

     // Particle.cs
    private void BitonicSort() {

        int logCount = (int)Mathf.Log(m_instanceCount, 2);

        for(int i = 0; i < logCount; ++i) {

            for(int j = 0; j <= i; j++) {

                int swapBit = 1 << (i - j);
                int upperBit = 2 << i;

                m_particleCalculator.SetInt(m_swapBitPropId, swapBit);
                m_particleCalculator.SetInt(m_upperBitPropId, upperBit);
                m_particleCalculator.Dispatch(m_bitonicSortKernel,
                                             m_particleGroupSize.x,
                                             m_particleGroupSize.y,
                                             m_particleGroupSize.z);
            }
        }
    }


// ParticleCalculator.compute
int _SwapBit;
int _UpperBit;

// パーティクルを格子のインデックスでソート
[numthreads(64, 1, 1)]
void ParallelBitonic(uint id : SV_DISPATCHTHREADID) {

    int low = id & (_SwapBit - 1);
    int swapPosX = (id << 1) - low;
    int swapPosY = swapPosX | _SwapBit;

    bool isUpper = (swapPosX & _UpperBit) == 0;

    if((_ParticleField[swapPosX].y > _ParticleField[swapPosY].y) == isUpper) {
        int2 temp = _ParticleField[swapPosX];
        _ParticleField[swapPosX] = _ParticleField[swapPosY];
        _ParticleField[swapPosY] = temp;
    }

}



各格子の始点と終点を求める

初めに格子の始点と終点を初期化します。 この処理だけGridの総数分行います。

// インデックス情報の初期化
[numthreads(64, 1, 1)]
void InitializeGridInfo(uint id : SV_DISPATCHTHREADID) {

    if(id >= _GridTotalRange)
        return;

    _GridInfo[id] = int2(0xffffff, 0xffffff);

}


次にソートされたパーティクルIDと格子indexのペアを調べます。あるidに入ってる格子indexと前後の格子indexを比較して、前のものと異なればbegin = id、後のものと異なれば end = id+1とします。

[numthreads(64, 1, 1)]
void BuildGridIndexRange(uint id : SV_DISPATCHTHREADID) {

    uint prevId = (id == 0) ? _ParticleCount - 1 : id - 1;
    uint nextId = (id == _ParticleCount - 1) ? 0 : id + 1;

    int index = _ParticleField[id].y;
    int prevIndex = _ParticleField[prevId].y;
    int nextIndex = _ParticleField[nextId].y;

    // 前の格子indexと異なればidがbeginとなる
    if(index != prevIndex) {
        _GridInfo[index].x = id;
    }

    // 後の格子indexと異なればid + 1がendとなる
    if(index != nextIndex) {
        _GridInfo[index].y = id + 1;
    }

}



パーティクルを格子のインデックス順に並び替える

一度パーティクル情報を別のBufferに移し変えてから、元のBufferに入れ直してパーティクルを格子のインデックス順に並び替えています。

[numthreads(64, 1, 1)]
void RearrengeBuffer(uint id : SV_DISPATCHTHREADID) {

    uint sortedId = _ParticleField[id].x;
    _SortedParticle[id] = _Particle[sortedId];

}


[numthreads(64, 1, 1)]
void CopyBuffer(uint id : SV_DISPATCHTHREADID) {

    _Particle[id] = _SortedParticle[id];

}



近傍探索

パーティクルの座標から格子のインデックスx, y, zを求めて、それをもとに周りの格子を求めまて、先ほど求めたbegin, endをもとにその格子に属するパーティクルを調べます。

ここでは調べる格子indexに入っているものを赤、周りの格子に入っているものを青色にしました。

[numthreads(64, 1, 1)]
void SearchNeighbor(uint id : SV_DISPATCHTHREADID) {

    float3 pos = _Particle[id].position;
    int3 gridIndex3 = CalcGridIndex3(pos);
    int pX = gridIndex3.x;
    int pY = gridIndex3.y;
    int pZ = gridIndex3.z;

    int gridIndex = CalcIndex3ToIndex(pX, pY, pZ);

    _Particle[id].color = float4(0, 0, 0, 1);

    if(_SearchIndex != gridIndex)
        return;

    // 自身の格子と周りの格子を調べる
    for(int z = max(pZ - 1, 0); z <= min(pZ + 1, _GridRange - 1); ++z) {
        for(int y = max(pY - 1, 0); y <= min(pY + 1, _GridRange - 1); ++y) {
            for(int x = max(pX - 1, 0); x <= min(pX + 1, _GridRange - 1); ++x) {

                int index = CalcIndex3ToIndex(x, y, z);
                int begin = _GridInfo[index].x;
                int end = _GridInfo[index].y;

                /* 格子で色をつける */
                if(_SearchIndex == index) {

                    for(int i = begin; i < end; ++i) {
                        int pI = _ParticleField[i].x;
                        _Particle[pI].color = float4(1, 0, 0, 1);
                    }

                }
                else {

                    for(int i = begin; i < end; ++i) {
                        int pI = _ParticleField[i].x;
                        _Particle[pI].color = float4(0, 1, 1, 1);
                    }

                }

            }
        }
    }

}



現在の課題

近傍探索では周りの格子のみ青色になるはずですが、それ以外の箇所でも青色になるパーティクルがありました。恐らく、パーティクルの入れ替えが終わる前に近傍探索を行っていることが原因ではないかと考えています。

これの解決方法はまだ模索中ですのでご了承ください。

f:id:vxd-naoshi-19961205-maro:20210116182857p:plain



また、近傍探索の実装をすると1つのプログラムファイルが大きくなってしまったので参考記事のように別クラスを作成して処理を分割したいと考えています。



最後に

今回近傍探索の実装を行いましたが、それとともにComputeShaderのクラス設計なども課題であると感じました。

これからこの近傍探索を使ってBoidsアルゴリズムの実装をしていく予定ですが、少し時間が掛かりそうなので気長にやっていこうと思います。



参考

重複になりますが、こちらの記事を参考に実装させて頂きました。

qiita.com


また、一部実装や内容で以下の記事も参考にさせて頂きました。

blog.livedoor.jp

blog.livedoor.jp

C++とOpenMPによるバイトニックソート

始めに

今回はバイトニックソートについてまとめていこうと思います。

参考にした記事がUnityとComputeShaderやJavaScriptでの実装だったので同じ言語だとあまり面白くならないのでC++を使用しました。


今回は図が多く長いので目次を用意します。



バイトニックソートについて

バイトニックソートはソートの並列アルゴリズムの一つだそうです。(wikipediaより引用)

計算量は O(n) = n { \log^ 2 n }クイックソートなどの高速なアルゴリズムと比べて少し大きいですが、n の部分の計算が並列で行えることからかなり速くソートすることが出来ます。


以下の動画の4:53あたりからバイトニックソートの処理を視覚的に見ることが出来ます。(純粋にソートアルゴリズムの勉強にもおすすめです)

少しだけ処理の順番が違いますが、とても参考になると思います。


私なりの考えになりますが、バイトニックソートには2つの流れがあると考えています。


大きな流れ

始めにこのような長さ2^ 3 = 8の数列があるとします。

f:id:vxd-naoshi-19961205-maro:20201229231825p:plain



この数列を2つずつ昇順、降順を繰り返した数列に整列させます。 赤が昇順、青が降順です。

f:id:vxd-naoshi-19961205-maro:20201229231837p:plain



次にこの数列を4つずつ昇順、降順を繰り返した配列に整列させます。

f:id:vxd-naoshi-19961205-maro:20201229231847p:plain



最後に全体を昇順に並べてソートが完了します。

f:id:vxd-naoshi-19961205-maro:20201229231859p:plain



このように2, 4, 8, 16. . . .と2の累乗の長さで昇順と降順を繰り返してソートを行っていきます。なので配列の長さが2の累乗でなくてはいけません。

また、この昇順降順の切り替えは数列の長さをnとすると \log_{2} n回になります。

ここの流れはマージソートを想像するとわかりやすいかもしれません。



小さな流れ

ではこの昇順と降順を作っていく流れを解説していきます。

2つずつの昇順降順

2つずつ昇順と降順を作る際は隣同士(距離1)を比較していきます。 赤が昇順、青が降順です。

f:id:vxd-naoshi-19961205-maro:20201229231937p:plain



4つずつの昇順降順

では4つずつ昇順と降順を作っていきます。始めに0番目と2番目、1番目と3番目と昇順、降順に並び替えます。すなわち距離2で要素を比較します。

f:id:vxd-naoshi-19961205-maro:20201229231948p:plain


次に隣同士(距離1)を比較して4つずつ昇順と降順を繰り返した数列になります。

f:id:vxd-naoshi-19961205-maro:20201229231958p:plain



8つずつの昇順降順

最後に2つの昇順と降順を1つにします。このとき長さ8の昇順を作っていくので距離4で並びかえを行います。なので、0番目と4番目、1番目と5番目と比較を行っていきます。

f:id:vxd-naoshi-19961205-maro:20201229232010p:plain


次に距離2で0番目と2番目、1番目と3番目とを比較していきます。

f:id:vxd-naoshi-19961205-maro:20201229232020p:plain


最後に距離1で比較してソートが完了します。

f:id:vxd-naoshi-19961205-maro:20201229232033p:plain


以上から2つずつの昇順降順を作る際は距離1で、4つずつの場合は距離2, 1で、8つずつの場合は4, 2, 1となっていることがわかります。

つまり、 2^ nずつの昇順降順を作る場合は距離を 2^ {n - 1}, 2^ {n - 2}, . . . , 2 ^0と距離を対数的に減らしながら整列させます。



bitで見るバイトニックソート

長さ8の配列があるとします。それのインデックスをbitで表現してバイトニックソートを考えていきます。


f:id:vxd-naoshi-19961205-maro:20201229224842p:plain


大きな流れ

始めに2つずつで昇順降順にしたとき、2bit目が0の場合に昇順、1の場合に降順になっていることがわかります。

4つずつの場合は3bit目が0か1で昇順降順を確認することが出来ます。

8つずつの場合は少し見づらいかもしれませんが、4bit目がすべて0であることから昇順になることがわかります。

f:id:vxd-naoshi-19961205-maro:20201229230321p:plain



よって、nこずつの昇順降順を確かめる場合は \log_{2} n + 1番目のbitが1か0かを見ればよいです。



小さな流れ

次に比較する際のbitを見ていきます。

距離1(001)の場合は1bit目以外が一致しているもの同士を比較します。

距離2(010)の場合は2bit目以外が一致しているもの同士を比較します。

距離4(100)の場合は3bit目以外が一致しているもの同士を比較します。

f:id:vxd-naoshi-19961205-maro:20201229235945p:plain


よって、距離が2^ nの場合は \log_{2} n + 1bit以外が一致しているものを比較していることがわかります。

また、Xとなっているbitを抜いてみると00, 01, 10, 11と一つずつ値が増えていることがわかります。



C++での実装

ではC++でバイトニックソートを実装していきます。

実装についてはこれらの記事を参考にさせて頂きました。

バイトニックソートの実装を理解する - e.blog

仮眠プログラマーのつぶやき : UnityでGPGPU応用編 バイトニックソートを高速化


以下、プログラムになります。

#include <cmath>
#include <iostream>
#include <vector>
#include <random>
using namespace std;

constexpr int shiftsize = 3;
constexpr int v_size = 1 << shiftsize;

vector<int> v(v_size);

void Randamize(int seed) {

    random_device rd;
    mt19937 mt(seed);

    for (int i = 0; i < v_size; ++i) 
        v[i] = mt();

}


void Swapping(int i, int j) {

    int swapDistance = 1 << (i - j);

    // 距離swapDistanceで並び替え
    for(int v_i = 0; v_i < v_size/2; ++v_i) {

        int low = v_i & (swapDistance - 1);
        int swapPos = (v_i << 1) - low;

        bool isUpper = (swapPos & (2 << i)) == 0;

        if((v[swapPos] > v[swapPos | swapDistance]) == isUpper)
            swap(v[swapPos], v[swapPos | swapDistance]);

    }

}

void BitonicSort() {

    int log_vsize = log2(v_size);

    // 大きな流れ
    for(int i = 0; i < log_vsize; ++i) {

        // 小さな流れ
        for(int j = 0; j < i + 1; ++j) {

            Swapping(i, j);

        }
    
    }
}

void PrintV() {

    for (int i = 0; i < v.size(); ++i)
        cout << v[i] << endl;
    cout << endl;

}

int main() {

    // 配列にランダムな値を入れる
    // 引数は乱数のseed
    Randamize(1);

    // ソート前の配列を確認
    cout << "before sort" << endl;
    PrintV();

    BitonicSort();

    // ソート後の配列を確認
    cout << "after sort" << endl;
    PrintV();

    return 0;
}



解説

まず、大きな流れの回数を求めます。

昇順と降順の長さは 2^ 1, 2^ 2, 2^ 3, . . .と2の累乗になっていることから、配列の長さを nとすると \log_{2} n回大きな流れをすることがわかります。

    int log_vsize = log2(v_size);

    // 大きな流れ
    for(int i = 0; i < log_vsize; ++i) {



次に、小さな流れの回数を求めます。

昇順と降順の長さが2の場合は入れ替える距離は1だけなので1回。

昇順と降順の長さが4の場合は入れ替える距離は2, 1なので2回。

昇順と降順の長さが8の場合は入れ替える距離は4, 2, 1なので3回。


よって、長さが 2^ 1のときは \log_{2} {2^ 1} = 1回、

長さが 2^ 2のときは \log_{2} {2^ 2} = 2回、

長さが 2^ 3のときは \log_{2} {2^ 3} = 3回となっていることがわかります。

なので並び替えを行う小さな流れの回数は1, 2, 3. . . と増えてきます。

    int log_vsize = log2(v_size);

    // 大きな流れ
    for(int i = 0; i < log_vsize; ++i) {

        // 小さな流れ
        for(int j = 0; j < i + 1; ++j) {

            Swapping(i, j);

        }
    
    }



次に、距離の求め方について解説します。

距離はn回目の大きな流れでは n bitから1bitまで距離を減らしながらソートするので1 << (大きな流れの回数 - 小さな流れの回数)となります。 プログラムでは i が大きな流れの回数、j が小さな流れの回数です。

    int swapDistance = 1 << (i - j);



では入れ替えを行う要素について解説していきます。

入れ替えを行う要素同士の個数は 配列の個数 / 2となります。もし、ここを配列の個数にしてしまうと余分にもう1回比較を行ってしまいます。

void Swapping(int i, int j) {

    int swapDistance = 1 << (i - j);

    // 距離swapDistanceで並び替え
    for(int v_i = 0; v_i < v_size/2; ++v_i) {


では、比較する要素を求めていきます。

始めに距離のbitが0の要素を求めます。v_iをbit列に直すと0, 1, 10, 11, 100, 101. . .となります。ここで距離のbit以下のv_iを取り出します。

        int low = v_i & (swapDistance - 1);


次にv_iを右に1つシフトしたbit列を作り、先ほど取り出した距離以下のbit列を引くことで距離のbitが0になったbit列を求めることができます。

        int swapPos = (v_i << 1) - low;


例として比較する距離 swapDistance が2 (010)だった場合、これらのbit列が求まります。この求められたbit列に距離をOR演算することで比較されるbit列を求めることができます。

swapDistance - 1 = 010 - 001 = 001

v_i = 000 :  
low = 000 & 001 = 000
swapPos = (v_i << 1) - low = 000 - 000 = 000
swapPos | swapDistance = 000 | 010 = 010

v_i = 001:
low = 001 & 001 = 001
swapPos = (v_i << 1) - low = 010 - 001 = 001
swapPos | swapDistance = 001 | 010 = 011

v_i = 010:
low = 010 & 001 = 000
swapPos = (v_i << 1) - low = 100 - 000 = 100
swapPos | swapDistance = 100 | 010 = 110

v_i = 011:
low = 011 & 001 = 001
swapPos = (v_i << 1) - low = 110 - 001 = 101 
swapPos | swapDistance = 101 | 010 = 111



最後に昇順か降順かの計算について解説します。

n回目の大きな流れでは(n + 1)bit目が1か0かで昇順降順を判断します。

プログラムでは i が大きな流れの回数なので 1 << (i + 1) = 2 << i を見れば、昇順か降順かを判断できます。 ここでは0の場合が昇順、1の場合が降順になります。

        bool isUpper = (swapPos & (2 << i)) == 0;


あとは、要素同士を比較してそれが昇順もしくは降順になっていなければ並び変えます。

        if((v[swapPos] > v[swapPos | swapDistance]) == isUpper)
            swap(v[swapPos], v[swapPos | swapDistance]);



実行結果

実行した結果は次の通りです。出力からソートされていることがわかります。

before sort
1791095845
-12091157
-1201197172
-289663928
491263
550290313
1298508491
-4120955

after sort
-1201197172
-289663928
-12091157
-4120955
491263
550290313
1298508491
1791095845



OpenMPによる並列処理

バイトニックソートで要素同士の比較、並び替えの処理は独立して行うことが可能です。 なので要素同士の比較はばらばらに実行しても結果は同じになります。

f:id:vxd-naoshi-19961205-maro:20201230165318p:plain


この処理をOpenMPを使って並列処理をすることが可能です。

    // 距離swapDistanceで並び替え
    // 並列処理
#pragma omp parallel for
    for(int v_i = 0; v_i < v_size/2; ++v_i) {

        int low = v_i & (swapDistance - 1);
        int swapPos = (v_i << 1) - low;

        bool isUpper = (swapPos & (2 << i)) == 0;

        if((v[swapPos] > v[swapPos | swapDistance]) == isUpper)
            swap(v[swapPos], v[swapPos | swapDistance]);

    }



計測、std::sortとの比較

一度バイトニックソートがどれだけ速いかをstd::sortと比較してみます。

計測用のプログラムは以下の通りです。 長さ 2^ {20}=1048576の配列を作り、ソートする前に一度ランダムな値を入れてからソートを開始します。

また、ソートするたびに掛かった時間を計測し最小値、最大値、平均値を出します。

C++プログラム

#include <cmath>
#include <iostream>
#include <random>
#include <vector>
#include <chrono>
#include <algorithm>
#include <limits>
#include <cfloat>
using namespace std;

constexpr int shiftsize = 20;
constexpr int v_size = 1 << shiftsize;

vector<int> v(v_size);

void Randamize(int x) {

    random_device rd;
    mt19937 mt(x);

    for (int i = 0; i < v_size; ++i) 
        v[i] = mt();

}

void Swapping(int i, int j) {
    int swapDistance = 1 << (i - j);
    
  // 距離swapDistanceで並び替え
  // 並列処理
#pragma omp parallel for
    for (int v_i = 0; v_i < v_size / 2; ++v_i) {

        int low = v_i & (swapDistance - 1);
        int swapPos = (v_i << 1) - low;

        bool isUpper = (swapPos & (2 << i)) == 0;

        if ((v[swapPos] > v[swapPos | swapDistance]) == isUpper)
            swap(v[swapPos], v[swapPos | swapDistance]);

    }
}

void BitonicSort() {

    int log_vsize = log2(v_size);

    for(int i = 0; i < log_vsize; ++i) {

        for(int j = 0; j < i + 1; ++j) {

            Swapping(i, j);

        }

    }

}

void PrintV() {

    for (int i = 0; i < v.size(); ++i)
        cout << v[i] << endl;
    cout << endl;

}

int main() {

    constexpr int count = 10;
    double minTime_b = DBL_MAX, maxTime_b = DBL_MIN, averageTime_b = 0;
    double minTime_s = DBL_MAX, maxTime_s = DBL_MIN, averageTime_s = 0;

    cout << "bitonicsort" << endl;
    // バイトニックソート計測
    for(int i = 0; i < count; ++i) {
        Randamize(i);

        auto start = chrono::system_clock::now();

        BitonicSort();

        auto end = chrono::system_clock::now();

        double elapsed = chrono::duration_cast<chrono::milliseconds>(end - start).count();

        cout << i << ":" << elapsed << endl;

        minTime_b = min(minTime_b, elapsed);
        maxTime_b = max(maxTime_b, elapsed);
        averageTime_b += elapsed;

    }

    cout << "---result---" << endl;
    cout << "min time:" << minTime_b << endl;
    cout << "max time:" << maxTime_b << endl;
    cout << "ave time:" << averageTime_b/count << endl;
    cout << endl;

    cout << "std::sort" << endl;
    // std::sort計測
    for (int i = 0; i < count; ++i) {
        Randamize(i);

        auto start = chrono::system_clock::now();

        sort(v.begin(), v.end());

        auto end = chrono::system_clock::now();

        double elapsed = chrono::duration_cast<chrono::milliseconds>(end - start).count();

        cout << i << ":" << elapsed << endl;

        minTime_s = min(minTime_s, elapsed);
        maxTime_s = max(maxTime_s, elapsed);
        averageTime_s += elapsed;
    }

    cout << "---result---" << endl;
    cout << "min time:" << minTime_s << endl;
    cout << "max time:" << maxTime_s << endl;
    cout << "ave time:" << averageTime_s / count << endl;
    cout << endl;

    return 0;
}


また、実行環境はwindowsでCPUはcorei7-9750H、コンパイラgccで実行時はなるべくタスクを実行しないようにしました。



結果

結果は以下の通りです。 バイトニックソートは少しばらつきはありますが、std::sortに比べて約1.2倍高速でした。

ただし、並列処理を行わなかった場合は実行時間が6倍ほど増えるので普通にソートする場合はstd::sortの方が高速です。

bitonicsort
0:256
1:260
2:273
3:293
4:265
5:276
6:271
7:261
8:216
9:219
---result---
min time:216
max time:293
ave time:259

std::sort
0:327
1:321
2:319
3:322
4:320
5:322
6:321
7:321
8:321
9:320
---result---
min time:319
max time:327
ave time:321.4



最後に

このソートを知ったきっかけとしてはniconicoの動画になりますが、そのときは「知らなくてもいいでしょ」と思っていました。

私の感想になりますが、ここまで理解するのに苦しんだソートアルゴリズムはなかったです。しかし、結果としてstd::sortよりも速くソート出来ることがわかりました。

これからになりますが、Boidsアルゴリズムの近傍探索などで使用してみたいです。



参考

重複になりますがこれらの記事を参考にさせて頂きました。

edom18.hateblo.jp

blog.livedoor.jp



また、OpenMPに関しては簡単ではありますがこの記事を参考にして使用しました。

www.sanko-shoko.net


OpenMPが実行できなかった際は英語になりますがこのサイトで調べて実行することが出来ました。

c++ - mingwでウィンドウ上でopenmpを使用する。 -lpthreadを見つけることができません



バーチャルでLTをしてみての感想(Online Game-Tech2020、stdout2020)

f:id:vxd-naoshi-19961205-maro:20201217224749p:plain

始めに

リモートでのLTイベントのOnline Game-Tech LT 2020が11/1に、stdout2020が12/13に開催されました。 この2つのイベントにバーチャル登壇をしてきました。その感想をまとめていこうと思います。


このバーチャル登壇で使用したものは以下の記事にまとめています。

shitakami.hatenablog.com



Online Game-Tech 2020について

これは@Shirotsuさんが開催したイベントです。

LTのテーマはゲームに関する技術でした。また、参加者の半分が一般、残りが学生で現場で働く方のお話が聞ける貴重なイベントとなりました。



Online Game-Tech LT 2020

私の発表は4:58:45~になります。



stdout2020について

こちらのイベントはカズ之助さんマヤミトさんが開催したLTイベントになります。Joken(日大工学部)とZli(会津大学)による合同LT会で学生メインのイベントなりました。が、登壇者は様々な大学の方がいらっしゃいました。

LTのテーマは技術系であればなんでもいい!でかなりゆるいものでした。

このようなイベントの参加が初めてだったので技術の幅もかなり広く、とても楽しいLTイベントでした。



#stdout2020 学生LT会

私の発表は42:20~になります。

リモートのLTイベントに参加した理由

率直に言いますと、作成したリモートプレゼンツール「Kawaii-Jyuku」を使って色んな人からの感想が聞きたいっていうのが参加理由です。

元々は勉強会で使う予定でそこまで公開する予定はありませんでした。

しかし、もうちょっと面白いものにしたいと考えて、色々の人の意見を聞いて必要な機能や修正すべき箇所を見つけるために参加を決意致しました。



LTの内容について

Online Game-Tech LT 2020ではGPUで遊ぶLTをしました。もともと夏休みの時期にUnityでGPUを使った機能の勉強をしていたのでそれをまとめてLTの内容にしました。

また、GPUを使った機能は派手なものが多くLT向けじゃないかなというのも理由の一つです。


次にstdout2020ですが、「LTの裏ではどんなふうになっているか知りたい」などの感想を頂いたのでいっそのことリモートプレゼンツールの話をしようと決めました。

ただし、裏側を全部みせてしまうとキャラクター性やバーチャル感がなくなってしまうと思ったので、概要的な内容に決めました。


最後に「すげぇ!」って言われたかったのでBoidsアルゴリズムを実装してアバターと共存させる機能を作成して、発表のトリにしました。

https://cdn-ak.f.st-hatena.com/images/fotolife/v/vxd-naoshi-19961205-maro/20201112/20201112181741.png



LT会に向けてやったこと

リモートプレゼンツールは部活の勉強会で1か月半ほど使用したので正しく動作することは確認済みでした。

ただし、コメント機能が未検証だったので何度かテストをしました。

テスト方法はTwitterを開き、トレンドを調べてVTuber関係のタグがあればそれを設定してタグがついたTweetをコメントで流してみるという流れでした。

テストをしようとした際に丁度ホロライブの星街すいせいの50万人記念ライブが開催されていました。


【3DLIVE】STARt IN THE SCREEEN!【#星街すいせい50万人記念ライブ】

この機会をありがたく使わせていただいて、コメント機能の動作を確認するとなんと3fps!となりました。

100以上のtweetが1秒間に発信されていたためそもそも動作しませんでした。


そこで、tweet数を制限し画面に10個までしか表示しないようにしました(実際にそれ以上のtweetは起こることはないと思っていましたが. . .)。


あとはBoidsアルゴリズムがかなりCPUとGPUのリソースを食うので出来るだけ処理を軽くするべく、表示するモニターの個数を減らしたり、プレゼンする際は別のアプリケーションを実行しないようにと気をつけました。



参加しての感想

正直、酷評やそもそもバーチャルキャラクターを受け付けないなどの感想がくることを覚悟していましたが、とても好評だと感じました。

自分が頑張って実装した部分すべてを良い!って言ってくださったのがとても嬉しかったです。 特にこだわって作成したUIもカッコいい!と言って頂けました。


参加してよかったこと

個人的に参加してよかったことはYoutubeに自分のプレゼンがアーカイブとして残ったことです。

このリモートプレゼンツールを使ってて問題だと感じた部分が「私が見てる風景」と「視聴者が見ている私」の間にギャップがあるということでした。

私がアバターを操作してプレゼンをしていますが、実際はちゃんと動いていないかもしれない、はなしと挙動がずれているかもしれないとずっと思っておりました。

そしてこの問題を解決する方法が見て頂いた方から感想を頂くしかないと思っていました。

しかし、アーカイブで自分のプレゼンを見ることでちゃんと正しく動作していることを確認することが出来ました。加えて、どこを修正すべきかどんな機能があればもっと面白くなるかを考える材料になりました。



今後について

LTを通してこれから修正したい部分は以下の通りです。 これらを少しずつやっていこうと思います。

  • スライドを別アプリから取得できるようにする
  • アプリをビルドする
  • 背景やステージを作る
  • ホワイトボードの改修とスライドにも書き込めるようにする
  • ちょっとしたリファクタリング


そして、気になる次の予定ですが残念ながら数か月間は使用することはないと思います。

LT会に登壇する際にいくつか必要な機材を部活からお借りしていましたが、それももう返さなくてはいけません。また、来年は社会人になる予定なのでどれだけ時間を作れるかが未定となっております。

ただし、これらのイベントを通してもう一度やってみたいという気持ちが強くなったのでこれから使用しないことはないでしょう。


これからについては全くわかりませんが、どこかしらでまた何かをしたいと思っております。



最後に

イベントを開催してくだった方々や見てくださった方々、本当にありがとうございました!

多次元配列を1次元配列に変換する計算

始めに

過去記事の内容で3次元配列を1次元配列にする解説が少しわかりにくいとなり、補足記事を書こうと思います。

shitakami.hatenablog.com

2次元配列を1次元配列にする

まずは2次元配列から1次元配列に変換することから始めます。

例として高さ3、幅5の2次配列があるとします。 横方向がx軸、縦方向がy軸として要素は0から始まるものとします。

下の図は各マス目をyとxで表したものです。

f:id:vxd-naoshi-19961205-maro:20201113190517p:plain


上の図から、(y, x) = (2, 1)とすれば上から3番目の左から2番目の要素を表します。



この表をyとxで表しましたが、次はx'だけで表したいと思います。要素の数え方は一番上の左から右方向へ数えて、一番右まで数えきったら下に移動して同様に左から右へ数えていきます。

下の図は各マス目をx'で表したものです。

f:id:vxd-naoshi-19961205-maro:20201113192345p:plain



ではここから(y, x)を使ってx'を求める計算をしていこうと思います。

x'の表からわかることをいくつか列挙していきます。

あるマスから右隣に移動するとx' + 1となり、下に移動するとx' + 5となります。

f:id:vxd-naoshi-19961205-maro:20201113194915p:plain


よってxが1増えるとx'+1となり、yが1増えるとx'+5となるので  x' = 5y + xとなります。

例えば、(y, x) = (2, 3)だった場合は x' = 5 \times 2 + 3= 13で(y, x)の表の位置とx'の表の位置が一致していることがわかります。



ここで、x' = 5y + xの5yの部分に注目すると5が幅の大きさに一致していることがわかります。

よって、高さheight、幅widthの2次元配列(y, x)を1次元配列x'にする計算は次のようになります。


x' = y \times width + x



3次元配列を1次元配列にする

例として高さ3、幅5、奥行3の3次元配列があるとします。

下の表(z, y, x)は手前がz=0で一つ奥の表に行くごとにzが増えていきます。

f:id:vxd-naoshi-19961205-maro:20201113225446p:plain


同様にこれもx'の表に直していきます。 1つの表が埋まれば次の表に移って数え続けます。

f:id:vxd-naoshi-19961205-maro:20201113230344p:plain



ここから(z, y, x)を使ってx'を求める計算を考えていきます。 xとyについては2次元配列と同様なので省きます。

2次元配列の表が次のものに移るとx'の値はx'+15されていることがわかります。

f:id:vxd-naoshi-19961205-maro:20201113230903p:plain



よって、zが1増えるとx'は15増えるので x' = 15z + 5y + xとなります。

例えば(z, y, x) = (2, 0, 3)とすると x' = 15 \times 2 + 5 \times 0 + 3 = 33となり(z, y, x)の表とx'の表の位置が一致していることが確かめられます。


ここで15zの15に注目すると高さ×幅=3×5=15と一致していることがわかります。

よって、奥行depth、高さheight、幅widthとすると3次元配列(z, y, x)を1次元配列x'に変換する数式は次のようになります。


x' = z(height \times width) + y \times width +x



多次元配列を1次元配列にする

ここまでこれば何となく法則性がわかってきたと思います。

例えば5次元配列があったとしてそれぞれの次元の大きさが(s_0, s_1, s_2, s_3, s_4)とすると5次元配列(i_0, i_1, i_2, i_3, i_4)を1次元配列i'に変換する数式は次のようになります。


i' = i_0(s_1 \times s_2 \times s_3 \times s_4) + i_1(s_2 \times s_3 \times s_4) + i_2(s_3 \times s_4) + i_3 \times s_4 + i_4



この数式を一般化したいと思います。(数学詳しくないので正しくないかもしれません!)

n次元配列(i_0, i_1, i_2, . . . i_n)がありそれぞれの次元の大きさを(s_0, s_1, s_2, . . . . s_n)とするとn次元配列を1次元配列i'にすると次のような数式になります。


i' = i_0 \prod_{x=1}^{N} s_x + i_1 \prod_{x=2}^{N} s_x + . . . . i_{n-1} \times s_n + i_n



すべての次元の大きさが同じ場合はn進数になる

ここでn次元配列のすべての次元の大きさがXで同じ場合を考えたいと思います。

このときn次元配列(i_0, i_1, i_2, . . . i_n)を1次元配列i'にすると次にようになります。


i' = i_0 X^{n - 1} + i_1 X^{n - 2} + . . . . + i_{n - 2} X^2 + i_{n - 1} X + i_n


これは実はn進数を10進数に変換する計算と一致します。 こう覚えておくと何となく便利かなと思ってます。

参考:【基本】n進法への変換(整数) | なかけんの数学ノート



最後に

少しだけ雑談しますと、前のブログの解説について後輩と友達に相談して結構議論になりました。ある人はこの内容が積分だと解釈しある人は行列と解釈し、話していて様々な見方が出来る内容なのかなと感じました。

私個人としてはこの内容をn進数と見て考えていまして、そういう風に今回の内容をまとめました。

少し競技プログラミングっぽいような感じの内容かもしれませんが、ときによってはゲームプログラミングに役立つかもしれません。

ベクトル場を使った衝突回避Boidsアルゴリズム

f:id:vxd-naoshi-19961205-maro:20201112181741p:plain

始めに

前回のブログで作成したベクトル場を使って衝突回避をするBoidアルゴリズムを実装したいと思います。

今回のサンプルです。

github.com



過去記事です。

shitakami.hatenablog.com

shitakami.hatenablog.com



やったこと

始めにベクトル場を計算します。 このベクトル場はモデルの法線ベクトルを使用しているので、モデルから離れるようなベクトルになります。

f:id:vxd-naoshi-19961205-maro:20201104112333g:plain



Boidsアルゴリズムにベクトル場の力を加えます。 計算は結合・分離・整列の3つの力にベクトル場の力を加えます。こうすることで個体がモデルに近づくと離れる力が働いて逃げるような動きができると考えました。

f:id:vxd-naoshi-19961205-maro:20201109152826p:plain



実装

かなりの行数があるのでプログラムを開く際は注意してください。

ComputeShader

#pragma kernel UpdateVectorField
#pragma kernel InitializeVectorField
#pragma kernel UpdateBoids

#include "Assets/cgincFiles/rotation.cginc"

#define COUNT 3

////////// ベクトル場用変数 ///////////////////////////////////////////

struct VectorFieldSquare {
    float3 position;
    float3 direction;
};


RWStructuredBuffer<VectorFieldSquare> _VectorFieldBuffer;
StructuredBuffer<float3> _PositionBuffer;
StructuredBuffer<float3> _NormalBuffer;

float _SquareScale;
int _VectorFieldLength;

float3 _MinField;
float3 _MaxField;

float3 _ModelPosition[COUNT];
float3 _ModelEulerAngle[COUNT];
float3 _ModelScale[COUNT];

int _VertexCountsThresholds[COUNT];

/////////////////////////////////////////////////////////////////////

//////////////////// Boidsアルゴリズム用変数 //////////////////////////
struct BoidsData {
    float3 position;
    float3 velocity;
};

RWStructuredBuffer<BoidsData> _BoidsData;

// Boidsパラメータ
float _CohesionForce;
float _SeparationForce;
float _AlignmentForce;

float _CohesionDistance;
float _SeparationDistance;
float _AlignmentDistance;

float _CohesionAngle;
float _SeparationAngle;
float _AlignmentAngle;

float _BoundaryForce;
float _BoundaryRange;
float3 _BoundaryCenter;

float _AvoidForce;

float _MaxVelocity;
float _MinVelocity;

int _InstanceCount;

float _MaxForce;

float _DeltaTime;
float _Epsilon;

/////////////////////////////////////////////////////////////////////

// 速度と座標から角度を求める
float CalcAngle(float3 velocity, float3 posX, float3 posY) {

    float3 vec = posY - posX;

    return acos(dot(normalize(velocity), normalize(vec)));

}

// 距離の2乗を求める
float CalcSqrDistance(float3 posX, float3 posY) {

    float3 vec = posY - posX;

    return dot(vec, vec);

}

// ベクトルの大きさを制限する
float3 limit(float3 vec, float max)
{
    float length = sqrt(dot(vec, vec)); // 大きさ
    return (length > max && length > 0) ? vec.xyz * (max / length) : vec.xyz;
}

[numthreads(8, 1, 1)]
void UpdateVectorField(uint id : SV_DISPATCHTHREADID) {

    int transformIndex = 0;

    for(int i = 0; (int)id < _VertexCountsThresholds[i]; ++i)
        transformIndex = i;

    float3 pos = _PositionBuffer[id];
    float3 normal = _NormalBuffer[id];
    float4x4 rotMat = EulerAnglesToMatrix(_ModelEulerAngle[transformIndex]);    // 回転行列を求める

    pos = mul(rotMat, pos);

    pos.x *= _ModelScale[transformIndex].x;
    pos.y *= _ModelScale[transformIndex].y;
    pos.z *= _ModelScale[transformIndex].z;

    pos += _ModelPosition[transformIndex];

    // 範囲外は計算しない
    if(pos.x < _MinField.x || _MaxField.x <= pos.x ||
       pos.y < _MinField.y || _MaxField.y <= pos.y ||
       pos.z < _MinField.z || _MaxField.z <= pos.z)
        return;

    // 頂点座標がどのマス目にいるか計算してindexを取得
    int3 index = (pos - _MinField) / _SquareScale;

    normal = mul(rotMat, normal);

    _VectorFieldBuffer[(index.x * _VectorFieldLength * _VectorFieldLength +
                        index.y * _VectorFieldLength +
                        index.z)].direction += normal;

    // ベクトルの要素の符号を取得
    int xd = normal.x / abs(normal.x + _Epsilon);
    int yd = normal.y / abs(normal.y + _Epsilon);
    int zd = normal.z / abs(normal.z + _Epsilon);

    if(0 <= index.x + xd && index.x + xd < _VectorFieldLength)
        _VectorFieldBuffer[(index.x + xd) * _VectorFieldLength * _VectorFieldLength +
                           index.y * _VectorFieldLength +
                           index.z].direction += normal * abs(normal.x);
    
    if(0 <= index.y + yd && index.y + yd < _VectorFieldLength)
        _VectorFieldBuffer[index.x * _VectorFieldLength * _VectorFieldLength +
                           (index.y + yd) * _VectorFieldLength +
                           index.z].direction += normal * abs(normal.y);
    
    if(0 <= index.z + zd && index.z + zd < _VectorFieldLength)
        _VectorFieldBuffer[index.x * _VectorFieldLength * _VectorFieldLength +
                           index.y * _VectorFieldLength +
                           index.z + zd].direction += normal * abs(normal.z);

}


[numthreads(8, 1, 1)]
void InitializeVectorField(uint id : SV_DISPATCHTHREADID) {

    _VectorFieldBuffer[id].direction = float3(0, 0, 0);

}

[numthreads(256, 1, 1)]
void UpdateBoids(uint id : SV_DISPATCHTHREADID) {

        float3 posX = _BoidsData[id.x].position;
    float3 velX = _BoidsData[id.x].velocity;

    float3 cohesionPositionSum = float3(0, 0, 0);
    float3 separationPositionSum = float3(0, 0, 0);
    float3 alignmentVelocitySum = float3(0, 0, 0);

    int cohesionCount = 0;
    int alignmentCount = 0;

    for(uint i = 0; i < (uint)_InstanceCount; ++i) {

        // 自身の計算は行わない
        if(i == id.x)
            continue;

        float3 posY = _BoidsData[i].position;
        float3 velY = _BoidsData[i].velocity;

        float sqrDistance = CalcSqrDistance(posX, posY);
        float angle = CalcAngle(velX, posX, posY);

        // 結合
        if(sqrDistance < _CohesionDistance && angle < _CohesionAngle) {
            cohesionPositionSum += posY;
            cohesionCount++;
        }

        // 分離
        if(sqrDistance < _SeparationDistance && angle < _SeparationAngle) {
            separationPositionSum += normalize(posX - posY) / sqrt(sqrDistance);        
        }

        // 整列
        if(sqrDistance < _AlignmentDistance && angle < _AlignmentAngle) {
            alignmentVelocitySum += velY;
            alignmentCount++;
        }

    }

    float3 cohesion = float3(0, 0, 0);
    float3 separation = separationPositionSum;
    float3 alignment = float3(0, 0, 0);
    float3 boundary = float3(0, 0, 0);
    float3 avoid = float3(0, 0, 0);

    if(cohesionCount != 0)
        cohesion = (cohesionPositionSum / (float)cohesionCount - posX) * _CohesionForce;
    
    if(alignmentCount != 0) 
        alignment = (alignmentVelocitySum / (float)alignmentCount - velX) * _AlignmentForce;
    
    separation *= _SeparationForce;

    // 範囲外から出た個体は範囲内に戻る力を加える
    
    float sqrDistFromCenter = dot(_BoundaryCenter - posX, _BoundaryCenter - posX);
    if(sqrDistFromCenter > _BoundaryRange)
        boundary = -_BoundaryForce * (posX - _BoundaryCenter) * (sqrDistFromCenter - _BoundaryRange) / sqrDistFromCenter;

    // ベクトル場内にいる際は、ベクトルの値を取得してその方向に逃げるようにする
    if(_MinField.x <= posX.x && posX.x < _MaxField.x &&
       _MinField.y <= posX.y && posX.y < _MaxField.y &&
       _MinField.z <= posX.z && posX.z < _MaxField.z) {
        int3 index = (posX - _MinField) / _SquareScale;
        avoid = _VectorFieldBuffer[(index.x * _VectorFieldLength * _VectorFieldLength +
                        index.y * _VectorFieldLength +
                        index.z)].direction * _AvoidForce;
       }

    // 結合、分離、整列の力を制限
    cohesion = limit(cohesion, _MaxForce);
    separation = limit(separation, _MaxForce);
    alignment = limit(alignment, _MaxForce);

    velX += (cohesion + separation + alignment + boundary + avoid) * _DeltaTime;

    float velXScale = length(velX);

    // 速度を制限
    if(velXScale < _MinVelocity) {
        velX = _MinVelocity * normalize(velX);
    }
    else if (velXScale > _MaxVelocity) {
        velX = _MaxVelocity * normalize(velX);
    }

    _BoidsData[id.x].velocity = velX;
    _BoidsData[id.x].position += velX;


}


AvoidanceBoids.cs(C#)

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using UnityEngine.Rendering;
using System.Runtime.InteropServices;

public class AvoidanceBoids : MonoBehaviour
{
    [SerializeField]
    private ComputeShader m_avoidanceBoidsSimulator;

    [SerializeField]
    private SkinnedMeshRenderer[] m_skinnedMeshRenderer;
    private Transform[] m_skinnedMeshTransform; 

    [Header("ベクトル場のマス目数")]
    [SerializeField]
    private int m_vectorFieldLength;

    [Header("ベクトル場のマスの大きさ")]
    [SerializeField]
    private float m_squareScale;

    [Header("ベクトル場の中心座標")]
    [SerializeField]
    private Vector3 m_vectorFieldCenter;

    private Mesh m_bakeMesh;
    private List<Vector3> m_vertices = new List<Vector3>();
    private List<Vector3> m_normals = new List<Vector3>();

    private Vector4[] m_meshPositions;
    private Vector4[] m_meshEulerAngles;
    private Vector4[] m_meshScales;
    private int[] m_vertexCountThresholds;

    #region ComputeShaderKernelAndGroupSize
    private int m_initializeVectorFieldKernel;
    private Vector3Int m_initializeVectorFieldGroupSize;

    private int m_updateVectorFieldKernel;
    private Vector3Int m_updateVectorFieldGroupSize;
    #endregion

    #region ComputeBuffers
    private ComputeBuffer m_vectorFieldBuffer;
    private ComputeBuffer m_positionBuffer;
    private ComputeBuffer m_normalBuffer;
    #endregion

    #region Shader_PropertyID
    private readonly int m_positionPropID = Shader.PropertyToID("_PositionBuffer");
    private readonly int m_normalPropID = Shader.PropertyToID("_NormalBuffer");
    private readonly int m_modelPositionPropID = Shader.PropertyToID("_ModelPosition");
    private readonly int m_modelEulerAnglePropID = Shader.PropertyToID("_ModelEulerAngle");
    private readonly int m_modelScalePropID = Shader.PropertyToID("_ModelScale");
    #endregion

    #region Properties
    public ComputeBuffer VectorField { get { return m_vectorFieldBuffer; } }
    public float SquareScale { get { return m_squareScale; } }
    public int SquareCount { get { return m_vectorFieldLength*m_vectorFieldLength*m_vectorFieldLength; } }
    
    public ComputeBuffer BoidsDataBuffer { get { return m_boidsDataBuffer; } }
    public int BoidsInstanceCount { get { return m_instanceCount; } }
    #endregion

    // ベクトル場のマス情報
    private struct VectorFieldSquare {
        public Vector3 position;
        public Vector3 direction;
    }

    [Space(20)]
    [SerializeField]
    private int m_instanceCount;

    [Header("力の強さ")]
    [Header("Boidsモデルのデータ")]
    [SerializeField]
    private float m_cohesionForce;
    [SerializeField]
    private float m_separationForce;
    [SerializeField]
    private float m_alignmentForce;

    [Space(5)]
    [Header("力の働く距離")]
    [SerializeField]
    private float m_cohesionDistance;
    [SerializeField]
    private float m_separationDistance;
    [SerializeField]
    private float m_alignmentDistance;

    [Space(5)]
    [Header("力の働く角度")]
    [SerializeField]
    private float m_cohesionAngle;
    [SerializeField]
    private float m_separationAngle;
    [SerializeField]
    private float m_alignmentAngle;

    [Space(5)]
    [SerializeField]
    private float m_boundaryForce;
    [SerializeField]
    private float m_boundaryRange;
    [SerializeField]
    private Vector3 m_boundaryCenter;

    [SerializeField]
    private float m_avoidForce;

    [Space(5)]
    [SerializeField]
    private float m_maxVelocity;

    [SerializeField]
    private float m_minVelocity;

    [SerializeField]    
    private float m_maxForce;

    private struct BoidsData {
        public Vector3 position;
        public Vector3 velocity;
    }

    private ComputeBuffer m_boidsDataBuffer;

    private int m_updateBoidsKernel;

    private int m_deltaTimeID = Shader.PropertyToID("_DeltaTime");

    private Vector3Int m_updateBoidsGroupSize;

    [Header("レンダラー")]
    [Space(20)]
    [SerializeField]
    private BoxRenderer m_boxRenderer;
    [SerializeField]
    private BoidsRenderer m_boidsRenderer;

    // Start is called before the first frame update
    void Start()
    {
        
        m_meshPositions = new Vector4[m_skinnedMeshRenderer.Length];
        m_meshEulerAngles = new Vector4[m_skinnedMeshRenderer.Length];
        m_meshScales = new Vector4[m_skinnedMeshRenderer.Length];
        m_vertexCountThresholds = new int[m_skinnedMeshRenderer.Length];

        InitializeAvoidanceBoids();
        InitializeBoidsSimulator();

        // GC.Allocを防ぐため
        m_bakeMesh = new Mesh();

        m_skinnedMeshTransform = new Transform[m_skinnedMeshRenderer.Length];
        for(int i = 0; i < m_skinnedMeshRenderer.Length; ++i)
            m_skinnedMeshTransform[i] = m_skinnedMeshRenderer[i].transform;
        // m_skinnedMeshTransform = m_skinnedMeshRenderer.transform;

        // AvoidanceBoidsの初期化が完了したらRendererを初期化する
        m_boxRenderer?.Initialize(this);
        m_boidsRenderer?.Initialize(this);

    }

    // Update is called once per frame
    void Update()
    {
        UpdateVectorField();
        UpdateBoids();
    }

    private void InitializeAvoidanceBoids() {

        // スレッドサイズ取得に使用
        uint x, y, z;

        // すべてのSkinnedMeshRendererの頂点数を求める
        // また頂点数の累積和を求めてComputeShaderの計算で使用する
        int vertexCount = 0;
        for(int i = 0; i < m_skinnedMeshRenderer.Length; ++i) {
            vertexCount += m_skinnedMeshRenderer[i].sharedMesh.vertexCount;
            m_vertexCountThresholds[i] = vertexCount;
        }
        // int vertexCount = m_skinnedMeshRenderer.sharedMesh.vertexCount;

        InitializeVectorFieldBuffer();
        InitializePositionNormalBuffer(vertexCount);

        m_initializeVectorFieldKernel = m_avoidanceBoidsSimulator.FindKernel("InitializeVectorField");
        m_avoidanceBoidsSimulator.GetKernelThreadGroupSizes(m_initializeVectorFieldKernel, out x, out y, out z);
        int squaresCount = m_vectorFieldLength*m_vectorFieldLength*m_vectorFieldLength;
        m_initializeVectorFieldGroupSize = new Vector3Int(squaresCount/(int)x, (int)y, (int)z);
        m_avoidanceBoidsSimulator.SetBuffer(m_initializeVectorFieldKernel, "_VectorFieldBuffer", m_vectorFieldBuffer);

        m_updateVectorFieldKernel = m_avoidanceBoidsSimulator.FindKernel("UpdateVectorField");
        m_avoidanceBoidsSimulator.GetKernelThreadGroupSizes(m_updateVectorFieldKernel, out x, out y, out z);
        
        m_updateVectorFieldGroupSize = new Vector3Int(vertexCount/(int)x, (int)y, (int)z);

        m_avoidanceBoidsSimulator.SetBuffer(m_updateVectorFieldKernel, "_PositionBuffer", m_positionBuffer);
        m_avoidanceBoidsSimulator.SetBuffer(m_updateVectorFieldKernel, "_NormalBuffer", m_normalBuffer);
        m_avoidanceBoidsSimulator.SetBuffer(m_updateVectorFieldKernel, "_VectorFieldBuffer", m_vectorFieldBuffer);
        m_avoidanceBoidsSimulator.SetFloat("_SquareScale", m_squareScale);
        m_avoidanceBoidsSimulator.SetInt("_VectorFieldLength", m_vectorFieldLength);

        float halfLength = m_squareScale * m_vectorFieldLength/2.0f;
        m_avoidanceBoidsSimulator.SetVector("_MinField", m_vectorFieldCenter - new Vector3(halfLength, halfLength, halfLength));
        m_avoidanceBoidsSimulator.SetVector("_MaxField", m_vectorFieldCenter + new Vector3(halfLength, halfLength, halfLength));

    }


    private void InitializeVectorFieldBuffer() {

        int squaresCount = m_vectorFieldLength*m_vectorFieldLength*m_vectorFieldLength;

        VectorFieldSquare[] squares = new VectorFieldSquare[squaresCount];

        float startPosX = m_vectorFieldCenter.x - m_vectorFieldLength * m_squareScale / 2.0f + m_squareScale / 2.0f;
        float startPosY = m_vectorFieldCenter.y - m_vectorFieldLength * m_squareScale / 2.0f + m_squareScale / 2.0f;
        float startPosZ = m_vectorFieldCenter.z - m_vectorFieldLength * m_squareScale / 2.0f + m_squareScale / 2.0f;

        for(int i = 0; i < m_vectorFieldLength; ++i) {
            for(int j = 0; j < m_vectorFieldLength; ++j) {
                for(int k = 0; k < m_vectorFieldLength; ++k) {
                    int index = i*m_vectorFieldLength*m_vectorFieldLength + j*m_vectorFieldLength + k;
                    squares[index].position = new Vector3(startPosX + i * m_squareScale, startPosY + j * m_squareScale, startPosZ + k * m_squareScale);
                    squares[index].direction = Vector3.zero;
                }
            }
        }

        m_vectorFieldBuffer = new ComputeBuffer(squaresCount, Marshal.SizeOf(typeof(VectorFieldSquare)));
        m_vectorFieldBuffer.SetData(squares);

    }

    private void InitializePositionNormalBuffer(int vertexCount) {

        m_positionBuffer = new ComputeBuffer(vertexCount, Marshal.SizeOf(typeof(Vector3)));
        m_normalBuffer = new ComputeBuffer(vertexCount, Marshal.SizeOf(typeof(Vector3)));

    }

    private void InitializeBoidsDataBuffer() {

        var boidsData = new BoidsData[m_instanceCount];

        for(int i = 0; i < m_instanceCount; ++i) {
            boidsData[i].position = Random.insideUnitSphere * m_boundaryRange;
            var velocity = new Vector3(Random.Range(-1f, 1f), Random.Range(-1f, 1f), Random.Range(-1f, 1f));
            boidsData[i].velocity = velocity.normalized * m_minVelocity; 
        }

        m_boidsDataBuffer = new ComputeBuffer(m_instanceCount, Marshal.SizeOf(typeof(BoidsData)));
        m_boidsDataBuffer.SetData(boidsData);


    }

    private void InitializeBoidsSimulator() {

        // 生成する個体の数を2の累乗にする(計算しやすくするため)
        m_instanceCount = Mathf.ClosestPowerOfTwo(m_instanceCount);

        InitializeBoidsDataBuffer();

        m_updateBoidsKernel = m_avoidanceBoidsSimulator.FindKernel("UpdateBoids");

        m_avoidanceBoidsSimulator.GetKernelThreadGroupSizes(m_updateBoidsKernel, out uint x, out uint y, out uint z);
        m_updateBoidsGroupSize = new Vector3Int(m_instanceCount / (int)x, (int)y, (int)z);

        m_avoidanceBoidsSimulator.SetFloat("_CohesionForce", m_cohesionForce);
        m_avoidanceBoidsSimulator.SetFloat("_SeparationForce", m_separationForce);
        m_avoidanceBoidsSimulator.SetFloat("_AlignmentForce", m_alignmentForce);

        // ComputeShader内の距離判定で2乗の値を使用しているので合わせる
        m_avoidanceBoidsSimulator.SetFloat("_CohesionDistance", m_cohesionDistance);
        m_avoidanceBoidsSimulator.SetFloat("_SeparationDistance", m_separationDistance);
        m_avoidanceBoidsSimulator.SetFloat("_AlignmentDistance", m_alignmentDistance);

        // ComputeShader内ではラジアンで判定するので度数方からラジアンに変更する
        m_avoidanceBoidsSimulator.SetFloat("_CohesionAngle", m_cohesionAngle * Mathf.Deg2Rad);
        m_avoidanceBoidsSimulator.SetFloat("_SeparationAngle", m_separationAngle * Mathf.Deg2Rad);
        m_avoidanceBoidsSimulator.SetFloat("_AlignmentAngle", m_alignmentAngle * Mathf.Deg2Rad);

        m_avoidanceBoidsSimulator.SetFloat("_BoundaryForce", m_boundaryForce);
        m_avoidanceBoidsSimulator.SetFloat("_BoundaryRange", m_boundaryRange * m_boundaryRange);
        m_avoidanceBoidsSimulator.SetVector("_BoundaryCenter", m_boundaryCenter);
        m_avoidanceBoidsSimulator.SetFloat("_AvoidForce", m_avoidForce);

        m_avoidanceBoidsSimulator.SetFloat("_MinVelocity", m_minVelocity);
        m_avoidanceBoidsSimulator.SetFloat("_MaxVelocity", m_maxVelocity);

        m_avoidanceBoidsSimulator.SetInt("_InstanceCount", m_instanceCount);

        m_avoidanceBoidsSimulator.SetFloat("_MaxForce", m_maxForce);
        
        m_avoidanceBoidsSimulator.SetBuffer(m_updateBoidsKernel, "_BoidsData", m_boidsDataBuffer);
        m_avoidanceBoidsSimulator.SetBuffer(m_updateBoidsKernel, "_VectorFieldBuffer", m_vectorFieldBuffer);

        m_avoidanceBoidsSimulator.SetFloat("_Epsilon", Mathf.Epsilon);
    }

    /// <summary>
    /// ベクトル場の更新
    /// </summary>
    private void UpdateVectorField() {
        // ベクトル場の初期化
        m_avoidanceBoidsSimulator.Dispatch(m_initializeVectorFieldKernel,
                                           m_initializeVectorFieldGroupSize.x,
                                           m_initializeVectorFieldGroupSize.y,
                                           m_initializeVectorFieldGroupSize.z);
        
        int bufferIndex = 0;
        for(int i = 0; i < m_skinnedMeshRenderer.Length; ++i) {
            m_meshPositions[i] = m_skinnedMeshTransform[i].position;
            m_meshEulerAngles[i] = m_skinnedMeshTransform[i].eulerAngles;
            m_meshScales[i] = m_skinnedMeshTransform[i].localScale;
            m_skinnedMeshRenderer[i].BakeMesh(m_bakeMesh);
            m_bakeMesh.GetVertices(m_vertices);
            m_bakeMesh.GetNormals(m_normals);
            int vertexCount = m_skinnedMeshRenderer[i].sharedMesh.vertexCount;
            m_positionBuffer.SetData(m_vertices, 0, bufferIndex, vertexCount);
            m_normalBuffer.SetData(m_normals, 0, bufferIndex, vertexCount);
            bufferIndex += vertexCount;
        }

        m_avoidanceBoidsSimulator.SetVectorArray(m_modelPositionPropID, m_meshPositions);
        m_avoidanceBoidsSimulator.SetVectorArray(m_modelEulerAnglePropID, m_meshEulerAngles);
        m_avoidanceBoidsSimulator.SetVectorArray(m_modelScalePropID, m_meshScales);

        // m_avoidanceBoidsSimulator.SetVector(m_modelPositionPropID, m_skinnedMeshTransform.position);
        // m_avoidanceBoidsSimulator.SetVector(m_modelEulerAnglePropID, m_skinnedMeshTransform.eulerAngles);
        // m_avoidanceBoidsSimulator.SetVector(m_modelScalePropID, m_skinnedMeshTransform.localScale);
        m_avoidanceBoidsSimulator.Dispatch(m_updateVectorFieldKernel,
                                           m_updateVectorFieldGroupSize.x,
                                           m_updateVectorFieldGroupSize.y,
                                           m_updateVectorFieldGroupSize.z);

    }

    /// <summary>
    /// Boidsアルゴリズムの更新
    /// </summary>
    private void UpdateBoids() {

        m_avoidanceBoidsSimulator.SetFloat(m_deltaTimeID, Time.deltaTime);
        m_avoidanceBoidsSimulator.Dispatch(m_updateBoidsKernel,
                                           m_updateBoidsGroupSize.x,
                                           m_updateBoidsGroupSize.y,
                                           m_updateBoidsGroupSize.z);

    }

    void OnDestroy() {

        m_vectorFieldBuffer?.Release();
        m_positionBuffer?.Release();
        m_normalBuffer?.Release();
        m_boidsDataBuffer?.Release();

    }

}


BoidsRenderer.cs(C#)

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

public class BoidsRenderer : MonoBehaviour
{
    [Header("DrawMeshInstancedIndirectのパラメータ")]
    [SerializeField]
    private Mesh m_mesh;

    [SerializeField]
    private Material m_material;

    [SerializeField]
    private Bounds m_bounds;

    [SerializeField]
    private ShadowCastingMode m_shadowCastingMode;

    [SerializeField]
    private bool m_receiveShadows;

    private AvoidanceBoids m_avoidanceBoids;

    private ComputeBuffer m_argsBuffer;


    public void Initialize(AvoidanceBoids avoidanceBoids) {
        m_avoidanceBoids = avoidanceBoids;
        InitializeArgsBuffer();

        m_material.SetBuffer("boidsDataBuffer", m_avoidanceBoids.BoidsDataBuffer);

    }

    void LateUpdate() {

        Graphics.DrawMeshInstancedIndirect(
            m_mesh,
            0,
            m_material,
            m_bounds,
            m_argsBuffer,
            0,
            null,
            m_shadowCastingMode,
            m_receiveShadows
        );
    }


    private void InitializeArgsBuffer() {

        var args = new uint[] { 0, 0, 0, 0, 0 };

        args[0] = m_mesh.GetIndexCount(0);
        args[1] = (uint)m_avoidanceBoids.BoidsInstanceCount;

        m_argsBuffer = new ComputeBuffer(1, 4 * args.Length, ComputeBufferType.IndirectArguments);

        m_argsBuffer.SetData(args);

    }


    private void OnDestroy() {

        m_argsBuffer?.Release();

    }

}


BoxRenderer.cs(C#)

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

public class BoxRenderer : MonoBehaviour
{
    [Header("DrawMeshInstancedInDirectの項目")]
    [Space(20)]
    [SerializeField]
    private Mesh _mesh;

    [SerializeField]
    private Material _material;

    [SerializeField]
    private Bounds _bounds;

    private ComputeBuffer _argsBuffer;

    private AvoidanceBoids m_avoidanceBoids;


    public void Initialize(AvoidanceBoids avoidanceBoids) {
        m_avoidanceBoids = avoidanceBoids;
        InitializeArgsBuffer();
        _material.SetFloat("_Scale", m_avoidanceBoids.SquareScale);
        _material.SetBuffer("_BoxDataBuffer", m_avoidanceBoids.VectorField);

    }

    void LateUpdate() {

            Graphics.DrawMeshInstancedIndirect(
            _mesh,
            0,
            _material,
            _bounds,
            _argsBuffer,
            0,
            null,
            ShadowCastingMode.Off,
            false
        );
    }

    private void InitializeArgsBuffer() {

        uint[] args = new uint[5] {0, 0, 0, 0, 0};

        args[0] = _mesh.GetIndexCount(0);
        int length = m_avoidanceBoids.SquareCount;
        args[1] = (uint)(length);

        _argsBuffer = new ComputeBuffer(1, args.Length * sizeof(uint), ComputeBufferType.IndirectArguments);
        _argsBuffer.SetData(args);
    }

    void OnDestroy() {

        _argsBuffer?.Release();

    }

}


Shader

Shader "Custom/BoidsShader"
{
    Properties
    {
        _Color ("Color", Color) = (1,1,1,1)
        _MainTex ("Albedo (RGB)", 2D) = "white" {}
        _Glossiness ("Smoothness", Range(0,1)) = 0.5
        _Metallic ("Metallic", Range(0,1)) = 0.0

        _ScaleX("ScaleX", float) = 1
        _ScaleY("ScaleY", float) = 1
        _ScaleZ("ScaleZ", float) = 1
    }
    SubShader
    {
        Tags { "RenderType"="Opaque" }
        LOD 200

        CGPROGRAM
        // Physically based Standard lighting model, and enable shadows on all light types
        #pragma surface surf Standard addshadow fullforwardshadows // 影を描画するためにはaddshadowが必要
        #pragma multi_compile_instancing    // GPU Instancingを可能にする
        #pragma instancing_options procedural:setup // setup関数を呼び出す


        // Use shader model 3.0 target, to get nicer looking lighting
        #pragma target 3.0

        sampler2D _MainTex;

        struct Input
        {
            float2 uv_MainTex;
        };

        struct BoidsData {
            float3 position;
            float3 velocity;
        };

        #ifdef UNITY_PROCEDURAL_INSTANCING_ENABLED
        StructuredBuffer<BoidsData> boidsDataBuffer;
        #endif

        float4x4 eulerAnglesToRotationMatrix(float3 angles) {

            float cx = cos(angles.x); float sx = sin(angles.x);
            float cy = cos(angles.z); float sy = sin(angles.z);
            float cz = cos(angles.y); float sz = sin(angles.y);

            return float4x4(
                cz*cy + sz*sx*sy, -cz*sy + sz*sx*cy, sz*cx, 0,
                sy*cx, cy*cx, -sx, 0,
                -sz*cy + cz*sx*sy, sy*sz + cz*sx*cy, cz*cx, 0,
                0, 0, 0, 1);

        }

        float4x4 CalcInverseMatrix(float3 position, float3 angle, float3 scale) {

            float4x4 inversScaleeMatrix = float4x4(
                1/scale.x, 0, 0, -position.x,
                0, 1/scale.y, 0, -position.y,
                0, 0, 1/scale.z, -position.z,
                0, 0, 0, 1);

            float4x4 mat = float4x4(
                1, 0, 0, 0,
                0, 1, 0, 0,
                0, 0, 1, 0,
                0, 0, 0, 1);

            float4x4 rotMatrix = mul(eulerAnglesToRotationMatrix(angle), mat);

            float4x4 inverseRotMatrix = float4x4(
                rotMatrix._11, rotMatrix._21, rotMatrix._31, 0,
                rotMatrix._12, rotMatrix._22, rotMatrix._32, 0,
                rotMatrix._13, rotMatrix._23, rotMatrix._33, 0,
                0, 0, 0, 1);

            return mul(inversScaleeMatrix, inverseRotMatrix);

        }

        fixed _ScaleX;
        fixed _ScaleY;
        fixed _ScaleZ;

        void setup() {

        #ifdef UNITY_PROCEDURAL_INSTANCING_ENABLED
            float3 position = boidsDataBuffer[unity_InstanceID].position;
            float3 velocity = boidsDataBuffer[unity_InstanceID].velocity;

            // スケーリング
            unity_ObjectToWorld._11_21_31_41 = float4(_ScaleX, 0, 0, 0);
            unity_ObjectToWorld._12_22_32_42 = float4(0, _ScaleY, 0, 0);
            unity_ObjectToWorld._13_23_33_43 = float4(0, 0, _ScaleZ, 0);

            // 速度から回転を求める
            float3 angle = float3(
                -asin(velocity.y/(length(velocity.xyz) + 1e-8)), // 0除算防止
                atan2(velocity.x, velocity.z),
                0);

            // 回転
            unity_ObjectToWorld = mul(eulerAnglesToRotationMatrix(angle), unity_ObjectToWorld);

            // 座標
            unity_ObjectToWorld._14_24_34_44 = float4(position, 1);

            // モデル行列を求める(間違っているかも. . .)
            // 参考:https://qiita.com/yuji_yasuhara/items/8d63455d1d277af4c270
            // 参考:http://gamemakerlearning.blog.fc2.com/blog-entry-196.html
            unity_WorldToObject = CalcInverseMatrix(position, angle, float3(_ScaleX, _ScaleY, _ScaleZ));


        #endif

        }


        half _Glossiness;
        half _Metallic;
        fixed4 _Color;

        // Add instancing support for this shader. You need to check 'Enable Instancing' on materials that use the shader.
        // See https://docs.unity3d.com/Manual/GPUInstancing.html for more information about instancing.
        // #pragma instancing_options assumeuniformscaling
        UNITY_INSTANCING_BUFFER_START(Props)
            // put more per-instance properties here
        UNITY_INSTANCING_BUFFER_END(Props)

        void surf (Input IN, inout SurfaceOutputStandard o)
        {
            // Albedo comes from a texture tinted by color
            fixed4 c = tex2D (_MainTex, IN.uv_MainTex) * _Color;
            o.Albedo = c.rgb;
            // Metallic and smoothness come from slider variables
            o.Metallic = _Metallic;
            o.Smoothness = _Glossiness;
            o.Alpha = c.a;
        }
        ENDCG
    }
    FallBack "Diffuse"
}



解説

複数のSkinnedMeshRendererに対応させる

使用するモデルによっては複数のSkinnedMeshRendererを使用しているので、そのようなモデルでも問題なく動作するようにしました。

ComputeBuffer.SetDataでは値を設定し始めるindexを指定することが出来るので、複数のSkinnedMeshRendererの頂点座標を1つのComputeBufferにいれて計算します。 法線ベクトルも同様に1つのComputeBufferに入れています。

        int bufferIndex = 0;
        for(int i = 0; i < m_skinnedMeshRenderer.Length; ++i) {
   . . . . . . .
   . . . . . . . 
            int vertexCount = m_skinnedMeshRenderer[i].sharedMesh.vertexCount;
            m_positionBuffer.SetData(m_vertices, 0, bufferIndex, vertexCount);
            m_normalBuffer.SetData(m_normals, 0, bufferIndex, vertexCount);
            bufferIndex += vertexCount;
        }


ただし、頂点座標や法線ベクトルを計算する際は各SkinnedMeshRendererのTransform情報が必要なのでそれらもComputeShaderに設定します。計算の際は頂点数の累積和を使用してどのTransform情報を使用するか判断しています。

    int transformIndex = 0;

    // 累積和とidからどのTransform情報を使用するか判断する
    for(int i = 0; (int)id < _VertexCountsThresholds[i]; ++i)
        transformIndex = i;



Boidsアルゴリズムにベクトル場の力を加える

Boidsの結合・分離・整列の計算をした後にベクトル場の力も加えます。ただし、個体がベクトル場の範囲内にいるかを確認してから計算します。

恐らく、ここら辺については計算方法を変えても良いと思います。

    float3 avoid = float3(0, 0, 0);
    // ベクトル場内にいる際は、ベクトルの値を取得してその方向に逃げるようにする
    if(_MinField.x <= posX.x && posX.x < _MaxField.x &&
       _MinField.y <= posX.y && posX.y < _MaxField.y &&
       _MinField.z <= posX.z && posX.z < _MaxField.z) {
        int3 index = (posX - _MinField) / _SquareScale;
        avoid = _VectorFieldBuffer[(index.x * _VectorFieldLength * _VectorFieldLength +
                        index.y * _VectorFieldLength +
                        index.z)].direction * _AvoidForce;
                        
       }

. . . . . . . .
. . . . . . . .
    velX += (cohesion + separation + alignment + boundary + avoid) * _DeltaTime;



レンダラーを分ける

DrawMeshInstancedIndirectで表示する機能も一つのクラスに入れると管理が難しくなると感じたので、Boidsを表示するクラスとベクトル場を表示するクラスを分けました。

注意すべき点として、ComputeBufferが初期化される前にMaterial.SetBufferをしてもマテリアルに設定されないのでBoidsアルゴリズムの初期化が終わった後にレンダラーの初期化を行います。

    void Start()
    {
   . . . . . .
   . . . . . . 
   . . . . . .
        // AvoidanceBoidsの初期化が完了したらRendererを初期化する
        m_boxRenderer?.Initialize(this);
        m_boidsRenderer?.Initialize(this);

    }



結果

VRIKとMapleちゃんを使ったデモです。


ベクトル場を使用した衝突回避Boidsアルゴリズム



まとめ

ある程度ちゃんと動きましたが、もうちょっといい動きが作れるんじゃないかなと感じました。

これからも改善の余地が見つかれば改良していこうと思います。



使用したもの

booth.pm

www.cgtrader.com

人型モデルの周りにベクトル場を生成する

f:id:vxd-naoshi-19961205-maro:20201102190840p:plain

始めに

Boidsアルゴリズムで衝突回避の実装について考えた結果、周りにベクトル場を生成すればいいのではと思いました。

そこで今回はSkinnedMeshRendererでアニメーションをするオブジェクトの周りにベクトル場を生成してみようと思います。


今回のサンプルリポジトリです。モデルデータやアニメーションは入っていないので、記事の最後に貼ってあるアセットなどを入れて試してみてください。

github.com



処理の流れ

モデルの法線を取得する

SkinnedMeshRenderer.BakeMeshを使用することでアニメーションしているモデルのメッシュ情報を取得できます。その情報をもとに法線方向を取得します。

f:id:vxd-naoshi-19961205-maro:20201102191513p:plain



空間を分割する

3次元空間をマス目に区切ります。

f:id:vxd-naoshi-19961205-maro:20201102220150p:plain



マス目に法線情報を書き込む

マス目に法線情報を書き込みます。書き込む際は加算して、使用する際はベクトルを正規化するようにします。

f:id:vxd-naoshi-19961205-maro:20201102220643p:plain


法線ベクトルを色として出力した結果が次の画像です。 法線が設定されていないマス目は描画していません。

f:id:vxd-naoshi-19961205-maro:20201102220758p:plain



実装

C#プログラム

using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System.Runtime.InteropServices;
using UnityEngine.Rendering;
public class VectorFieldCalculator : MonoBehaviour
{

    [SerializeField]
    private ComputeShader m_vectorGridCalculator;

    private Vector3Int m_calcVectorGridGroupSize;
    private Vector3Int m_initializeGridGroupSize;
    private int m_calcVectorGridKernel;
    private int m_initializeGridKernel;

    [SerializeField]
    private SkinnedMeshRenderer m_skinnedMeshRenderer;

    [Space(20)]
    [SerializeField]
    private float m_boxScale;

    [SerializeField]
    private int m_boxLength;

    [SerializeField]
    private Vector3 m_boxCenter;

    [Header("DrawMeshInstancedInDirectの項目")]
    [Space(20)]
    [SerializeField]
    private Mesh m_mesh;

    [SerializeField]
    private Material m_material;

    [SerializeField]
    private Bounds m_bounds;

    private ComputeBuffer m_argsBuffer;
    private ComputeBuffer m_boxDataBuffer;
    private ComputeBuffer m_normalsBuffer;
    private ComputeBuffer m_positionsBuffer;
    private List<Vector3> m_vertices = new List<Vector3>(); // GC.Collectを防ぐため
    private List<Vector3> m_normals = new List<Vector3>();  // GC.Collectを防ぐため
    private Mesh m_bakeMesh;
    private Transform m_modelTransform;
    private readonly int m_normalsPropID = Shader.PropertyToID("normals");
    private readonly int m_positionsPropID = Shader.PropertyToID("positions");
    private readonly int m_modelRotationPropID = Shader.PropertyToID("modelRotation");
    private readonly int m_modelEulerAnglePropID = Shader.PropertyToID("modelEulerAngle");
    private readonly int m_modelLocalScale = Shader.PropertyToID("modelLocalScale");
    private readonly int m_modelPositionPropID = Shader.PropertyToID("modelPosition");

    
    struct BoxData {
        public Vector3 position;
        public Vector3 direction;
    }


    // Start is called before the first frame update
    void Start()
    {
        InitializeArgsBuffer();
        InitializeVectorGridCalculator();
        
        m_material.SetFloat("_Scale", m_boxScale);

        m_modelTransform = m_skinnedMeshRenderer.transform;

        // GC.Allocを防ぐため
        m_bakeMesh = new Mesh();
    }

    // Update is called once per frame
    void Update()
    {

        UpdateVectorField();

    }

    void LateUpdate() {

        Graphics.DrawMeshInstancedIndirect(
            m_mesh,
            0,
            m_material,
            m_bounds,
            m_argsBuffer,
            0,
            null,
            ShadowCastingMode.Off,
            false
        );

    }

    public void UpdateVectorField() {

        // ベクトル場の初期化
        m_vectorGridCalculator.Dispatch(m_initializeGridKernel, 
                                       m_initializeGridGroupSize.x, 
                                       m_initializeGridGroupSize.y, 
                                       m_initializeGridGroupSize.z);
        
        m_skinnedMeshRenderer.BakeMesh(m_bakeMesh);
        m_bakeMesh.GetVertices(m_vertices);
        m_bakeMesh.GetNormals(m_normals);

        m_normalsBuffer.SetData(m_normals);
        m_positionsBuffer.SetData(m_vertices);

        m_vectorGridCalculator.SetVector(m_modelLocalScale, m_modelTransform.localScale);
        m_vectorGridCalculator.SetVector(m_modelEulerAnglePropID, m_modelTransform.eulerAngles);
        m_vectorGridCalculator.SetVector(m_modelPositionPropID, m_modelTransform.position);
        m_vectorGridCalculator.Dispatch(m_calcVectorGridKernel,
                                       m_calcVectorGridGroupSize.x,
                                       m_calcVectorGridGroupSize.y,
                                       m_calcVectorGridGroupSize.z);

    }

    private void InitializeArgsBuffer() {

        uint[] args = new uint[5] {0, 0, 0, 0, 0};

        args[0] = m_mesh.GetIndexCount(0);
        args[1] = (uint)(m_boxLength * m_boxLength * m_boxLength);

        m_argsBuffer = new ComputeBuffer(1, args.Length * sizeof(uint), ComputeBufferType.IndirectArguments);
        m_argsBuffer.SetData(args);
    }

    private void InitializeBoxDataBuffer() {

        int boxesCount = m_boxLength * m_boxLength * m_boxLength;
        BoxData[] boxDatas = new BoxData[boxesCount];
        float startPosX = m_boxCenter.x - m_boxLength * m_boxScale / 2.0f + m_boxScale/2.0f;
        float startPosY = m_boxCenter.y - m_boxLength * m_boxScale / 2.0f + m_boxScale/2.0f;
        float startPosZ = m_boxCenter.z - m_boxLength * m_boxScale / 2.0f + m_boxScale/2.0f;

        for(int i = 0; i < m_boxLength; ++i) {
            for(int j = 0; j < m_boxLength; ++j) {
                for(int k = 0; k < m_boxLength; ++k) {
                    int index = i * m_boxLength * m_boxLength + j * m_boxLength + k;
                    boxDatas[index].position = new Vector3(startPosX + i * m_boxScale, startPosY + j * m_boxScale, startPosZ + k * m_boxScale);
                    boxDatas[index].direction = Vector3.zero;
                }
            }
        }

        m_boxDataBuffer = new ComputeBuffer(boxesCount, Marshal.SizeOf(typeof(BoxData)));
        m_boxDataBuffer.SetData(boxDatas);
        m_material.SetBuffer("_BoxDataBuffer", m_boxDataBuffer);
    }

    private void InitializeVectorGridCalculator() {

        InitializeBoxDataBuffer();

        uint x, y, z;

        m_initializeGridKernel = m_vectorGridCalculator.FindKernel("InitializeGrid");
        m_vectorGridCalculator.GetKernelThreadGroupSizes(m_initializeGridKernel, out x, out y, out z);        
        m_initializeGridGroupSize = new Vector3Int(m_boxLength*m_boxLength*m_boxLength/(int)x, (int)y, (int)z);
        m_vectorGridCalculator.SetBuffer(m_initializeGridKernel, "boxData", m_boxDataBuffer);

        m_calcVectorGridKernel = m_vectorGridCalculator.FindKernel("CalcVectorGrid");

        m_vectorGridCalculator.GetKernelThreadGroupSizes(m_calcVectorGridKernel, out x, out y, out z);
        m_calcVectorGridGroupSize = new Vector3Int(m_skinnedMeshRenderer.sharedMesh.vertexCount / (int)x, (int)y, (int)z);

        float halfLength = m_boxScale * m_boxLength/2.0f;

        Vector3 minField = m_boxCenter - new Vector3(halfLength, halfLength, halfLength);
        Vector3 maxField = m_boxCenter + new Vector3(halfLength, halfLength, halfLength);

        m_vectorGridCalculator.SetVector("minField", minField);
        m_vectorGridCalculator.SetVector("maxField", maxField);

        m_vectorGridCalculator.SetFloat("boxScale", m_boxScale);
        m_vectorGridCalculator.SetInt("boxLength", m_boxLength); 

        m_vectorGridCalculator.SetBuffer(m_calcVectorGridKernel,
                                        "boxData",
                                        m_boxDataBuffer);

        m_normalsBuffer = new ComputeBuffer(m_skinnedMeshRenderer.sharedMesh.vertexCount, Marshal.SizeOf(typeof(Vector3)));
        m_positionsBuffer = new ComputeBuffer(m_skinnedMeshRenderer.sharedMesh.vertexCount, Marshal.SizeOf(typeof(Vector3)));
        m_vectorGridCalculator.SetBuffer(m_calcVectorGridKernel, m_normalsPropID, m_normalsBuffer);
        m_vectorGridCalculator.SetBuffer(m_calcVectorGridKernel, m_positionsPropID, m_positionsBuffer);
    }

    void OnDestroy() {
        m_argsBuffer?.Release();
        m_boxDataBuffer?.Release();
        m_normalsBuffer?.Release();
        m_positionsBuffer?.Release();
    }

}


ComputeShader

#pragma kernel CalcVectorGrid
#pragma kernel InitializeGrid

#include "Assets/cgincFiles/rotation.cginc"

struct BoxData {
    float3 position;
    float3 direction;
};

RWStructuredBuffer<BoxData> boxData;
StructuredBuffer<float3> normals;
StructuredBuffer<float3> positions;

float3 minField;
float3 maxField;

float boxScale;
int boxLength;

float4 modelRotation;
float3 modelPosition;
float3 modelLocalScale;
float3 modelEulerAngle;

[numthreads(8,1,1)]
void CalcVectorGrid (uint id : SV_DispatchThreadID)
{

    float3 pos = positions[id];
    float3 normal = normals[id];
    float4x4 rotMat = EulerAnglesToMatrix(modelEulerAngle); // 回転行列を求める

    pos = mul(rotMat, pos);
   
    pos.x *= modelLocalScale.x;
    pos.y *= modelLocalScale.y;
    pos.z *= modelLocalScale.z;

    pos += modelPosition;

    // 範囲外は計算しない
    if(pos.x < minField.x || maxField.x <= pos.x ||
       pos.y < minField.y || maxField.y <= pos.y ||
       pos.z < minField.z || maxField.z <= pos.z)
        return;

    int3 index = (pos - minField) / boxScale;

    normal = mul(rotMat, normal);

    boxData[(index.x * boxLength * boxLength + index.y * boxLength + index.z)].direction += normal;

}

[numthreads(8, 1, 1)]
void InitializeGrid(uint id : SV_DISPATCHTHREADID) {

    boxData[id].direction = float3(0, 0, 0);

}


Shader

Shader "Hidden/BoxShader"
{
    Properties
    {
        _Color("Color", Color) = (1, 1, 1, 1)
        _MainTex("Texture", 2D) = "white"
    }
    SubShader
    {
        Tags { "RenderType"="Transparent" "Queue" = "Transparent" }
        LOD 100
    //    Blend SrcAlpha OneMinusSrcAlpha
        Pass
        {
            CGPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            // make fog work
            #pragma multi_compile_fog

            #include "UnityCG.cginc"

            struct appdata
            {
                float4 vertex : POSITION;
                float2 uv : TEXCOORD0;
            };

            struct v2f
            {
                float2 uv : TEXCOORD;
                UNITY_FOG_COORDS(1)
                float4 vertex : SV_POSITION;
                float3 normal : TEXCOORD1;
            };

            struct BoxData {
                float3 position;
                float3 normal;
            };

            StructuredBuffer<BoxData> _BoxDataBuffer;
            float4 _Color;
            sampler2D _MainTex;
            half _Scale;

            v2f vert (appdata v, uint instancedId : SV_INSTANCEID)
            {
                v2f o;

               float3 worldPos = _BoxDataBuffer[instancedId].position;
                
                float4x4 scale = float4x4(
                    float4(_Scale, 0, 0, 0),
                    float4(0, _Scale, 0, 0),
                    float4(0, 0, _Scale, 0),
                    float4(0, 0, 0, 1)
                );

                v.vertex = mul(v.vertex, scale);
                o.vertex = UnityObjectToClipPos(v.vertex + worldPos);
                o.uv = v.uv;
                o.normal = _BoxDataBuffer[instancedId].normal;
                UNITY_TRANSFER_FOG(o,o.vertex);
                return o;
            }

            fixed4 frag (v2f i) : SV_Target
            {
                float4 col = tex2D(_MainTex, i.uv) * float4(i.normal, 1);
                if(length(i.normal))
                    col.a = 1;
                else {
                    discard;
                    // 全Cubeを出力する場合は以下の処理をする
                    // col = tex2D(_MainTex, i.uv) * float4(1, 1, 1, 1); 
                }
                return col;
            }
            ENDCG
        }
    }
}



解説

3次元配列を1次元配列にする

マス目自体は3次元ですが、1次元配列でデータを管理します。

一片のマス目の個数を Length とすると1次元配列の長さはLength3になります。

また、マス目の座標を(x, y, z)とするとこの3つの変数を使って1次元配列のインデックスは x * Length2 + y * Length + z と書き換えることが出来ます。 (n進数的な考え方がわかりやすいかもしれません)

参考 : 【JavaScript】二次元配列を一次元配列にする、とはどういうこと?要するにこういうこと。 - オドフラン ~いつもどこかに「なるほど」を~


追記:解説記事を書きました。 shitakami.hatenablog.com

基本的にはこの計算方法でマス目を求めていきます。



ベクトル場の更新

ベクトル場を更新する際は始めにベクトル場の初期化を行います。

初期化処理は単純にすべてのマス目に(0, 0, 0)を代入するだけです。ただ、量が多いのでComputeShaderで実行します。

        // ベクトル場の初期化
        m_vectorGridCalculator.Dispatch(m_initializeGridKernel, 
                                       m_initializeGridGroupSize.x, 
                                       m_initializeGridGroupSize.y, 
                                       m_initializeGridGroupSize.z);
[numthreads(8, 1, 1)]
void InitializeGrid(uint id : SV_DISPATCHTHREADID) {

    boxData[id].direction = float3(0, 0, 0);

}



次に、現在のモデルのMeshを取得します。SkinnedMeshRenderer.BakeMeshを使用することで、現在のMeshを取得できます。 (BakeMeshについては公式を参照してください: SkinnedMeshRenderer-BakeMesh - Unity スクリプトリファレンス

取得したMeshから頂点座標と法線ベクトルを取得して、ComputeShaderに設定します。

        m_skinnedMeshRenderer.BakeMesh(m_bakeMesh);
        m_bakeMesh.GetVertices(m_vertices);
        m_bakeMesh.GetNormals(m_normals);

        m_normalsBuffer.SetData(m_normals);
        m_positionsBuffer.SetData(m_vertices);



最後にSkinnedMeshRendererのTransform情報もComputeShaderに設定して、ComputeShaderを実行します。

 m_vectorGridCalculator.SetVector(m_modelLocalScale, m_modelTransform.localScale);
        m_vectorGridCalculator.SetVector(m_modelEulerAnglePropID, m_modelTransform.eulerAngles);
        m_vectorGridCalculator.SetVector(m_modelPositionPropID, m_modelTransform.position);
        m_vectorGridCalculator.Dispatch(m_calcVectorGridKernel,
                                       m_calcVectorGridGroupSize.x,
                                       m_calcVectorGridGroupSize.y,
                                       m_calcVectorGridGroupSize.z);



ComputeShader内での法線更新処理についても解説します。

BakeMeshで受け取った頂点座標はローカル座標なので、モデルの座標、回転、スケールを使ってワールド座標に変換します。

    float3 pos = positions[id];
    float3 normal = normals[id];
    float4x4 rotMat = EulerAnglesToMatrix(modelEulerAngle); // 回転行列を求める

    pos = mul(rotMat, pos);
   
    pos.x *= modelLocalScale.x;
    pos.y *= modelLocalScale.y;
    pos.z *= modelLocalScale.z;

    pos += modelPosition;



次に頂点座標がベクトル場内にあるかを調べます。

ベクトル場の範囲は ベクトル場の中心座標 ± ベクトル場の一片の長さ / 2 になります。C#プログラムでは次のように計算しています。

        float halfLength = m_boxScale * m_boxLength/2.0f;

        Vector3 minField = m_boxCenter - new Vector3(halfLength, halfLength, halfLength);
        Vector3 maxField = m_boxCenter + new Vector3(halfLength, halfLength, halfLength);



この範囲を使用して先ほどの頂点座標が範囲内にあるか確かめて、範囲内であれば法線ベクトルもワールド空間のベクトルに直してベクトル場に加算しています。

また、頂点座標からインデックスを取得する際は (頂点座標 - 境界値) / マスの大きさ で求まります。

 // 範囲外は計算しない
    if(pos.x < minField.x || maxField.x <= pos.x ||
       pos.y < minField.y || maxField.y <= pos.y ||
       pos.z < minField.z || maxField.z <= pos.z)
        return;

    int3 index = (pos - minField) / boxScale;

    normal = mul(rotMat, normal);

    boxData[(index.x * boxLength * boxLength + index.y * boxLength + index.z)].direction += normal;



GC.Colloectを防ぐ

ベクトル場を更新する際にMeshクラスとListクラスを使用するので、毎回newするのではなく変数として保持して使いまわすことでGC.Collectを防ぎます。

また、Mesh.verticesを使用するとその都度Vector3[]がnewされるみたいなのでMesh.GetVerticesを使って頂点座標を取得しています。

    private List<Vector3> m_vertices = new List<Vector3>(); // GC.Collectを防ぐため
    private List<Vector3> m_normals = new List<Vector3>();  // GC.Collectを防ぐため
. . . . . . . . . . 

    // Start is called before the first frame update
    void Start()
    {
        . . . . . . . . . .

        // GC.Allocを防ぐため
        m_bakeMesh = new Mesh();
    }
. . . . . . . . . . . .
. . . . . . . . . . . .

        m_skinnedMeshRenderer.BakeMesh(m_bakeMesh);
        m_bakeMesh.GetVertices(m_vertices);
        m_bakeMesh.GetNormals(m_normals);



ベクトル場の情報を色で出力する

vert関数で法線ベクトルをStructuredBufferから取り出してfrag関数で法線ベクトルをそのまま色で出力しています。もし、ベクトルの大きさが 0 である場合は出力しません。

fixed4 frag (v2f i) : SV_Target
            {
                float4 col = tex2D(_MainTex, i.uv) * float4(i.normal, 1);
                if(length(i.normal))
                    col.a = 1;
                else {
                    discard;
                    // 全Cubeを出力する場合は以下の処理をする
                    // col = tex2D(_MainTex, i.uv) * float4(1, 1, 1, 1); 
                }
                return col;
            }



追記

ベクトル場の拡大

ベクトル場の範囲が少し狭いと感じたのでベクトル場を広げようと思います。

ベクトル場の拡大は次のような計算を行いました。

f:id:vxd-naoshi-19961205-maro:20201109223736p:plain



プログラムではComputeShaderで次の計算を加えました。

. . . . . .
    boxData[(index.x * boxLength * boxLength + index.y * boxLength + index.z)].direction += normal;

    
    // ベクトルの要素の符号を取得、Epsilonは0除算回避
    int xd = normal.x / abs(normal.x + _Epsilon);
    int yd = normal.y / abs(normal.y + _Epsilon);
    int zd = normal.z / abs(normal.z + _Epsilon);

    if(0 <= index.x + xd && index.x + xd < boxLength)
        boxData[(index.x + xd) * boxLength * boxLength +
                           index.y * boxLength +
                           index.z].direction += normal * abs(normal.x);
    
    if(0 <= index.y + yd && index.y + yd < boxLength)
        boxData[index.x * boxLength * boxLength +
                           (index.y + yd) * boxLength +
                           index.z].direction += normal * abs(normal.y);
    
    if(0 <= index.z + zd && index.z + zd < boxLength)
        boxData[index.x * boxLength * boxLength +
                           index.y * boxLength +
                           index.z + zd].direction += normal * abs(normal.z);



結果は次のようにベクトル場の範囲が大きくなりました。

f:id:vxd-naoshi-19961205-maro:20201109230126p:plain

結果

以下のモデルのアニメーションからベクトル場を生成しました。

f:id:vxd-naoshi-19961205-maro:20201104113215g:plain

f:id:vxd-naoshi-19961205-maro:20201104112333g:plain

参考

法線ベクトルを出力する際にお借りしました。

qiita.com



使用したAsset

assetstore.unity.com

assetstore.unity.com

ML-Agentsの訓練設定ファイルを生成するEditor拡張の作成

始めに

ML-Agentsの訓練設定ファイルをUnityで編集できるエディタ拡張を作成致しました。

github.com



※まだ開発途中です!今後更新又は修正していく予定です

TrainingSettingFileEditor

f:id:vxd-naoshi-19961205-maro:20201011151856p:plain



機能

各パラメータを設定して.yamlファイルを生成することが出来ます。

開き方

Winodw/ML-Agents/TrainingSettingFileEditorからエディタを開くことが出来ます。 f:id:vxd-naoshi-19961205-maro:20201011222238p:plain



Curiosityやセルフプレイなどのオプション項目

Curiosityやセルフプレイなどを使って学習したい場合はチェックボックスをチェックすることで設定項目が出てきます。

f:id:vxd-naoshi-19961205-maro:20201011162912p:plain



保存

設定が完了したら一番下の "Generate TrainingSettingFile"を押すと保存先を指定するウィンドウが開きます。

そこでファイル名と保存先をしていすることで、訓練設定ファイルを生成することが出来ます。

ファイル名はデフォルトで "BehaviorName" + ppr or sacとなります。

f:id:vxd-naoshi-19961205-maro:20201011170604p:plain

f:id:vxd-naoshi-19961205-maro:20201011170728p:plain

例として生成されたファイルは次のようになります。

TestBehaivor_ppo.yaml

behaviors:
  TestBehavior:
    trainer_type: ppo
    max_steps: 0
    time_horizon: 0
    summary_freq: 0
    keep_checkpoints: 5
    checkpoint_interval: 500000
    threaded: true

    hyperparameters:
      batch_size: 0
      buffer_size: 0
      learning_rate: 0
      learning_rate_schedule: linear
      beta: 0
      epsilon: 0
      lambd: 0
      num_epoch: 0

    network_settings:
      vis_encoder_type: simple
      normalize: false
      hidden_units: 0
      num_layers: 0

    reward_signals:
      extrinsic:
        strength: 0
        gamma: 0

      curiosity:
        gamma: 0
        strength: 0
        encoding_size: 64
        learning_rate: 0.0003

    self_play:
      save_steps: 0
      team_change: 0
      swap_steps: 0
      window: 0
      play_against_latest_model_ratio: 0
      initial_elo: 0

これから追加予定の機能

tooltip

今のところすべて "tooltip" と出力されますが、これから各パラメータの説明文を追加する予定です。

f:id:vxd-naoshi-19961205-maro:20201011161812p:plain



複数のbehaviorの設定

公式のML-Agentsサンプルでのサッカーでは2つのbehaviorの設定をしています。 以下にその訓練設定ファイルを載せます。

このようなファイルの生成がまだ出来ないのでこれからその設定項目を作成する予定です。

StrikersVsGoalie.yaml

behaviors:
  Goalie:
    trainer_type: ppo
    hyperparameters:
      batch_size: 2048
      buffer_size: 20480
      learning_rate: 0.0003
      beta: 0.005
      epsilon: 0.2
      lambd: 0.95
      num_epoch: 3
      learning_rate_schedule: constant
    network_settings:
      normalize: false
      hidden_units: 512
      num_layers: 2
      vis_encode_type: simple
    reward_signals:
      extrinsic:
        gamma: 0.99
        strength: 1.0
    keep_checkpoints: 5
    max_steps: 50000000
    time_horizon: 1000
    summary_freq: 10000
    threaded: true
    self_play:
      save_steps: 50000
      team_change: 200000
      swap_steps: 1000
      window: 10
      play_against_latest_model_ratio: 0.5
      initial_elo: 1200.0
  Striker:
    trainer_type: ppo
    hyperparameters:
      batch_size: 2048
      buffer_size: 20480
      learning_rate: 0.0003
      beta: 0.005
      epsilon: 0.2
      lambd: 0.95
      num_epoch: 3
      learning_rate_schedule: constant
    network_settings:
      normalize: false
      hidden_units: 512
      num_layers: 2
      vis_encode_type: simple
    reward_signals:
      extrinsic:
        gamma: 0.99
        strength: 1.0
    keep_checkpoints: 5
    max_steps: 50000000
    time_horizon: 1000
    summary_freq: 10000
    threaded: true
    self_play:
      save_steps: 50000
      team_change: 200000
      swap_steps: 4000
      window: 10
      play_against_latest_model_ratio: 0.5
      initial_elo: 1200.0



カリキュラム学習、環境パラメータのランダム化

これらの設定も出来るようにしたかったのですが、ブログを上げるまでに間に合いませんでした。



参考

EditorWindowでもtooltipを表示させる方法を調べた際に参考にさせて頂きました。

answers.unity.com


ファイルを保存する際に保存先を指定する方法で参考にさせて頂きました。 baba-s.hatenablog.com


訓練設定ファイルの項目については以下の本を参考にさせて頂きました。