【Python 项目】类鸟群:仿真鸟群

QuantumStack 2024-07-24 14:05:02 阅读 77

类鸟群:仿真鸟群

仔细观察一群鸟或一群鱼,你会发现,虽然群体由个体生物组成,但该群体作为一个整体似乎有它自己的生命。鸟群中的鸟在移动、飞越和绕过障碍物时,彼此之间相互定位。受到打扰或惊吓时会破坏编队,但随后重新集结,仿佛被某种更大的力量控制。

1986年,Craig Reynolds创造鸟类群体行为的一种逼真模拟,称为“类鸟群(Boids)”模型。关于类鸟群模型,值得注意的是,只有 3 个简单的规则控制着群体中个体间的相互作用,但该模型产生的行为类似于真正的鸟群。类鸟群模型被广泛研究,甚至被用来制作群体的计算机动画,如电影“蝙蝠侠归来(1992)”中的行军企鹅。

本项目将利用 Reynolds 的 3 个规则,创建一个类鸟群,模拟 N 只鸟的群体行为,并画出随着时间的推移,它们的位置和运动方向。我们还会提供一个方法,向鸟群中添加一只鸟,以及一种驱散效果,可以用于研究群体的局部干扰效果。类鸟群被称为“N 体模拟”,因为它模拟了 N 个粒子的动态系统,彼此之间施加作用力。

工作原理

模拟类鸟群的三大核心规则如下:

分离:保持类鸟个体之间的最小距离;

列队:让每个类鸟个体指向其局部同伴的平均移动方向;

内聚:让每个类鸟个体朝其局部同伴的质量中心移动。

类鸟群模拟也可以添加其他规则,如避开障碍物,或受到打扰时驱散鸟群,在随后的小节中我们将会探讨这些。这个版本的类鸟群在模拟的每一步中,实现了这些核心规则。

对于群体中的所有类鸟个体,做以下几件事:

应用三大核心规则;应用所有附加规则;应用所有边界条件。更新类鸟个体的位置和速度。绘制新的位置和速度。

如你所见,这些简单的规则创造了一个鸟群,它具有演变的复杂行为。

所需模块

下面是该模拟要用到的 Python 工具:

numpy 数组,用于保存类鸟群的位置和速度;matplotlib 库,用于生成类鸟群动画;argparse,用于处理命令行选项;scipy.spatial.distance 模块,包含一些非常简洁的方法,计算点之间的距离。

代码

首先,要计算类鸟群的位置和速度。接下来,要为模拟设置边界条件,看看如何绘制类鸟群,并实现前面讨论的类鸟群模拟规则。最后,我们会为模拟添加一些有趣的事件,即添加一些类鸟个体和驱散类鸟群。

计算类鸟群的位置和速度

类鸟群仿真需要从 numpy 数组取得信息,计算每一步中类鸟群个体的位置和速度。模拟开始时,将所有类鸟群个体大致放在屏幕中央,速度设置为随机的方向。

<code>import math

import numpy as np

width, height = 640, 480

pos = [width/2.0, height/2.0] + 10*np.random.rand(2*N).reshape(N, 2)

angles = 2*math.pi*np.random.rand(N)

vel = np.array(list(zip(np.sin(angles), np.cos(angles))))

开始在第一行导入 math 模块,用于接下来的计算。在第二行,将 numpy 库导入为 np(少一些录入)。然后,设置屏幕上模拟窗口的宽度和高度。在第四行,创建一个 numpy 数组 pos,对窗口中心加上 10 个单位以内的随机偏移。代码np.random.rand(2 * N)创建了一个一维数组,包含范围在[0,1]的 2N 个随机数。然后 reshape()调用将它转换成二维数组的形状(N,2),它将用于保存类鸟群个体的位置。也要注意,numpy 的广播规则在这里生效:1×2 的数组加到 N×2 的数组的每个元素上。

接下来,用以下方法,创建随机单位速度矢量数组(这些都是模为 1.0 的矢量,指向随机的方向):给定一个角度 t,数字对(cos(t), sin(t))位于半径为 1.0 的圆上,中心在原点(0, 0)。如果从原点到圆上的一点画一条线,就得到一个单位矢量,它取决于角度 A。如果随机选择角度 A,就得到一个随机速度矢量。下图展示了这个方案。

在这里插入图片描述

在第五行,生成一个数组,包含 N 个随机角度,范围在[0, 2pi],在第六行,用前面讨论的随机向量方法生成一个数组,并用内置的 zip()方法将坐标分组。

设置边界条件

鸟儿飞翔在无际的天空,但类鸟群必须在有限的空间中运动。要创建这个空间,就要创建边界条件,在这个例子中,我们采用“平铺小块边界条件”。

