Boids集群模拟 - Sand and Foam - FlowingCrescent's Blog

前言

时光荏苒，不知不觉到了中秋国庆佳节，距离上一次的更新又过了一个季度。
其实期间有一次也有个小技术能够分享，不过姑且没时间写，之后有时间就再补上吧。

那么进入正题，这次在公司闲来无事的时候学习了下同事的boids算法，虽然代码没有照抄但是与其思路基本相同，也很感谢同事的分享，让我学习到了之前尚有些微了解但未曾直接动手尝试的boids算法。

本文需要读者拥有一定的ComputeShader相关知识，如果对ComputeShader不甚了解，可以先阅读：
Unity 3D : ComputeShader 全面詳解
 Compute Shaders - Rendering One Million Cubes

算法介绍

Boids算法是由Craig Reynolds于1986年提出的集群模拟算法，“boid"则是”"bird-oid object"的缩写，"oid"便是“某一类”的意思。
那么简单而言，这个算法对于一个集群中的单个单位的运动趋势进行了“受力分解”，将其复杂的运动趋势分解为Separation（分离力）、Cohesion（凝聚力）、Alignment（对齐力）这简单的三种力，并且只要该单位受这三种力的驱动，它们相互影响之下就可以实现如集群般的模拟效果。
B77(IPBUINRB7TMEIJ48E1.png

Separation（分离力）

分离力便是指每个boid单位都会对临近的单位进行排斥，这也就让它们保持了一定的距离而不至于全部撞在一起，通常而言是两个单位间距离越近则分离力越强。也可以将这个力用于躲避障碍物。

8W$PX8N~}J6VV6UC{B2N(A5.png

Alignment（对齐力）

对齐力指每个boid单位都会有一个与临近单位进行同样方向运动的趋势，大伙都往前走，总不能就你突然想往右或者左吧，这就是所谓大趋势了。

@4_}17))UKZN4LUYLPS5K$S.png

Cohesion（凝聚力）

凝聚力指每个单位都会有往相邻单位的中心点移动的趋势，跟分离力恰恰相反，非常有趣。

更具体的介绍可以参照wiki或者读这篇雷火的鱼群模拟分享。

简单CPU版本

毕竟CPU版并不是本文重点，如果要如以前的文章一样详实地写下每一步编写过程就有些冗繁了，这里我就直接放出CPU代码并带上一定的注释和解释。
CPU版代码很简单，我也没有细致地整理过代码，主打一个能跑就行。总之就是将每个boid都作为Prefab生成后再进行boid逻辑执行。
只需有一个Controller进行一些全局设置，然后各个单位有个类各自为战即可模拟出大略效果。

using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using Random = UnityEngine.Random;

public class BoidsController : MonoBehaviour
{
    [Range(1, 1000)]
    public int boidNumber = 100;//boid数量
    public Vector3 spawnRange = new Vector3(10, 10, 10);//出生时随机的范围大小
    public GameObject boidPrefab;//boid prefab
    public float targetRange = 50f;//目标点随机位置范围
    public float moveSpeed = 10f;//boid移动速度
    public float rotateSpeed = 1f;//boid转向速度
    public float neighborDistance = 3f;//临近距离，只有在这距离之内的boid才会被认为是临近单位
    public float separationStrength = 500f;//分离力强度，试下来发现这个参数应该要比较大
    public float alignmentStrength  = 5f;//对齐力强度
    public float cohesionStrength  = 5f;//凝聚力强度


    public List<Boid> allBoids;
    public Vector3 targetPos = Vector3.one;

    public static BoidsController _instance;
    private void Awake()
    {
        _instance = this;
        allBoids = new List<Boid>();
        InitBoids();
        InvokeRepeating("RandomTarget", 0f, 5);//每5秒会改变目标点
    }

    public void InitBoids()
    {
        for (int i = 0; i < boidNumber; i++)
        {
            float randomX = Random.Range(0, spawnRange.x);
            float randomY = Random.Range(0, spawnRange.y);
            float randomZ = Random.Range(0, spawnRange.z);
            var randomPosition = new Vector3(randomX, randomY, randomZ);
            var boid = Instantiate(boidPrefab, randomPosition, Quaternion.identity);

            boid.transform.SetParent(this.transform);
            allBoids.Add(boid.GetComponent<Boid>());

        }

    }

