Touchdesigner中實現群體模擬boids behaviours simulation

 

前段時間在touchdesigner中學着derivative forum中的大神寫了一個集群的script data。

其實集群效果（boids behavious)本來就是一個非常經典的入門級模擬算法，很多國外的計算機圖形課程都有這個作業，尤其是專攻特效的一些CS課程更是少不了這個。之前在這篇博文中也專門講了怎么在houdini里面用sop solver實現這個模擬效果。另外在touchdesigner中實現這一效果的過程中，我又發現了另一個非常有效直接解釋這個集群模擬的程序模型（偽代碼）：http://www.kfish.org/boids/pseudocode.html

這是實現的效果：

整個算法這次是在python中實現的，在貼出源代碼之前，要講一講touchdesigner(TD)的一個局限。在Houdini中我們處理大量點雲的過程中很多時候會使用到pcopen，這是Houdini內置的非常高效的以索引點為中心搜尋該點一定半徑范圍內的點的功能函數，大方向來講應該就是把KD樹的算法集成優化給我們使用了。但是在TD中沒有這個方法，竟然沒有TvT！！！所以這次的集群模擬只能局限在很少的點數內（幾百個撐死），而且每個個體對於集群的行為是對除了它本身以外的所有點進行采樣計算出來的，這樣的話基本上就看不到小群體行為。因為如果每個點只采樣一定半徑內鄰居的行為的話，就能夠形成大群體有一定規律，同時有些小部分團體也有自己的行為規律，這樣更接近真實的集群效果。

 

另外就是在python中因為沒有像vex中現成的pcopen，所以對其他點的采樣使用了numpy這個函數庫，把采集的所有鄰居點一次性全部堆到一個array中做統一操作，這樣稍稍能把速度提高一些。由於numpy是python中用於科學計算的函數庫，內置的一些操作矩陣的方法效率還是非常高的，但是首先需要對矩陣的變換和計算要有一定基礎，也就是所謂的線性代數不能丟。感覺自己的丟得也差不多了 TvT

 

最后python腳本不是一個點一個點的進行加載運算的。他的方法是一次性把輸入端的幾何體全部拿過來就開始自己的計算和操作。而輸入端的每一個點相當集群里面的一個個體，所以要使用for each循環遍歷每一個點，而每一個點下面又有一堆矩陣的計算，所以計算量duang的一下就上去了。

下面就是代碼了：

# me is this DAT.
# 
# scriptOP is the OP which is cooking.

def cook(scriptOP):
    import numpy as np

    button = op(scriptOP.par.string0.eval())

    #get controller's parm
    controller = op(scriptOP.par.string1.eval())

    cohesionScale = float(controller['cohesionScale',1])
    separationScale = float(controller['separationScale',1])
    alignmentScale = float(controller['alignmentScale',1])
    targetScale = float(controller['targetScale',1])
    boidMinDist = float(controller['boidMinDist',1])
    speedLimit = float(controller['speedLimit',1])

    boundaryScale = float(controller['boundaryScale',1])
    boundarySize = float(controller['boundarySize',1])

    # get the target position
    targetOP = op(scriptOP.par.string2.eval())
    targetPos = np.array([targetOP[0].eval(), targetOP[1].eval(), targetOP[2].eval()])

    #use button to refresh the screen
    if button[0].eval() != 0:
        scriptOP.clear()

    #get all points from the first input
    points = scriptOP.points
    if not points:
        scriptOP.copy(scriptOP.inputs[0])
        points = scriptOP.points

    #create the array for position and velocity
    oldPosition = np.empty([len(points),3])
    oldVelocity = np.empty([len(points),3])
    
    #input the position and vel data to array
    for i in range(len(points)):
        oldPosition[i] = [points[i].x, points[i].y, points[i].z]
        oldVelocity[i] = [points[i].v[0], points[i].v[1], points[i].v[2]]

    #iterate each point
    for i in range(len(points)):
        currentPos = oldPosition[i]
        currentVel = oldVelocity[i]
        restPos = np.delete(oldPosition, i, 0)
        restVel = np.delete(oldPosition, i, 0)
        
        #cohesion rule
        cohesionPos = np.mean(restPos, axis = 0)
        cohensionForce = (cohesionPos - currentPos) / cohesionScale

        #separation rule
        distance = (np.sqrt((restPos-currentPos)**2)).sum(axis=1)
        np.place(distance, distance>boidMinDist, 0)
        np.place(distance, distance>0, separationScale)
        separationForce = -((restPos - currentPos)*distance[:,np.newaxis]).sum(axis = 0)

        #alignment rule
        alignVel = np.mean(restVel, axis = 0)
        alignmentForce = (alignVel - currentVel) / alignmentScale

        #aim to the target
        targetForce = (targetPos - currentPos) / targetScale

        #set the boundary for particles
        limits = np.array([[boundarySize,boundarySize,boundarySize],[-boundarySize,-boundarySize,-boundarySize]])
        boundForce = np.array([0.0,0.0,0.0])
        if currentPos[0] > limits[0,0]:
            boundForce[0] = -boundaryScale
        elif currentPos[0] < limits[1,0]:
            boundForce[0] = boundaryScale
        if currentPos[1] > limits[0,1]:
            boundForce[1] = -boundaryScale
        elif currentPos[1] < limits[1,1]:
            boundForce[1] = boundaryScale
        if currentPos[2] > limits[0,2]:
            boundForce[2] = -boundaryScale
        elif currentPos[2] < limits[1,2]:
            boundForce[2] = boundaryScale

        ##combine all the forces
        forces = currentVel + cohensionForce + separationForce + alignmentForce + targetForce + boundForce

        #limit speed
        lengthSpeed = np.linalg.norm(forces)
        if lengthSpeed > speedLimit:
            forces = forces / lengthSpeed * speedLimit

        #update P, v and N
        # update point location and save velocity into Z
        points[i].x += forces[0]
        points[i].y += forces[1]
        points[i].z += forces[2]
        points[i].v[0] = forces[0]
        points[i].v[1] = forces[1]
        points[i].v[2] = forces[2]
        points[i].N[0] = forces[0]
        points[i].N[1] = forces[1]
        points[i].N[2] = forces[2]


    return

 

最后捎帶講一講做這個例子時自己寫的一個scatter方法，TD里面很遺憾也沒有scatter這個功能，作為Houdini玩家這個絕對不能忍。簡單參考了一下scatter的思路，用另一個python腳本做出了一個scatter原型：

# me is this DAT.
# 
# scriptOP is the OP which is cooking.

def cook(scriptOP):
    scriptOP.clear()

    seed = tdu.rand(2325)

    geo = op(scriptOP.inputs[0])
    prims = geo.prims

    for i in range(200): #"200" is the number of points we create
        for prim in prims:
            u = tdu.rand(seed)
            seed = seed + 1
            v = tdu.rand(seed)
    
            point = scriptOP.appendPoint()
            point.P = prim.eval(u,v)
    return