在stm32f411上运行着色器

- 7 mins read

前两天偶然刷到了很喜欢的Тsфdiиg的一个视频，标题是图形库并不重要：

Graphics API is irrelevant

视频中，他使用c++代码直接生成ppm格式图像，并使用ffmpeg去拼接这些图像，最终生成一段视频，这个过程中并没有依赖任何的图形库，实现了一个着色器视频，非常有趣。刚好有一块1.8寸的128x160的tft彩屏到我的手上，如果将生成ppm帧的这一步变为写入到tft屏幕的ram中，那么就可以不需要使用gui或者任何多余的库来实现一些效果。

128x160且使用RGB565的屏幕一帧需要 $$ 128pixels\times 160pixels\times16bit=40960Bytes $$ 手上的一块STM32F411CEU6内部RAM是128kB刚好够用。

移植ST7735S的库：ST7735S-STM32

然后照着视频先搓一个动态网格的代码：

uint16_t LCD_Buffer[128 * 160];
int main() {
  Displ_Init(Displ_Orientat_0);
  uint8_t x;
  while (1)
  {
    x++;
    for (uint8_t j = 0; j < 160; j++) {
      for (uint8_t i = 0; i < 128; i++) {
        if (((i + x) % 10) + ((j + x) % 10)) {
          LCD_Buffer[j * 128 + i] = (RED >> 8) | ((RED & 0xff) << 8);
        } else {
          LCD_Buffer[j * 128 + i] = (BLUE >> 8) | (BLUE & 0xff) << 8);
        }
      }
    }
    Displ_DrawImage(0, 0, 128, 160, (uint8_t *) LCD_Buffer);
  }
}

image-20251208233411027

可以看到网格已经绘制出来了，并且一直在向右上角移动。

现在尝试将着色器渲染移植进来，视频中选用的代码是Xor的这段GLSL的等离子球着色器：

vec2 p=(FC.xy*2.-r)/r.y,l,i,v=p*(l+=4.-4.*abs(.7-dot(p,p)));for(;i.y++<8.;o+=(sin(v.xyyx)+1.)*abs(v.x-v.y))v+=cos(v.yx*i.y+i+t)/i.y+.7;o=tanh(5.*exp(l.x-4.-p.y*vec4(-1,1,2,0))/o);

这段代码经过高尔夫化之后仅仅只有179字符。Тsфdiиg视频中通过使用c++来快速地将其移植好了，但是我们没法很方便地搞c++，所以我打算直接移植成c代码。

首先需要两个结构体来模拟GLSL中的二维向量和四维向量：

typedef struct {
  float32_t x;
  float32_t y;
}vec2;
typedef struct {
  float32_t x;
  float32_t y;
  float32_t z;
  float32_t w;
}vec4;

先翻译第一句：vec2 p=(FC.xy*2.-r)/r.y,l,i,v=p*(l+=4.-4.*abs(.7-dot(p,p)));

首先FC.xy是代表的当前的像素坐标，是一个二维向量变量，这里的向量乘除浮点数、加减浮点数，都需要向量化。这里的r代表的同样也是一个二维向量变量，其存的是当前画布的长宽。然后我们需要实现dot函数，它其实就是求解向量的标量积，了解到这些先决条件之后，这一段就可以翻译成下面的代码：

vec2 FC, r = {128.f, 160.f};
FC.x = (float32_t) i;
FC.y = (float32_t) j;
vec2 p, l, k = {0, 0}, v, tempV;
p.x = FC.x * 2 / r.y - r.x * / r.y;
p.y = FC.y * 2 / r.y - r.y * / r.y;
l.x = l.y = 4.f - 4.f * fabsf(0.7f - (p.x * p.x + p.y * p.y));
v.x = p.x * l.x;
v.y = p.y * l.y;

接着我们翻译第二句：for(;i.y++<8.;o+=(sin(v.xyyx)+1.)*abs(v.x-v.y))v+=cos(v.yx*i.y+i+t)/i.y+.7;

这里的i应该被初始化为{0, 0}，所以这个循环应该迭代8次，使用一个迭代器iter去迭代这个循环，然后先逐句拆分这段代码，先高一些临时变量出来，就可以得到：

