SoftGLRender源码：SIMD技术

SIMD技术
- 介绍
工作原理
- 使用__m128
- 常用_mm_函数
SoftGLRender的SIMD

SIMD技术

介绍

SIMD（Single Instruction, Multiple Data），意为单指令多数据，是一种让CPU用一条指令同时处理多个数据的技术，广泛用于加速图形处理、物理模拟、音视频编解码、机器学习等。

CPU有以下几种工作方式：

模式	描述
SISD	Single Instruction, Single Data（普通 CPU）
SIMD	Single Instruction, Multiple Data（并行处理数组）
MISD	多指令同数据（罕见）
MIMD	多指令多数据（多核/多线程）

SIMD 是现代 CPU 的标准配置，用来加速数据并行任务，而不是任务并行.

目前，不同平台的技术实现：

平台/厂商	指令集	向量宽度	数据类型支持
Intel/AMD	SSE（128位）	4 × float	`__m128`
Intel/AMD	AVX（256位）	8 × float	`__m256`
ARM	NEON	4–8 × T	各种整型/浮点
Apple	Accelerate	多种	高级封装的 SIMD

工作原理

以普通加法 vs SIMD为例

普通加法处理4个数，代码可能这样：

// 加法处理 4 个数：
for (int i = 0; i < 4; i++) {
    result[i] = a[i] + b[i];
}

需要执行4次加法.

如果使用SIMD技术：

__m128 va = _mm_loadu_ps(a);  // 加载 4 个 float
__m128 vb = _mm_loadu_ps(b);
__m128 vr = _mm_add_ps(va, vb);  // 一次完成 4 个加法

只需要1条指令，一次完成4次浮点数加法.

__m128 是一个 SIMD（Single Instruction, Multiple Data）向量类型，它由 Intel 的 SSE（Streaming SIMD Extensions）指令集提供的一种 128 位数据结构，专门用于在一个 CPU 指令中并行处理 4 个 32 位浮点数（float）.

简单说，__m128是一个装着4个float的128bit浮点数寄存器，一次能进行4个float计算.

__m128 (128 bits) 
 ┌────────┬────────┬────────┬────────┐
 | float3 | float2 | float1 | float0 |   ← 4 个 32-bit 浮点数
 └────────┴────────┴────────┴────────┘

逻辑上，__m128等价于装有4个float的数组

__m128 v = [x, y, z, w];

使用__m128

#include <xmmintrin.h>  // 包含 SSE 的头文件

__m128 a = _mm_set_ps(4.0f, 3.0f, 2.0f, 1.0f);  // 设置向量 [1,2,3,4]
__m128 b = _mm_set1_ps(2.0f);                  // 所有元素都设置为 2.0

__m128 result = _mm_add_ps(a, b);              // 每个元素相加：a[i] + b[i]

相当于执行了：

result = [1+2, 2+2, 3+2, 4+2] = [3, 4, 5, 6]

一次完成4个浮点数加法，很高效.

常用_mm_函数

函数名	含义
`_mm_set_ps`	设置向量（从高到低）
`_mm_add_ps`	4 元素浮点加法
`_mm_sub_ps`	减法
`_mm_mul_ps`	乘法
`_mm_div_ps`	除法
`_mm_loadu_ps`	从内存加载 4 个 float
`_mm_storeu_ps`	把向量写回内存

使用__m128实施SIMD技术时，需要注意:

__m128只能存放float. 如果想存放int，需要使用__m128i；存double，需要使用__m128d；

2）需要include <xmmintrin.h>；

3）编译时，需要开启SSE（现代编译器通常默认开启）.

_mm_load*系列函数

_mm_load*系列函数都是SSE SIMD指令，用于从内存加载数据到SIMD寄存器（e.g. __m128，__m128d）. 用途、数据类型、速度、内存要求不同.

简要对比：

指令	内存要求	速度（通常）	说明
`_mm_load_ps`	16字节对齐	更快（理想情况）	快速加载，要求地址对齐到 16 字节边界
`_mm_loadu_ps`	可以不对齐	稍慢（通用性强）	“Unaligned” 加载，不需要对齐
`_mm_load_ss`	任意对齐	通用	只加载1个 float到最低位，其它为 0
`_mm_load1_ps`	任意对齐	通用	加载 1 个 float 并复制到所有通道
`_mm_load_pd`	16字节对齐	更快	加载 2 个 double 值

例：

使用 _mm_load_ps（对齐）

// 内存必须 16 字节对齐
float data[4] __attribute__((aligned(16))) = {1.0f, 2.0f, 3.0f, 4.0f};
__m128 vec = _mm_load_ps(data);  // OK

若 data 没有对齐到16字节边界，可能会崩溃或性能极差！

使用 _mm_loadu_ps（不对齐）（SSE）

