SSE图像算法优化系列十七：多个图像处理中常用函数的SSE实现。

SSE图像算法优化系列⼗七：多个图像处理中常⽤函数的SSE实现。

在做图像处理的SSE优化时，也会经常遇到⼀些⼩的过程、数值优化等代码，本⽂分享⼀些个⼈收藏或实现的代码⽚段给⼤家。

⼀、快速求对数运算
对数运算在图像处理中也是个经常会遇到的过程，特备是在⼀些数据压缩和空间转换时常常会⽤到，⽽且是个⽐较耗时的函数，标准的SSE库⾥并没有提供该函数的实现，如果需要⾼精度的SSE版本，⽹络上已经有了，参考：，这个的精度和标准库的精度基本⼀致了，稍作整理后的代码如下：
// 对数函数的SSE实现，⾼精度版
inline __m128 _mm_log_ps(__m128 x)
{
static const __declspec(align(16)) int _ps_min_norm_pos[4] = { 0x00800000, 0x00800000, 0x00800000, 0x00800000 };
static const __declspec(align(16)) int _ps_inv_mant_mask[4] = { ~0x7f800000, ~0x7f800000, ~0x7f800000, ~0x7f800000 };
static const __declspec(align(16)) int _pi32_0x7f[4] = { 0x7f, 0x7f, 0x7f, 0x7f };
static const __declspec(align(16)) float _ps_1[4] = { 1.0f, 1.0f, 1.0f, 1.0f };
static const __declspec(align(16)) float _ps_0p5[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
static const __declspec(align(16)) float _ps_sqrthf[4] = { 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f };
static const __declspec(align(16)) float _ps_log_p0[4] = { 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f };
static const __declspec(align(16)) float _ps_log_p1[4] = { -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f };
static const __declspec(align(16)) float _ps_log_p2[4] = { 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f };
static const __declspec(align(16)) float _ps_log_p3[4] = { -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f };
static const __declspec(align(16)) float _ps_log_p4[4] = { 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f };
static const __declspec(align(16)) float _ps_log_p5[4] = { -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f };
static const __declspec(align(16)) float _ps_log_p6[4] = { 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f };
static const __declspec(align(16)) float _ps_log_p7[4] = { -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f };
static const __declspec(align(16)) float _ps_log_p8[4] = { 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f };
static const __declspec(align(16)) float _ps_log_q1[4] = { -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f };
static const __declspec(align(16)) float _ps_log_q2[4] = { 0.693359375f, 0.693359375f, 0.693359375f, 0.693359375f };
__m128 one = *(__m128*)_ps_1;
__m128 invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
/* cut off denormalized stuff */
x = _mm_max_ps(x, *(__m128*)_ps_min_norm_pos);
__m128i emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
/* keep only the fractional part */
x = _mm_and_ps(x, *(__m128*)_ps_inv_mant_mask);
x = _mm_or_ps(x, _mm_set1_ps(0.5f));
emm0 = _mm_sub_epi32(emm0, *(__m128i *)_pi32_0x7f);
__m128 e = _mm_cvtepi32_ps(emm0);
e = _mm_add_ps(e, one);
__m128 mask = _mm_cmplt_ps(x, *(__m128*)_ps_sqrthf);
__m128 tmp = _mm_and_ps(x, mask);
x = _mm_sub_ps(x, one);
e = _mm_sub_ps(e, _mm_and_ps(one, mask));
x = _mm_add_ps(x, tmp);
__m128 z = _mm_mul_ps(x, x);
__m128 y = *(__m128*)_ps_log_p0;
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p1);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p2);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p3);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p4);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p5);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p6);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p7);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p8);
y = _mm_mul_ps(y, x);
y = _mm_mul_ps(y, z);
tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q1);
y = _mm_add_ps(y, tmp);
tmp = _mm_mul_ps(z, *(__m128*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q2);
x = _mm_add_ps(x, y);
x = _mm_add_ps(x, tmp);
x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN
return x;
}
看上去有⼀⼤堆代码，不过实测这个的速度越是标准库（本⽂是指启动增强指令集选项设置为：未设置，设计上编译器在此种情况下会⾃动设置为SSE2增强，这可以从反编译logf函数看到，因此，这⾥的速度⽐较还不是和纯Fpu实现的⽐较）的2倍，如果稍微降低点精度，⽐如_ps_log_p5到_ps_log_p8之间的代码，还能提⾼点速度。

