OpenCV 实现图像结构相似度算法 (SSIM 算法)

　　看完原论文后（英语不太好，只看懂了个大概），大致上明白了它是做什么以及如何实现的等，于是决定写一篇博客，所以该篇文章简单介绍一下 SSIM（有能力的话，看原论文或维基比我来介绍更好），并给出我用 opencv 在 C++ 中的 SSIM 算法实现。

　　首先什么是 SSIM 算法，该算法主要用于检测两张尺寸相同的图像的相似度，但注意到论文标题的中的 Structural，所以实际上它主要通过分别比较两个图像的亮度(l)、对比度(c)、结构(s)，然后对这三个要素加权并乘积表示，而在论文中这三个要素用下面公式来表示：

$$l(x, y) = frac{(2mu _xmu _y + C_1)}{(mu _{x}^{2} + mu _{y}^{2} + C_1)}$$

$$c(x, y) = frac{(2sigma _{x}sigma _{y} + C_2)}{(sigma _{x}^{2} + sigma _{y}^{2} + C_2)}$$

$$s(x, y) = frac{(sigma _{xy} + C_3)}{(sigma _{x} sigma _{y} + C_3)}$$

　　这里 $mu _x$ 为均值，$sigma _{x}$ 为方差，$sigma _{xy}$ 表示协方差。这里 $C_1$、$C_2$、$C_3$ 是为了避免当分母为 0 时造成的不稳定问题（所以写算法的时候可以放心，一定不会出现除 0 的情况）

　　而 SSIM 的一般方程为：

$$ssim(x, y) = [l(x, y)^alpha cdot c(x, y)^eta cdot s(x, y)^gamma ]$$

　　这里一般 $alpha$，$eta$，$gamma$ 取 $1$，并且令 $C_3 = frac{C_2}{2}$，这样就得到简化的 SSIM 公式：

$$ssim(x, y) = frac{(2mu _xmu _y + C_1)(sigma _{xy} + C_2)}{(mu _{x}^{2} + mu _{y}^{2} + C_1)(sigma _{x}^{2} + sigma _{y}^{2} + C_2)}$$

　　论文中还指出该公式满足下面三个条件：

　　　　1.对称性，即 $ssim(x, y) = ssim(y, x)$.

　　　　2.有界性，即 $ssim(x, y) <= 1$.

　　　　3.有唯一最大值，即 $ssim(x, y) <= 1$，这里当且仅当 $x = y$ 时取等并且这里还有一点是说：两个矩阵的每个矩阵元素都一一对应相等.

　　上面三个条件很容易理解，对称性保证了当两张图交换位置带入公式时，显然测量指标不应该发生变化，有界性保证了归一化到 0 ~ 1（话说我觉得应该不可能出现小于 0 的情况），评估可以更直观的按百分比来说明，最后一个条件保证了我们使用相同的图像作为输入，一定能得到相等的结果。

　　实际上公式是十分显然的，我们只要分别求出两张图像的均值以及方差，再对两张图像求协方差就可以带入到 SSIM 公式进行计算，这里计算结果其实得到的是一张图像，它显示了两张图的混叠，借用一下python 以及两张 opencv 的示例图：

　　(pic5.png)

　　(pic6.png)

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv
from skimage import data, img_as_float
from skimage.metrics import structural_similarity as ssim


img1 = cv.imread('C:/opencv/samples/data/pic5.png', cv.CV_8U)
img2 = cv.imread('C:/opencv/samples/data/pic6.png', cv.CV_8U)


def mse(x, y):
    return np.linalg.norm(x - y)


fig, axes = plt.subplots()

mse_none = mse(img1, img2)
mssim, s = ssim(img1, img2, data_range=img1.max() - img1.min(), full=True)

label = 'MSE: {:.2f}, SSIM: {:.2f}'

axes.imshow(s, cmap=plt.cm.gray, vmin=0, vmax=1)
axes.set_xlabel(label.format(mse_none, mssim))
axes.set_title('SSIM image')

plt.tight_layout()
plt.show()

　　会看到 ssim 输出：

　　而评估算法的性能指标公式是：

$$mssim(x, y) = frac{1}{M}sum_{j = 1}^{M} ssim(x_j, y_j)$$

　　也就是像素取均值。

　　该算法的优点在于对图像的结构信息敏感，所以相似图像的识别率很高，但缺点也有：

　　1.如果将原图经过旋转（宽度和高度不同的话试试旋转 180 度）、平移等变换或加噪点后再和原图进行比较就显得吃力了：