    public void RandomTarget()
    {
        float randomX = Random.Range(0, targetRange);
        float randomY = Random.Range(0, targetRange);
        float randomZ = Random.Range(0, targetRange);

        targetPos.x = randomX;
        targetPos.y = randomY;
        targetPos.z = randomZ;
    }


    private void OnDrawGizmos()
    {
        Gizmos.DrawWireCube(targetPos, Vector3.one);
    }
}

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

public class Boid : MonoBehaviour
{
    public int nearbyCount = 0;
    private Vector3 targetPos = Vector3.zero;
    private Vector3 force = Vector3.zero;
    private Vector3 velocity = Vector3.zero;
    private Vector3 separation = Vector3.zero;
    private Vector3 alignment = Vector3.zero;
    private Vector3 cohesion = Vector3.zero;
    
    private void Start()
    {
        velocity = transform.forward;
    }

    void Update()
    {
        GetTarget();//获取目标点
        GetForce();//受力
        UpdateVelocity();//更新方向
        MoveAndRotate();//进行位移
    }

    public void GetTarget()
    {
        targetPos = BoidsController._instance.targetPos;
    }

    public void GetForce()
    {
        nearbyCount = 0;
        separation = Vector3.zero;
        alignment = Vector3.zero;
        cohesion = Vector3.zero;
        
        force = targetPos - transform.position; // target force
        float neighborDistance = BoidsController._instance.neighborDistance;

        //遍历当前所有的boid，找到所有的临近boid然后进行受力
        for (int i = 0; i < BoidsController._instance.allBoids.Count; i++)
        {
            var otherBoid = BoidsController._instance.allBoids[i];
            if(otherBoid == this)
                continue;

            var distanceVector = otherBoid.transform.position - transform.position;
            float distance = distanceVector.magnitude;
            
            if (distance < neighborDistance)
            {
                //排斥力
                separation += -distanceVector.normalized * (neighborDistance - distance) / neighborDistance;
                alignment += otherBoid.transform.forward;//对齐力
                cohesion += otherBoid.transform.position;//凝聚力
                
                nearbyCount++;
            }
        }

        float avg = 1.0f / Mathf.Max(nearbyCount, 1);
        separation *= avg * BoidsController._instance.separationStrength;
        alignment *= avg * BoidsController._instance.alignmentStrength;
        cohesion *= avg;

        cohesion = (cohesion - transform.position) * BoidsController._instance.cohesionStrength;;


        force += separation + alignment + cohesion;

    }

    public void UpdateVelocity()
    {
        velocity = transform.forward;;
        velocity += force;
        velocity = Vector3.Normalize(velocity);
    }

    public void MoveAndRotate()
    {
        var rotateSpeed = BoidsController._instance.rotateSpeed;
        var speed = BoidsController._instance.moveSpeed;
        
        transform.rotation = Quaternion.Lerp(transform.rotation, Quaternion.LookRotation(velocity), rotateSpeed * Time.deltaTime);
        transform.Translate(transform.forward * Time.deltaTime * speed, Space.World);
    }
    
}

于是在场景中设置一个BoidController后对其进行一些设置，再运行工程即可。

此时的开销主要在boid之间相互遍历的O(n2)算法以及大量的Boid MonoBehavior物体上，渲染方面性能大量盈余，于是我给Prefab加了TrailRenderer玩了下（虽说生成Trail网格还是CPU的活233），顺便观察每条鱼的运行轨迹。
这种直接使用GameObject，让每个Boid遍历其他所有Boid的方案，性能显然是最差的，这里我的3060+i7 12700H，跑个400单位就已经只能25fps运行，1000单位则直接卡到个位数fps。

如果要绘制大规模的boids集群，那么很显然每个boid单位不可能使用MonoBehavior，直接使用Instance进行绘制即可。
那么关键就在于“优化boid获取临近单位的算法”了，通常而言思路就是进行空间划分：

对空间进行八叉树划分，使每个单位方便查找到临近格子中包含的单位，这种方法需要持续更新每个格子中储存的boid信息。
或者是对空间进行等比例划分后根据划分好的空间格子ID对其中的boid进行排序，这样一来“查找临近的单位”就变成在排序后的数组中直接“访问特定区间里的单位”，因为临近的单位在根据Grid ID排序后必然连续。

这次我尝试的是第二种方法，因为同事教了我以前不太熟悉的Counting Sort（笑）

简单GPU版本

那么GPU版本的方案就是在ComputeShader中进行Boids的模拟，然后调用Graphics.DrawMeshInstancedProcedural进行Instance绘制，在shader中直接取用对应index的位置和方向信息即可。

