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的一下就上去了。

下面就是代码了：

  1 # me is this DAT.
  2 # 
  3 # scriptOP is the OP which is cooking.
  4 
  5 def cook(scriptOP):
  6     import numpy as np
  7 
  8     button = op(scriptOP.par.string0.eval())
  9 
 10     #get controller's parm
 11     controller = op(scriptOP.par.string1.eval())
 12 
 13     cohesionScale = float(controller['cohesionScale',1])
 14     separationScale = float(controller['separationScale',1])
 15     alignmentScale = float(controller['alignmentScale',1])
 16     targetScale = float(controller['targetScale',1])
 17     boidMinDist = float(controller['boidMinDist',1])
 18     speedLimit = float(controller['speedLimit',1])
 19 
 20     boundaryScale = float(controller['boundaryScale',1])
 21     boundarySize = float(controller['boundarySize',1])
 22 
 23     # get the target position
 24     targetOP = op(scriptOP.par.string2.eval())
 25     targetPos = np.array([targetOP[0].eval(), targetOP[1].eval(), targetOP[2].eval()])
 26 
 27     #use button to refresh the screen
 28     if button[0].eval() != 0:
 29         scriptOP.clear()
 30 
 31     #get all points from the first input
 32     points = scriptOP.points
 33     if not points:
 34         scriptOP.copy(scriptOP.inputs[0])
 35         points = scriptOP.points
 36 
 37     #create the array for position and velocity
 38     oldPosition = np.empty([len(points),3])
 39     oldVelocity = np.empty([len(points),3])
 40     
 41     #input the position and vel data to array
 42     for i in range(len(points)):
 43         oldPosition[i] = [points[i].x, points[i].y, points[i].z]
 44         oldVelocity[i] = [points[i].v[0], points[i].v[1], points[i].v[2]]
 45 
 46     #iterate each point
 47     for i in range(len(points)):
 48         currentPos = oldPosition[i]
 49         currentVel = oldVelocity[i]
 50         restPos = np.delete(oldPosition, i, 0)
 51         restVel = np.delete(oldPosition, i, 0)
 52         
 53         #cohesion rule
 54         cohesionPos = np.mean(restPos, axis = 0)
 55         cohensionForce = (cohesionPos - currentPos) / cohesionScale
 56 
 57         #separation rule
 58         distance = (np.sqrt((restPos-currentPos)**2)).sum(axis=1)
 59         np.place(distance, distance>boidMinDist, 0)
 60         np.place(distance, distance>0, separationScale)
 61         separationForce = -((restPos - currentPos)*distance[:,np.newaxis]).sum(axis = 0)
 62 
 63         #alignment rule
 64         alignVel = np.mean(restVel, axis = 0)
 65         alignmentForce = (alignVel - currentVel) / alignmentScale
 66 
 67         #aim to the target
 68         targetForce = (targetPos - currentPos) / targetScale
 69 
 70         #set the boundary for particles
 71         limits = np.array([[boundarySize,boundarySize,boundarySize],[-boundarySize,-boundarySize,-boundarySize]])
 72         boundForce = np.array([0.0,0.0,0.0])
 73         if currentPos[0] > limits[0,0]:
 74             boundForce[0] = -boundaryScale
 75         elif currentPos[0] < limits[1,0]:
 76             boundForce[0] = boundaryScale
 77         if currentPos[1] > limits[0,1]:
 78             boundForce[1] = -boundaryScale
 79         elif currentPos[1] < limits[1,1]:
 80             boundForce[1] = boundaryScale
 81         if currentPos[2] > limits[0,2]:
 82             boundForce[2] = -boundaryScale
 83         elif currentPos[2] < limits[1,2]:
 84             boundForce[2] = boundaryScale
 85 
 86         ##combine all the forces
 87         forces = currentVel + cohensionForce + separationForce + alignmentForce + targetForce + boundForce
 88 
 89         #limit speed
 90         lengthSpeed = np.linalg.norm(forces)
 91         if lengthSpeed > speedLimit:
 92             forces = forces / lengthSpeed * speedLimit
 93 
 94         #update P, v and N
 95         # update point location and save velocity into Z
 96         points[i].x += forces[0]
 97         points[i].y += forces[1]
 98         points[i].z += forces[2]
 99         points[i].v[0] = forces[0]
100         points[i].v[1] = forces[1]
101         points[i].v[2] = forces[2]
102         points[i].N[0] = forces[0]
103         points[i].N[1] = forces[1]
104         points[i].N[2] = forces[2]
105 
106 
107     return

 

最后捎带讲一讲做这个例子时自己写的一个scatter方法，TD里面很遗憾也没有scatter这个功能，作为Houdini玩家这个绝对不能忍。简单参考了一下scatter的思路，用另一个python脚本做出了一个scatter原型：

 1 # me is this DAT.
 2 # 
 3 # scriptOP is the OP which is cooking.
 4 
 5 def cook(scriptOP):
 6     scriptOP.clear()
 7 
 8     seed = tdu.rand(2325)
 9 
10     geo = op(scriptOP.inputs[0])
11     prims = geo.prims
12 
13     for i in range(200): #"200" is the number of points we create
14         for prim in prims:
15             u = tdu.rand(seed)
16             seed = seed + 1
17             v = tdu.rand(seed)
18     
19             point = scriptOP.appendPoint()
20             point.P = prim.eval(u,v)
21     return

 

posted @ 2015-03-26 00:27  嘉豪一点点  阅读(3602)  评论(1)  编辑  收藏  举报