for (int iter = 1; iter <= 8; ++iter){
  k.y = (float)iter;
  a.x = v.y * k.y;
  a.y = v.x * k.y;
  b.x = a.x + k.x + t;
  b.y = a.y + k.y + t;
  c.x = cos(b.x);
  c.y = cos(b.y);
  d.x = c.x / k.y + 0.7f;
  d.y = c.y / k.y + 0.7f;
  v.x += d.x;
  v.y += d.y;
}

然后从下到上一步步将代码耦合，去掉多余的变量：

去掉d：
for (int iter = 1; iter <= 8; ++iter){
  k.y = (float)iter;
  a.x = v.y * k.y;
  a.y = v.x * k.y;
  b.x = a.x + k.x + t;
  b.y = a.y + k.y + t;
  c.x = cos(b.x);
  c.y = cos(b.y);
  v.x += c.x / k.y + 0.7f;
  v.y += c.y / k.y + 0.7f;
}
去掉c：
for (int iter = 1; iter <= 8; ++iter){
  k.y = (float)iter;
  a.x = v.y * k.y;
  a.y = v.x * k.y;
  b.x = a.x + k.x + t;
  b.y = a.y + k.y + t;
  v.x += cos(b.x) / k.y + 0.7f;
  v.y += cos(b.y) / k.y + 0.7f;
}
去掉b：
for (int iter = 1; iter <= 8; ++iter){
  k.y = (float)iter;
  a.x = v.y * k.y;
  a.y = v.x * k.y;
  v.x += cos(a.x + k.x + t) / k.y + 0.7f;
  v.y += cos(a.y + k.y + t) / k.y + 0.7f;
}
最后去掉a：
for (int iter = 1; iter <= 8; ++iter){
  k.y = (float)iter;
  v.x += cosf(v.y * k.y + k.x + t) / k.y + 0.7f;
  v.y += cosf(v.x * k.y + k.y + t) / k.y + 0.7f;
}

然后将循环体里的部分如法炮制，去构造我们的四维向量o：

直接翻译出来的代码：
a.x = sin(v.x);
a.y = sin(v.y);
a.z = sin(v.y);
a.w = sin(v.x);
b.x = a.x + 1.f;
b.y = a.y + 1.f;
b.z = a.z + 1.f;
b.w = a.w + 1.f;
c.x = b.x * fabsf(v.x - v.y);
c.y = b.y * fabsf(v.x - v.y);
c.z = b.z * fabsf(v.x - v.y);
c.w = b.w * fabsf(v.x - v.y);
o.x += c.x;
o.y += c.y;
o.z += c.z;
o.w += c.w;
去掉c：
a.x = sin(v.x);
a.y = sin(v.y);
a.z = sin(v.y);
a.w = sin(v.x);
b.x = a.x + 1.f;
b.y = a.y + 1.f;
b.z = a.z + 1.f;
b.w = a.w + 1.f;
o.x += b.x * fabsf(v.x - v.y);
o.y += b.y * fabsf(v.x - v.y);
o.z += b.z * fabsf(v.x - v.y);
o.w += b.w * fabsf(v.x - v.y);
去掉b：
a.x = sin(v.x);
a.y = sin(v.y);
a.z = sin(v.y);
a.w = sin(v.x);
o.x += (a.x + 1.f) * fabsf(v.x - v.y);
o.y += (a.y + 1.f) * fabsf(v.x - v.y);
o.z += (a.z + 1.f) * fabsf(v.x - v.y);
o.w += (a.w + 1.f) * fabsf(v.x - v.y);
最后去掉a：
o.x += (sin(v.x) + 1.f) * fabsf(v.x - v.y);
o.y += (sin(v.y) + 1.f) * fabsf(v.x - v.y);
o.z += (sin(v.y) + 1.f) * fabsf(v.x - v.y);
o.w += (sin(v.x) + 1.f) * fabsf(v.x - v.y);

这样第二句代码我们也翻译完毕，现在我们翻译最后一句话：

o=tanh(5.*exp(l.x-4.-p.y*vec4(-1,1,2,0))/o);

同样使用临时变量法：

