截屏2022-04-01 下午5.54.46.png

在B站上偶然看到一个这样的鱼群模拟的视频，很有意思。大自然很多东西都是类似分形、群体涌现的模式。我们可以用OpenGL简单地模拟一下，先实现一个二维模型，再尝试扩展到三维。

相关资料：
https://blog.csdn.net/liweizhao/article/details/82106886
介绍了简单的集群算法原理
https://www.jianshu.com/p/9881630fb542
用计算着色器实现的集群算法，非常强大

原理分析：每条鱼可以先用点代替，在二维平面上生成一些随机点，它们有坐标和速度，每个点会寻找周围的点，跟随群体移动。关键在于寻找周围的点，每个点不能去遍历所有点的位置，那样时间复杂度是O(N2)。需要对空间进行划分，有四叉树、希尔伯特曲线等方法。感觉在这个问题上将空间均匀划分更好搜索相邻的点，可以划分成8x8或16x16的矩形块，每次搜索3x3的区块，这样时间复杂度就下降到了O(N)。当然极端情况下如果所有点都集中到了3x3的空间里，时间复杂度就会回到O(N2)，这是需要优化的地方。

1 模拟算法

1.1 初始化所有点

用一个数组保存所有点的位置，数组按[x0,y0,x1,y1...]方式保存；另一个数组保存所有点的速度，数组按[vx0,vy0,vx1,vy1...]方式保存；再用一个数组保存所有点的临时新位置，因为要先用老位置计算，全部计算完才能刷新位置，不能用新老位置夹杂计算。点的速度控制在一定范围内。

    float[] pos = new float[MAX_COUNT * 2];
    float[] vel = new float[MAX_COUNT * 2];

    float[] tempPos = new float[pos.length];

    void controlSpeed(int i) {
        float vx = vel[i];
        float vy = vel[i + 1];
        float vv = (float) Math.sqrt(vx * vx + vy * vy);
        float mul = 1f;
        if (vv >= vmax) {
            mul = vmax / vv;
        } else if (vv < vmin) {
            mul = vmin / vv;
        }
        vel[i] *= mul;
        vel[i + 1] *= mul;
    }

1.2 划分区块

划分成16x16个空间区块。用Part表示区块，实际上里面是一个HashSet，保存了区块里面的点；再用一个HashMap保存所有的区块，key是区块的索引，也就是从0~16x16-1顺序排列。每个区块第一次用到时才创建，防止初始化太慢。

    // PART index -> Part

    static int PARTS = 16; // 8 * 8
    static int partW = BOUND / PARTS;
    HashMap<Integer, Part> partMap = new HashMap<>();

    static class Part {

        public Part() {
            set = new HashSet<>();
        }

        HashSet<Integer> set;

    }

查找一个点所在区块

    Part findPart(int posIndex) {
        float x = pos[posIndex];
        float y = pos[posIndex + 1];
        int px = (int)x / partW;
        int py = (int)y / partW;
        int partIndex = px + py * PARTS;
        Part p = partMap.get(partIndex);
        if (p == null) {
            p = new Part();
            partMap.put(partIndex, p);
            Log.d("chao", "create part " + partIndex);
        }
        return p;
    }

查找3x3相邻区块

    int[] neighborsDir = {
            -1, -1, 0, -1, 1, -1,
            -1,  0, 0,  0, 1,  0,
            -1,  1, 0,  1, 1,  1,
    };
    void findNeighborParts(int posIndex) {
        float x = pos[posIndex];
        float y = pos[posIndex + 1];
        int px = (int)x / partW;
        int py = (int)y / partW;
        int iNeighbor = 0;
        for (int i = 0; i < neighborsDir.length - 1; i+=2) {
            int rx = px + neighborsDir[i];
            int ry = py + neighborsDir[i + 1];
            if (rx >= 0 && rx < PARTS && ry >= 0 && ry < PARTS) {
                int partIndex = rx + ry * PARTS;
                Part p = partMap.get(partIndex);
                if (p != null) {
                    for (Integer iPos: p.set) {
                        neighbors[iNeighbor++] = iPos;
                    }
                }
            }
        }
        neighbors[iNeighbor] = -1; // 结尾标志
    }

1.3 核心算法：更新位置

