类鸟群：仿真鸟群

仿真鸟群实现与解析

1、思路解析

• 在画布上随机生成一系列boids（鸟类个体）
• 对于所有的鸟类个体，做以下几件事：
  • 应用三大核心规则（分离，列队，内聚）
  • 应用所有附加规则（速度最大值）
  • 应用所有边界条件
• 更新鸟类个体的位置和速度
• 绘制新的位置和速度

2、难点是理解三大核心规则：

规则1：分离

#规则1：分离
        D1 = self.distMatrix < 50.0#若阈值太大，则后面速度vel有可能过大，导致鸟群分崩离析,若阈值太小，则鸟群分离范围探测过小，使个体探测群体范围过小
        # print(f"distMatrix={self.distMatrix}\n")
        pos_f = D1.dot(self.pos) 
        pos_multiply = self.pos*D1.sum(axis=1).reshape(self.N,1)
        vel = pos_multiply - pos_f #通过点的位置生成速度vel向量
        #print(f"vel={vel}\n") #测试
        self.limit(vel,self.maxRuleVel)

规则2：列队

#定义队列的距离阈值，使鸟群沿着特定方向走
        D = self.distMatrix < 10.0 #值越大，探测距离越大，即鸟头开始改变的时刻，即为探测到
        #规则2：列队(alignment)
        vel_2 = D.dot(self.vel)
        self.limit(vel_2,self.maxRuleVel)
        vel += vel_2

规则3：聚集

#规则3：聚集,若没有聚集规则，则鸟的个体之间的位置之间就没有约束
        vel_3 = D.dot(self.pos) - self.pos
        self.limit(vel_3,self.maxRuleVel)
        vel += vel_3

规则2，3都是在vel的基础上增加的，想象一下向量的加法意味什么

3、用范式来对速度vel进行“正则化”？

$$\frac{v_i\ast \epsilon }{norm(v)}=\frac{v_i\ast \epsilon }{\sqrt{{v_1}^2+{v_2}^2+...+{v_n}^2}}(n为v的维度)$$

    def limitVec(self, vec, maxVal):
        """limit magnitide of 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):
        for vec in X:
            self.limitVec(vec,maxVal)

4、边界值的判断

采用平铺小块边界条件

    def applyBC(self):
        """应用边界条件"""
        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

5、完整代码

import numpy as np
import math
import sys,argparse
from scipy.spatial.distance import squareform,pdist,cdist
from numpy.linalg import norm
import matplotlib.pyplot as plt
import matplotlib.animation as animation

width, height = 640, 480

class Boids:
    """仿真鸟群类"""
    def __init__(self,N):
        """初始化鸟群仿真"""
        self.pos = [width/2.0,height/2.0] + 10*np.random.rand(N*2).reshape(N,2)
        # ax1.scatter(self.pos[:,0],self.pos[:,1],c="red")
        angles = 2*math.pi*np.random.rand(N)
        self.vel =np.array(list(zip(np.sin(angles),np.cos(angles))))
        self.N = N
        #最小距离设置为25
        self.minDist = 25.0
        #规定的速度最大增幅
        self.maxRuleVel = 0.03
        #限制最终速度(每一帧距离的变化值)
        self.maxVel = 2.0

    def tick(self, frameNum, pts, beak):
        """Update the simulation by one time step."""
        # get pairwise distances
        self.distMatrix = squareform(pdist(self.pos))
        #应用规则，获取速度vel
        self.vel += self.applyRules()
        self.limit(self.vel,self.maxVel)
        self.pos += self.vel
        self.applyBC()
        #通过vel更新鸟的躯干和头部
        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 magnitide of 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):
        for vec in X:
            self.limitVec(vec,maxVal)

    def applyRules(self):
        #规则1：分离
        D1 = self.distMatrix < 50.0#若阈值太大，则后面速度vel有可能过大，导致鸟群分崩离析,若阈值太小，则鸟群分离范围探测过小，使个体探测群体范围过小
        # print(f"distMatrix={self.distMatrix}\n")
        pos_f = D1.dot(self.pos)
        pos_multiply = self.pos*D1.sum(axis=1).reshape(self.N,1)
        vel = pos_multiply - pos_f #通过点的位置生成速度vel向量
        print(f"vel={vel}\n")
        self.limit(vel,self.maxRuleVel)

        #定义队列的距离阈值，使鸟群沿着特定方向走
        D = self.distMatrix < 10.0 #值越大，探测距离越大，即鸟头开始改变的时刻，即为探测到
        #规则2：列队(alignment)
        vel_2 = D.dot(self.vel)
        self.limit(vel_2,self.maxRuleVel)
        vel += vel_2

        #规则3：聚集,若没有聚集规则，则鸟的个体之间的位置之间就没有约束
        vel_3 = D.dot(self.pos) - self.pos
        self.limit(vel_3,self.maxRuleVel)
        vel += vel_3

        #规则2，3都是在vel的基础上增加的，想象一下向量的加法意味什么

        return vel

    def applyBC(self):
        """应用边界条件"""
        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 buttonPress(self,event):
        #左键添加boid
        if event.button == 1:
            self.pos = np.concatenate((self.pos,np.array([[event.xdata,event.ydata]])),axis=0)
            #随机生成速度
            print(f"add:{event.xdata},{event.ydata}")
            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
        #右键分散boids
        elif event.button == 3:
            #给boids提供一个速度
            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("-------------开始鸟群仿真------------")
    parser = argparse.ArgumentParser(description="实现Craig Reynold鸟群仿真...")
    parser.add_argument('--num-boids',dest='N',required=False)
    args = parser.parse_args()

    N = 3
    if args.N:
        N=int(args.N)

    boids = Boids(N)


    fig,ax1 = plt.subplots()
    ax1 = plt.axes(xlim=(0, width), ylim=(0, height))

    pts, = ax1.plot([],[],markersize=10,c='b',marker='o',ls='None')
    beak, = ax1.plot([],[],markersize=4,c='r',marker='o',ls='None')
    anim = animation.FuncAnimation(fig,tick,fargs=(pts,beak,boids),interval=50)

    #添加一个”点击事件“
    cid = fig.canvas.mpl_connect('button_press_event',boids.buttonPress)
    plt.show()

#call main
if __name__ == "__main__":
    main()