GPU版的简易版本大致如下，ComputeShader中：

[numthreads(512,1,1)]
void UpdateMoving (uint3 id : SV_DispatchThreadID)
{
    float3 pos = _PosBuffer[id.x];
    float3 dir = _RotationBuffer[id.x];
    int nearbyCount = 0;
    float3 separation = 0;
    float3 alignment = 0;
    float3 cohesion = 0;
    
    float3 force = _TargetPos.xyz - pos.xyz; // target force

    for (int i = 0; i < _Count; i++) //常规for循环查询
    {
        if(i == id.x)
            continue;
        
        float4 otherBoid = _PosBuffer[i];
        float3 distanceVector = otherBoid.xyz - pos.xyz;
        float distance = length(distanceVector);
            
        if (distance < _NeighborDistance)
        {
            separation += -normalize(distanceVector) * (_NeighborDistance - distance) / _NeighborDistance;
            alignment += _RotationBuffer[i].xyz;
            cohesion += otherBoid.xyz;
                
            nearbyCount++;
        }
    }

    float avg = 1.0 / max(nearbyCount, 1);
    separation *= avg * _SeparationStrength;
    alignment *= avg * _AlignmentStrength;
    cohesion *= avg;
    
    cohesion = (cohesion - pos) * _CohesionStrength;;
    
    force *= _TargerStrength;
    force += separation + alignment + cohesion;
    force = normalize(force);
    float3 newDir = normalize(lerp(dir, force, _RotateSpeed * _TimeDelta));
    newDir.y = max(-0.5, min(0.5, newDir.y));
    newDir = normalize(newDir);
    
    _RotationBuffer[id.x].xyz = newDir;
    _PosBuffer[id.x].xyz = pos + newDir * _MoveSpeed * _TimeDelta;
}

在BoidController中进行ComputeShader以及DrawInstance的调用

   void InitBuffer()
   {
       boidPosBuffer = new ComputeBuffer(count, sizeof(float)*4); //float4
       boidRotationBuffer = new ComputeBuffer(count, sizeof(float)*4); //float4 
   }


void UpdateMoving()
   {
       int kernel = computeShader.FindKernel("UpdateMoving");
       computeShader.SetBuffer(kernel, _PosBuffer, boidPosBuffer);
       computeShader.SetBuffer(kernel, _RotationBuffer, boidRotationBuffer);
       computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);
       computeShader.SetInt(_Count, count);
       computeShader.SetFloat(_MoveSpeed, moveSpeed);
       computeShader.SetFloat(_RotateSpeed, rotateSpeed);
       computeShader.SetFloat(_NeighborDistance, neighborDistance);
       computeShader.SetFloat(_TargerStrength, targerStrength);
       computeShader.SetFloat(_SeparationStrength, separationStrength);
       computeShader.SetFloat(_AlignmentStrength, alignmentStrength);
       computeShader.SetFloat(_CohesionStrength, cohesionStrength);
       computeShader.SetVector(_TargetPos, targetPos);
       computeShader.SetFloat(_TimeDelta, Time.deltaTime);
       
       computeShader.SetBuffer(kernel,_NeighborCountBuffer, neighborCountBuffer);
       
       computeShader.Dispatch(kernel, count/512, 1, 1);

   }


   void Update()
   {
       if (boidPosBuffer.count != count)
           InitBuffer();

       UpdateMoving();

       material.SetBuffer(_PosBuffer, boidPosBuffer);
       material.SetBuffer(_RotationBuffer, boidRotationBuffer);
       Graphics.DrawMeshInstancedProcedural(mesh, 0, material, maxBounds, count, prop, castShadow, false, gameObject.layer);
       
   }

Shader中直接取用Buffer即可。

