NekoLab

用p5.js模拟一串粒子的二维运动

代码见这里哦,可以用键盘上下左右交互的网页在这里

这里想模拟的一个效果是一串粒子,只由一个“头部”的牵引发生运动,需要有自然流畅的动态和一定随机性。最终想在三维中做成类似陈星汉Flower中玩家操控的一串花瓣一样的运动,参考这个视频。细想之,并非一个简单的动画效果可以达到,所以想到先做一个二维的原型来试验,之后再用Unity做三维的就容易得多了。

1. p5.js简介与编程踩过的坑

p5.js是Processing语言移植到js语言的一个库,Processing是专门为艺术与交互设计创造出的语言,不需太多技术性的设置,编代码只需要专注于绘图和交互,大大降低了交互图形绘制的门槛。之前看Youtube上的Coding Challenge就用p5.js编写很多图形和可视化算法,给人不少图形学方面的灵感。

做这个不太复杂的原型用Unity C#觉得太重量级了一点,于是借此机会学个新技能,下载p5.min.js只有300多k的源文件,打开文本编辑器就能开始写了。

绘图的核心就是两个函数setup()和draw(),setup只在程序加载之初执行一次,做一些初始化的工作,draw会在每一帧都调用一次,用来更新画布内容和接收用户输入。每一帧时间并不是固定的,而我们后面的动力学模拟需要用到每帧的时间差,所以设置个全局变量计算时间。
var time = 0;
var deltaTime = 0;
function getTime() {
return millis() / 1000;
}
function draw() {
deltaTime = getTime() - time;
time += deltaTime;
// take input,
// and update everything else
}

写习惯了C#, C++这种强类型的面向对象语言,再去写脚本语言真的很不习惯。首先,一切变量的类型都是不确定的,函数参数没有类型检查(不是很熟悉这个,但就算有,也要多些一些代码,还是省省吧),所以类型是一切随缘,出错后果自负。

js的类型系统设计的很奇葩,它有一个运算符===,三个等号代表严格相等,那要问如果不严格相等是什么情况呢?比如一个布尔型false可以等于一个空数组,也等于整数0,还等于只有一个元素0的数组[0]。有下面一张类型比较表,是不是很令人崩溃?

红色:=== 橙色:== 黄色:<= 和 >= 同时成立,== 不成立 蓝色:只有 >= 绿色:只有 <=

语言这么设计貌似是因为早期互联网要防止丢包数据出错的考虑,但我们现在做数值模拟,不需要那么多数据类型,慎重起见所有的等号都用===来判断严格相等吧。