将类鸟群模拟想象成发生在一个平铺的空间:如果类鸟群个体离开一个小块,它将从相反的方向进入到相同的小块。环形边界条件和小块边界条件之间的主要区别是,类鸟群模拟不会发生在离散的网格上,而是在一个连续区域移动。下图展示了这些小块边界条件的样子。请看下图中间的小块。飞出右侧的鸟儿正进入右边的小块,但该边界条件确保它们实际上通过平铺在左边的小块,又回到了中心的小块。在顶部和底部的小块,可以看到同样的事情发生。

在这里插入图片描述

下面是如何为类鸟群模拟实现平铺小块边界条件:

<code>def applyBC(self):

"""apply boundary conditions"""

deltaR = 2.0

for coord in self.pos:

if coord[0] > width + deltaR:

coord[0] = - deltaR

if coord[0] < - deltaR:

coord[0] = width + deltaR

if coord[1] > height + deltaR:

coord[1] = - deltaR

if coord[1] < - deltaR:

coord[1] = height + deltaR

在第五行,如果 x 坐标比小块的宽度大,则将它设置回小块的左侧边缘。该行中的 deltaR 提供了一个微小的缓冲区,它允许类鸟群个体开始从相反方向回来之前,稍稍移出小块之外一点,从而产生更好的视觉效果。在小块的左侧、顶部和底部边缘执行类似的检查。

绘制类鸟群

要生成动画,需要知道类鸟群个体的位置和速度,并有办法在每个时间步骤中表示位置和运动方向。

绘制类鸟群个体的身体和头部

为了生成类鸟群动画,我们用 matplotlib 和一点小技巧来绘制位置和速度。将每个类鸟群个体画成两个圆,如下图所示。较大的圆代表身体,较小的圆表示头部。点 P 是身体的中心,H 是头部的中心。根据公式 H = P + k × V 来计算 H 的位置,其中 V 是类鸟群个体的速度,k 是常数。在任何给定时间,类鸟群个体的头指向运动的方向。这指明了类鸟群个体的移动方向,比只画身体更好。

在这里插入图片描述

在下面的代码片段中,利用 matplotlib,用圆形标记画出类鸟群个体的身体。

<code>fig = plt.figure()

ax = plt.axes(xlim=(0, width), ylim=(0, height))

pts, = ax.plot([], [], markersize=10, c='k', marker='o', ls='None')code>

beak, = ax.plot([], [], markersize=4, c='r', marker='o', ls='None')code>

anim = animation.FuncAnimation(fig, tick, fargs=(pts, beak, boids),interval=50)

在第三和第四行分别为类鸟群个体的身体(pts)和头部(beak)标记设置大小和形状。在第五行为动画窗口添加鼠标按钮事件。既然知道了如何绘制身体和喙,让我们看看如何更新它们的位置。

更新类鸟群个体的位置

动画开始后,需要更新身体和头的位置,它指明了类鸟群个体移动的方向。用以下代码来实现:

vec = self.pos + 10*self.vel/self.maxVel

beak.set_sdata(vec.rehape(2*self.N)[::2], vec.reshape(2*self.N)[1::2])

在第一行,计算头部的位置,即在速度(vel)的方向上增加 10 个单位的位移。该位移确定了喙和身体之间的距离。在第二行,用头部位置的新值来更新(reshape)matplotlib 的轴(set_data)。[::2]从速度列表中选出偶数元素(x 轴的值),[1::2]选出奇数元素(Y 轴的值)。

应用类鸟群规则

现在,要在 Python 中实现类鸟群的 3 个规则。我们用“numpy 的方式”来完成这件事,避免循环并利用高度优化的 numpy 方法。

import numpy as np

from scipy.spatial.distance import squareform, pdist, cdist

def test2(pos, radius):

# get distance matrix

distMatrix = squareform(pdist(pos))

# apply threshold

D = distMatrix < radius

# compute velocity

vel = pos*D.sum(axis=1).reshape(N, 1) - D.dot(pos)

return vel

在第五行,用 squareform()和 pdist()方法(在 scipy 库中定义),来计算一组点之间两两的距离(从数组中任意取两点,计算距离,然后针对所有可能的两点对这么做)。

squareform()方法给出一个 3×3 矩阵,其中项 Mij 给出了点 Pi和 Pj 之间的距离。接下来,在第七行,基于距离筛选这个矩阵。

在第九行的方法有点复杂。D.sum()方法按列对矩阵中的 True 值求和。reshape 是必需的,因为和是 N 个值的一维数组(形如(N,)),而你希望它形如(N,1),这样它就能够与位置数组相乘。D.dot()就是矩阵和位置矢量的点积(乘法)。

下面的方法利用前面讨论的 numpy 技术,应用类鸟群的 3 个规则:

def applyRules(self):

# apply rule #1: Separation

D = distMatrix < 25.0

