Monogame C#中反应扩散算法性能优化：600x600网格仅10FPS

我按照Karl Sims的教程实现了反应扩散算法，功能正常，但在600x600网格下仅能达到10FPS，帧率太低没法用到我的应用里。想请教怎么优化实现，或者是不是我误用了Monogame的机制？

现有Update实现 #

public void Update()
{
    for (int x = 1; x < N-1; x++)
    {
        for (int y = 1; y < M-1; y++)
        {
            float a = _a[Index(x, y)];
            float b = _b[Index(x, y)];

            _bufferA[Index(x, y)] = a + DT * (DA * laplacian(x, y, _a) - a * b * b + Feed * (1 - a));
            _bufferB[Index(x, y)] = b + DT * (DB * laplacian(x, y, _b) + a * b * b - (Kill + Feed) * b);
        }
    }

    // Update the grid to the new values.
    (_a, _bufferA) = (_bufferA, _a);
    (_b, _bufferB) = (_bufferB, _b);
}

// Convolution that calculates the difference between chemical concentrations of neighbouring cells.
private float laplacian(int x, int y, float[] arr)
{
    float sum = 0;
    for (int i = -1; i <= 1; i++)
        for (int j = -1; j <= 1; j++)
            switch (Math.Abs(i) + Math.Abs(j))
            {
                case 0: sum += -1 * arr[Index(x, y)]; break;
                case 1: sum += 0.2f * arr[Index(x + i, y + j)]; break;
                case 2: sum += 0.05f * arr[Index(x + i, y + j)]; break;
            }
    return sum;
}

现有Draw实现 #

public void Draw(Texture2D canvas)
{
    Color[] data = new Color[Const.Width * Const.Height];
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < M; j++)
        {
            data[i + j * N] = new Color((int)(255 * _a[Index(i, j)]), 0, (int)(255 * _b[Index(i, j)]));
        }
    }
    canvas.SetData(data);
}

优化建议 #

1. 内联拉普拉斯计算，消除函数调用开销

每次调用laplacian都会带来函数栈开销，且重复计算邻域索引。直接把拉普拉斯计算逻辑嵌入Update循环，能显著减少开销：

public void Update()
{
    int maxX = N - 1;
    int maxY = M - 1;
    for (int x = 1; x < maxX; x++)
    {
        for (int y = 1; y < maxY; y++)
        {
            int idx = Index(x, y);
            float a = _a[idx];
            float b = _b[idx];

            // 计算A的拉普拉斯
            float lapA = -1 * a;
            lapA += 0.2f * (_a[Index(x-1, y)] + _a[Index(x+1, y)] + _a[Index(x, y-1)] + _a[Index(x, y+1)]);
            lapA += 0.05f * (_a[Index(x-1, y-1)] + _a[Index(x-1, y+1)] + _a[Index(x+1, y-1)] + _a[Index(x+1, y+1)]);

            // 计算B的拉普拉斯
            float lapB = -1 * b;
            lapB += 0.2f * (_b[Index(x-1, y)] + _b[Index(x+1, y)] + _b[Index(x, y-1)] + _b[Index(x, y+1)]);
            lapB += 0.05f * (_b[Index(x-1, y-1)] + _b[Index(x-1, y+1)] + _b[Index(x+1, y-1)] + _b[Index(x+1, y+1)]);

            _bufferA[idx] = a + DT * (DA * lapA - a * b * b + Feed * (1 - a));
            _bufferB[idx] = b + DT * (DB * lapB + a * b * b - (Kill + Feed) * b);
        }
    }

    (_a, _bufferA) = (_bufferA, _a);
    (_b, _bufferB) = (_bufferB, _b);
}

2. 预计算邻域索引，减少重复计算

Index(x,y)在循环中被多次调用，提前计算当前点和所有邻域的索引，避免重复计算：