遍历一遍所有点，寻找每个点在相邻区块内的相邻点，保存到neighbors数组里，计算出它们的平均速度，应用到当前点，用该速度计算出新位置，保存到临时列表中。
再次遍历所有点，应用新位置。并比较新老位置，如果新老位置在不同区块中，那么更新这两个区块。

        long time = System.currentTimeMillis();

        for (int i = 0; i < pos.length - 1; i+=2) {

            findNeighborParts(i);

            // 计算群速度：可见的别的点的速度和
            float vxSum = 0;
            float vySum = 0;
            float count = 0;
            for (Integer j : neighbors) {
                if (j == -1) {
                    break;
                }
                if (j == i) {
                    continue;
                }
                float dx = pos[i] - pos[j];
                float dy = pos[i + 1] - pos[j + 1];
                float dis = dx * dx + dy * dy;
                if (dis < visiDis) {
                    vxSum += -dx;
                    vySum += -dy;
                    count++;
                } else if (dis < closeDis) {
                    vxSum -= dx;
                    vySum -= dy;
                    count++;
                }
            }

            // 改变速度
            if (count > 0) {
                float vxAve = vxSum / count;
                float vyAve = vySum / count;
                vel[i] = vel[i] * 0.8f + vxAve * 0.2f;
                vel[i + 1] = vel[i + 1] * 0.8f + vyAve * 0.2f;
                controlSpeed(i);
            }

            // 更新位置，保存到temp里
            float x1 = pos[i] + dt * vel[i];
            float y1 = pos[i + 1] + dt * vel[i + 1];
            boolean out = false;
            if (x1 < 0 || x1 >= BOUND) {
                vel[i] = -vel[i];
                out = true;
            }
            if (y1 < 0 || y1 >= BOUND) {
                vel[i + 1] = -vel[i + 1];
                out = true;
            }
            if (out) {
                tempPos[i] = pos[i];
                tempPos[i + 1] = pos[i + 1];
            } else {
                tempPos[i] = x1;
                tempPos[i + 1] = y1;
            }
        }

        // 更新Part，写入新速度
        for (int i = 0; i < pos.length - 1; i+=2) {
            Part p0 = findPart(i);
            pos[i] = tempPos[i];
            pos[i + 1] = tempPos[i + 1];
            Part p1 = findPart(i);
            if (p0 == p1) {
                continue;
            }
            p0.set.remove(i);
            p1.set.add(i);
        }

        posBuffer.position(0);
        posBuffer.asFloatBuffer().put(pos);

        time = System.currentTimeMillis() - time;
        Log.d("chao", "updae time " + time);
    }

2 OpenGL绘制

2.1 创建顶点数组

顶点数组ByteBuffer posBuffer大小跟位置数组相同，

    int[] vao;
    int[] vbo;

    ByteBuffer posBuffer;

    @Override
    public void onSurfaceCreated(GL10 gl, EGLConfig config) {

//        initPos();
        initAll();

        program = ShaderUtils.loadProgramGroup();

//        //分配内存空间,每个浮点型占4字节空间
        posBuffer = ByteBuffer
                .allocateDirect(2 * 4 * MAX_COUNT)
                .order(ByteOrder.nativeOrder());

        vao = new int[1];
        glGenVertexArrays(1, vao, 0);
        glBindVertexArray(vao[0]);

        vbo = new int[1];
        glGenBuffers(1, vbo, 0);
        glBindBuffer(GL_ARRAY_BUFFER, vbo[0]);
//        glBufferData(GL_ARRAY_BUFFER, vertices.length * 4, vertexBuffer, GL_STATIC_DRAW);
        glBufferData(GL_ARRAY_BUFFER, MAX_COUNT * 4 * 2, posBuffer, GL_STREAM_DRAW);

        // Load the vertex data
        glVertexAttribPointer(0, 2, GL_FLOAT, false, 2 * 4, 0);
        glEnableVertexAttribArray(0);

        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glBindVertexArray(0);

        glClearColor(0.0f, 0.0f, 0.0f, 1.0f);

//        glCullFace(GL_BACK);
//        glEnable(GL_CULL_FACE);
    }

2.2 刷新和绘制

计算完位置后将位置数组传给posBuffer，再将posBuffer传给vbo，用glDrawArrays(GL_POINTS, 0, MAX_COUNT)方式绘制点精灵。

    void updateAll() {
        ...
        posBuffer.position(0);
        posBuffer.asFloatBuffer().put(pos);
    }

    @Override
    public void onDrawFrame(GL10 gl) {
        super.onDrawFrame(gl);

//        update();
        updateAll();

        // Clear the color buffer
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        // 刷新vbo数据
        glBindBuffer(GL_ARRAY_BUFFER, vbo[0]);
        glBufferSubData(GL_ARRAY_BUFFER, 0, MAX_COUNT * 4 * 2,  posBuffer);
        glBindBuffer(GL_ARRAY_BUFFER, 0);


        // Use the program object
        glUseProgram(program);
        glBindVertexArray(vao[0]);

        glDrawArrays(GL_POINTS, 0, MAX_COUNT);
    }

2.3 着色器

顶点着色器将传入的位置转换一下，设置点的大小
shader_group_v.glsl

#version 300 es
layout (location = 0) in vec2 vPosition;

out vec2 vPos;

void main() {
     vPos = (vPosition / 1024.0f - 0.5f) * 2.0f;
     gl_PointSize = 10.0f;
     gl_Position  = vec4(vPos, 0.0f, 1.0f);
}

片段着色器根据位置选择不同颜色

#version 300 es
precision mediump float;

in vec2 vPos;

out vec4 fragColor;

void main() {
     fragColor = vec4((vPos + 1.0f) * 0.5f, 0.5f, 1.0f);
}

这样就完成了！

3 扩展到三维

