前段时间在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