Shader "Custom/Boids" {
    Properties {
        _BaseMap ("Example Texture", 2D) = "white" {}
        _BaseColor ("Example Colour", Color) = (0, 0.66, 0.73, 1)
    }
    SubShader {
        Tags { "RenderType"="Opaque" "RenderPipeline"="UniversalRenderPipeline" }
        
        HLSLINCLUDE
        #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Core.hlsl"
        #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Lighting.hlsl"

        
        CBUFFER_START(UnityPerMaterial)
        float4 _BaseMap_ST;
        float4 _BaseColor;
        CBUFFER_END
        ENDHLSL
        
        Pass {
            Tags { "LightMode"="UniversalForward" }
            
            Cull back
            Zwrite on
            Ztest LEqual
            
            HLSLPROGRAM
            #pragma vertex vert
            #pragma fragment frag
            #pragma multi_compile_fog
            
            struct a2v {
                float4 positionOS   : POSITION;
                float2 uv           : TEXCOORD0;
                float3 normal       : NORMAL;
                float4 color        : COLOR;
            };
            
            struct v2f {
                float4 positionCS  : SV_POSITION;
                float2 uv           : TEXCOORD0;
                float4 color        : TEXCOORD1;
                float3 normalWS     : TEXCOORD2;
                float3 positionWS     : TEXCOORD3;
            };
            
            TEXTURE2D(_BaseMap);
            SAMPLER(sampler_BaseMap);

            StructuredBuffer<float4> _PosBuffer;
            StructuredBuffer<float4> _RotationBuffer;
            
            v2f vert(a2v v, uint instID : SV_InstanceID) {
                v2f o;
                o.color = v.positionOS.y > 0;

                float3 forward = normalize(_RotationBuffer[instID].xyz);
                float3 right = normalize(cross(forward, float3(0, 1, 0)));
                float3 up = normalize(cross(right, forward));
                right = normalize(cross(forward, up));
                
                v.positionOS.xyz = mul(v.positionOS.xyz, float3x3(right, up, forward));
                
                float3 centerWS = _PosBuffer[instID].xyz;
                float3 posWS = v.positionOS + centerWS;
                o.positionCS = TransformWorldToHClip(posWS);

                o.positionWS = posWS;
                o.uv = TRANSFORM_TEX(v.uv, _BaseMap);
                o.normalWS = mul(v.normal, float3x3(right, up, forward));

                return o;
            }
            
            half4 frag(v2f i) : SV_Target {
                half4 baseMap = SAMPLE_TEXTURE2D(_BaseMap, sampler_BaseMap, i.uv);

                half3 normalWS = normalize(i.normalWS);
                InputData inputData = (InputData)0;
                inputData.normalWS = normalWS;
                inputData.positionWS = i.positionWS;
                inputData.viewDirectionWS = normalize(GetWorldSpaceViewDir(i.positionWS));
                // inputData.fogCoord = ComputeFogFactor(i.positionCS.z);
                inputData.bakedGI = SampleSH(inputData.normalWS);

                SurfaceData surfaceData = (SurfaceData)0;
                surfaceData.albedo = baseMap * _BaseColor;
                surfaceData.alpha = 1;
                surfaceData.metallic = 0;
                surfaceData.smoothness = 0;
                surfaceData.occlusion = 1;
                half4 color = UniversalFragmentPBR(inputData, surfaceData);

                color.rgb = MixFog(color.rgb, ComputeFogFactor(i.positionCS.z));

                // color.rgb = i.color;
                return color;
            }
            ENDHLSL
        }
        
    }
}

此时已经可以流畅绘制2w个boid
以上只是把上文的CPU版本照搬到GPU，并没有进行算法层面的优化，但即便如此，这样的方案已经能流畅地画上万个boid了，毕竟使用了GPU多线程，也节省了MonoBehavior的开销。但如果达到4万左右的数量，便会开始产生卡顿。

优化GPU版本

那么进一步优化的关键就在于ComputeShader中模拟Boids的部分，正如之前所说，重点在于对Boids根据Grid的顺序进行排序，算法思路来自于Fast Fixed-Radius Nearest Neighbor Search on the GPU

Counting Sort

介绍

NIVIDA的这张PDF展示了Counting Sort的全过程，在此干脆就对这张图进行翻译：

我们拥有一组Input数据，是每个boid单位的index和坐标以及它们对应的Grid ID（通过坐标计算可得）
计算每个Grid中有多少boid单位，以及每个boid是该Grid中的“第几个”。比如图上的5号单位是50号grid中的第1个，而8号单位是50号grid中的第2个。
对每个Grid中的boid数量进行前缀和（Prefix Sum）。“前缀和”指这个数组元素之前（也有将当前元素包含进去的算法）的所有元素值之和，比如上图中第50号grid之前没有元素，那么它的前缀和为0，但51号grid的前缀和值由于50号grid中boid数量为2，则51号grid的前缀值就为2，以此类推，60号grid的前缀和为3。
结合第一步计算的Grid ID进行第三步中计算的Grid前缀和查询（记为sum），然后查询第二步计算的“它是每个Grid中的第几个”（记为offset），我们便可以直接计算出每个boid在整个数组中排序后的位置（sum + offset）。