vel = self.pos*D.sum(axis=1).reshape(self.N, 1) - D.dot(self.pos)#应用分离规则时,每个个体都被“推离”相邻个体一定距离

#计算出的速度被限制在某个最大值以内

self.limit(vel, self.maxRuleVel)

# distance threshold for alignment (different from separation)

D = distMatrix < 50.0

# apply rule #2: Alignment

vel2 = D.dot(self.vel)#应用列队规则时,50 个单位的半径内,所有相邻个体的速度之和限制为一个最大值

self.limit(vel2, self.maxRuleVel)

vel += vel2;

# apply rule #3: Cohesion

vel3 = D.dot(self.pos) - self.pos#为每个个体增加一个速度矢量,它指向一定半径内相邻个体的重心或几何中心

self.limit(vel3, self.maxRuleVel)

vel += vel3

return vel

添加个体

类鸟群模拟的核心规则会导致类鸟群展示出群聚行为。但是,让我们在模拟过程中添加一个个体,看看表现如何,让事情变得更有趣。

下面的代码创建一个鼠标事件,让你点击鼠标左键添加一个个体。个体将出现在光标的位置,具有随机指定的速度

# 用 mpl_connect()方法向 matplotlib 画布添加一个按钮按下事件

cid = fig.canvas.mpl_connect('button_press_event', buttonPress)

现在,为了处理鼠标事件,实际创建类鸟群个体,添加以下代码:

def buttonPress(self, event):

"""event handler for matplotlib button presses"""

#确保鼠标事件是左键点击

if event.button is 1:

#将(event.xdata,event.ydata)给出的鼠标位置添加到类鸟群的位置数组

self.pos = np.concatenate((self.pos, np.array([[event.xdata, event.ydata]])), axis=0)

# 将一个随机速度矢量添加到类鸟群的速度数组,并将类鸟群的计数增加 1

angles = 2*math.pi*np.random.rand(1)

v = np.array(list(zip(np.sin(angles), np.cos(angles))))

self.vel = np.concatenate((self.vel, v), axis=0)

self.N += 1

驱散类鸟群

3 个模拟规则保持类鸟群在移动时成为一个群体。但是,群体受到惊扰时,会发生什么?为了模拟这种情况,可以引入一种“驱散”效果:如果在用户界面(UI)窗口中单击右键,群体就会分散。你可以认为这是群体面对突然出现的捕食者的反应,或突然出现一声巨响惊吓了鸟群。下面是实现该效果的一种方式,它作为buttonPress()方法的延续:

# 检查鼠标按键是否是右键单击事件

elif event.button is 3:

# 改变每个个体的速度,在干扰出现的点(即点击鼠标的位置)的相反的方向上增加一个分量

self.vel += 0.1*(self.pos - np.array([[event.xdata, event.ydata]]))

最初,类鸟群将飞离该点,但你会看到,3 个规则胜出,类鸟群将作为群体再次会聚。

命令行参数

下面是类鸟群程序如何处理命令行参数:

parser = argparse.ArgumentParser(description="Implementing CraigReynolds's Boids...")code>

# add arguments

parser.add_argument('--num-boids', dest='N', required=False)code>

args = parser.parse_args()

# set the initial number of boids

N = 100

if args.N:

N = int(args.N)

# create boids

boids = Boids(N)

Boids 类

接下来看看 Boids 类,它代表了模拟。

class Boids:

"""class that represents Boids simulation"""

def __init__(self, N):

"""initialize the Boid simulation"""

# 初始化位置和速度数组

self.pos = [width/2.0, height/2.0] + 10*np.random.rand(2*N).reshape(N, 2)

# normalized random velocities

angles = 2*math.pi*np.random.rand(N)

self.vel = np.array(list(zip(np.sin(angles), np.cos(angles))))

self.N = N

# minimum distance of approach

self.minDist = 25.0

# maximum magnitude of velocities calculated by "rules"

self.maxRuleVel = 0.03

# maximum magnitude of the final velocity

self.maxVel = 2.0

Boid 类处理初始化,更新动画,并应用规则。

boids.tick()在每个时间步骤被调用,以便更新动画,如下所示:

def tick(frameNum, pts, beak, boids):

#print frameNum

"""update function for animation"""

boids.tick(frameNum, pts, beak)

return pts, beak

我们还需要一种方法来限制某些矢量的值。否则,速度将在每个时间步骤无限制地增加,模拟将崩溃。

def limitVec(self, vec, maxVal):

"""limit the magnitude of the 2D vector"""

mag = norm(vec)

if mag > maxVal:

vec[0], vec[1] = vec[0]*maxVal/mag, vec[1]*maxVal/mag

#限制了数组中的值,采用模拟规则计算出的值

def limit(self, X, maxVal):

"""limit the magnitude of 2D vectors in array X to maxValue"""

for vec in X:

self.limitVec(vec, maxVal)

完整代码

