x265中图像质量评估算法SSIM( 结构相似性 Structural Similarity )的实现
x265中的SSIM是在calculateSSIM函数中实现的。源代码如下,其主要的步骤:将像素块以4x4的尺寸进行划分,分别计算其基本计算单元,包括∑∑a(i,j)\sum\sum a(i,j)∑∑a(i,j)、∑∑b(i,j)\sum\sum b(i,j)∑∑b(i,j)、∑∑[b(i,j)2+a(i,j)2]\sum\sum [b(i,j)^2+a(i,j)^2]∑∑[b(i,j)2+a..
x265中的SSIM是在calculateSSIM函数中实现的。源代码如下,其主要的步骤:
- 将像素块以4x4的尺寸进行划分,分别计算其基本计算单元,包括 ∑ ∑ a ( i , j ) \sum\sum a(i,j) ∑∑a(i,j)、 ∑ ∑ b ( i , j ) \sum\sum b(i,j) ∑∑b(i,j)、 ∑ ∑ [ b ( i , j ) 2 + a ( i , j ) 2 ] \sum\sum [b(i,j)^2+a(i,j)^2] ∑∑[b(i,j)2+a(i,j)2]、 ∑ ∑ a ( i , j ) b ( i , j ) \sum\sum a(i,j) b(i,j) ∑∑a(i,j)b(i,j)。其中a、b分别表示将要比较的像素块a和b。该过程在ssim_4x4x2_core函数中完成。
- 然后将基本计算单元进行组合,得到以8x8为尺寸的计算单元。由ssim_end_4函数完成。
- 最后在ssim_end_4函数中,将计算单元送进ssim_end_1函数,由该函数完成最后的组合得到一个8x8块的ssim值。
/* Function to calculate SSIM for each row */
//SSIM is done for each row in blocks of 4x4
static float calculateSSIM(pixel *pix1, intptr_t stride1, pixel *pix2, intptr_t stride2, uint32_t width, uint32_t height, void *buf, uint32_t& cnt)
{
uint32_t z = 0;
float ssim = 0.0;
int(*sum0)[4] = (int(*)[4])buf; //指向长度为4的数组的指针
int(*sum1)[4] = sum0 + (width >> 2) + 3;
width >>= 2;
height >>= 2; //除以4,因为SSIM的计算是以4x4为基本单位的
for (uint32_t y = 1; y < height; y++) //处理列
{
for (; z <= y; z++) //处理列,保证当前有两行处于处理过程中
{
std::swap(sum0, sum1);
for (uint32_t x = 0; x < width; x += 2) //处理行
primitives.ssim_4x4x2_core(&pix1[4 * (x + (z * stride1))], stride1, &pix2[4 * (x + (z * stride2))], stride2, &sum0[x]);
//计算两个4x4块的基本单元,并将结果存储在sum0或sum1中
}
for (uint32_t x = 0; x < width - 1; x += 4)//计算当前两行中的所有4x4块
ssim += primitives.ssim_end_4(sum0 + x, sum1 + x, X265_MIN(4, width - x - 1));
//将上述计算的基本单元进行组合得到ssim值,装换为以8x8为基本单元
}
cnt = (height - 1) * (width - 1);
return ssim;
}
然后看ssim_4x4x2_core函数:
/* structural similarity metric */
//结构的相似性度量
static void ssim_4x4x2_core(const pixel* pix1, intptr_t stride1, const pixel* pix2, intptr_t stride2, int sums[2][4])
{
for (int z = 0; z < 2; z++)
{
uint32_t s1 = 0, s2 = 0, ss = 0, s12 = 0;
for (int y = 0; y < 4; y++)
{
for (int x = 0; x < 4; x++)
{
int a = pix1[x + y * stride1];
int b = pix2[x + y * stride2];
s1 += a;
s2 += b;
ss += a * a;
ss += b * b;
s12 += a * b;
}
}
sums[z][0] = s1;
sums[z][1] = s2;
sums[z][2] = ss;
sums[z][3] = s12;
pix1 += 4;
pix2 += 4;
}
}
这里就非常明确
- s1表示 ∑ ∑ a ( i , j ) \sum\sum a(i,j) ∑∑a(i,j).
- s2表示 ∑ ∑ b ( i , j ) \sum\sum b(i,j) ∑∑b(i,j).
- ss表示 ∑ ∑ [ b ( i , j ) 2 + a ( i , j ) 2 ] \sum\sum [b(i,j)^2+a(i,j)^2] ∑∑[b(i,j)2+a(i,j)2].
- s12表示 ∑ ∑ a ( i , j ) b ( i , j ) \sum\sum a(i,j) b(i,j) ∑∑a(i,j)b(i,j).
并且该函数中计算了一行中两个4x4块,分别将结果存储在sums[0][n]、sums[1][n]中。
然后是ssim_end_4函数,在执行该函数的过程中必定有两行的4x4块已经完成了基本单元的计算,也就是当前处理的块的尺寸是8 x width (像素块原宽)。
//将原本以4x4为基本单位扩展到8x8为基本单位
static float ssim_end_4(int sum0[5][4], int sum1[5][4], int width)
{
float ssim = 0.0;
for (int i = 0; i < width; i++)
{
ssim += ssim_end_1(sum0[i][0] + sum0[i + 1][0] + sum1[i][0] + sum1[i + 1][0],
sum0[i][1] + sum0[i + 1][1] + sum1[i][1] + sum1[i + 1][1],
sum0[i][2] + sum0[i + 1][2] + sum1[i][2] + sum1[i + 1][2],
sum0[i][3] + sum0[i + 1][3] + sum1[i][3] + sum1[i + 1][3]); //以8x8为基本单位
}
return ssim;
}
其中的8x8块是存在重叠的。
最后是ssim_end_1函数。
static float ssim_end_1(int s1, int s2, int ss, int s12)
{
/* Maximum value for 10-bit is: ss*64 = (2^10-1)^2*16*4*64 = 4286582784, which will overflow in some cases.
* s1*s1, s2*s2, and s1*s2 also obtain this value for edge cases: ((2^10-1)*16*4)^2 = 4286582784.
* Maximum value for 9-bit is: ss*64 = (2^9-1)^2*16*4*64 = 1069551616, which will not overflow. */
#if HIGH_BIT_DEPTH
X265_CHECK((X265_DEPTH == 10) || (X265_DEPTH == 12), "ssim invalid depth\n");
#define type float
static const float ssim_c1 = (float)(.01 * .01 * PIXEL_MAX * PIXEL_MAX * 64);
static const float ssim_c2 = (float)(.03 * .03 * PIXEL_MAX * PIXEL_MAX * 64 * 63);
#else
X265_CHECK(X265_DEPTH == 8, "ssim invalid depth\n");
#define type int
static const int ssim_c1 = (int)(.01 * .01 * PIXEL_MAX * PIXEL_MAX * 64 + .5); //C1 =
static const int ssim_c2 = (int)(.03 * .03 * PIXEL_MAX * PIXEL_MAX * 64 * 63 + .5); //k1 = 0.01; k2 = 0.03; L = PIXEL_MAX
#endif
type fs1 = (type)s1; //a
type fs2 = (type)s2; //b
type fss = (type)ss; //a * a + b * b
type fs12 = (type)s12; //a * b
type vars = (type)(fss * 64 - fs1 * fs1 - fs2 * fs2); //
type covar = (type)(fs12 * 64 - fs1 * fs2);//cov(x,y)
return (float)(2 * fs1 * fs2 + ssim_c1) * (float)(2 * covar + ssim_c2)
/ ((float)(fs1 * fs1 + fs2 * fs2 + ssim_c1) * (float)(vars + ssim_c2));
//分子分母第一项结果乘以系数N--64; 第二项结果乘以系数N*(N-1)--64*63
#undef type
#undef PIXEL_MAX
}
这里就需要去了解SSIM的计算式了,中间过程这里就不去论述了,直接给出表达式:
S
S
I
M
(
a
,
b
)
=
(
2
μ
a
μ
b
+
C
1
)
(
2
σ
a
b
+
C
2
)
(
μ
a
2
+
μ
b
2
+
C
1
)
(
σ
a
2
+
σ
b
2
+
C
2
)
SSIM(a,b)=\frac{(2\mu_a\mu_b+C_1)(2\sigma_{ab}+C_2)}{({\mu_a}^2+{\mu_b}^2+C_1)({\sigma_a}^2+{\sigma_b}^2+C_2)}
SSIM(a,b)=(μa2+μb2+C1)(σa2+σb2+C2)(2μaμb+C1)(2σab+C2)
其中
C
1
C_1
C1取
(
K
1
∗
L
)
2
(K1*L)^2
(K1∗L)2;
C
2
C_2
C2取
(
K
2
∗
L
)
2
(K2*L)^2
(K2∗L)2。K1=0.01, K2=0.03, L=PIXEL_MAX (是像素值的动态范围)。
按照代码中的计算方法就可以得到该表达式,但其中多个地方存在 ( ∗ 64 ) (*64) (∗64)和 ( ∗ 63 ) (*63) (∗63)。是因为计算过程中没有加入系数N和(N-1)。
比如均值、方差、协方差的计算公式如下。
μ
a
=
1
64
∑
∑
a
(
i
,
j
)
\mu_a=\frac{1}{64}\sum\sum a(i,j)
μa=641∑∑a(i,j)
σ
a
2
=
1
64
−
1
∑
∑
(
a
(
i
,
j
)
−
μ
a
)
2
{\sigma_a}^2=\frac{1}{64-1}\sum\sum (a(i,j)-\mu_a)^2
σa2=64−11∑∑(a(i,j)−μa)2
σ
a
b
=
1
64
−
1
∑
∑
(
a
(
i
,
j
)
−
μ
a
)
(
b
(
i
,
j
)
−
μ
b
)
\sigma_{ab}=\frac{1}{64-1}\sum\sum(a(i,j)-\mu_a)(b(i,j)-\mu_b)
σab=64−11∑∑(a(i,j)−μa)(b(i,j)−μb)
要用到的公式就是这些,将其展开得到的就是代码中的表现形式。但分子分母中的各项都会多出一个系数,可以将部分提取出来:分子分母第一项可以提取出64、第二项可以提取出64*63。可以消除。但最后
2
μ
a
μ
b
2\mu_a\mu_b
2μaμb和
μ
a
2
+
μ
b
2
{\mu_a}^2+{\mu_b}^2
μa2+μb2都多出一个系数64无法消除。可能这一系数是作为亮度系数的权值保留的。
![Logo](https://devpress.csdnimg.cn/79de2bf0b7994defa4242ef90d5513fa.jpg)
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)