整个算法开销为O(n+k)，n为boid单位数而k为grid数量。

那么接下来的主要内容就是如何在Compute Shader中实现Counting Sort。

计算Grid ID并统计grid中boid数据

由于我们每帧都要从0开始重新统计每个Grid中的boid数据，因此每帧我们都要先Clear一下Buffer，
C#中：

void ClearBuffer()
{
    int totalGridCount = gridDivde * gridDivde * gridDivde;
    int kernel = computeShader.FindKernel("ClearBuffer");
    //int2[totalGridCount] buffer, x存储当前grid中boid数量，y存储该grid之前的boid数
    computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);

    computeShader.Dispatch(kernel, totalGridCount/64, 1, 1);
}

void UpdateGrid()
{
    int kernel = computeShader.FindKernel("UpdateGrid");
    //int2[boidCount] buffer, x存储boid所属grid，y存储boid在所属grid中的index
    computeShader.SetBuffer(kernel, _GridIdIndexBuffer, gridIDBuffer);
    computeShader.SetBuffer(kernel, _PosBuffer, boidPosBuffer);
    computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);

    computeShader.SetInt(_GridDivide, gridDivde);
    computeShader.SetVector(_BoundCenter, transform.position);
    computeShader.SetFloat(_BoundSize, boundSize);
    computeShader.Dispatch(kernel, count/64, 1, 1);
}

ComputeShader中：

[numthreads(64,1,1)]
void ClearBuffer (uint3 id : SV_DispatchThreadID)
{
    _GridCountBuffer[id.x].xy = 0;
}

int getGridID(float3 pos)
{
    float3 relativePos = pos - (_BoundCenter.xyz - _BoundSize.xxx * 0.5);
    //将计算后的相对坐标限制在Bound之内，用于保证ID正确
    relativePos = min(max(0.0001, relativePos), _BoundSize.xxx - 0.0001);
    float gridSize = _BoundSize / _GridDivide;
    int gridID = floor(relativePos.x / gridSize) + floor(relativePos.z / gridSize) * _GridDivide + floor(relativePos.y / gridSize) * _GridDivide * _GridDivide;

    return gridID;
}

[numthreads(64,1,1)]
void UpdateGrid (uint3 id : SV_DispatchThreadID)
{
    int gridID = getGridID(_PosBuffer[id.x]);
    _GridIdIndexBuffer[id.x].x = gridID;//计算每个boid所属的grid id

    //https://learn.microsoft.com/en-us/windows/win32/direct3dhlsl/interlockedadd
    InterlockedAdd(_GridCountBuffer[gridID].x, 1, _GridIdIndexBuffer[id.x].y);//计算每个grid中有多少boid，以及boid在该grid中的index
}

其中我们用到了InterlockedAdd，它是HLSL提供的函数，能够使多个线程同时访问共享变量时，以原子方式执行加法操作，确保线程之间不会发生竞态条件或数据不一致的情况。
而InterlockedAdd(dest, value, original_value)这一写法能够时value加到dest上的同时，将dest加之前的原值copy给original_value。
因此在这里的效果则是每个boid都会给对应的_GridCountBuffer[gridID].x + 1，然后加之前的值写入_GridIdIndexBuffer[id.x].y（正是boid在该grid中的index值）。可以说相当之妙。

对grid进行前缀和

由于我们需要排序的数组长度过长（我这里使用646464的grid划分，因此grid数高达262144），需要对这个长数组进行切分后各个线程组内部对各个片段分别进行求前缀和。
而使用GPU（或者说并行计算）求前缀和同样也有比较高超的算法，依然主要参考NIVIDA的文档：https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda
说来非常简单，记循环次数为c（从0开始），每次循环让array[i] += array[i - 2c]，若i - 2c小于0则保持原值不变。说来是感觉简单得令人难以理解，但是如果看以下例图则可豁然开朗：
)}D)EJ19DHL9D5_CN4QXDC7.png

//C#中：
void UpdateGridPrefixSum()
{
    int totalGridCount = gridDivde * gridDivde * gridDivde;

    int kernel = computeShader.FindKernel("UpdatePrefixSum");
    computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);
    computeShader.SetBuffer(kernel, _PrefixSumEndBuffer, prefixSumEndBuffer);
    computeShader.SetInt(_GridDivide, gridDivde);
    computeShader.Dispatch(kernel, totalGridCount/64, 1, 1); //各线程组内求前缀和

    //...
}