另外，在很多场合我们还可以使⽤另外⼀种低精度的log函数，其C代码如下所⽰：
//https:///questions/9411823/fast-log2float-x-implementation-c
inline float IM_Flog(float val)
{
union
{
float val;
int x;
} u = { val };
float log_2 = (float)(((u.x >> 23) & 255) - 128);
u.x &= ~(255 << 23);
u.x += (127 << 23);
log_2 += ((-0.34484843f) * u.val + 2.02466578f) * u.val - 0.67487759f;
return log_2 * 0.69314718f;
}
这个函数⼤概有⼩数点后2位精度。

上述代码⼤约也是标准函数的2倍速度左右。

但是上述函数是可以向量化的，我们来尝试实现。

我们⾸先来看联合体，其实这个东西就是两个东西占同⼀个内存空间，然后外部⽤不同的规则去读取他，在SSE⾥，有着丰富的cast函数，他也是⼲这个事情的，⽐如这⾥的联合体就可以⽤_mm_castps_si128来转换，⽽实际上这个Intrinsic并不会产⽣任何的汇编语句。

那么后⾯的那些移位、或运算、⾮运算、加减乘除之类的就是直接翻译了，毫⽆难处，完整的代码如下所⽰：
inline __m128 _mm_flog_ps(__m128 x)
{
__m128i I = _mm_castps_si128(x);
__m128 log_2 = _mm_cvtepi32_ps(_mm_sub_epi32(_mm_and_si128(_mm_srli_epi32(I, 23), _mm_set1_epi32(255)), _mm_set1_epi32(128)));
I = _mm_and_si128(I, _mm_set1_epi32(-2139095041)); // 255 << 23
I = _mm_add_epi32(I, _mm_set1_epi32(1065353216)); // 127 << 23
__m128 F = _mm_castsi128_ps(I);
__m128 T = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(-0.34484843f), F), _mm_set1_ps(2.02466578f));
T = _mm_sub_ps(_mm_mul_ps(T, F), _mm_set1_ps(0.67487759f));
return _mm_mul_ps(_mm_add_ps(log_2, T), _mm_set1_ps(0.69314718f));
}
经过实测，这个速度可以达到标准库的7到8倍的优势。

⼆、快速求幂运算
⼀般图像编程中有log出现的地⽅就会有exp出现，因此exp的优化也尤为重要，同样在中也有exp的优化（还有sin,cos的SSE优化语句呢），我这⾥就不贴那个的代码了，我们同样关注下⽤联合体实现的近似快速算法，其C代码如下所⽰：
inline float IM_Fexp(float Y)
{
union
{
double Value;
int X[2];
} V;
V.X[1] = (int)(Y * 1512775 + 1072632447 + 0.5F);
V.X[0] = 0;
return (float)V.Value;
}
测试这个和标准的exp库函数速度居然差不多，不晓得为啥，但我们来试下他的SSE优化版本了。

　V.X[1] = (int)(Y * 1512775 + 1072632447 + 0.5F);这句话没啥难度，直接翻译就可以了，注意⼏个强制类型转化就可以了，如下所⽰：
__m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps(1512775)), _mm_set1_ps(1072632447)));
由于我们想⼀次性处理4个float类型的数据，因此也就需要4个union的空间，这样就需要2个__m128i变量来保存数据，每个XMM寄存器的数据应该分别为：
T1 0 T0 0 + T3 0 T2 0 （⾼位----》低位）
这个可以使⽤unpack来实现，具体如下：
__m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
__m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
最后我们认为__m128i⾥的数据是double数据，直接⼀个cast就可以了，然后因为我们只需要单精度的数据，再使⽤_mm_cvtpd_ps将double转换为float类型，注意这个时候还需要将他们连接再⼀起形成⼀个完整的__m128变量，最终的代码如下：
inline __m128 _mm_fexp_ps(__m128 Y)
{
__m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps(1512775)), _mm_set1_ps(1072632447)));
__m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
__m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
return _mm_movelh_ps(_mm_cvtpd_ps(_mm_castsi128_pd(TL)), _mm_cvtpd_ps(_mm_castsi128_pd(TH)));
}
实测这个的提速⼤概有10倍。

如果要求double的exp，其SSE代码你会了吗？
三、pow函数的优化。

⼀种常⽤的近似算法如下所⽰：
inline float IM_Fpow(float a, float b)
{
union
{
double Value;
int X[2];
} V;
V.X[1] = (int)(b * (V.X[1] - 1072632447) + 1072632447);
V.X[0] = 0;
return (float)V.Value;
}
和exp很类似，留给有兴趣的⼈⾃⼰实现。