　　2.算法只能处理相同形状的两个图像，但可以多通道。

　　最后，给出我的 C++ 实现（不支持多通道图像，以后或许会考虑实现多通道）：

/********************************************************************************************************************************************************
* 该部分主要功能实现 SSIM 算法
*
* 关于 SSIM 算法参考文献：https://ece.uwaterloo.ca/~z70wang/publications/ssim.pdf
* 公式：$$ssim(x, y) = frac{(2mu _xmu _y + C_1)(sigma _{xy} + C_2)}{(mu _{x}^{2} + mu _{y}^{2} + C_1)(sigma _{x}^{2} + sigma _{y}^{2} + C_2)}$$
*
* 参数：
*   im1: 图像 1
*   im2: 图像 2
*   window: 滑动窗口大小，用于卷积滤波
*   k1: 可调节常数，默认 k1 = 0.01
*   k2: 可调节常数，默认 k2 = 0.03
*   L: 单通道灰度图像像素值范围，默认 L = 255.0
*
* 返回值：
*   相似度指标 (类型：double)
*********************************************************************************************************************************************************/
float ssim(Mat im1, Mat im2, int window = 7, float k1 = 0.01f, float k2 = 0.03f, float L = 255.f)
{
    CV_Assert(im1.size() == im2.size());
    
    int ndim = im1.dims;
    float NP = std::powf(window, ndim);
    float cov_norm = NP / (NP - 1);
    float C1 = (k1 * L) * (k1 * L);
    float C2 = (k2 * L) * (k2 * L);
    
    Mat ux, uy;
    Mat uxx = im1.mul(im1);
    Mat uyy = im2.mul(im2);
    Mat uxy = im1.mul(im2);
    
    blur(im1, ux, Size(window, window), Point(-1, -1));
    blur(im2, uy, Size(window, window), Point(-1, -1));
    
    blur(uxx, uxx, Size(window, window), Point(-1, -1));
    blur(uyy, uyy, Size(window, window), Point(-1, -1));
    blur(uxy, uxy, Size(window, window), Point(-1, -1));
    
    Mat ux_sq = ux.mul(ux);
    Mat uy_sq = uy.mul(uy);
    Mat uxy_m = ux.mul(uy);

    Mat vx = cov_norm * (uxx - ux_sq);
    Mat vy = cov_norm * (uyy - uy_sq);
    Mat vxy = cov_norm * (uxy - uxy_m);

    Mat A1 = 2 * uxy_m;
    Mat A2 = 2 * vxy;
    Mat B1 = ux_sq + uy_sq;
    Mat B2 = vx + vy;
    
    Mat ssim_map = (A1 + C1).mul(A2 + C2) / (B1 + C1).mul(B2 + C2);
    
    Scalar mssim = mean(ssim_map);
    ssim_map.convertTo(ssim_map, CV_8UC1, 255, 0);
    
    imshow("ssim", ssim_map);

    return mssim[0];
}

　　用下面代码即可测试算法：

#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/videoio.hpp"
#include <iostream>

using namespace std;
using namespace cv;
// ssim()
int main(int argc, const char** argv)
{
    Mat orginal_im;
    if (argc > 1)
        orginal_im = imread(argv[1]);
    else
        orginal_im = imread("C:/opencv/samples/data/lena.jpg");

    Mat frame = imread("C:/opencv/samples/data/lena_tmpl.jpg");
    Mat im1f, im2f;
    orginal_im.convertTo(im1f, CV_32FC1);
    frame.convertTo(im2f, CV_32FC1);

    ostringstream text;
    text << "mssim:" << ssim(im1f, im2f);

    putText(frame, text.str(), Point(30, 30),
        FONT_HERSHEY_COMPLEX_SMALL, 0.8, Scalar(255, 200, 0), 1, 8);
    imshow("Original image", orginal_im);
    imshow("Frame", frame);
    waitKey(0);

    return 0;
}

　　运行结果：

　　参考文献：

　　　　1.https://ece.uwaterloo.ca/~z70wang/publications/ssim.pdf

　　　　2.https://en.wikipedia.org/wiki/Structural_similarity

　　算法实现参考：

　　　　1.https://www.cns.nyu.edu/~lcv/ssim/ssim_index.m

　　　　2.https://github.com/PAHdb/ssim/blob/master/ssim.pro

　　　　3.skimage.metrics.structural_similarity 源码