//Compute Shader中：
RWStructuredBuffer<int> _PrefixSumEndBuffer;
groupshared int groupCache[512];
[numthreads(512,1,1)]
void UpdatePrefixSum(uint3 id : SV_DispatchThreadID,//在所有线程中的id
                     uint3 group_tid : SV_GroupThreadID,//在该线程组中的id
                     uint3 group_id :  SV_GroupID)//该线程组的id
{
    groupCache[group_tid.x] = _GridCountBuffer[id.x].x;

    GroupMemoryBarrierWithGroupSync();
    for(int offset = 1; offset < 512; offset *= 2)
    {
        int offsetedIndex = group_tid.x - offset;
        if(offsetedIndex > 0)
            groupCache[group_tid.x] += groupCache[offsetedIndex];

        GroupMemoryBarrierWithGroupSync();
    }

    _GridCountBuffer[id.x].y = group_tid.x == 0 ? 0 : groupCache[group_tid.x];
    
    if(group_tid.x == 511)//储存最后一个前缀和的值
        _PrefixSumEndBuffer[group_id.x] = groupCache[group_tid.x];
    
}

在这里我们用到了512个（646464/512）线程组，此时一个线程组中处理512个单位的前缀和。
此处我们用到了GroupMemoryBarrierWithGroupSync()，该指令能够强制线程组中的线程在执行到这一命令时先暂停（或称阻塞，即Block），直到该线程组中的所有线程都达到了该指令，简单来说就是一个线程同步指令。具体可以参考:https://learn.microsoft.com/en-us/windows/win32/direct3dhlsl/groupmemorybarrierwithgroupsync

那么来调试一下这个算法是否正确，在C#中：

int[] testArray = new int[64*64*64*2];//因为buffer为int2类型所以要乘2
int totalCount = 0;
for (int i = 0; i < testArray.Length; i++)
{
    testArray[i] = (i + 1) % 2; // 1, 0, 1, 0, 1, 0,....  每个x设为1，y为0
}
gridCountBuffer.SetData(testArray);

int kernel = computeShader.FindKernel("UpdatePrefixSum");
computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);
computeShader.SetBuffer(kernel, _PrefixSumEndBuffer, prefixSumEndBuffer);
computeShader.SetInt(_GridDivide, gridDivde);
computeShader.Dispatch(kernel, totalGridCount/512, 1, 1); //各线程组内求前缀和

int[] testArray2 = new int[Mathf.RoundToInt(testArray.Length/2.0f/512)];
prefixSumEndBuffer.GetData(testArray2);
for (int i = 0; i < testArray2.Length; i++)
{
    Debug.Log("output:" + i + "  =  "  + testArray2[i]);//输出所有线程组中前缀和最后的值
}

理想结果而言它的log全都应该是511，但是现实很骨感：

很遗憾，还是得遵从NIVIDA的文档使用两个buffer进行来回copy，否则以上单buffer的写法很可能导致在这次循环中线程a正想读array[i]的值，但这个array[i]却被线程b刚刚修改过导致结果错误。

因此在Compute Shader中要改成这样

//双buffer版
RWStructuredBuffer<int> _PrefixSumEndBuffer;
groupshared int groupCache[512];
groupshared int groupCache2[512];
[numthreads(512,1,1)]
void UpdatePrefixSum(uint3 id : SV_DispatchThreadID,//在所有线程中的id
                     uint3 group_tid : SV_GroupThreadID,//在该线程组中的id
                     uint3 group_id :  SV_GroupID)//该线程组的id
{
    groupCache[group_tid.x] = _GridCountBuffer[id.x].x;
    groupCache2[group_tid.x] = _GridCountBuffer[id.x].x;
    bool flag = false;
    GroupMemoryBarrierWithGroupSync();
    for(int offset = 1; offset < 512; offset *= 2)
    {
        int offsetedIndex = group_tid.x - offset;

        if(flag)
        {
            if(offsetedIndex > 0)
                groupCache[group_tid.x] = groupCache2[group_tid.x] + groupCache2[offsetedIndex];
            else
                groupCache[group_tid.x] = groupCache2[group_tid.x];

        }
        else
        {
            if(offsetedIndex > 0)
                groupCache2[group_tid.x] = groupCache[group_tid.x] + groupCache[offsetedIndex];
            else
                groupCache2[group_tid.x] = groupCache[group_tid.x];
        }

        flag = !flag;
        GroupMemoryBarrierWithGroupSync();
    }

    _GridCountBuffer[id.x].y = group_tid.x == 0 ? 0 : groupCache2[group_tid.x];
    
    if(group_tid.x == 511)
        _PrefixSumEndBuffer[group_id.x] = groupCache2[group_tid.x];
    
}