float data[4] = {1.0f, 2.0f, 3.0f, 4.0f};
__m128 vec = _mm_loadu_ps(data);  // 安全，即使 data 没对齐

不要求data对齐. 适用于结构体、动态分配数组等内存布局不确定的情形.

_mm_set_ps函数

_mm_set_ps设置__m128寄存器的指令，从4个float数创建一个128bit SIMD向量.

__m128 _mm_set_ps(float z, float y, float x, float w);

// index-value对应关系：
index: [3] [2] [1] [0]
value:  z   y   x   w

注：_mm_set_ps参数顺序是逆序，参数从右到左放进SIMD 寄存器.

如果要正常的参数顺序，用_mm_setr_ps：

__m128 a = _mm_set_ps(4, 3, 2, 1);   // a = [1, 2, 3, 4]
__m128 b = _mm_setr_ps(1, 2, 3, 4);  // b = [1, 2, 3, 4]

_mm_shuffle_ps与_MM_SHUFFLE

_mm_shuffle_ps, _MM_SHUFFLE：强大SSE指令组合，重排（shuffle）1 or 2个 __m128 SIMD 寄存器中的 float 元素.

_mm_shuffle_ps用于重排浮点向量中的元素：

__m128 _mm_shuffle_ps(__m128 a, __m128 b, const int imm8);

从a中选2个元素（用于输出的低2个元素）
从b中选2个元素（用于输出的高2个元素）
选择哪2个，由imm8（可通过_MM_SHUFFLE生成）控制

_MM_SHUFFLE 构造imm8控制字：

#define _MM_SHUFFLE(z, y, x, w) (((z)<<6) | ((y)<<4) | ((x)<<2) | ((w)))

参数	取自哪里	用于输出哪个位置
`z`	第2个参数 `b`	输出位置 [3]（最高位）
`y`	第2个参数 `b`	输出位置 [2]
`x`	第1个参数 `a`	输出位置 [1]
`w`	第1个参数 `a`	输出位置 [0]（最低位）

例如，

__m128 a = _mm_set_ps(4.0, 3.0, 2.0, 1.0); // [4, 3, 2, 1]
__m128 b = _mm_set_ps(8.0, 7.0, 6.0, 5.0); // [8, 7, 6, 5]

// 重排
__m128 result = _mm_shuffle_ps(a, b, _MM_SHUFFLE(2, 1, 2, 1));

拆解_MM_SHUFFLE(2, 1, 2, 1)：

z: 输出[3] = b[2] = 6
y: 输出[2] = b[1] = 7
x: 输出[1] = a[2] = 2
w: 输出[0] = a[1] = 3

最终结果：result=[3, 2, 7, 6]

_mm_mul_ps乘法

_mm_mul_ps float逐元素乘法.

__m128 a = _mm_set_ps(4.0f, 3.0f, 2.0f, 1.0f); // [1, 2, 3, 4]
__m128 b = _mm_set_ps(0.5f, 1.5f, 2.0f, 4.0f); // [4, 2, 1.5, 0.5]

// 逐元素相乘
__m128 result = _mm_mul_ps(a, b); // result = [4, 4, 4.5, 2]

_mm_fmsub_ps融合乘法

__m128 _mm_fmsub_ps(__m128 a, __m128 b, __m128 c);

等价于

result = (a * b) - c;

不过，这是一次性完成乘法、减法操作，在一次指令中执行完成.

__m128 a = _mm_set_ps(4.0f, 3.0f, 2.0f, 1.0f);  // a = [1, 2, 3, 4]
__m128 b = _mm_set1_ps(2.0f);                  // b = [2, 2, 2, 2]
__m128 c = _mm_set1_ps(1.0f);                  // c = [1, 1, 1, 1]

__m128 result = _mm_fmsub_ps(a, b, c);         // result = (a * b) - c

执行结果：

result = [1*2 - 1, 2*2 - 1, 3*2 - 1, 4*2 - 1] = [1, 3, 5, 7]

SoftGLRender的SIMD

SIMD模块

文件：Base/SIMD.h

//SIMD.h

#ifdef _MSC_VER
#define MM_F32(v, i) v.m128_f32[i]
#else
#define MM_F32(v, i) v[i]
#endif // _MSC_VER

#define PTR_ADD(p) ((size_t)p)

这段代码定义了2个宏，MM_F32 和 PTR_ADD，用于跨平台的SIMD向量访问和指针转换.

宏MM_F32(v, i)，访问SIMD向量中第i个浮点数，要求v是__m128类型.

不同编译器对__m128的访问支持不同：

1）MSVC编译器提供数组成员.m128_f32；
2）GCC/Clang则直接像访问数组一样 v[i].

宏PTR_ADD，将指针强制转换为整数地址，常用于：

1）指针偏移；
2）地址比较或计算哈希值；
3）调试时，内存打印.