public void Update()
{
    int maxX = N - 1;
    int maxY = M - 1;
    for (int x = 1; x < maxX; x++)
    {
        for (int y = 1; y < maxY; y++)
        {
            int idx = x + y * N; // 直接内联Index逻辑，避免函数调用
            int idxLeft = (x-1) + y * N;
            int idxRight = (x+1) + y * N;
            int idxUp = x + (y-1) * N;
            int idxDown = x + (y+1) * N;
            int idxUpLeft = (x-1) + (y-1) * N;
            int idxUpRight = (x+1) + (y-1) * N;
            int idxDownLeft = (x-1) + (y+1) * N;
            int idxDownRight = (x+1) + (y+1) * N;

            float a = _a[idx];
            float b = _b[idx];

            float lapA = -1 * a 
                + 0.2f * (_a[idxLeft] + _a[idxRight] + _a[idxUp] + _a[idxDown])
                + 0.05f * (_a[idxUpLeft] + _a[idxUpRight] + _a[idxDownLeft] + _a[idxDownRight]);

            float lapB = -1 * b 
                + 0.2f * (_b[idxLeft] + _b[idxRight] + _b[idxUp] + _b[idxDown])
                + 0.05f * (_b[idxUpLeft] + _b[idxUpRight] + _b[idxDownLeft] + _b[idxDownRight]);

            _bufferA[idx] = a + DT * (DA * lapA - a * b * b + Feed * (1 - a));
            _bufferB[idx] = b + DT * (DB * lapB + a * b * b - (Kill + Feed) * b);
        }
    }

    (_a, _bufferA) = (_bufferA, _a);
    (_b, _bufferB) = (_bufferB, _b);
}

3. 优化Draw方法，避免重复内存分配

每次Draw都创建新的Color[]会触发GC，导致帧率波动。把Color数组作为类成员初始化一次，每次只更新内容：

private Color[] _drawData;

// 在类构造函数中初始化
public YourClass()
{
    _drawData = new Color[Const.Width * Const.Height];
}

public void Draw(Texture2D canvas)
{
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < M; j++)
        {
            int idx = i + j * N;
            _drawData[idx] = new Color(
                (byte)(255 * _a[idx]), // 用byte代替int，减少类型转换开销
                0, 
                (byte)(255 * _b[idx])
            );
        }
    }
    canvas.SetData(_drawData);
}

4. 用SIMD指令并行加速计算

利用.NET的System.Numerics.Vector<float>一次处理多个float值，大幅提升计算效率：

using System.Numerics;

public void Update()
{
    int maxX = N - 1;
    int maxY = M - 1;
    int vectorSize = Vector<float>.Count;
    for (int y = 1; y < maxY; y++)
    {
        // 处理能被向量长度整除的批量像素
        for (int x = 1; x <= maxX - vectorSize; x += vectorSize)
        {
            int baseIdx = x + y * N;
            // 批量加载当前A、B值
            Vector<float> aVec = new Vector<float>(_a, baseIdx);
            Vector<float> bVec = new Vector<float>(_b, baseIdx);

            // 批量加载邻域值
            Vector<float> aLeft = new Vector<float>(_a, baseIdx - 1);
            Vector<float> aRight = new Vector<float>(_a, baseIdx + 1);
            Vector<float> aUp = new Vector<float>(_a, baseIdx - N);
            Vector<float> aDown = new Vector<float>(_a, baseIdx + N);

            Vector<float> bLeft = new Vector<float>(_b, baseIdx - 1);
            Vector<float> bRight = new Vector<float>(_b, baseIdx + 1);
            Vector<float> bUp = new Vector<float>(_b, baseIdx - N);
            Vector<float> bDown = new Vector<float>(_b, baseIdx + N);

            // 计算拉普拉斯
            Vector<float> lapAVec = -1 * aVec 
                + 0.2f * (aLeft + aRight + aUp + aDown)
                + 0.05f * (
                    new Vector<float>(_a, baseIdx - N -1) + new Vector<float>(_a, baseIdx - N +1)
                    + new Vector<float>(_a, baseIdx + N -1) + new Vector<float>(_a, baseIdx + N +1)
                );

            Vector<float> lapBVec = -1 * bVec 
                + 0.2f * (bLeft + bRight + bUp + bDown)
                + 0.05f * (
                    new Vector<float>(_b, baseIdx - N -1) + new Vector<float>(_b, baseIdx - N +1)
                    + new Vector<float>(_b, baseIdx + N -1) + new Vector<float>(_b, baseIdx + N +1)
                );

            // 计算更新值
            Vector<float> bSquared = bVec * bVec;
            Vector<float> aUpdate = aVec + DT * (DA * lapAVec - aVec * bSquared + Feed * (Vector<float>.One - aVec));
            Vector<float> bUpdate = bVec + DT * (DB * lapBVec + aVec * bSquared - (Kill + Feed) * bVec);

            // 保存到缓冲区
            aUpdate.CopyTo(_bufferA, baseIdx);
            bUpdate.CopyTo(_bufferB, baseIdx);
        }

        // 处理剩余的单个像素
        for (int x = maxX - (maxX % vectorSize); x < maxX; x++)
        {
            int idx = x + y * N;
            float a = _a[idx];
            float b = _b[idx];
            
            float lapA = -1 * a 
                + 0.2f * (_a[idx-1] + _a[idx+1] + _a[idx-N] + _a[idx+N])
                + 0.05f * (_a[idx-N-1] + _a[idx-N+1] + _a[idx+N-1] + _a[idx+N+1]);
            
            float lapB = -1 * b 
                + 0.2f * (_b[idx-1] + _b[idx+1] + _b[idx-N] + _b[idx+N])
                + 0.05f * (_b[idx-N-1] + _b[idx-N+1] + _b[idx+N-1] + _b[idx+N+1]);
            
            _bufferA[idx] = a + DT * (DA * lapA - a * b * b + Feed * (1 - a));
            _bufferB[idx] = b + DT * (DB * lapB + a * b * b - (Kill + Feed) * b);
        }
    }

    (_a, _bufferA) = (_bufferA, _a);
    (_b, _bufferB) = (_bufferB, _b);
}