改成这样的写法后再进行调试，每次结果都能是完美的511：

那么接下来需要将被划分的前缀和数组们进行合并，算法也同样很简单，如果仔细看上面的代码可以发现我在线程组中的最后把每个前缀和数组中的最后一个值保存到了另一个buffer中，正是要它用来连接各个前缀和数组。

比如以下是4个刚被划分出来尚未进行前缀和计算的数组：

	[0]	[1]	[2]	[3]	end
array1	5	2	1	3	/
array2	3	1	4	6	/
array3	6	5	7	2	/
array4	2	1	4	3	/

那么我们对它们分别进行前缀和后结果如下：

	[1]	[2]	[3]	end
array1	5	7	8	11
array2	3	4	8	14
array3	6	11	18	20
array4	2	3	7	10

要让array2的前缀和结果与array1的接上非常简单，只需要让array2的所有值都加上array1的"end"值即可，而array3想要与array2和array1的接上，则要让它的所有值加上array1和array2的end……
是不是很熟悉？没错，又是求前缀和。
我们对{11, 14, 20, 10}求前缀和，可以获得end[]{0, 11, 25, 45}，只要让每个array中的所有值分别加上其array的index对应的end[index]即可。

那么我们继续实现，C#中：

kernel = computeShader.FindKernel("UpdateEndPrefixSum");
computeShader.SetInt(_GridDivide, gridDivde);
computeShader.SetBuffer(kernel, _PrefixSumEndBuffer, prefixSumEndBuffer);
computeShader.Dispatch(kernel, 1, 1, 1); //各线程组最后的值求前缀和

kernel = computeShader.FindKernel("UpdateFullPrefixSum");
computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);
computeShader.SetBuffer(kernel, _PrefixSumEndBuffer, prefixSumEndBuffer);
computeShader.Dispatch(kernel, totalGridCount/64, 1, 1); //获得完整的前缀和列表

ComputeShader中：

[numthreads(512,1,1)]
void UpdateEndPrefixSum(uint3 id : SV_DispatchThreadID)
{
    groupCache[id.x] = _PrefixSumEndBuffer[id.x].x;
    groupCache2[id.x] = _PrefixSumEndBuffer[id.x].x;

    bool flag = false;
    GroupMemoryBarrierWithGroupSync();
    for(int offset = 1; offset < 512; offset *= 2)
    {
        int offsetedIndex = id.x - offset;

        if(flag)
        {
            if(offsetedIndex > 0)
                groupCache[id.x] = groupCache2[id.x] + groupCache2[offsetedIndex];
            else
                groupCache[id.x] = groupCache2[id.x];

        }
        else
        {
            if(offsetedIndex > 0)
                groupCache2[id.x] = groupCache[id.x] + groupCache[offsetedIndex];
            else
                groupCache2[id.x] = groupCache[id.x];

        }

        flag = !flag;
        GroupMemoryBarrierWithGroupSync();
    }

    _PrefixSumEndBuffer[id.x].x = id.x == 0 ? 0 : groupCache2[id.x];
    
}

[numthreads(512,1,1)]
void UpdateFullPrefixSum(uint3 id : SV_DispatchThreadID)
{
    int offsetIndex = id.x / 512;
	_GridCountBuffer[id.x].y += _PrefixSumEndBuffer[offsetIndex].x;
}

那么让我们验证一下，我这里输出了合并后的最后100个值

没有问题，最后也的确是262143（512*512-1），少了1是因为不包含自己（不包含当前grid内的boid）。

对boid根据grid进行排序

到了这一步就相当简单了，C#中：

//C#
kernel = computeShader.FindKernel("SortBoidIndexByGrid");
computeShader.SetBuffer(kernel, _BoidIndexSortByGridBuffer, boidIndexSortByGridBuffer);
computeShader.SetBuffer(kernel, _GridCountBuffer, gridCountBuffer);
computeShader.SetBuffer(kernel, _GridIdIndexBuffer, gridIDBuffer);
computeShader.Dispatch(kernel, count/512, 1, 1); //将Boid Index依照Grid ID排序