e.g. 下面例子将指针ptr转换为size_t后，打印

float* ptr = some_array;
size_t address = PTR_ADD(ptr);
printf("地址是：%zu\n", address);

SIMD加速求重心坐标

软件渲染中，有一个重要步骤是求三角形的重心坐标.

原理参见计算机图形：三角形及重心空间

这段代码通过2个版本计算三角形的重心坐标：
1）SIMD优化的版本，需要开启SOFTGL_SIMD_OPT；
2）普通的glm库方法计算.

// Render/Software/RendererSoft.cpp

// @param vert[0] {vertPos[2].x, vertPos[1].x, vertPos[0].x, 0.f}
//        vert[1] {vertPos[2].y, vertPos[1].y, vertPos[0].y, 0.f}
//        vert[2] {vertPos[0].z, vertPos[1].z, vertPos[2].z, 0.f}
//        vert[3] {vertPos[0].w, vertPos[1].w, vertPos[2].w, 0.f}
// @param v0 A顶点坐标. 
//        v1 B顶点坐标
//        v2 C顶点坐标
// @param p 三角形ABC中, 当前采样坐标
// @param [out] bc重心坐标
bool RendererSoft::barycentric(glm::aligned_vec4 *vert, glm::aligned_vec4 &v0, glm::aligned_vec4 &p,
                               glm::aligned_vec4 &bc) {
#ifdef SOFTGL_SIMD_OPT
  // Ref: https://geometrian.com/programming/tutorials/cross-product/index.php
  __m128 vec0 = _mm_sub_ps(_mm_load_ps(&vert[0].x), _mm_set_ps(0, p.x, v0.x, v0.x));
  __m128 vec1 = _mm_sub_ps(_mm_load_ps(&vert[1].x), _mm_set_ps(0, p.y, v0.y, v0.y));

  __m128 tmp0 = _mm_shuffle_ps(vec0, vec0, _MM_SHUFFLE(3, 0, 2, 1));
  __m128 tmp1 = _mm_shuffle_ps(vec1, vec1, _MM_SHUFFLE(3, 1, 0, 2));
  __m128 tmp2 = _mm_mul_ps(tmp0, vec1);
  __m128 tmp3 = _mm_shuffle_ps(tmp2, tmp2, _MM_SHUFFLE(3, 0, 2, 1));
  __m128 u = _mm_sub_ps(_mm_mul_ps(tmp0, tmp1), tmp3);

  if (std::abs(MM_F32(u, 2)) < FLT_EPSILON) {
    return false;
  }

  u = _mm_div_ps(u, _mm_set1_ps(MM_F32(u, 2)));
  bc = {1.f - (MM_F32(u, 0) + MM_F32(u, 1)), MM_F32(u, 1), MM_F32(u, 0), 0.f};
#else
  glm::vec3 u = glm::cross(glm::vec3(vert[0]) - glm::vec3(v0.x, v0.x, p.x),
                           glm::vec3(vert[1]) - glm::vec3(v0.y, v0.y, p.y));
  if (std::abs(u.z) < FLT_EPSILON) {
    return false;
  }

  u /= u.z;
  bc = {1.f - (u.x + u.y), u.y, u.x, 0.f};
#endif

  if (bc.x < 0 || bc.y < 0 || bc.z < 0) {
    return false;
  }

  return true;
}

为了简化参数记忆，建立有名称的顶点与顶点数组之间的关系.

\(△ABC\)中，顶点对应关系：

A: vertPos[0], v0
B: vertPos[1], v1
C: vertPos[2], v2

于是，

vert[0] = {Cx, Bx, Ax, 0}
vert[1] = {Cy, By, Ay, 0}
vert[2] = {Az, Bz, Cz, 0}
vert[3] = {Aw, Bw, Cw, 0}
v0 = {Ax, Ay, Az, Aw}
v1 = {Bx, By, Bz, Bw}
v2 = {Cx, Cy, Cz, Cw}

第2个版本（普通glm库计算）

先看barycentric第2个版本实现，即没有使用SIMD的版本：

...

  glm::vec3 u = glm::cross(glm::vec3(vert[0]) - glm::vec3(v0.x, v0.x, p.x),
                           glm::vec3(vert[1]) - glm::vec3(v0.y, v0.y, p.y));
  if (std::abs(u.z) < FLT_EPSILON) {
    return false;
  }

  u /= u.z;
  bc = {1.f - (u.x + u.y), u.y, u.x, 0.f};
#endif

  if (bc.x < 0 || bc.y < 0 || bc.z < 0) { // P不在三角形ABC内
    return false;
  }

  return true;

这段代码是不是又回到了之前的重心坐标插值？
参见myrender our_gl.cpp

根据顶点A、B、C坐标，求重心坐标的原理：