5. 转移计算到GPU（Shader实现）

这是大网格场景下最有效的优化方式，利用GPU并行计算能力：

• 创建两个RenderTarget2D作为A、B的纹理载体，交替作为输入输出
• 编写HLSL像素着色器实现反应扩散逻辑
• 每帧在Monogame中渲染到RenderTarget完成GPU计算

示例HLSL像素着色器（简化版）：

texture InputA;
texture InputB;
sampler InputSamplerA = sampler_state
{
    Texture = <InputA>;
    Filter = POINT;
    AddressU = CLAMP;
    AddressV = CLAMP;
};
sampler InputSamplerB = sampler_state
{
    Texture = <InputB>;
    Filter = POINT;
    AddressU = CLAMP;
    AddressV = CLAMP;
};

float DT;
float DA;
float DB;
float Feed;
float Kill;
float2 GridSize;

float4 PixelShaderFunction(float2 coords : TEXCOORD0) : COLOR0
{
    float2 pixelSize = 1.0f / GridSize;
    float a = tex2D(InputSamplerA, coords).r;
    float b = tex2D(InputSamplerB, coords).g;

    // 计算拉普拉斯A
    float lapA = -1 * a;
    lapA += 0.2f * (
        tex2D(InputSamplerA, coords + float2(-pixelSize.x, 0)).r +
        tex2D(InputSamplerA, coords + float2(pixelSize.x, 0)).r +
        tex2D(InputSamplerA, coords + float2(0, -pixelSize.y)).r +
        tex2D(InputSamplerA, coords + float2(0, pixelSize.y)).r
    );
    lapA += 0.05f * (
        tex2D(InputSamplerA, coords + float2(-pixelSize.x, -pixelSize.y)).r +
        tex2D(InputSamplerA, coords + float2(-pixelSize.x, pixelSize.y)).r +
        tex2D(InputSamplerA, coords + float2(pixelSize.x, -pixelSize.y)).r +
        tex2D(InputSamplerA, coords + float2(pixelSize.x, pixelSize.y)).r
    );

    // 计算拉普拉斯B
    float lapB = -1 * b;
    lapB += 0.2f * (
        tex2D(InputSamplerB, coords + float2(-pixelSize.x, 0)).g +
        tex2D(InputSamplerB, coords + float2(pixelSize.x, 0)).g +
        tex2D(InputSamplerB, coords + float2(0, -pixelSize.y)).g +
        tex2D(InputSamplerB, coords + float2(0, pixelSize.y)).g
    );
    lapB += 0.05f * (
        tex2D(InputSamplerB, coords + float2(-pixelSize.x, -pixelSize.y)).g +
        tex2D(InputSamplerB, coords + float2(-pixelSize.x, pixelSize.y)).g +
        tex2D(InputSamplerB, coords + float2(pixelSize.x, -pixelSize.y)).g +
        tex2D(InputSamplerB, coords + float2(pixelSize.x, pixelSize.y)).g
    );

    float newA = a + DT * (DA * lapA - a * b * b + Feed * (1 - a));
    float newB = b + DT * (DB * lapB + a * b * b - (Kill + Feed) * b);

    return float4(newA, newB, 0, 1);
}

technique Technique1
{
    pass Pass1
    {
        PixelShader = compile ps_2_0 PixelShaderFunction();
    }
}

在Monogame中设置RenderTarget、传递Shader参数并渲染四边形，即可完成GPU端计算，600x600网格能轻松达到60FPS以上。

内容的提问来源于stack exchange，提问作者Fractal Salamander