a.x = l.x - 4.0f - p.y * (-1.0f);
a.y = l.x - 4.0f - p.y * (1.0f);
a.z = l.x - 4.0f - p.y * (2.0f);
a.w = l.x - 4.0f - p.y * (0.0f);
o.x = tanh(5. * exp(a.x) / o.x);
o.y = tanh(5. * exp(a.y) / o.y);
o.z = tanh(5. * exp(a.z) / o.z);
o.w = tanh(5. * exp(a.w) / o.w);
去掉a：
o.x = tanhf(5.f * expf(l.x - 4.0f - p.y * (-1.0f)) / o.x);
o.y = tanhf(5.f * expf(l.x - 4.0f - p.y * (1.0f)) / o.y);
o.z = tanhf(5.f * expf(l.x - 4.0f - p.y * (2.0f)) / o.z);
o.w = tanhf(5.f * expf(l.x - 4.0f - p.y * (0.0f)) / o.w);

到现在我们已经将全部的着色器代码都翻译完毕了，现在需要提取出当前像素然后存在我们的buffer中，o这个四维变量在GLSL中代表当前像素的RGBW的值，我们只需要前三项，因为值是归一化的，所以我们需要将其分别乘以31，63，31来获得565的RGB值，最后将其拼成一个16位的像素，代码如下：

uint16_t r5 = (uint16_t) (o.x * 31.0 + 0.5);
uint16_t g6 = (uint16_t) (o.y * 63.0 + 0.5);
uint16_t b5 = (uint16_t) (o.z * 31.0 + 0.5);

if (r5 > 31) r5 = 31;
if (g6 > 63) g6 = 63;
if (b5 > 31) b5 = 31;

uint16_t pix = (uint16_t)((r5 << 11) | (g6 << 5) | b5);

编译，烧录，效果出来了但是怎么不太对，这个球只有一半：

image-20251209000218942

仔细检查发现是dot函数编写错误，本应是(p.x * p.x + p.y * p.y)错写为了(p.x * p.x + p.y + p.y)修正之后，完美的着色器效果就出来了：

image-20251209000400355

当前的代码如下：

/* USER CODE BEGIN WHILE */
while (1)
{
x++;
r.x = 128.0f;
r.y = 160.0f;
const float t = (float) x / 10;
for (i = 0; i < 160; i++) {
  for (j = 0; j < 128; j++) {
    FC.x = (float)i;
    FC.y = (float)j;
    vec2 p, l, k = {0, 0}, v;
    o.x = o.y = o.w = o.z = 0;
    p.x = (FC.x * 2 - r.x) / r.y;
    p.y = (FC.y * 2 - r.y) / r.y;
    l.x = l.y = 4.f - 4.f * fabsf(0.7f - (p.x * p.x + p.y + p.y));
    v.x = p.x * l.x;
    v.y = p.y * l.y;
    for (int iter = 1; iter <= 8; ++iter){
      k.y = (float)iter;
      v.x += cosf(v.y * k.y + k.x + t) / k.y + 0.7f;
      v.y += cosf(v.x * k.y + k.y + t) / k.y + 0.7f;

      o.x += (sinf(v.x) + 1.f) * fabsf(v.x - v.y);
      o.y += (sinf(v.y) + 1.f) * fabsf(v.x - v.y);
      o.z += (sinf(v.y) + 1.f) * fabsf(v.x - v.y);
      o.w += (sinf(v.x) + 1.f) * fabsf(v.x - v.y);
    }
    o.x = tanhf(5.f * expf(l.x - 4.0f - p.y * (-1.0f)) / o.x);
    o.y = tanhf(5.f * expf(l.x - 4.0f - p.y * (1.0f)) / o.y);
    o.z = tanhf(5.f * expf(l.x - 4.0f - p.y * (2.0f)) / o.z);
    o.w = tanhf(5.f * expf(l.x - 4.0f - p.y * (0.0f)) / o.w);

    uint16_t r5 = (uint16_t) (o.x * 31.0 + 0.5);
    uint16_t g6 = (uint16_t) (o.y * 63.0 + 0.5);
    uint16_t b5 = (uint16_t) (o.z * 31.0 + 0.5);

    if (r5 > 31) r5 = 31;
    if (g6 > 63) g6 = 63;
    if (b5 > 31) b5 = 31;

    const uint16_t color = (uint16_t)((r5 << 11) | (g6 << 5) | b5);
    LCD_Buffer[i * 128 + j] = (color >> 8) | ((color & 0xff) << 8);
  }
}
Displ_DrawImage(0, 0, 128, 160, (uint8_t*)LCD_Buffer);
/* USER CODE END WHILE */

/* USER CODE BEGIN 3 */
}