s1 = (Cx-Ax, Bx-Ax, Ax-Px)
s2 = (Cy-Ay, By-Ay, Ay-Py)
u = s1 x s2
重心坐标[alpha, beta, gamma]满足P=alpha*A + beta*B + gamma*C, 
其中，
alpha=1-(u.x+u.y)/u.z, 
beta=u.y/u.z, 
gamma=u.x/u.z

综上，我们知道barycentric用于：

1）判断采样坐标P是否为\(△ABC\)内一点（含边）；
2）输出重心坐标\([\alpha, \beta, \gamma]=[1-(u.x+u.y)/u.z,u.y/u.z, u.x/u.z]\)到坐标bc

第1个版本（SIMD优化）

    // vec0 = [Cx, Bx, Ax, 0] - [Ax, Bx, Px, 0]
  __m128 vec0 = _mm_sub_ps(_mm_load_ps(&vert[0].x), _mm_set_ps(0, p.x, v0.x, v0.x));
   // vec1 = [Cy, By, Ay, 0] - [Ay, By, Py, 0]
  __m128 vec1 = _mm_sub_ps(_mm_load_ps(&vert[1].x), _mm_set_ps(0, p.y, v0.y, v0.y));

得到2个__m128类型向量vec0, vec1

下面这段代码，用于求vec0, vec1的叉积：

\[x\times y= \begin{pmatrix} x_1y_2-x_2y_1\\x_2y_0-x_0y_2\\x_0y_1-x_1y_0 \end{pmatrix} \]

下面代码主要用SSE指令求向量叉积，而barycentric中求\(\bm{u}\)代码类似：

// Ref: https://geometrian.com/resources/cross_product/

// vec0=[x0,x1,x2,x3]
// vec1=[y0,y1,y2,y3]
//Method 5: Fewer swizzles, swizzle before subtraction
[[nodiscard]] inline static __m128 cross_product(
	__m128 const& vec0, __m128 const& vec1
) noexcept {
    // tmp0=[vec0[1], vec0[2], vec0[0], vec0[3]]
    //     =[x1,x2,x0,x3]
	__m128 tmp0 = _mm_shuffle_ps( vec0,vec0, _MM_SHUFFLE(3,0,2,1) );

    // tmp1=[vec1[2], vec1[0], vec1[1], vec1[3]]
    //     =[y2,y0,y1,y3]
	__m128 tmp1 = _mm_shuffle_ps( vec1,vec1, _MM_SHUFFLE(3,1,0,2) );

    // tmp2 = tmp0 * vec1
    //      = [x1y0, x2y1, x0y2,x3y3]
	__m128 tmp2 = _mm_mul_ps( tmp0, vec1 );

    // tmp3 = tmp0 * tmp1
    //      = [x1y2, x2y0, x0y1, x3y3]
	__m128 tmp3 = _mm_mul_ps( tmp0, tmp1 );

    // tmp4=[tmp2[1], tmp2[2], tmp2[0], tmp2[3]]
    //     =[x2y1, x0y2, x1y0, x3y3]
	__m128 tmp4 = _mm_shuffle_ps( tmp2,tmp2, _MM_SHUFFLE(3,0,2,1) );

    // tmp3 - tmp4 = [x1y2-x2y1, x2y0-x0y2, x0y1-x1y0, 0]
	return _mm_sub_ps( tmp3, tmp4 );
}

接下来，根据\(\bm{u}\)求出重心坐标bc=(1-(u.x+u.y)/u.z, u.y/u.z, u.x/u.z, 0)

  // |u.z| < FLT_EPSILON
  if (std::abs(MM_F32(u, 2)) < FLT_EPSILON) {
    return false;
  }

  // u = u/u.z
  u = _mm_div_ps(u, _mm_set1_ps(MM_F32(u, 2)));
  // bc = (1-(u.x+u.y), u.y, u.x, 0)
  bc = {1.f - (MM_F32(u, 0) + MM_F32(u, 1)), MM_F32(u, 1), MM_F32(u, 0), 0.f};

其他

SOFTGL_SIMD_OPT

使用SIMD优化版本的宏定义开关SOFTGL_SIMD_OPT，可在根目录CMakeLists.txt中开启.

# enable SIMD
add_definitions("-DSOFTGL_SIMD_OPT")

FLT_EPSILON

FLT_EPSILON 是标准库在<cfloat>/<float.h>中定义的宏，表示单精度浮点数 (float) 的机器 epsilon（增量）.

既是float能表示的精度，也是float能表示的最小增量.

#include <cfloat>  // C++
#include <float.h> // C

#define FLT_EPSILON 1.192092896e-07F  // 典型值（32-bit float）

表示最小的正浮点数ε，使得 1.0 + ε ≠ 1.0（float运算中）

posted @ 2025-05-18 16:54 明明1109 阅读(84) 评论(0) 收藏举报

刷新页面返回顶部

明明1109