下面是类鸟群模拟的完整程序:

import sys, argparse

import math

import numpy as np

import matplotlib.pyplot as plt

import matplotlib.animation as animation

from scipy.spatial.distance import squareform, pdist, cdist

from numpy.linalg import norm

width, height = 640, 480

class Boids:

"""class that represents Boids simulation"""

def __init__(self, N):

"""initialize the Boid simulation"""

# 初始化位置和速度数组

self.pos = [width/2.0, height/2.0] + 10*np.random.rand(2*N).reshape(N, 2)

# normalized random velocities

angles = 2*math.pi*np.random.rand(N)

self.vel = np.array(list(zip(np.sin(angles), np.cos(angles))))

self.N = N

# minimum distance of approach

self.minDist = 25.0

# maximum magnitude of velocities calculated by "rules"

self.maxRuleVel = 0.03

# maximum magnitude of the final velocity

self.maxVel = 2.0

def tick(self, frameNum, pts, beak):

#print frameNum

"""update function for animation"""

self.distMatrix = squareform(pdist(self.pos))

self.vel += self.applyRules()

self.limit(self.vel, self.maxVel)

self.pos += self.vel

self.applyBC()

pts.set_data(self.pos.reshape(2*self.N)[::2],

self.pos.reshape(2*self.N)[1::2])

vec = self.pos + 10*self.vel/self.maxVel

beak.set_data(vec.reshape(2*self.N)[::2],

vec.reshape(2*self.N)[1::2])

def limitVec(self, vec, maxVal):

"""limit the magnitude of the 2D vector"""

mag = norm(vec)

if mag > maxVal:

vec[0], vec[1] = vec[0]*maxVal/mag, vec[1]*maxVal/mag

#限制了数组中的值,采用模拟规则计算出的值

def limit(self, X, maxVal):

"""limit the magnitude of 2D vectors in array X to maxValue"""

for vec in X:

self.limitVec(vec, maxVal)

def applyBC(self):

"""apply boundary conditions"""

deltaR = 2.0

for coord in self.pos:

if coord[0] > width + deltaR:

coord[0] = - deltaR

if coord[0] < - deltaR:

coord[0] = width + deltaR

if coord[1] > height + deltaR:

coord[1] = - deltaR

if coord[1] < - deltaR:

coord[1] = height + deltaR

def applyRules(self):

# apply rule #1: Separation

D = self.distMatrix < 25.0

vel = self.pos*D.sum(axis=1).reshape(self.N, 1) - D.dot(self.pos)

self.limit(vel, self.maxRuleVel)

# distance threshold for alignment (different from separation)

D = self.distMatrix < 50.0

# apply rule #2: Alignment

vel2 = D.dot(self.vel)

self.limit(vel2, self.maxRuleVel)

vel += vel2;

# apply rule #3: Cohesion

vel3 = D.dot(self.pos) - self.pos

self.limit(vel3, self.maxRuleVel)

vel += vel3

return vel

def buttonPress(self, event):

"""event handler for matplotlib button presses"""

# left-click to add a boid

if event.button is 1:

self.pos = np.concatenate((self.pos, np.array([[event.xdata, event.ydata]])), axis=0)

# generate a random velocity

angles = 2*math.pi*np.random.rand(1)

v = np.array(list(zip(np.sin(angles), np.cos(angles))))

self.vel = np.concatenate((self.vel, v), axis=0)

self.N += 1

# right-click to scatter boids

elif event.button is 3:

# add scattering velocity

self.vel += 0.1*(self.pos - np.array([[event.xdata, event.ydata]]))

def tick(frameNum, pts, beak, boids):

boids.tick(frameNum, pts, beak)

return pts, beak

def main():

print('starting boids...')

parser = argparse.ArgumentParser(description="Implementing CraigReynolds's Boids...")code>

# add arguments

parser.add_argument('--num-boids', dest='N', required=False)code>

args = parser.parse_args()

# set the initial number of boids

N = 100

if args.N:

N = int(args.N)

# create boids

boids = Boids(N)

fig = plt.figure()

ax = plt.axes(xlim=(0, width), ylim=(0, height))

pts, = ax.plot([], [], markersize=10, c='k', marker='o', ls='None')code>

beak, = ax.plot([], [], markersize=4, c='r', marker='o', ls='None')code>

anim = animation.FuncAnimation(fig, tick, fargs=(pts, beak, boids),interval=50)

# 用 mpl_connect()方法向 matplotlib 画布添加一个按钮按下事件

cid = fig.canvas.mpl_connect('button_press_event', boids.buttonPress)

plt.show()

if __name__ == '__main__':

main()



声明

本文内容仅代表作者观点,或转载于其他网站,本站不以此文作为商业用途
如有涉及侵权,请联系本站进行删除
转载本站原创文章,请注明来源及作者。