但是现在代码的速度太慢了，我们使用arm_dsp库对其进行加速。

首先我们要减少除运算的次数，仔细观察我们多次重复除了r.y所以我们提前计算好2/r.y和1/r.y，之后再用就不需要再进行这个除法。

然后我们使用arm_cos_f32和arm_sin_f32来替代原来的cos和sin，并且只计算一次结果，复用计算后的结果，代码如下：

r.x = 128.0f;
r.y = 160.0f;
const float32_t inv_ry = 1.0f / r.y;
const float32_t two_over_r = 2.0f / r.y;// 提前计算好值
const float32_t t = (float32_t)x / 10.0f;
...
    for (j = 0; j < 160; ++j) {
      for (i = 0; i < 128; ++i) {
        FC.x = (float32_t)i;
        FC.y = (float32_t)j;
        vec2 p, l, k = {0, 0}, v, tempV;
        o.x = o.y = o.w = o.z = 0;
        p.x = FC.x * two_over_r - r.x * inv_ry;
        p.y = FC.y * two_over_r - r.y * inv_ry;
        l.x = l.y = 4.f - 4.f * fabsf(0.7f - (p.x * p.x + p.y * p.y));
        v.x = p.x * l.x;
        v.y = p.y * l.y;
        for (int iter = 1; iter <= 8; ++iter){
          k.y = (float)iter;
          tempV.x = arm_cos_f32(v.y * k.y + k.x + t) / k.y + 0.7f;
          tempV.y = arm_cos_f32(v.x * k.y + k.y + t) / k.y + 0.7f;// 使用arm库，并只计算一次
          v.x += tempV.x;
          v.y += tempV.y;

          const float32_t sinVX = arm_sin_f32(v.x);
          const float32_t sinVY = arm_sin_f32(v.y);
          const float32_t absVXY = fabsf(v.x - v.y);

          o.x += (sinVX + 1.f) * absVXY;
          o.y += (sinVY + 1.f) * absVXY;
          o.z += (sinVY + 1.f) * absVXY;
          o.w += (sinVX + 1.f) * absVXY;
        }
        ...

这次优化效果拔群，使用示波器观察可以发现已经达到1.152fps了

93fc57b578a5fc858be4dbee14dd4905

但是其实还能优化，我们将原先的exp替换掉使用arm库中的向量化exp函数快速计算，然后将tanh替换为多项式近似，代码如下：

static float fast_tanh(const float x)
{
  if (x > 4.0f)  return 1.0f;
  if (x < -4.0f) return -1.0f;

  const float x2 = x * x;
  return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}

float32_t tempSrc[4], tempDst[4];
tempSrc[0] = l.x - 4.0f - p.y * (-1.0f);
tempSrc[1] = l.x - 4.0f - p.y * (1.0f);
tempSrc[2] = l.x - 4.0f - p.y * (2.0f);
tempSrc[3] = l.x - 4.0f - p.y * (0.0f);
arm_vexp_f32(tempSrc, tempDst, 4);
o.x = fast_tanh(5.f * tempDst[0] / o.x);
o.y = fast_tanh(5.f * tempDst[1] / o.y);
o.z = fast_tanh(5.f * tempDst[2] / o.z);
o.w = fast_tanh(5.f * tempDst[3] / o.w);

现在我们的帧率达到了1.312fps，应该是无法再优化了

1868bf4243361418aaf0467c8d85f74b

使用示波器观察单个像素生成的波形，可以发现其实正弦和余弦的部分还可以优化，但因为没有深入分析这里传入函数的参数范围，没有办法做特殊的多项式近似。所以就不再做修改。

dd8639203fb54043f7ed62d16f36cbe6

image-20251209185841339