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