這是我偶然發現的一個程式化生成的小demo，並且意外的發現這個Demo可以使用了幾個非常典型的GPU特效製作的技巧，改一改都可以出一道TA面試題了。所以在這裡推薦給大佬們練練手，給萌新們入入門。
簡介
由圓開始的演化
初次接觸到這個案例是在觀看一個Houdini系列教學影片的時候（可以看這裡），這個影片描述了一個叫Sage Jenson的藝術家做的程式化藝術作品（原作者Sage的部落格在這裡），而這位藝術家的靈感又是源自於一篇介紹多頭絨泡菌的生物學論文（好吧我承認我是在套娃）。在這裡會使用Unity的Compute Shader 相關的功能對這個效果進行實現。
我的Unity開原始碼公佈如下：
最終效果展示
老規矩先上效果圖
從《星空》開始進行2D演化
https：//www。zhihu。com/video/1220162023923052544
從《星空》開始進行2D演化的幾個階段
從圓球開始進行3D演化
https：//www。zhihu。com/video/1220523510684090368
從圓球開始進行3D演化 靜幀I
從圓球開始進行3D演化 靜幀II
【演算法原理】
生物學模型在原Paper裡有介紹，但是Sage的部落格裡介紹的更加清楚，這裡借用部落格的圖來解釋一下。
整個細菌的演化系統可以看作是一個巨大的粒子系統，每個細菌是一個粒子。所以上面看到的小白點就代表了在運動中的細菌粒子。為方便計算，假定每個粒子的運動速度恆定，只是運動方向在不斷變化。
那麼在圖中畫出的有顏色的部分又是什麼呢？每個粒子在每一幀都會透過蒐集空間的資訊來調整運動的方向，同時也會在空間內留下自己的資訊。這裡的生物資訊用一個一維向量表示，可看做是資訊素濃度。在上述圖中有顏色的部分其實就代表了資訊素的濃度。我們把空間中的生物資訊儲存在一個二維貼圖裡，我們把它稱之為Trail，並且我們人為的定義所研究的二維空間四方連續（即當粒子走出+x邊界時，會從y座標相同的-x邊界點進入空間，其他三個方向類似）。
系統的更新分為以下所述的六個步驟：
細菌網路演化迭代演算法示意
SENSE： 每個粒子蒐集附近的資訊濃度，實際上只需要取樣三個點，如上圖左上所示。
ROTATE： 根據蒐集到的資訊濃度，進行運動方向的調整，入上圖右上所示。
MOVE： 更新粒子的位置資訊。
DEPOSIT：每個粒子留下各自的資訊，並存儲在Trail裡。
DIFFUSE： 更新Trail，每個格點中的資訊向四周擴散，可以理解為流體的流動。
DECAY： 更新Trail，每個格點中的資訊衰減，一般為百分比衰減，可以理解為揮發。
上述六個步驟，其中1-3是粒子相關，4-6是Trail相關。其中1和4是粒子和Trail互動的步驟，1是需要從Trail讀入資料，並更新粒子的資訊，4則是需要讀取粒子的資訊並更新Trail內的資料。
這裡有一個小細節，就是在1和4互動步驟裡遍歷誰的問題。這個問題在CPU上很好解決，但是GPU程式設計需要考慮架構的問題，即對同一個貼圖無法同時進行讀寫操作，因此在這種資料互動的問題上會比較容易遇到問題。在下面的文章中我會展開講述這個問題，在此先不表。這個方案僅僅是對於這個案例可行的，未必是最優解，如果有更好的方法歡迎留言評論。
瞭解了演算法之後，是不是很想上手敲程式碼了呢。感覺這麼簡明的核心演算法，二三十行的程式碼也就能寫完了呀，甚至數學原理比什麼快排來說可簡單多了。然而悲慘的事實是，就我個人的體驗而言，
真正調演算法的時間也就只佔10%（o(╥﹏╥)o）
，大部分時間還是花在了整體架構的搭建和最後摳細節上面。下面我會把一個我認為比較科學的製作流程給大家分享一下。
【資料結構】
在對核心演算法有初步的瞭解之後，就可以著手設計資料結構了。在設計資料結構的時候，還需要同時考慮資料的生命週期和維護等一系列問題。資料結構的具體細節在開發過程中是會經歷無數次迭代的，甚至可能會經歷若干次大重構。
所以實際的開發過程中，資料結構的大方向找準即可，在前期不需要花太多時間去精雕細琢
。
參考Sage部落格中的內容，一幅圖沒有幾十萬個粒子是做不出效果來的（別問我看著2000粒子結果的表情），所以只能考慮全GPU實現。我採用的是GPU粒子+RT貼圖繪製的方案。
粒子部分，需要每幀更新每個粒子的位置，動力學部分粒子的速度方向需要被記錄，視覺部分可能需要做粒子引數的隨機，可能會考慮把標量速度，標量角速度，感測器引數加入進來。儲存的容器考慮為一個一維陣列，基本上ComputeBuffer可以完成。
Trail部分，是一個二維的網格，很明顯就是RenderTexture（以下簡稱為RT）。由於只要記錄一個浮點數，即生物資訊濃度，考慮只開r通道的浮點數RT。並且考慮到生物濃度資訊的絕對數值是沒有意義的，在演算法裡只需要進行相對數值的計算，所以儘量把儲存的數值往1靠保持精度（好吧，其實只是後面為了做視覺效果需要保持在0-1 （￢_￢））
資料結構示意
粒子
這裡我只選擇了最基礎的資訊進行儲存：
struct
ParticleInfo
{
float3
pos
；
float3
vel
；
}；
在C#裡定義如下（Physarum2D。cs L50）
private
ComputeBuffer
particleInfo
；
在compute shader裡定義如下（Physarum2DUpdate。comput L28）
RWStructuredBuffer
在。shader裡定義如下（BillboardParticles。shader L45）
StructuredBuffer
<
ParticleInfo
>
ParticleInfoBuffer
；
當然上述的ParticleInfo在C#端和Shader端都需要宣告，並且位數需要對齊。我用了一個單獨的ComputeShader來進行資料結構的初始化。
C#端程式碼（Physarum2D。cs L156）
particleInfo
=
new
ComputeBuffer
（
particleCount
，
Marshal
。
SizeOf
（
typeof
（
ParticleInfo
）））；
InitParticle
。
SetBuffer
（
InitParticleHandle
，
ParticleInfoID
，
particleInfo
）；
InitParticle
。
Dispatch
（
InitParticleHandle
，
particleCount
/
threadNum
，
1
，
1
）；
shader端程式碼（Physarum2DInit。compute L20）
［
numthreads
（
THREAD_NUM
，
1
，
1
）］
void
InitParticle
（
uint3
id
：
SV_DispatchThreadID
）
{
ParticleInfoBuffer
［
id
。
x
］。
pos
=
float3
（
rand2
（
id
。
x
）
*
2。0
*
_Size
-
float2
（
1
，
1
）
*
_Size
，
0
）；
ParticleInfoBuffer
［
id
。
x
］。
vel
=
normalize
（
float3
（
rand2
（
id
。
x
*
2
）
*
2
-
float2
（
1
，
1
）
，
0
））；
}
Trail
實際上用了3張RT，因為同時對一張RT進行讀和寫是會產生錯誤的，所以在迭代Trail的時候，需要用兩張RT進行左右互博。並且在步驟4 Deposit的時候，為了方便，我會把粒子的資訊先採集到一個Deposit貼圖裡，然後再和原來的進行merge，所以這個Merge的步驟需要3張貼圖。
關於兩張貼圖左右互博，有一個
小技巧
（Physarum2D。cs L179）
public
void
SwitchRT
（）
{
var
tem
=
trailRT
［
0
］；
trailRT
［
0
］
=
trailRT
［
1
］；
trailRT
［
1
］
=
tem
；
}
public
void
UpdateShader
（）
{
。。。
PhysarumUpdate
。
SetTexture
（
DiffuseTrailHandle
，
TrailReadID
，
trailRT
［
READ
］）；
PhysarumUpdate
。
SetTexture
（
DiffuseTrailHandle
，
TrailWriteID
，
trailRT
［
WRITE
］）；
PhysarumUpdate
。
Dispatch
（
DiffuseTrailHandle
，
texGroupCount
，
texGroupCount
，
1
）；
SwitchRT
（）；
PhysarumUpdate
。
SetTexture
（
DecayTrailHandle
，
TrailReadID
，
trailRT
［
READ
］）；
PhysarumUpdate
。
SetTexture
（
DecayTrailHandle
，
TrailWriteID
，
trailRT
［
WRITE
］）；
PhysarumUpdate
。
Dispatch
（
DecayTrailHandle
，
texGroupCount
，
texGroupCount
，
1
）；
SwitchRT
（）；
。。。
}
嗯，對，就是每次Dispatch完，就交換兩個RT的引用，這樣就每次讀和寫的引用都是固定的，就不用去算哪個RT是輸入，哪個RT是輸出了。
RT的生命週期管理可以參考Physarum2D。cs的SetupRT（）和ReleaseBuffers（）兩個函式。這種管理方法雖然簡單但是可能會過時，至少在SRP裡已經不再使用，所以這裡就不貼程式碼了。
【視覺化】
在製作一個效果的時候，最先會被考慮的往往是效果的演算法，效率之類的問題。事實上視覺化也是非常重要的一環，在實際生產中視覺化水平的高低甚至能很大程度上決定開發的效率（說的就是你Unity）。直接進行演算法的實現非常容易暴斃，特別是GPU程式設計，compute shader的編寫宛若黑盒。在寫完一大段程式碼之後，信心滿滿的按下play鍵~~~誒，黑了耶，等等，剛剛這裡寫錯了，馬上好，馬上好~~~~哇哦，紫了~。所以，這裡建議在進行比較大的專案的編寫時，可以先做好輸出結果的視覺化。從零開始一步一個腳印的進行實現，慢慢寫，比較快。
以這個專案為例，在粒子部分，我們可以“假定”所有粒子開始時隨機分佈在空間中，以勻速運動。先實現粒子的渲染，渲染製作完成後，再把粒子運動的6個迭代步驟加進來。
GPU粒子
關於GPU粒子，有一個很好的開源專案可以參考。想要深入瞭解的可以去看看這個專案，這裡只講一下我在這個專案中進行的最基礎的實現。
由於粒子的資訊被儲存在ComputeBuffer裡，所以可以很簡單的在shader裡拿到粒子的位置資訊，同時需要使用instance的feature，進行粒子的遍歷，在vert shader中完成粒子位置的重新組織。看一下程式碼就懂了：
StructuredBuffer
<
ParticleInfo
>
particles
；
StructuredBuffer
<
float3
>
quad
；
struct
v2f
{
float4
pos
：
POSITION
；
float2
uv
：
TEXCOORD0
；
float4
col
：
COLOR
；
}；
v2f
vert
（
uint
id
：
SV_VertexID
，
uint
inst
：
SV_InstanceID
）
{
v2f
o
；
float3
q
=
quad
［
id
］；
// 把粒子位置投影到View空間，進行偏移，然後再投影到螢幕空間，保證Billboard始終面朝相機
o
。
pos
=
mul
（
UNITY_MATRIX_P
，
mul
（
UNITY_MATRIX_V
，
float4
（
particles
［
inst
］。
pos
，
1。0
f
））
+
float4
（
q
，
0。0
f
）
*
_Size
）；
o
。
uv
=
q
+
0。5
f
；
o
。
col
=
_Color
；
return
o
；
}
fixed4
frag
（
v2f
i
）
：
COLOR
{
return
tex2D
（
_MainTex
，
i
。
uv
）
*
i
。
col
*
i
。
col
。
a
；
}
這個quad陣列是什麼呢，實際上它是一個常量陣列，表示了一個billboard面片的6個點（兩個三角形，有兩個點需要計算兩次，實際上模型是4個點）的座標資訊。
Quad陣列內資訊對應正方形的四個頂點
在C#裡初始化以及渲染的程式碼如下（Physarum2D。cs）
public
void
SetupBuffer
（）
{
。。。
quad
。
SetData
（
new
［］
{
new
Vector3
（-
0。5f
，
0。5f
），
new
Vector3
（
0。5f
，
0。5f
），
new
Vector3
（
0。5f
，-
0。5f
），
new
Vector3
（
0。5f
，-
0。5f
），
new
Vector3
（-
0。5f
，-
0。5f
），
new
Vector3
（-
0。5f
，
0。5f
）
}）；
}
void
OnRenderObject
（）
{
renderMaterial
。
SetBuffer
（
“particles”
，
particleInfo
）；
//傳遞粒子位置資訊
renderMaterial
。
SetBuffer
（
“quad”
，
quad
）；
//傳遞Billboard模型資訊
renderMaterial
。
SetPass
（
0
）；
Graphics
。
DrawProceduralNow
（
MeshTopology
。
Triangles
，
6
，
particleCount
）；
//進行渲染
}
用勻速運動模型來檢查GPU粒子是否正確完成
Trail視覺化
由於Trail是用RT進行儲存，所以渲染Shader裡直接讀RT就可以了。在這基礎上，我根據濃度進行了進行了一個LUT處理。
進行LUT處理前後的Trail渲染結果
Trail渲染LUT實現程式碼（VisualizeTrail。shader L48）
fixed4
frag
（
v2f
i
）
：
SV_Target
{
// sample the texture
fixed4
trail
=
tex2D
（
_MainTex
，
i
。
uv
）；
fixed4
col
=
tex2D
（
_LUT
，
float2
（
clamp
（
trail
。
r
，
0。001
，
0。995
）
，
0
））；
return
col
；
}
【演算法實現】（以及遇到的坑）
墨跡了這麼久，終於可以開始做核心演算法的實現了。其實只要上述部分打通，核心演算法其實實現起來非常的快。
首先核心演算法有6個步驟，分為粒子和Trail的更新兩個部分。我們先看粒子
SENSE： 每個粒子蒐集附近的資訊濃度，實際上只需要取樣三個點。
ROTATE： 根據蒐集到的資訊濃度，進行運動方向的調整。
MOVE： 更新粒子的位置資訊。
基本就是一個簡單的粒子更新嘛。由於每個粒子都是獨立的，互相之間不透過粒子模擬的資料影響，所以可以把這三個步驟合併到一起（Physarum2DUpdate。compute L37）
［
numthreads
（
THREAD_NUM
，
1
，
1
）］
void
UpdateParticle
（
uint3
id
：
SV_DispatchThreadID
）
{
。。。
// Step I ： sense the density
float
densityForward
=
TrailRead
［
IDForward
。
xy
］。
r
；
float
densityLeft
=
TrailRead
［
IDLeft
。
xy
］。
r
；
float
densityRight
=
TrailRead
［
IDRight
。
xy
］。
r
；
// Step II ： Rotate according to the density
if
（
densityForward
<
densityLeft
&&
densityForward
<
densityRight
）
{
// turn randomly
vel
=
rand
（
pos
+
_Time
。
yyy
*
5
）
<
0。5
？
mul
（
_TurnLeftMat
，
vel
）
：
mul
（
_TurnRightMat
，
vel
）；
}
else
if
（
densityLeft
<
densityForward
&&
densityForward
<
densityRight
）
{
// turn right
vel
=
mul
（
_TurnRightMat
，
vel
）；
}
else
if
（
densityLeft
>
densityForward
&&
densityForward
>
densityRight
）
{
// turn left
vel
=
mul
（
_TurnLeftMat
，
vel
）；
}
vel
=
normalize
（
vel
）；
// Step III ： Move the Particles
pos
+=
vel
*
_Speed
*
_DeltaTime
；
pos
=
RepeatPosition
（
pos
，
_Size
）；
ParticleInfoBuffer
［
id
。
x
］。
pos
=
pos
；
ParticleInfoBuffer
［
id
。
x
］。
vel
=
vel
；
}
有兩個小細節需要注意，一個是粒子的運動用旋轉矩陣，為了提高效率，矩陣需要在C#中預算好。二是粒子的運動是四方連續的，簡單來說就是從右邊出去的粒子會再從左邊進來，所以需要套一個RepeatPosition來讓粒子保持在固定範圍內運動。
然後說說5，6步，由於第5步涉及相鄰網格的資料讀取的問題，所以需要兩張Render Texture 來完成這個計算，一張用於儲存上一幀的資訊，一張用來儲存更新過後的資訊。由於GPU計算式並行的，如果只使用一張貼圖的話會觸發死鎖。
程式碼如下（Physarum2DUpdate。compute L103）
// Step V ： Diffuse To the Neighbour
［
numthreads
（
THREAD_NUM
，
THREAD_NUM
，
1
）］
void
DiffuseTrail
（
uint3
id
：
SV_DispatchThreadID
）
{
float4
den
=
TrailRead
［
id
。
xy
］；
float4
depositDen
=
DepositTex
［
id
。
xy
］；
//get the deposit from last step
den
+=
depositDen
；
den
*=
（
1
-
_DiffuseRate
*
9
）；
for
（
int
i
=
-
1
；
i
<=
1
；
++
i
）
{
for
（
int
j
=
-
1
；
j
<=
1
；
++
j
）
{
uint3
target
=
id
；
target
。
x
=
（
target
。
x
+
i
+
_TrailResolution
）
%
_TrailResolution
；
target
。
y
=
（
target
。
y
+
j
+
_TrailResolution
）
%
_TrailResolution
；
den
+=
TrailRead
［
target
。
xy
］
*
_DiffuseRate
；
}
}
TrailWrite
［
id
。
xy
］
=
den
；
}
// Step VI ： Decay the Trail
［
numthreads
（
THREAD_NUM
，
THREAD_NUM
，
1
）］
void
DecayTrail
（
uint3
id
：
SV_DispatchThreadID
）
{
float4
den
=
TrailRead
［
id
。
xy
］；
den
*=
lerp
（
1。0
，
_DecayRate
，
_DeltaTime
）
；
TrailWrite
［
id
。
xy
］
=
den
；
}
最後說說比較煩人的步驟4吧。按照一般做法，是遍歷每一個粒子，然後把對應粒子的對應的Trail格點找出來，在對應格點上增加一個_DepositRate。但是，這時候就需要考慮GPU程式設計的並行問題了，在並行執行的情況下，是沒有辦法同時對一個RT進行讀和寫的：
Trial
［
GetTextureID
（
ParticleID
）
］
+=
_DepositRate
；
//無法對同時對一個RT進行讀和寫
對於這個問題的第一個解決方法，就是逆向思維：既然無法用粒子去改變RT，那麼可以反過來在RT端接收來自粒子的影響就好了。所以這個方法是，不對粒子進行遍歷，而是對先在網格內進行遍歷，每個網格再對所有粒子進行遍歷。即對於RT上的每個畫素，遍歷粒子，考察粒子是否在該畫素上方，若是，則積累一個_DepositRate。
這種方法是能夠精確的計算出來模擬的結果的，然而代價就是執行的複雜度很高，達到了O（nm^2），n為粒子數，m為貼圖尺寸。這會產生什麼問題呢，在一個2000粒子的1024*1024貼圖規模的環境下，就已經無法保持實時了。用RenderDoc捕捉得到的渲染時間如下：
該方法在每幀渲染時間（單位為微秒）
那麼2000粒子的規模能夠達到怎麼樣的效果呢，如圖所示：
不~許~笑！這2000個粒子真的已經很努力了！（其實什麼也幹不了。。。
雖然可以用多級cache之類的方法降複雜度，但是工作量會增加的非常大，並且指標不治本。
所以我們需要回到一開始的方案，能不能找到一個方法，既解決同時讀寫的問題，又把只遍歷一遍粒子的方法變成可能？答案是肯定的。
這又要回到一個做渲染時經常會遇到的問題，即效率和精度之間的平衡。學CS的同學經常會碰到空間和時間的取捨。但是在渲染領域，要是能碰到只通過用空間換時間就能解決的問題那可真的是太幸運了。事實上大部分情況下，渲染效率問題還需要透過平衡精度來解決——畢竟只要肉眼看上去正確，數學上的那點區別就不必在意了。於是在這個案例中，我們選擇了需要犧牲部分精度來換取渲染的效率。
我們認為，在Trail的一個格點裡，無論堆積了多少個粒子，都會產生同樣的Deposit積累。也就是說，我們直接把資訊素積累的+=變成=，這樣RT就從讀寫變成只寫了。實際上這種近似是合理的，因為比如10萬個粒子，被1024*1024個網格一分，實際上每個格子裡平均也就0。1個，我們可以認為這個案例裡的粒子足夠稀疏（好吧，其實也沒有什麼理論的根據，但是從結果來看這個近似是可行的）。但是這只是把新增的部分抽離出來，實際上隨著時變的資訊素濃度還是需要讀寫操作的，所以這裡的做法是額外再開一個DepositRT，把粒子對Trail的改變都記錄到這個DepositRT上，然後再把這個RT和上一幀的TrailRT進行合併。
最終程式碼如下（Physarum2DUpdate。compute L89）
// Step IV ： Deposit from the particle to the trail
［
numthreads
（
THREAD_NUM
，
1
，
1
）］
void
Deposit
（
uint3
id
：
SV_DispatchThreadID
）
{
float3
pos
=
ParticleInfoBuffer
［
id
。
x
］。
pos
；
uint3
pID
=
GetTextureID
（
pos
，
_Size
，
_TrailResolution
）；
DepositTex
［
pID
。
xy
］
=
_DepositRate
*
_DeltaTime
；
}
這種做法的複雜度為o（n），在本人的工作臺式機上跑50萬個粒子也是槓槓的
使用新方法在每幀渲染時間（單位為微秒）
小結
上一個圖展示一下不同sense angle下形成菌落形態
和論文中的實驗結果略有出入，不過就不細摳了，畢竟我們只是做視覺效果的。
2D實現的部分基本已經已經講完了。這個教程以嘮嗑聊心得為主，關於具體的演算法細節還是鼓勵去github裡瞭解，我就不再多講了。3D實現與2D類似，但是會遇到一些新的問題，比如3D的視覺化就沒有2D的那麼輕鬆，同時3D的粒子的演算法在數學上也會更加麻煩一些。另外2D部分的初始化現在是隨機的，如果初始化需要根據一個貼圖來做該怎麼辦呢。這些內容我會在下篇和大家分享。
夾帶一些私貨
如果有人問起，TA是什麼，或許會有回答Technical Artist，幫助技術和美術實現功能對接，幫助生產規範化，提高生產效率的職位。
僅僅如此嗎？
事實上，我相信在每一個TA的心目中都會有一片星辰大海，或是那些讓畫面表現超越了硬體桎梏的精妙演算法，或是那些完美還原了自然界電磁波形態的寫實渲染，又或者，像今天介紹的這個專案的原作者Sage Jenson一樣，是那些用美妙的數學公式搭建出來的，只存在於虛擬世界中的瑰麗圖景。
Procedure Art 也正是在這樣一群對虛擬視覺效果有著熱愛的人們中發展起來的。我在大學的圖形學課時，也曾沉迷過用噪聲製作出各種奇奇怪怪的圖案，那時候換來的只有助教同學們的不解。幸好後來讀研究生時，我逐漸瞭解了這一種藝術流派，並且學習了各種強大的工具，讓自己有機會能夠繼續在數字視覺的海洋中繼續自我放飛。作為一點私心，希望這篇文章能夠拋轉引玉，能夠吸引更多的同學參與到Procedural Art的創作中來。
引用
［1］ Sage Jenson，Personal Blog
［2］Physarum， Entagma，
［3］Characteristics of pattern formation and evolution in approximations of Physarum transport networks， Jeff Jones， 2011
［4］GPU Particle， Robert-K
