始めに
趣味でJobSystemを使ったBoidsシミュレーションを作成しています。そのときにNativeMultiHashMapがとても有用だったので備忘録として記事を書きます。
前提
今回の記事ではBoidsとJobSystemについては一切解説しません。
Boidsについては、Unity Graphics Programming vol.1 の第3章を参照してください。こちらの内容をもとに実装をしています。
github.com
JobSystemはUnity公式の解説動画やネット上にたくさん記事が挙げられているのでそれらを参照してください。
www.youtube.com
環境
- Unity 2021.3.10f1
- Unity.Burst 1.6.6
- Unity.Collections 1.2.4
- Unity.Mathematics
- CPU core i7 8700k
プロジェクトはこちらです。
github.com
Boids実装概要
この記事では全探索と近傍探索それぞれのBoidsを実装します。
その両方で共通となる実装について簡単に触れていきます。
Boidsの設定データ
Boidsでは3つの力(結合、整列、分離)があり、それぞれ重みと影響範囲を持ちます。
また、加えて個体の最大速度やシミュレーション範囲などを設定する必要があります。
これらをInspectorで管理するのは大変なので、ScriptableObjectとして作成しております。
以下のコードは後で解説する全探索用の設定データになります。
AllSearchBoidsSettings
using Unity.Mathematics;
using UnityEngine;
namespace Boids.Settings
{
[CreateAssetMenu(fileName = "AllSearchBoidsSetting", menuName = "Boids/AllSearchSetting")]
public class AllSearchBoidsSetting : ScriptableObject
{
[Header("結合")]
[SerializeField] private float _cohesionWeight;
[SerializeField] private float _cohesionAffectedRadius;
[Header("分離")]
[SerializeField] private float _separationWeight;
[SerializeField] private float _separationAffectedRadius;
[Header("整列")]
[SerializeField] private float _alignmentWeight;
[SerializeField] private float _alignmentAffectedRadius;
[Header("シミュレーション空間")]
[SerializeField] private Vector3 _simulationAreaCenter;
[SerializeField] private Vector3 _simulationAreaScale;
[SerializeField] private float _avoidSimulationAreaWeight;
[Space(20)]
[SerializeField] private float _maxSpeed;
[SerializeField] private float _maxSteerForce;
[Header("個体のスケール")]
[SerializeField] private float3 _instanceScale;
public float CohesionWeight => _cohesionWeight;
public float CohesionAffectedRadiusSqr => _cohesionAffectedRadius * _cohesionAffectedRadius;
public float SeparateWeight => _separationWeight;
public float SeparateAffectedRadiusSqr => _separationAffectedRadius * _separationAffectedRadius;
public float AlignmentWeight => _alignmentWeight;
public float AlignmentAffectedRadiusSqr => _alignmentAffectedRadius * _alignmentAffectedRadius;
public Vector3 SimulationAreaCenter => _simulationAreaCenter;
public Vector3 SimulationAreaScale => _simulationAreaScale / 2;
public float AvoidSimulationAreaWeight => _avoidSimulationAreaWeight;
public float MaxSpeed => _maxSpeed;
public float MaxSteerForce => _maxSteerForce;
public float3 InstanceScale => _instanceScale;
}
}
Boids実行クラス
Boidsのシミュレーションを呼び出して、結果を描画するクラスを作ります。
今回描画する個体数は1000以上なので、描画がボトルネックにならないよう RenderMeshInstanced
を使用します。
今回はInspectorで設定されたenum値に応じて、全探索と近傍探索のシミュレーションを分岐させます。
BoidsSimulator
using System;
using Boids;
using Boids.Settings;
using Unity.Collections;
using UnityEngine;
using UnityEngine.Rendering;
namespace BoidsSimulator
{
public class BoidsSimulator : MonoBehaviour
{
private enum BoidsSimulationType
{
AllSearch,
NeighborSearch,
}
[Header("シミュレーション法")]
[SerializeField] private BoidsSimulationType _boidsSimulationType;
[SerializeField] private int _instanceCount;
[Header("描画関係")]
[SerializeField] private Mesh _mesh;
[SerializeField] private Material _material;
private BoidsData[] _boidsDatas;
private RenderParams _renderParams;
[Header("AllSearchSetting")]
[SerializeField] private AllSearchBoidsSetting _allSearchBoidsSetting;
[Header("NeighborSearchSetting")]
[SerializeField] private NeighborSearchBoidsSetting _neighborSearchBoidsSetting;
private void Start()
{
InitializeBoidsInstance();
_renderParams = new RenderParams(_material) { receiveShadows = true, shadowCastingMode = ShadowCastingMode.On };
}
private void Update()
{
switch (_boidsSimulationType)
{
case BoidsSimulationType.AllSearch:
UpdateBoidsByAllSearch();
break;
case BoidsSimulationType.NeighborSearch:
UpdateBoidsByNeighborSearch();
break;
default:
throw new ArgumentOutOfRangeException();
}
}
private void UpdateBoidsByAllSearch()
{
var allSearchBoidsSimulator = new AllSearchBoidsSimulator(_allSearchBoidsSetting);
var matricesArray = new NativeArray<Matrix4x4>(_instanceCount, Allocator.TempJob);
allSearchBoidsSimulator.Calculate(matricesArray, _boidsDatas);
RendererUtility.InstanceRenderUtility.DrawAll(_mesh, _renderParams, matricesArray);
matricesArray.Dispose();
}
private void UpdateBoidsByNeighborSearch()
{
var matricesArray = new NativeArray<Matrix4x4>(_instanceCount, Allocator.TempJob);
var neighborSearchBoidsSimulator = new NeighborSearchBoidsSimulator(_neighborSearchBoidsSetting);
neighborSearchBoidsSimulator.Calculate(matricesArray, _boidsDatas);
RendererUtility.InstanceRenderUtility.DrawAll(_mesh, _renderParams, matricesArray);
matricesArray.Dispose();
}
private void InitializeBoidsInstance()
{
_boidsDatas = new BoidsData[_instanceCount];
BoidsInitializer.In(_boidsDatas, _allSearchBoidsSetting.SimulationAreaCenter, _allSearchBoidsSetting.SimulationAreaScale, 0.1f);
}
}
}
RendererUtility
using Unity.Collections;
using UnityEngine;
namespace RendererUtility
{
public static class InstanceRenderUtility
{
public static void DrawAll(Mesh mesh, RenderParams renderParams, NativeArray<Matrix4x4> matricesArray)
{
const int instanceCountPerDraw = 1023;
var instanceCount = matricesArray.Length;
for (int i = 0; i < instanceCount; i += instanceCountPerDraw)
{
var length = Mathf.Min(instanceCountPerDraw, instanceCount - i);
Graphics.RenderMeshInstanced(renderParams, mesh, 0, matricesArray, length, i);
}
}
}
}
Boids計算クラス
設定データをもとにBoidsの計算をして、描画用Matrixを求めるクラスを定義します。
このクラス内で計算に必要となるNativeArrayを生成してBoidsの計算を行うJobSystemを実行します。
計算が完了したら、NativeArrayの破棄と計算結果の反映を行います。(コードは全探索から)
AllSearchBoidsSimulator
using Boids.Job;
using Boids.Settings;
using Unity.Collections;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEngine;
namespace Boids
{
public class AllSearchBoidsSimulator
{
private readonly AllSearchBoidsSetting _allSearchBoidsSetting;
public AllSearchBoidsSimulator(AllSearchBoidsSetting boidsSetting)
{
_allSearchBoidsSetting = boidsSetting;
}
public void Calculate(NativeArray<Matrix4x4> instanceMatricesArray, BoidsData[] boidsDatas)
{
var instanceCount = boidsDatas.Length;
var boidsDataArray = new NativeArray<BoidsData>(instanceCount, Allocator.TempJob);
var boidsSteerArray = new NativeArray<float3>(instanceCount, Allocator.TempJob);
boidsDataArray.CopyFrom(boidsDatas);
var boidsJob = new AllSearchBoidsSimulatorJob
(
_allSearchBoidsSetting.CohesionWeight,
_allSearchBoidsSetting.CohesionAffectedRadiusSqr,
_allSearchBoidsSetting.SeparateWeight,
_allSearchBoidsSetting.SeparateAffectedRadiusSqr,
_allSearchBoidsSetting.AlignmentWeight,
_allSearchBoidsSetting.AlignmentAffectedRadiusSqr,
_allSearchBoidsSetting.MaxSpeed,
_allSearchBoidsSetting.MaxSteerForce,
boidsDataArray,
boidsSteerArray
);
var boidsJobHandler = boidsJob.Schedule(instanceCount, 0);
boidsJobHandler.Complete();
var applySteerForce = new ApplySteerForceJob
(
boidsDataArray,
boidsSteerArray,
instanceMatricesArray,
_allSearchBoidsSetting.SimulationAreaCenter,
_allSearchBoidsSetting.SimulationAreaScale,
_allSearchBoidsSetting.AvoidSimulationAreaWeight,
Time.deltaTime,
_allSearchBoidsSetting.MaxSpeed,
_allSearchBoidsSetting.InstanceScale
);
var applySteerForceHandler = applySteerForce.Schedule(instanceCount, 0);
applySteerForceHandler.Complete();
boidsDataArray.CopyTo(boidsDatas);
boidsDataArray.Dispose();
boidsSteerArray.Dispose();
}
}
}
JobSystem
Boidsの計算を2段階に分けて計算します。
- 各個体の操舵力を求める(全探索と近傍探索の2種類)
- 操舵力を各個体に反映、描画用Matrixの計算
1.の各個体の操舵力を求める計算は次のようになります
BoidsSimulatorJob.Execute
public void Execute(int ownIndex)
{
var ownPosition = _boidsDatasRead[ownIndex].Position;
var ownVelocity = _boidsDatasRead[ownIndex].Velocity;
var cohesionPositionSum = new float3();
var cohesionTargetCount = 0;
var separateRepluseSum = new float3();
var separateTargetCount = 0;
var alignmentVelocitySum = new float3();
var alignmentTargetCount = 0;
{
if (ownIndex == targetIndex)
{
continue;
}
var targetPosition = _boidsDatasRead[targetIndex].Position;
var targetVelocity = _boidsDatasRead[targetIndex].Velocity;
var diff = ownPosition - targetPosition;
var distanceSqr = math.dot(diff, diff);
if (distanceSqr <= _cohesionAffectedRadiusSqr)
{
cohesionPositionSum += targetPosition;
cohesionTargetCount++;
}
if (distanceSqr <= _separateAffectedRadiusSqr)
{
separateRepluseSum += math.normalize(diff) / math.sqrt(distanceSqr);
separateTargetCount++;
}
if (distanceSqr <= _alignmentAffectedRadiusSqr)
{
alignmentVelocitySum += targetVelocity;
alignmentTargetCount++;
}
}
var cohesionSteer = new float3();
if (cohesionTargetCount > 0)
{
var cohesionPositionAverage = cohesionPositionSum / cohesionTargetCount;
var cohesionDirection = cohesionPositionAverage - ownPosition;
var cohesionVelocity = math.normalize(cohesionDirection) * _maxSpeed;
cohesionSteer = MathematicsUtility.Limit(cohesionVelocity - ownVelocity, _maxForceSteer);
}
var separateSteer = new float3();
if (separateTargetCount > 0)
{
var separateRepulseAverage = separateRepluseSum / separateTargetCount;
var separateVelocity = math.normalize(separateRepulseAverage) * _maxSpeed;
separateSteer = MathematicsUtility.Limit(separateVelocity - ownVelocity, _maxForceSteer);
}
var alignmentSteer = new float3();
if (alignmentTargetCount > 0)
{
var alignmentVelocityAverage = alignmentVelocitySum / alignmentTargetCount;
var alignmentVelocity = math.normalize(alignmentVelocityAverage) * _maxSpeed;
alignmentSteer = MathematicsUtility.Limit(alignmentVelocity - ownVelocity, _maxForceSteer);
}
_boidsSteerWrite[ownIndex] =
cohesionSteer * _cohesionWeight +
separateSteer * _separateWeight +
alignmentSteer * _alignmentWeight;
}
2.では1.で求めた操舵力を各個体に反映します。こちらは全探索と近傍探索の両方で共通です。
また、ここで個体の描画用Matrixを求めます。(本当は別で分けた方が良さそう)
ApplySteerForceJob
using Boids.Mathematics;
using Unity.Burst;
using Unity.Collections;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEngine;
namespace Boids.Job
{
[BurstCompile]
internal struct ApplySteerForceJob : IJobParallelFor
{
private NativeArray<BoidsData> _boidsDatasWrite;
[ReadOnly] private readonly NativeArray<float3> _boidsForceRead;
[WriteOnly] private NativeArray<Matrix4x4> _instanceMatrices;
[ReadOnly] private readonly float3 _simulationAreaCenter;
[ReadOnly] private readonly float3 _simulationAreaScale;
[ReadOnly] private readonly float _avoidWallWeight;
[ReadOnly] private readonly float _deltaTime;
[ReadOnly] private readonly float _maxSpeed;
[ReadOnly] private readonly float3 _instanceScale;
public ApplySteerForceJob(
NativeArray<BoidsData> boidsDatasWrite,
NativeArray<float3> boidsForceRead,
NativeArray<Matrix4x4> instanceMatrices,
float3 simulationAreaCenter,
float3 simulationAreaScale,
float avoidWallWeight,
float deltaTime,
float maxSpeed,
float3 instanceScale
)
{
_boidsDatasWrite = boidsDatasWrite;
_boidsForceRead = boidsForceRead;
_instanceMatrices = instanceMatrices;
_simulationAreaCenter = simulationAreaCenter;
_simulationAreaScale = simulationAreaScale;
_avoidWallWeight = avoidWallWeight;
_deltaTime = deltaTime;
_maxSpeed = maxSpeed;
_instanceScale = instanceScale;
}
public void Execute(int ownIndex)
{
var boidsData = _boidsDatasWrite[ownIndex];
var force = _boidsForceRead[ownIndex];
force += AvoidAreaEdge(boidsData.Position, _simulationAreaCenter, _simulationAreaScale) * _avoidWallWeight;
var velocity = boidsData.Velocity + (force * _deltaTime);
boidsData.Velocity = MathematicsUtility.Limit(velocity, _maxSpeed);
boidsData.Position += velocity * _deltaTime;
_boidsDatasWrite[ownIndex] = boidsData;
var rotationY = math.atan2(boidsData.Velocity.x, boidsData.Velocity.z);
var rotationX = (float) -math.asin(boidsData.Velocity.y / (math.length(boidsData.Velocity.xyz) + 1e-8));
var rotation = quaternion.Euler(rotationX, rotationY, 0);
_instanceMatrices[ownIndex] = float4x4.TRS(boidsData.Position, rotation, _instanceScale);
}
private static float3 AvoidAreaEdge(float3 position, float3 simulationAreaCenter, float3 simulationAreaScale)
{
var acc = new float3();
acc.x = position.x < simulationAreaCenter.x - simulationAreaScale.x
? acc.x + 1.0f
: acc.x;
acc.x = position.x > simulationAreaCenter.x + simulationAreaScale.x
? acc.x - 1.0f
: acc.x;
acc.y = position.y < simulationAreaCenter.y - simulationAreaScale.y
? acc.y + 1.0f
: acc.y;
acc.y = position.y > simulationAreaCenter.y + simulationAreaScale.y
? acc.y - 1.0f
: acc.y;
acc.z = position.z < simulationAreaCenter.z - simulationAreaScale.z
? acc.z + 1.0f
: acc.z;
acc.z = position.z > simulationAreaCenter.z + simulationAreaScale.z
? acc.z - 1.0f
: acc.z;
return acc;
}
}
}
全探索Boids
実装
全探索なので、各個体は全ての個体を見て自身に与える力を求めます。計算量は O(n2) になります。
AllSearchBoidsSimulatorJob
using Boids.Mathematics;
using Unity.Burst;
using Unity.Collections;
using Unity.Jobs;
using Unity.Mathematics;
using UnityEditor.U2D;
namespace Boids.Job
{
[BurstCompile]
internal struct AllSearchBoidsSimulatorJob : IJobParallelFor
{
[ReadOnly] private readonly float _cohesionWeight;
[ReadOnly] private readonly float _cohesionAffectedRadiusSqr;
[ReadOnly] private readonly float _separateWeight;
[ReadOnly] private readonly float _separateAffectedRadiusSqr;
[ReadOnly] private readonly float _alignmentWeight;
[ReadOnly] private readonly float _alignmentAffectedRadiusSqr;
[ReadOnly] private readonly float _maxSpeed;
[ReadOnly] private readonly float _maxForceSteer;
[ReadOnly] private readonly NativeArray<BoidsData> _boidsDatasRead;
[WriteOnly] private NativeArray<float3> _boidsSteerWrite;
public AllSearchBoidsSimulatorJob(
float cohesionWeight,
float cohesionAffectedRadiusSqr,
float separateWeight,
float separateAffectedRadiusSqr,
float alignmentWeight,
float alignmentAffectedRadiusSqr,
float maxSpeed,
float maxForceSteer,
NativeArray<BoidsData> boidsDatasRead,
NativeArray<float3> boidsSteerWrite
)
{
_cohesionWeight = cohesionWeight;
_cohesionAffectedRadiusSqr = cohesionAffectedRadiusSqr;
_separateWeight = separateWeight;
_separateAffectedRadiusSqr = separateAffectedRadiusSqr;
_alignmentWeight = alignmentWeight;
_alignmentAffectedRadiusSqr = alignmentAffectedRadiusSqr;
_maxSpeed = maxSpeed;
_maxForceSteer = maxForceSteer;
_boidsDatasRead = boidsDatasRead;
_boidsSteerWrite = boidsSteerWrite;
}
public void Execute(int ownIndex)
{
var ownPosition = _boidsDatasRead[ownIndex].Position;
var ownVelocity = _boidsDatasRead[ownIndex].Velocity;
var cohesionPositionSum = new float3();
var cohesionTargetCount = 0;
var separateRepluseSum = new float3();
var separateTargetCount = 0;
var alignmentVelocitySum = new float3();
var alignmentTargetCount = 0;
for (int targetIndex = 0; targetIndex < _boidsDatasRead.Length; ++targetIndex)
{
if (ownIndex == targetIndex)
{
continue;
}
var targetPosition = _boidsDatasRead[targetIndex].Position;
var targetVelocity = _boidsDatasRead[targetIndex].Velocity;
var diff = ownPosition - targetPosition;
var distanceSqr = math.dot(diff, diff);
if (distanceSqr <= _cohesionAffectedRadiusSqr)
{
cohesionPositionSum += targetPosition;
cohesionTargetCount++;
}
if (distanceSqr <= _separateAffectedRadiusSqr)
{
separateRepluseSum += math.normalize(diff) / math.sqrt(distanceSqr);
separateTargetCount++;
}
if (distanceSqr <= _alignmentAffectedRadiusSqr)
{
alignmentVelocitySum += targetVelocity;
alignmentTargetCount++;
}
}
var cohesionSteer = new float3();
if (cohesionTargetCount > 0)
{
var cohesionPositionAverage = cohesionPositionSum / cohesionTargetCount;
var cohesionDirection = cohesionPositionAverage - ownPosition;
var cohesionVelocity = math.normalize(cohesionDirection) * _maxSpeed;
cohesionSteer = MathematicsUtility.Limit(cohesionVelocity - ownVelocity, _maxForceSteer);
}
var separateSteer = new float3();
if (separateTargetCount > 0)
{
var separateRepulseAverage = separateRepluseSum / separateTargetCount;
var separateVelocity = math.normalize(separateRepulseAverage) * _maxSpeed;
separateSteer = MathematicsUtility.Limit(separateVelocity - ownVelocity, _maxForceSteer);
}
var alignmentSteer = new float3();
if (alignmentTargetCount > 0)
{
var alignmentVelocityAverage = alignmentVelocitySum / alignmentTargetCount;
var alignmentVelocity = math.normalize(alignmentVelocityAverage) * _maxSpeed;
alignmentSteer = MathematicsUtility.Limit(alignmentVelocity - ownVelocity, _maxForceSteer);
}
_boidsSteerWrite[ownIndex] =
cohesionSteer * _cohesionWeight +
separateSteer * _separateWeight +
alignmentSteer * _alignmentWeight;
}
}
}
結果
5000匹と10000匹を動かした結果になります。
5000匹はそこまで重くはないですが、10000匹まで行くと25fpsまで落ちてしまいます。
Profilerを見ると AllSearchBoidsSimulatorJob
がボトルネックになっています。
個体数 |
5000匹 |
10000匹 |
動作 |
|
|
Profiler |
|
|
fps |
約80fps |
約25fps |
Updateの時間 |
9.13ms |
36.82ms |
近傍探索Boids
NativeMultiHashMapについて
Struct NativeMultiHashMap<TKey, TValue> | Package Manager UI website
NativeContainerの一つで、キーに対して複数の値を持つことができます。
NativeMultiHashMapのコンストラクタで格納するデータのサイズと AllocatorManager.AllocatorHandle
を指定します。(例のコードでは Allocator.TemJob
)
new NativeMultiHashMap<TKey, TValue>(instanceCount, Allocator.TempJob);
要素の追加は Add
が使えますが、Job内で要素を追加する場合は NativeMultiHashMap<TKey, TValue>.ParallelWriter
で行う必要があります。
var job = new HogeJob { nativeMultiHashMap.AsParallelWriter() };
.....
nativeMultiHasmMapParallelWriter.Add(key, value);
要素へのアクセス方法ですが、NativeMultiHashMapIterator<TKey>
を使います。
TryGetFirstValue
でイテレータと最初の要素を取得して、TryGetNextValue
で次の要素を取得していきます。要素の取得に失敗した場合は false が返されます。
for (var success = nativeMultiHashMap.TryGetFirstValue(key, out var value, out var iterator);
success;
success = nativeMultiHashMap.TryGetNextValue(out value, ref iterator))
{
......
}
実装
始めにシミュレーション空間を格子上に分割して、各個体が所属する格子に個体のIndexを NativeMultiHashMap
保存します。
格子の座標計算ではシミュレーション空間内に限定し、範囲外も空間内に納めています。
NeighborSearchBoidsSetting.Calculate
public void Calculate(NativeArray<Matrix4x4> instanceMatricesArray, BoidsData[] boidsDatas)
{
var instanceCount = boidsDatas.Length;
var boidsDataArray = new NativeArray<BoidsData>(instanceCount, Allocator.TempJob);
var gridMultiHashMap = new NativeMultiHashMap<int3, int>(instanceCount, Allocator.TempJob);
boidsDataArray.CopyFrom(boidsDatas);
var registerInstanceToGridJob = new RegisterNeighborSearchGridJob
(
gridMultiHashMap.AsParallelWriter(),
boidsDataArray,
_neighborSearchBoidsSetting.MinNeighborSearchGridPoint,
_neighborSearchBoidsSetting.NeighborSearchGridScale,
_neighborSearchBoidsSetting.NeighborSearchGridCount
);
var registerInstanceToGridHandler = registerInstanceToGridJob.Schedule(instanceCount, 0);
registerInstanceToGridHandler.Complete();
.........
.........
RegisterNeighborSearchGridJob
using Boids.Mathematics;
using Unity.Collections;
using Unity.Jobs;
using Unity.Mathematics;
namespace Boids.Job
{
internal struct RegisterNeighborSearchGridJob : IJobParallelFor
{
[WriteOnly] private NativeMultiHashMap<int3, int>.ParallelWriter _gridWriter;
[ReadOnly] private readonly NativeArray<BoidsData> _boidsDatasRead;
[ReadOnly] private readonly float3 _minGridPoint;
[ReadOnly] private readonly float _gridScale;
[ReadOnly] private readonly int3 _gridCount;
public RegisterNeighborSearchGridJob(
NativeMultiHashMap<int3, int>.ParallelWriter gridWriter,
NativeArray<BoidsData> boidsDatasRead,
float3 minGridPoint,
float gridScale,
int3 gridCount)
{
_gridWriter = gridWriter;
_boidsDatasRead = boidsDatasRead;
_minGridPoint = minGridPoint;
_gridScale = gridScale;
_gridCount = gridCount;
}
public void Execute(int index)
{
var boidsDataPosition = _boidsDatasRead[index].Position;
var gridIndex = MathematicsUtility.CalculateGridIndex(boidsDataPosition, _minGridPoint, _gridScale, _gridCount);
_gridWriter.Add(gridIndex, index);
}
}
}
MatematicsUtility
internal static int3 CalculateGridIndex(float3 position, float3 minGridPoint, float gridScale, int3 gridCount)
{
return new int3(
(int) math.clamp((position.x - minGridPoint.x) / gridScale, 0, gridCount.x),
(int) math.clamp((position.y - minGridPoint.y) / gridScale, 0, gridCount.y),
(int) math.clamp((position.z - minGridPoint.z) / gridScale, 0, gridCount.z)
);
}
次にBoidsの操舵力を求める処理についてです。
個体が所属する格子を求めて、周囲の格子から個体のIndexを取得します。これにより、周囲の個体のみで操舵力を求められます。
それ以外の処理は基本全探索と同じです。(共通部分は省略)
NeighborSearchBoidsSimulatorJob
using Boids.Mathematics;
using Unity.Burst;
using Unity.Collections;
using Unity.Jobs;
using Unity.Mathematics;
namespace Boids.Job
{
[BurstCompile]
internal struct NeighborSearchBoidsSimulatorJob : IJobParallelFor
{
[ReadOnly] private readonly float _cohesionWeight;
[ReadOnly] private readonly float _cohesionAffectedRadiusSqr;
[ReadOnly] private readonly float _separateWeight;
[ReadOnly] private readonly float _separateAffectedRadiusSqr;
[ReadOnly] private readonly float _alignmentWeight;
[ReadOnly] private readonly float _alignmentAffectedRadiusSqr;
[ReadOnly] private readonly float _maxSpeed;
[ReadOnly] private readonly float _maxForceSteer;
[ReadOnly] private NativeMultiHashMap<int3, int> _grid;
[ReadOnly] private readonly float _gridScale;
[ReadOnly] private readonly int3 _gridCount;
[ReadOnly] private readonly float3 _minGridPoint;
[ReadOnly] private readonly NativeArray<BoidsData> _boidsDatasRead;
[WriteOnly] private NativeArray<float3> _boidsSteerWrite;
public NeighborSearchBoidsSimulatorJob(
float cohesionWeight,
float cohesionAffectedRadiusSqr,
float separateWeight,
float separateAffectedRadiusSqr,
float alignmentWeight,
float alignmentAffectedRadiusSqr,
float maxSpeed,
float maxForceSteer,
NativeMultiHashMap<int3, int> grid,
float gridScale,
int3 gridCount,
float3 minGridPoint,
NativeArray<BoidsData> boidsDatasRead,
NativeArray<float3> boidsSteerWrite
)
{
_cohesionWeight = cohesionWeight;
_cohesionAffectedRadiusSqr = cohesionAffectedRadiusSqr;
_separateWeight = separateWeight;
_separateAffectedRadiusSqr = separateAffectedRadiusSqr;
_alignmentWeight = alignmentWeight;
_alignmentAffectedRadiusSqr = alignmentAffectedRadiusSqr;
_maxSpeed = maxSpeed;
_maxForceSteer = maxForceSteer;
_grid = grid;
_gridScale = gridScale;
_gridCount = gridCount;
_minGridPoint = minGridPoint;
_boidsDatasRead = boidsDatasRead;
_boidsSteerWrite = boidsSteerWrite;
}
public void Execute(int ownIndex)
{
var gridIndex = MathematicsUtility.CalculateGridIndex(ownPosition, _minGridPoint, _gridScale, _gridCount);
int minX = gridIndex.x - 1 < 0 ? 0 : gridIndex.x - 1;
int minY = gridIndex.y - 1 < 0 ? 0 : gridIndex.y - 1;
int minZ = gridIndex.z - 1 < 0 ? 0 : gridIndex.z - 1;
int maxX = gridIndex.x + 1 >= _gridCount.x ? gridIndex.x : gridIndex.x + 1;
int maxY = gridIndex.y + 1 >= _gridCount.y ? gridIndex.y : gridIndex.y + 1;
int maxZ = gridIndex.z + 1 >= _gridCount.z ? gridIndex.z : gridIndex.z + 1;
for (int x = minX; x <= maxX; ++x)
for (int y = minY; y <= maxY; ++y)
for (int z = minZ; z <= maxZ; ++z)
{
var key = new int3(x, y, z);
for (var success = _grid.TryGetFirstValue(key, out var targetIndex, out var iterator);
success;
success = _grid.TryGetNextValue(out targetIndex, ref iterator))
{
if (ownIndex == targetIndex)
{
continue;
}
var targetPosition = _boidsDatasRead[targetIndex].Position;
var targetVelocity = _boidsDatasRead[targetIndex].Velocity;
}
}
}
}
結果
同じく5000匹と10000匹の実行結果になります。
全探索に比べてかなり高速になったことがわかります。また、5000匹の場合は操舵力の計算を求めるのに1ms以下に抑えられております。
個体数 |
5000匹 |
10000匹 |
動作 |
|
|
Profiler |
|
|
fps |
約240fps |
約130fps |
Updateの時間 |
1.55ms |
4.64ms |
比較
全探索と近傍探索で操舵力を求めるのにかかった時間を比較してみます。
(時間は添付した画像から、近傍探索は 格子に登録する時間+操舵力を計算する時間 )
|
5000匹 |
10000匹 |
全探索 |
8.6ms |
36.1ms |
近傍探索 |
0.4ms+0.7ms=1.1ms |
0.7ms+3.2ms=3.9ms |
倍率 |
7.8倍 |
9倍 |
表からわかる通り、10倍近く高速になりました。
全探索では計算量 O(n2) だったのに対して、近傍探索では計算量がおおよそ O(n) になるので大幅に高速化できました。
まとめ
NativeMultiHashMapを使うことでここまで最適化できたことに驚きでした。
これまでBoidsシミュレーションをコンテンツで扱うのは少し難しいと感じていましたが、この結果からそこまで躊躇しなくても良いのではないかと考えております。
今後Boidsで何かしらコンテンツを作っていくので、また何か知見を得られたらまとめます。
始めに
並列処理が出来るJobSystemのサンプルを作成して、どれほどの性能なのかを計ってみようと思います。
何も考えずにサンプルを作って動作確認しているので、JobSystemについての詳細については一切ふれません。
もし、JobSystemについて詳しく知りたい方はUnity公式の動画がおすすめです。
www.youtube.com
www.youtube.com
作成するサンプル
指定された立方体内を反射し続けるCubeを大量に作って、以下の環境でどれだけの性能が出るかを計測しました。
- 愚直な実装 + MeshRenderer
- JobSystem + Burst + MeshRenderer
- JobSystem + Burst + RenderMeshInstanced
また、負荷を高めるために O(n2) の計算を加えています。
ここで載せるコードは雑に書いたものになるので、参考程度にして下さい。
動作環境
愚直な実装 + MeshRenderer
MonobehaviorのUpdateにそのままCubeの移動処理を書きます。
CubeMove.cs
using System.Collections.Generic;
using UnityEngine;
public class CubeMove : MonoBehaviour
{
[SerializeField] private int _instanceCount;
[SerializeField] private GameObject _prefab;
[SerializeField] private float _range;
[SerializeField] private float _velocity;
private readonly List<ObjData> _objDatas = new();
private readonly List<GameObject> _gameObjects = new();
private class ObjData
{
public Vector3 Position;
public Vector3 Velocity;
}
private void Start()
{
for (int i = 0; i < _instanceCount; ++i)
{
var position = RandomPositionInRange();
var instance = Instantiate(_prefab, position, Quaternion.identity);
_objDatas.Add(new ObjData { Position = position, Velocity = RandomDirectionVelocity() } );
_gameObjects.Add(instance);
}
}
private void Update()
{
UpdateObjects();
}
private void UpdateObjects()
{
for (int i = 0; i < _instanceCount; ++i)
{
SumVelocity();
var obj = _objDatas[i];
var deltaTime = Time.deltaTime;
var diff = obj.Velocity * deltaTime;
if (Mathf.Abs(obj.Position.x + diff.x) >= _range) obj.Velocity.x *= -1;
if (Mathf.Abs(obj.Position.y + diff.y) >= _range) obj.Velocity.y *= -1;
if (Mathf.Abs(obj.Position.z + diff.z) >= _range) obj.Velocity.z *= -1;
obj.Position += obj.Velocity * deltaTime;
_gameObjects[i].transform.position = obj.Position;
}
}
private void SumVelocity()
{
var sumVelocity = Vector3.zero;
var sumPosition = Vector3.zero;
for (int i = 0; i < _objDatas.Count; ++i)
{
sumVelocity += _objDatas[i].Velocity;
sumPosition += _objDatas[i].Position;
}
}
private Vector3 RandomPositionInRange()
{
return new Vector3(
Random.Range(-_range, _range),
Random.Range(-_range, _range),
Random.Range(-_range, _range)
);
}
private Vector3 RandomDirectionVelocity()
{
return Random.onUnitSphere * _velocity;
}
}
こちらの実行結果が次のようになります。
cube数が100個の場合は約100fps、1000個の場合は約13fpsになります。
n = 100 |
n = 1000 |
|
|
JobSystem + Burst + MeshRenderer
次にJobSystemとBurstを使ってサンプルを作ってみます。
CubeMoveByJobSystem.cs
using System.Collections.Generic;
using System.Linq;
using Unity.Burst;
using Unity.Collections;
using UnityEngine;
using UnityEngine.Jobs;
using Random = UnityEngine.Random;
public class CubeMoveByJobSystem : MonoBehaviour
{
[SerializeField] private int _instanceCount;
[SerializeField] private GameObject _prefab;
[SerializeField] private float _range;
[SerializeField] private float _velocity;
private ObjData[] _objDatas;
private readonly List<GameObject> _gameObjects = new();
private TransformAccessArray _transformAccessArray;
private struct ObjData
{
public Vector3 Position;
public Vector3 Velocity;
}
private void Start()
{
_objDatas = new ObjData[_instanceCount];
var objDataList = new List<ObjData>();
for (int i = 0; i < _instanceCount; ++i)
{
var position = RandomPositionInRange();
var instance = Instantiate(_prefab, position, Quaternion.identity);
objDataList.Add(new ObjData { Position = position, Velocity = RandomDirectionVelocity() } );
_gameObjects.Add(instance);
}
_objDatas = objDataList.ToArray();
_transformAccessArray = new TransformAccessArray(_gameObjects.Select(obj => obj.transform).ToArray());
}
private void Update()
{
UpdateByJobSystem();
}
private void UpdateByJobSystem()
{
var objDataArray = new NativeArray<ObjData>(_instanceCount, Allocator.TempJob);
var objDataReadArray = new NativeArray<ObjData>(_instanceCount, Allocator.TempJob);
objDataArray.CopyFrom(_objDatas);
objDataReadArray.CopyFrom(_objDatas);
var job = new MoveObjectJob
{
objDatas = objDataArray,
objDatasRead = objDataReadArray,
range = _range,
deltaTime = Time.deltaTime
};
var handler = job.Schedule(_transformAccessArray);
handler.Complete();
objDataArray.CopyTo(_objDatas);
objDataArray.Dispose();
objDataReadArray.Dispose();
}
private Vector3 RandomPositionInRange()
{
return new Vector3(
Random.Range(-_range, _range),
Random.Range(-_range, _range),
Random.Range(-_range, _range)
);
}
private Vector3 RandomDirectionVelocity()
{
return Random.onUnitSphere * _velocity;
}
private void OnDestroy()
{
_transformAccessArray.Dispose();
}
[BurstCompile]
private struct MoveObjectJob : IJobParallelForTransform
{
public NativeArray<ObjData> objDatas;
[ReadOnly] public NativeArray<ObjData> objDatasRead;
[ReadOnly] public float range;
[ReadOnly] public float deltaTime;
public void Execute(int index, TransformAccess transform)
{
SumVelocity();
var obj = objDatas[index];
var diff = obj.Velocity * deltaTime;
if (Mathf.Abs(obj.Position.x + diff.x) >= range) obj.Velocity.x *= -1;
if (Mathf.Abs(obj.Position.y + diff.y) >= range) obj.Velocity.y *= -1;
if (Mathf.Abs(obj.Position.z + diff.z) >= range) obj.Velocity.z *= -1;
obj.Position += obj.Velocity * deltaTime;
objDatas[index] = obj;
transform.localPosition = obj.Position;
}
private void SumVelocity()
{
var sumVelocity = Vector3.zero;
var sumPosition = Vector3.zero;
for (int i = 0; i < objDatasRead.Length; ++i)
{
sumVelocity += objDatasRead[i].Velocity;
sumPosition += objDatasRead[i].Position;
}
}
}
}
こちらの実行結果は次のようになります。
n = 1000の場合は安定して80fps近く出せてますが、n = 5000で40fps行かないぐらい、n = 10000で約15fpsまで落ちました。
n = 1000 |
n = 5000 |
n = 10000 |
|
|
|
JobSystem + Burst + RenderMeshInstanced
先程のサンプルの処理速度を見るに、計算処理ではなく描画処理がボトルネックになっているとわかりました。(計算負荷を上げている個所をコメントアウトしてもフレームレートが変わらない)
なので、MeshRendererではなく、Graphics.RenderMeshInstancedを使って描画をしてみます。(
Unity - Scripting API: Graphics.RenderMeshInstanced)
CubeMoveByJobySystemAndGpuInstancing.cs
using System.Collections.Generic;
using Unity.Burst;
using Unity.Collections;
using Unity.Jobs;
using UnityEngine;
using UnityEngine.Rendering;
using Random = UnityEngine.Random;
public class CubeMoveByJobySystemAndGpuInstancing : MonoBehaviour
{
[SerializeField] private int _instanceCount;
[SerializeField] private float _range;
[SerializeField] private Mesh _mesh;
[SerializeField] private Material _material;
[SerializeField] private float _velocity;
private ObjData[] _objDatas;
private RenderParams _renderParams;
private const int InstanceCountPerDraw = 1023;
private struct ObjData
{
public Vector3 Position;
public Vector3 Velocity;
}
private void Start()
{
_objDatas = new ObjData[_instanceCount];
var objDataList = new List<ObjData>();
for (int i = 0; i < _instanceCount; ++i)
{
var position = RandomPositionInRange();
objDataList.Add(new ObjData { Position = position, Velocity = RandomDirectionVelocity() } );
}
_objDatas = objDataList.ToArray();
_renderParams = new RenderParams(_material) { receiveShadows = true, shadowCastingMode = ShadowCastingMode.On };
}
private void Update()
{
UpdateByJobSystem();
}
private void UpdateByJobSystem()
{
var objDataArray = new NativeArray<ObjData>(_instanceCount, Allocator.TempJob);
var objDataReadArray = new NativeArray<ObjData>(_instanceCount, Allocator.TempJob);
var objMatrix = new NativeArray<Matrix4x4>(_instanceCount, Allocator.TempJob);
objDataArray.CopyFrom(_objDatas);
objDataReadArray.CopyFrom(_objDatas);
var job = new MoveObjectJob
{
objDatas = objDataArray,
objDatasRead = objDataReadArray,
objMatrix = objMatrix,
range = _range,
deltaTime = Time.deltaTime
};
var handler = job.Schedule(_instanceCount, 0);
handler.Complete();
objDataArray.CopyTo(_objDatas);
DrawAll(objMatrix);
objDataArray.Dispose();
objDataReadArray.Dispose();
objMatrix.Dispose();
}
private Vector3 RandomPositionInRange()
{
return new Vector3(
Random.Range(-_range, _range),
Random.Range(-_range, _range),
Random.Range(-_range, _range)
);
}
private Vector3 RandomDirectionVelocity()
{
return Random.onUnitSphere * _velocity;
}
private void DrawAll(NativeArray<Matrix4x4> matricesArray)
{
for (int i = 0; i < _instanceCount; i += InstanceCountPerDraw)
{
var length = Mathf.Min(InstanceCountPerDraw, _instanceCount - i);
Graphics.RenderMeshInstanced(_renderParams, _mesh, 0, matricesArray, length, i);
}
}
[BurstCompile]
private struct MoveObjectJob : IJobParallelFor
{
public NativeArray<ObjData> objDatas;
[ReadOnly] public NativeArray<ObjData> objDatasRead;
[ReadOnly] public float range;
[ReadOnly] public float deltaTime;
[WriteOnly] public NativeArray<Matrix4x4> objMatrix;
public void Execute(int index)
{
SumVelocity();
var obj = objDatas[index];
var diff = obj.Velocity * deltaTime;
if (Mathf.Abs(obj.Position.x + diff.x) >= range) obj.Velocity.x *= -1;
if (Mathf.Abs(obj.Position.y + diff.y) >= range) obj.Velocity.y *= -1;
if (Mathf.Abs(obj.Position.z + diff.z) >= range) obj.Velocity.z *= -1;
obj.Position += obj.Velocity * deltaTime;
objDatas[index] = obj;
objMatrix[index] = Matrix4x4.TRS(obj.Position, Quaternion.identity, Vector3.one);
}
private void SumVelocity()
{
var sumVelocity = Vector3.zero;
var sumPosition = Vector3.zero;
for (int i = 0; i < objDatasRead.Length; ++i)
{
sumVelocity += objDatasRead[i].Velocity;
sumPosition += objDatasRead[i].Position;
}
}
}
}
実行結果は次のとおりです。(マテリアルを設定する必要があるので、Cubeの色が変わっています)
n = 5000 |
n = 10000 |
n = 20000 |
|
|
|
n = 5000では安定して約90fps、n = 10000でも安定して50fps以上でました。
ただし、n = 20000まで来ると約20fpsになります。
計測結果 まとめ
|
n = 1000 |
n = 5000 |
n = 10000 |
愚直な実装 + MeshRenderer |
13fps |
- |
- |
JobSystem + Burst + MeshRenderer |
80fps |
40fps |
15fps |
JobSystem + Burst + RenderMeshInstanced |
140fps |
90fps |
50fps |
まとめ・感想
JobSystemで簡単なサンプルを作ってみましたが、とても取っつきやすい印象でした。実行速度もなかなか良い結果が出たので実用レベルと感じてみます。
また、Graphics.RenderMeshInstancedはNativeArrayをそのまま渡せるので、JobSystemと相性が良いんじゃないかなと思います。
Graphics.RenderMeshInstanced(_renderParams, _mesh, 0, matricesArray);
計測結果から、極限のパフォーマンスを目指すとき以外はJobSystemでも事足りると感じています。(JobSystemを使うときは極限のパフォーマンスを目指すときかもしれませんが....)
ComputeShaderと比べて、学習コストの低さやC#でコードが書けること、Profilerで実行速度を観測できる点でJobSystemの方が圧倒的に扱いやすいです。
参考
shibuya24.info
shibuya24.info
www.youtube.com
目次
始めに
過去のBoidsに関する記事です。Boidsに関してはここではあまり触れないのでご了承ください。
過去記事
shitakami.hateblo.jp
shitakami.hateblo.jp
shitakami.hateblo.jp
障害物から"逃げる"
始めに障害物から逃げる実装について簡単にお話します。
- 障害物の中心座標から個体へのベクトル を求めて、ベクトルの長さが障害物から逃げる範囲の半径 以下であるかをチェック
- 個体が範囲内である場合はベクトル を使って斥力を求めて、個体の速度ベクトルに加算する
結果
こちらを実装した結果が次のgifになります。見て分かる通り、この実装だけでは"避ける"ではなく"逃げる"となり、障害物へ向かってくる個体はそのままぶつかって跳ね返る挙動をします。
この方法でパラメータを調整したり、斥力を距離に反比例するようにしても避けるような挙動を作るのは難しかったです。
障害物を"回避する"
"逃げる"とは別に"回避する"をBoidsに実装します。
障害物を回避する方法は次のようになります。
- 障害物を回避する範囲内にいる、かつ逃げる範囲外にいるかチェック
- 前方向に障害物があるかをチェック
- 個体の速度ベクトル を延長して、障害物の中心点との距離 を求めて障害物から逃げる範囲の半径 以下であるかをチェック
- 1~3の条件を満たした個体の速度ベクトル と、個体から障害物の中心点へのベクトル との外積を計算して回転軸を求める。
- 回転軸と各速度を組み合わせて、速度ベクトルを回転させる。
1. 障害物を回避する範囲内にいる、かつ逃げる範囲外にいるか
こちらは個体から障害物への距離 を求めて、それが < <= であるかをチェックします。( : 逃げる範囲、 : 回避範囲 )範囲外の個体は回避処理は適用させません。
2. 前方向に障害物があるかをチェック
個体の速度ベクトル と個体から障害物へのベクトル との内積 を求めます。(内積を求める前に2つのベクトルを正規化する。) 内積が0以下の場合は と がなす角が鈍角になり、前方に障害物がない状態になります。なので、この場合は回避処理は行いません。
3. 個体の速度ベクトル を延長して、障害物の中心点との距離 を求めて障害物から逃げる範囲の半径 以下であるかをチェック
個体の速度ベクトルをまっすぐ延ばして直線を作ります。次に、障害物の中心点Oから直線に向けて垂直な線分OPを作って、その線分の距離 を求めます。
この距離 が逃げる範囲 以下であるかをチェックします。
こちらの計算は次のように行います。
- 個体の速度ベクトル と 個体から障害物の中心点へのベクトル との外積 を求める
- 求めた外積のベクトルの長さ を計算して、速度ベクトルの長さ で割って距離 を求める
- 距離 が 以下であるかチェックする
上記を1つの式にまとめますと次にようになります。
外積で求められるベクトルの大きさは2つのベクトルが生成する平行四辺形の面積になる性質があります。 (参考 : 「外積の長さ = 平行四辺形の面積」 証明 - 理数アラカルト -)
次に底辺をなすベクトル の大きさで割ることで平行四辺形の高さ = 中心点からの直線への距離を求められます。
4, 5. 回転軸を求めて速度ベクトルを回転させる
これまでの条件を満たした個体に対して障害物を避けるような回転を加えます。
回転軸は3.で求めた外積の単位ベクトル になります。この回転軸と各速度を組み合わせて回転を生成します。
回転の生成ではロドリゲスの回転公式を使っています。(参考 : 3D数学の復習と実践(ロドリゲスの回転公式) - なおしのこれまで、これから)
後は個体の速度ベクトルに回転を適用することで、障害物を避けるような挙動ができます。
コード
以上の実装はComputeShaderで行っております。 参考程度にして頂けると助かります。
inline float3x3 CalcAvoidObstacleTorque(const float3 position, const float3 velocity)
{
const float3 diff = _ObstaclePosition.xyz - position;
const float distance = sqrt(dot(diff, diff));
if(distance > _AvoidObstacleRadius || _EscapeObstacleRadius > distance)
{
return float3x3(float3(1, 0, 0), float3(0, 1, 0), float3(0, 0, 1));
}
const float3 target2ObstacleDirection = normalize(diff);
const float3 forward = normalize(velocity);
const float directionDot = dot(target2ObstacleDirection, forward);
if(directionDot < 0)
{
return float3x3(float3(1, 0, 0), float3(0, 1, 0), float3(0, 0, 1));
}
const float3 forward2DiffCross = cross(forward, diff);
const float forward2ObstacleDistance = length(forward2DiffCross);
if(forward2ObstacleDistance > _EscapeObstacleRadius)
{
return float3x3(float3(1, 0, 0), float3(0, 1, 0), float3(0, 0, 1));
}
const float crossElementSum = forward2DiffCross.x + forward2DiffCross.y + forward2DiffCross.z;
const float3 axis = crossElementSum != 0 ? normalize(forward2DiffCross) : float3(0, 1, 0);
return CalcRotateMatrixByAxis(axis, _AvoidObstacleAngularVelocity);
}
....
....
const float3x3 avoidTorque = CalcAvoidObstacleTorque(boidData.position, velocity);
const float3 avoidVelocity = mul(velocity, avoidTorque);
結果
見やすいようにBoidsの行動範囲を2次元にしています。 分かれずに球体に向かっていく個体は障害物を上下方向に避けようとして詰まっています。
応用
障害物を複数追加して動作するようにした結果は次にようになりました。 (はっきりした方法ではないですが、複数の障害物を回避する場合は回転軸を合成しています。)
障害物で輪っかを作ったら真ん中の穴を通ってくれました。
考察
問題点
この方法で障害物を組み合わせて壁や窪みを作った場合は上手く避けてくれません。(回避先については考慮されていない) この方法で避けられるのは線状のものに限られるでしょう。
壁や窪みが静的なオブジェクトであれば前回のような力場を一度だけ計算して、衝突判定させるのが良いかもしれません。
shitakami.hatenablog.com
複数の障害物に対しての厳密な回避を考える
複数の障害物が接触している場合は、1つの球体に合成して回避する方法もありかもしれません。(しかし、障害物の穴などを塞いでしまうかも)
一番良いと考えているのが、速度ベクトルをずらして複数のRayを飛ばして障害物に当たらないor障害物との交点が最も遠い方向に逃げるようにするのが良いかもしれません。
この方法についてはこちらの動画がとても参考になると思います。
www.youtube.com
最後に
今回は前回できなかったBoidsの回避について実装と考察を行いました。
成果物としては自分の求める要件を満たすものになったので、考察で述べた内容についてはまた遠い将来になると思われます。
これからはこの機能を作りたいコンテンツで使っていく予定です。
参考
今回のBoidsの実装は「Unity Gracphics Programming vol.1」を元にしています。
UnityGraphicsProgramming-vol1.pdf - Google ドライブ
回避のアルゴリズムを考えるうえでこちらの本の計算幾何学の内容を参考にしています。