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

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

近傍探索の実装とまとめ

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