js可以写面向对象的程序,但其中的对象机制其实是个假的面向对象,就使得面向对象的代码写出来特别奇怪,比如我们的ParticleStream类包含数据nodes和成员函数render,就要这么写:
function ParticleStream() {
this.nodes = [];
// and so on
}
ParticleStream.prototype.render = function() {// balabala}

prototype是个内置关键字,理解成定义成员函数的方法就好了,不深究。

接下来就是恶心的部分了。js没有类型重载(归根结底还是因为它不是真正面向对象的),在我们的模拟里面要计算几何位置,会大量用到p5.Vector这个二维矢量类,但是!你想写矢量相加c=a+b是不行的,会报错,应该写成c=p5.Vector.add(a, b),减同理。但是!还有一种常用的操作是a+=b,js就要写成a.add(b),执行这个语句自动会改变a的值,这两个add方法刚开始一不小心就会弄混。

如果不记得不能类型重载这回事的话,有些时候出问题是不会报错的,而是悄悄跳过了执行的语句,或者将变量变成0!比如a和b都是向量,想计算a=b*3,不好意思,这一步过后a就可能变成0了,因为乘法根本不能执行。正确的写法是a = p5.Vector.mult(b, 3)

函数传参默认是传引用,比如a是一个对象作为函数参数,函数里定义新变量var b=a,然后修改b的某个成员,函数返回之后发现原来那个a也变化了!这是因为a和b其实引用的是同一个内存数据,并没有想象中的用复制构造函数建立一个新对象……正确的方法要么是自己定义一个复制构造函数(有点麻烦),要么手动把一个个成员变量赋值过去。既然我们要快速编程算结果,不是为了解决计算机语言问题的,还是选择好理解的后者方案吧。

还有要注意的,类的成员变量要记得加this访问,比如this.velocity就是这个对象的速度,如果没有加this,并且函数范围内恰好有个叫velocity的变量的话……等着算出奇怪的结果然后一行行的debug吧。第一次遇到这种bug极难发现问题,因为C++的思维定势会让人觉得写的velocity就是对象成员变量的那个,根本想不到别的可能性。

说了这么多,都是因为第一次用js被坑得惨,写这个程序有一半时间都在找这些原因导致的诡异bug……

2. 骨架:弹簧质点系统

回到我们的问题上吧。

首先,这一串粒子的运动不能像贪吃蛇一样呆板,必须在运动的时候曲线自然伸展,而停下来的时候能柔和缓慢的收缩,所以就想着用弹簧连接一串质点来模拟骨架的运动。

运动模拟其实就是,在每一帧绘图的时候,知道上一帧时刻的位置与速度,知道这一帧物体受力引起的加速度,来求解当前更新过的位置与速度。在我们程序里面,每一帧draw()调用的时候就调用所有对象的update函数(自己定义的),在update函数里首先把所有力作用在受力物体上,让受力物体积累加速度(牛顿定律:a=F/m),再用数值积分的算法计算新的速度和位置。

积分算法的选取直接影响到显示的结果,一开始用欧拉法发现运动很生硬粗糙,总觉得像卡顿不流畅一样,但帧率是完全没有问题的。换用蛙跳法,结果算出来的运动就如丝流畅。蛙跳法计算量只比欧拉法大一点,精度和稳定性能提高很多,推荐这个方法。

开始采用的弹簧是能伸能缩,伸缩时都会产生按胡克定律F=-kx产生弹力,但这样一串弹簧质点会产生混沌的运动,质点疯狂抖动像一群磕大了的虫子。这应该是一个物理常识,大一时候学过双摆听说那种两根棍子两个枢轴的“简单”系统都会产生混沌运动,更何况一串弹簧,差点忘记了。最后改成了弹簧被拉伸时产生吸引力,被挤压(长度短于原始值)时就没有弹力的作用。

弹簧的另一种力会导致耗散效应,来源于弹簧自身各个部分的摩擦损耗,让弹性力导致的振动不能永久持续下去,而是幅度越来越小最终停止。这里假设耗散力与两端点的相对速度成正比,弹簧伸缩的速度越快,耗散就越多。

我们这个系统唯一的用户输入是头节点的位置,后面所有节点都随之而运动。这一考虑只是因为想要限制用户只通过上下左右操作。头节点是接受输入的节点,设定它运动不受任何力的影响,只跟随输入的坐标值发生平移。

加上这些效应之后,发现前进和停止都没问题,但是转弯的时候很生硬。一个纯物理的模型是有局限的,我们需要加上一些非物理的假想的作用。希望看到的结果是,头部往一个方向转弯时,曲线的身体就能自然往反方向摆动,所以想到可以加上一个与头部运动方向相反的“风力”。

风力应该有多大?在靠近链的中间位置比较大,靠近端点比较小,才能让这条“虫子”的“身体”圆润地转过弯。

风力持续多长时间?我的设定是,在刚开始转弯(也就是头部接收到输入坐标改变时)计时一两秒钟,期间风力大小随时间递减到零。也是为了让运动平滑,不发生速度的突变。

给每个模拟的对象写一个render函数,指定显示在屏幕上的形状。这里让节点、弹簧和节点的速度都绘制出来:

3. 粒子运动:集群模拟(Flock Simulation)

集群模拟是复杂性系统的一个经典问题,常用来解释自然界中鸟群、鱼群这种大规模具有自组织特征的运动。群落中的每个个体依照简单的规则运动,可以用来研究整体大尺度的运动特征。

p5.js官方示例就有对这个问题的模拟,代码直接参考了Nature of Code,是一本讲用代码模拟自然现象的书,可以免费看在线版也可以捐献任意金额获得pdf,epub,mobi的电子书。我捐了一美金下载了,还没来得及看。

集群中的每个个体称为一个Boid,会受到三种基本力的作用:

  • seperation: 靠的太近的时候会分开
  • alignment: 向近邻Boids的运动方向前进
  • cohesion: 向近邻Boids的平均位置(质心)靠近

如何让Boid靠近一个位置?定义了一个seek方法,函数输入值是想要靠近的目标位置,输出会给粒子一个作用力,这个力会把粒子往目标方向推动。

在我们这个问题里,为了让一群粒子跟随定位好的节点,将粒子分组,每个组的粒子跟随一个弹簧链上的节点。

Nature of Code的原始版本中,一个粒子会受到一定空间范围所有粒子的sep, ali, coh力的作用,只要在这个范围内,所有粒子的贡献都是均等的。但我们并不希望这么做,而是希望确定粒子在骨架链上的相对位置,只与在链上相对位置接近的粒子作用,而不是全局空间范围内的粒子,运动起来才会按照链的形态排成一列而不是在空间聚成一团。我目前用的方法是,让粒子只能和同组与前后邻近组的粒子相互作用,相互作用的规则还依照Nature of Code里面的基本算法。

但是写这篇文章的时候想到了另一个方案,可以给每个粒子规定代表在链上整体位置的变量,如同做一个空间坐标变换,比如这样做:设每组a个粒子,第x组的第y个粒子序号就是n=ax+y,具有不同n值的粒子相互作用强度随|n1-n2|差值递减,比如按反比规律。

我完全改变了Cohesion方法的实现,不让Boid靠近近邻的质心,而是让它的速度向同组节点(称之为guide吧)速度靠近。因为是靠近速度而不是靠近位置,加了一个seekVelocity方法,输入想要接近的目标速度,输出会施加一个力,使得粒子速度更接近目标速度。同时还给粒子加了一个空间约束,移动范围不能离guide太远,超出范围就施加一个指向guide位置的力。

4. 改进

前面写的这些基本就可以运行了,但是看起来很死板,粒子会集中在每一组的节点附近,运动起来是一节一节的。做了几个简单改进:

  • 想让粒子流呈现前密后疏的形态,越靠近尾部粒子之间距离越大。使用的方法有:1. 弹簧固有长度前短后长; 2. 弹簧的胡克系数前大后小(越小的系数则弹力越小); 3. Flock的Cohesion方法中约束范围越靠后越大

  • Cohesion方法使用“软约束”,粒子在约束范围距离d以内不受约束力,距离大于d的时候约束力从0开始增长到一个有限值,比如用一个反比例函数来实现。

  • 弹簧的长度、弹性系数与Cohesion方法里的约束范围、约束力相适应,这个多要调试

  • 调节了sep, ali, coh三种力的大小,这里发现相对大小变化一点,最后的运动模式会相差很大,过程中看到了不少有意思的运动状态。