四：两个求倒数函数的优化误区
SSE提供了连个快速求倒数的函数，_mm_rcp_ps，_mm_rsqrt_ps，他们都是近似值，只有12bit的精度，如果想通过他们得到精确的倒数值，需要⽜顿 - 拉弗森⽅法，⽐如利⽤_mm_rcp_ps 求精确倒数的代码如下：
__forceinline __m128 _mm_prcp_ps(__m128 a)
{
__m128 rcp = _mm_rcp_ps(a); // 此函数只有12bit的精度.
return _mm_sub_ps(_mm_add_ps(rcp, rcp), _mm_mul_ps(a, _mm_mul_ps(rcp, rcp))); // x1 = x0 * (2 - d * x0) = 2 * x0 - d * x0 * x0，使⽤⽜顿 - 拉弗森⽅法这种⽅法可以提⾼精度到23bit }
但是实测这个还不如直接⽤_mm_div_ps的速度，即使是下⾯的函数：
__forceinline __m128 _mm_fdiv_ps(__m128 a, __m128 b)
{
return _mm_mul_ps(a, _mm_rcp_ps(b));
}
似乎速度也不够好，⽽且精度还低了。

特别低，如果使⽤_mm_rcp_ps和_mm_rsqrt_ps联合求近似sqrt，即如下代码，速度好像还慢了，真搞不明⽩为什么。

__forceinline __m128 _mm_fsqrt_ps(__m128 a)
{
return _mm_rcp_ps(_mm_rsqrt_ps(a));
}
五：避免除数为0时的获得⽆效效果
在SSE指令中，没有提供整数的除法指令，不知道这是为什么，所以整数除法⼀般只能借⽤浮点版本的指令，同时，除法存在的⼀个问题就是如果除数为0，可能会触发异常，不过SSE在这种情况下不会抛出异常，但是我们应该避免，避免的⽅式有很多，⽐如判断如果除数为0，就做特殊处理，或者如果除数为0就除以⼀个很⼩的数，不过⼤部分的需求是，除数为0，则返回0，此时就可以使⽤下⾯的SSE指令代替_mm_div_ps:
// 四个浮点数的除法a/b，如果b中某个分量为0，则对应位置返回0值
inline __m128 _mm_divz_ps(__m128 a, __m128 b)
{
__m128 Mask = _mm_cmpeq_ps(b, _mm_setzero_ps());
return _mm_blendv_ps(_mm_div_ps(a, b), _mm_setzero_ps(), Mask);
}
即先把除数和0进⾏⽐较，然后在把_mm_div_ps的返回值中，除数为0的部分⽤0代替，当然，这会带来⼀定的性能下降。

实际上，利⽤位运算，上述代码还可以稍作优化如下：
inline __m128 _mm_divz_ps(__m128 a, __m128 b)
{
return _mm_and_ps(_mm_div_ps(a, b), _mm_cmpneq_ps(b, _mm_setzero_ps()));
}
六、将4个32位整数转换为字节数并保存
很多情况下，我们的算法计算需要将字节类型扩展到32位，计算完成后再保存的字节数据中，这个时候使⽤SSE的话是没有直接的指令的，不过SSE4提供了⼀条_mm_cvtsi128_si32指令，可以将XMM寄存器的4个int32数据转换为4个字节数据并保存到⼀个普通的寄存器中，因此，我们只要调⽤⼏个合适的pack语句就可以实现这个功能了，如下所⽰：
// 将4个32位整形变量数据打包到4个字节数据中
inline void _mm_storesi128_4char(unsigned char *Dest, __m128i P)
{
__m128i T = _mm_packs_epi32(P, P);
*((int *)Dest) = _mm_cvtsi128_si32(_mm_packus_epi16(T, T));
}
七、读取12个字节数到⼀个XMM寄存器中
XMM寄存器是16个字节⼤⼩的，⽽且SSE的很多计算是以4的整数倍字节位单位进⾏的，但是在图像处理中，70%情况下处理的是彩⾊的24位图像，即⼀个像素占⽤3个字节，如果直接使⽤load指令载⼊数据，⼀次性可载⼊5加1/3个像素，这对算法的处理是很不⽅便的，⼀般状况下都是加载4个像素，即12个字节，然后扩展成16个字节（给每个像素增加⼀个Alpha值），我们当然可以直接使⽤load加载16个字节，然后每次跳过12个字节在进⾏load加载，但是其实也可以使⽤下⾯的加载12个字节的函数：
// 从指针p处加载12个字节数据到XMM寄存器中，寄存器最⾼32位清0
inline __m128i _mm_loadu_epi96(const __m128i * p)
{
return _mm_unpacklo_epi64(_mm_loadl_epi64(p), _mm_cvtsi32_si128(((int *)p)[2]));
}
保存当然也可以只保存XMM的低12位：
// 将寄存器Q的低位12个字节数据写⼊到指针P中。

inline void _mm_storeu_epi96(__m128i *P, __m128i Q)
{
_mm_storel_epi64(P, Q);
((int *)P)[2] = _mm_cvtsi128_si32(_mm_srli_si128(Q, 8));
}
不过实际测试，可能还是直接使⽤_mm_loadu_si128和_mm_storeu_si128快点，但是要注意循环的结束为⽌。

⼋、整除255
整除255在图像处理是⾮常⾮常常见的操作，前⾯说了，SSE⾥没有整数除法的指令，如果转换到浮点在除那就慢的死了，⼀般情况下如果要求精度不⾼可以使⽤右移8位实现，如果⾮要精确值可以使⽤如下的C代码：
// 计算整数整除255的四舍五⼊结果。

inline int IM_Div255(int V)
{
return (((V >> 8) + V + 1) >> 8); // 似乎V可以是负数
}
翻译为SSE为：
// 返回16位⽆符号整形数据整除255后四舍五⼊的结果： x = ((x + 1) + (x >> 8)) >> 8:
inline __m128i _mm_div255_epu16(__m128i x)
{
return _mm_srli_epi16(_mm_adds_epu16(_mm_adds_epu16(x, _mm_set1_epi16(1)), _mm_srli_epi16(x, 8)), 8);
}
九、求XMM寄存器内所有元素的累加值
这也是个常见的需求，我们可能把某个结果重复的结果保存在寄存器中，最后结束时在把寄存器中的每个元素想加，你当然可以通过访问__m128i变量的内部的元素实现，但是据说这样会降低循环内的优化，⼀种⽅式是直接⽤SSE指令实现，⽐如对8个有符号的short类型的相加代码如下所⽰：
// 8个有符号的16位的数据相加的和。

// https:///questions/31382209/computing-the-inner-product-of-vectors-with-allowed-scalar-values-0-1-and-2-usi/31382878#31382878
inline int _mm_hsum_epi16(__m128i V) // V7 V6 V5 V4 V3 V2 V1 V0
{
// V = _mm_unpacklo_epi16(_mm_hadd_epi16(V, _mm_setzero_si128()), _mm_setzero_si128()); 也可以⽤这句，_mm_hadd_epi16似乎对计算结果超出32768能获得正确结果
__m128i T = _mm_madd_epi16(V, _mm_set1_epi16(1)); // V7+V6 V5+V4 V3+V2 V1+V0
T = _mm_add_epi32(T, _mm_srli_si128(T, 8)); // V7+V6+V3+V2 V5+V4+V1+V0 0 0
T = _mm_add_epi32(T, _mm_srli_si128(T, 4)); // V7+V6+V3+V2+V5+V4+V1+V0 V5+V4+V1+V0 0 0
return _mm_cvtsi128_si32(T); // 提取低位
}
对于epi32或者ps类型也是使⽤类似的过程的。

10、求16个字节的最⼩值
⽐如我们要求⼀个字节序列的最⼩值，我们肯定会使⽤_mm_min_epi8这样的函数保存每隔16个字节的最⼩值，这样最终我们得到16个字节的⼀个XMM寄存器，整个序列的最⼩值肯定在这个16个字节⾥⾯，这个时候我们可以巧妙的借⽤下⾯的SSE语句得到这16个字节的最⼩值：
// 求16个字节数据的最⼩值, 只能针对字节数据。

inline int _mm_hmin_epu8(__m128i a)
{
__m128i L = _mm_unpacklo_epi8(a, _mm_setzero_si128());
__m128i H = _mm_unpackhi_epi8(a, _mm_setzero_si128());
return _mm_extract_epi16(_mm_min_epu16(_mm_minpos_epu16(L), _mm_minpos_epu16(H)), 0);
}
SSE3提供了_mm_minpos_epu16函数，他能获取8个⽆符号数的的最⼩值及其最⼩值的索引，放置在寄存器的低16和低32位，我们把字节数据扩展到16位，然后在通过两次⽐较就可以获得相应的最⼩值了。

那如果是求最⼤值呢，可惜SSE没有提供_mm_maxpos_epu16函数，但是也⽆妨，稍微修改下上⾯的代码就可以了，如下所⽰：
// 求16个字节数据的最⼤值, 只能针对字节数据。

inline int _mm_hmax_epu8(__m128i a)
{
__m128i b = _mm_subs_epu8(_mm_set1_epi8(255), a);
__m128i L = _mm_unpacklo_epi8(b, _mm_setzero_si128());
__m128i H = _mm_unpackhi_epi8(b, _mm_setzero_si128());
return 255 - _mm_extract_epi16(_mm_min_epu16(_mm_minpos_epu16(L), _mm_minpos_epu16(H)), 0);
}
⼗⼀、其他⼀些优化技巧
在以及等⽹站上还有很多参考的资料，希望⼤家⾃⼰去学习下。