在二维的基础上修改，主要是位置和速度都增加了z方向，空间划分的区块成了16x16x16个，每个点需要寻找周围3x3x3=27个区块。其实算法的主要时间和空间的复杂度基本上是线性增长，变成原来的1.5倍，是完全没问题的。

    float[] pos = new float[MAX_COUNT * 3];
    float[] vel = new float[MAX_COUNT * 3];

    // 空间中27个方块的三维坐标偏移
    int[] neighborsDir = {-1, -1, -1, -1, -1, 0, -1, -1, 1, -1, 0, -1, -1, 0, 0, -1, 0, 1, -1, 1, -1, -1, 1, 0, -1, 1, 1, 0, -1, -1, 0, -1, 0, 0, -1, 1, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1, -1, 0, 1, 0, 0, 1, 1, 1, -1, -1, 1, -1, 0, 1, -1, 1, 1, 0, -1, 1, 0, 0, 1, 0, 1, 1, 1, -1, 1, 1, 0, 1, 1, 1};
    int[] neighbors = new int[pos.length / 3 + 1];

    Part findPart(int posIndex) {
        float x = pos[posIndex];
        float y = pos[posIndex + 1];
        float z = pos[posIndex + 2];
        int px = (int)x / partW;
        int py = (int)y / partW;
        int pz = (int)z / partW;
        int partIndex = px + py * PARTS + pz * PARTS * PARTS;
        Part p = partMap.get(partIndex);
        if (p == null) {
            p = new Part();
            partMap.put(partIndex, p);
            Log.d("chao", "create part " + partIndex);
        }
        return p;
    }

    void findNeighborParts(int posIndex) {
        float x = pos[posIndex];
        float y = pos[posIndex + 1];
        float z = pos[posIndex + 2];
        int px = (int)x / partW;
        int py = (int)y / partW;
        int pz = (int)z / partW;
        int iNeighbor = 0;
        for (int i = 0; i < neighborsDir.length - 2; i+=3) {
            int rx = px + neighborsDir[i];
            int ry = py + neighborsDir[i + 1];
            int rz = pz + neighborsDir[i + 2];
            if (rx >= 0 && rx < PARTS && ry >= 0 && ry < PARTS && rz >= 0 && rz < PARTS) {
                int partIndex = rx + ry * PARTS + rz * PARTS * PARTS;
                Part p = partMap.get(partIndex);
                if (p != null) {
                    for (Integer iPos: p.set) {
                        neighbors[iNeighbor++] = iPos;
                    }
                }
            }
        }
        neighbors[iNeighbor] = -1; // 结尾标志
    }

OpenGL绘制上有一些增加，变成三维后需要一些空间矩阵变换

    float[] modelMat = new float[16];
    float[] viewMat = new float[16];
    float[] projectionMat = new float[16];

    @Override
    public void onDrawFrame(GL10 gl) {
        super.onDrawFrame(gl);
//        update();
        updateAll();

        // Clear the color buffer
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

        // 刷新vbo数据
        glBindBuffer(GL_ARRAY_BUFFER, vbo[0]);
        glBufferSubData(GL_ARRAY_BUFFER, 0, MAX_COUNT * 4 * 3,  posBuffer);
        glBindBuffer(GL_ARRAY_BUFFER, 0);


        // Use the program object
        glUseProgram(program);

//        Matrix.setIdentityM(modelMat, 0);
        Matrix.setIdentityM(viewMat, 0);
        Matrix.setIdentityM(projectionMat, 0);

        Matrix.perspectiveM(projectionMat, 0, OurCamera.radians(ourCamera.Zoom), width / height, 0.1f, 100.0f);
        ourCamera.GetViewMatrix(viewMat);

        int loc1 = glGetUniformLocation(program, "view");
        glUniformMatrix4fv(loc1, 1, false, viewMat, 0);
        int loc2 = glGetUniformLocation(program, "projection");
        glUniformMatrix4fv(loc2, 1, false, projectionMat, 0);
        Matrix.setIdentityM(modelMat, 0);
        Matrix.translateM(modelMat, 0, 0, 0, 0);
        Matrix.rotateM(modelMat, 0, rx, 0.0f, 1.0f, 0.0f);
        Matrix.rotateM(modelMat, 0, ry, 1.0f, 0.0f, 0.0f);
        int loc3 = glGetUniformLocation(program, "model");
        glUniformMatrix4fv(loc3, 1, false, modelMat, 0);


        glBindVertexArray(vao[0]);
        glDrawArrays(GL_POINTS, 0, MAX_COUNT);

        // 先画不透明的物体
        // 设置标志位
        glDepthMask(false);
        // 再画半透明物体
        boxShader.draw(modelMat, viewMat, projectionMat);
        // 重置标志位
        glDepthMask(true);
    }

顶点着色器
shader_group3d_v.glsl

#version 300 es
layout (location = 0) in vec3 vPosition;

uniform mat4 model;
uniform mat4 view;
uniform mat4 projection;

out vec3 vPos;

void main() {
     vPos = (vPosition / 1024.0f - 0.5f) * 2.0f;
     gl_PointSize = 10.0f;
     gl_Position  = projection * view * model * vec4(vPos, 1.0f);
}