//Compute Shader
[numthreads(512,1,1)]
void SortBoidIndexByGrid(uint3 id : SV_DispatchThreadID)
{
    int indexOffset = _GridIdIndexBuffer[id.x].y;//boid所属grid中的index
    int gridID = _GridIdIndexBuffer[id.x].x;//boid所属的grid id
    int gridPrefixSum = _GridCountBuffer[gridID].y;//费了大功夫算的grid前缀和

    int index = gridPrefixSum + indexOffset;//在所有boid中的排序
    _BoidIndexSortByGridBuffer[index] = id.x;
}

至此终于是把Counting Sort实现完毕。

查询相邻格子中的Boid

那么终于到了检验的时刻，这里我直接用三层循环来遍历空间中共27（333)个grid，查询它们包含的boid。

int2 boidGridID = _GridIdIndexBuffer[id.x].xy; //x存储boid所属grid，y存储boid在grid中的index
//还原出原本的grid坐标
int gridX = boidGridID % (_GridDivide * _GridDivide) % _GridDivide;
int gridY = boidGridID % (_GridDivide * _GridDivide) / _GridDivide;
int gridZ = boidGridID / (_GridDivide * _GridDivide);

int indexMin;
int indexMax;
for(int x = -1; x <= 1; x++)
{
    for(int y = -1; y <= 1; y++)
    {
        for(int z = -1; z <= 1; z++)
        {
            int searchX = gridX + x;
            int searchY = gridY + y;
            int searchZ = gridZ + z;

            // grid合法
            if(searchX >= 0 && searchX < _GridDivide && searchY >= 0 && searchY < _GridDivide && searchZ >= 0 && searchZ < _GridDivide)
            {
                int searchGridID = searchX + searchZ * _GridDivide + searchY * _GridDivide * _GridDivide;

                indexMin = _GridCountBuffer[searchGridID].y;
                indexMax = indexMin + _GridCountBuffer[searchGridID].x;
                for(int i = indexMin; i < indexMax; i++)
                {
                    int otherBoidId = _BoidIndexSortByGridBuffer[i];

                    if(id.x == otherBoidId)
                        continue;

                    float4 otherBoid = _PosBuffer[otherBoidId];
                    float3 distanceVector = otherBoid.xyz - pos.xyz;
                    float distance = length(distanceVector);

                    if (distance < _NeighborDistance)
                    {
                        separation += -normalize(distanceVector) * (_NeighborDistance - distance) / _NeighborDistance;
                        alignment += _RotationBuffer[otherBoidId].xyz;
                        cohesion += otherBoid.xyz;

                        nearbyCount++;
                    }
                }
            }
        }
    }
}

大功告成，此时运行，能够流畅绘制20w boids，共约3240w vertices

最后随手录了个视频，一共102400只boid，在空间中每5秒会随机布16个点（橙色球），它们将会以离自己最近的点为目标点进行移动：

不过似乎已经被B站压得不知道变成什么东西了（笑）

结语

这次的核心在于实现大规模boid绘制，如果要真正地投入项目中使用，还应当实现鱼群/鸟群的vat，如果是鸟的话可能也需要一个状态机来进行多种动画的管理（振翅、滑翔、停留），这是下一步应当研究的。
如果想要与场景产生碰撞（躲避障碍物），需要烘焙场景的SDF，或者将场景以简单的几何体进行隐式表示（类似Collider）后传入ComputeShader，让每个单位都感知最近的障碍物位置与距离信息，并对障碍物方向上的速度进行衰减。
之后如果还有更多空余时间，我可能会把这次的工程上传github以供分享学习，这次国庆8天假期里有4天应酬，2天玩伊苏10，2天往返于上海与家乡，时间实在不算充足（迫真）。

那么本文到此结束，感谢阅读，若你能从此文中有所收获，便再好不过了。

感谢您阅读完本文，若您认为本文对您有所帮助，可以将其分享给其他人；若您发现文章对您的利益产生侵害，请联系作者进行删除。

前言

算法介绍

简单CPU版本

简单GPU版本

优化GPU版本

Counting Sort

介绍

计算Grid ID并统计grid中boid数据

对grid进行前缀和

对boid根据grid进行排序

查询相邻格子中的Boid

结语

FEATURED TAGS