目标

在 OpenCV 中的视频输入与相似度测量教程中，我已经介绍了 PSNR 和 SSIM 方法来检查两幅图像之间的相似度。正如你所见，执行过程需要相当长的时间，特别是 SSIM。但是，如果 CPU 版 OpenCV 实现的性能不尽如人意，而你的系统恰好有一个 NVIDIA CUDA GPU 设备，那么并非无计可施。你可以尝试将自己的算法移植或编写到显卡上。

本教程将帮助你很好地掌握如何通过使用 OpenCV 的 GPU 模块进行编码。作为先决条件，你应该已经了解如何处理 core、highgui 和 imgproc 模块。因此，我们的主要目标是：

与 CPU 相比有什么不同？
为 PSNR 和 SSIM 创建 GPU 代码
优化代码以获得最大性能

源代码

你也可以在 OpenCV 源代码库的 samples/cpp/tutorial_code/gpu/gpu-basics-similarity/gpu-basics-similarity 目录下找到源代码和视频文件，或者从这里下载。完整的源代码相当长（由于通过命令行参数和性能测量控制应用程序）。因此，为了避免在本节中充斥这些内容，这里只包含函数本身。

PSNR 返回一个浮点数，如果两个输入相似，则该数值在 30 到 50 之间（越高越好）。

double getPSNR(const Mat& I1, const Mat& I2)
{
    Mat s1;
    absdiff(I1, I2, s1); // |I1 - I2|
s1.convertTo(s1, CV_32F); // 无法对 8 位整数进行平方运算
s1 = s1.mul(s1); // |I1 - I2|^2
 
    Scalar s = sum(s1); // 对每个通道求和
 
    double sse = s.val[0] + s.val[1] + s.val[2]; // 对通道求和
 
    double mse = sse /(double)(I1.channels() * I1.total());
 
    // 对于非常小的 SSE，加上 epsilon 来近似无穷大的 PSNR (~361 dB)
    if( sse <= 1e-10) mse+= DBL_EPSILON;
 
    double psnr = 10.0*log10((255*255)/mse);
    return psnr;
}

double getPSNR_CUDA(const Mat& I1, const Mat& I2)
{
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
gI1.upload(I1);
gI2.upload(I2);
 
gI1.convertTo(t1, CV_32F);
gI2.convertTo(t2, CV_32F);
 
cuda::absdiff(t1.reshape(1), t2.reshape(1), gs);
cuda::multiply(gs, gs, gs);
 
    Scalar s = cuda::sum(gs);
    double sse = s.val[0] + s.val[1] + s.val[2];
 
    double mse = sse /(double)(gI1.channels() * I1.total());
 
    // 对于非常小的 SSE，加上 epsilon 来近似无穷大的 PSNR (~361 dB)
    if( sse <= 1e-10) mse+= DBL_EPSILON;
 
    double psnr = 10.0*log10((255*255)/mse);
    return psnr;
}

struct BufferPSNR // 优化的 CUDA 版本
{ // CUDA 上的数据分配成本很高。使用缓冲区来解决：一次分配，稍后重用。
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
    cuda::GpuMat buf;
};

double getPSNR_CUDA_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)
{
b.gI1.upload(I1);
b.gI2.upload(I2);
 
b.gI1.convertTo(b.t1, CV_32F);
b.gI2.convertTo(b.t2, CV_32F);
 
cuda::absdiff(b.t1.reshape(1), b.t2.reshape(1), b.gs);
cuda::multiply(b.gs, b.gs, b.gs);
 
    double sse = cuda::sum(b.gs, b.buf)[0];
 
    double mse = sse /(double)(I1.channels() * I1.total());
 
    // 对于非常小的 SSE，加上 epsilon 来近似无穷大的 PSNR (~361 dB)
    if (sse <= 1e-10) mse += DBL_EPSILON;
 
    double psnr = 10.0*log10((255*255)/mse);
    return psnr;
}

SSIM 返回图像的 MSSIM。这也是一个零到一之间的浮点数（越高越好），但是我们每个通道都有一个。因此，我们返回一个*Scalar* OpenCV 数据结构。

Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
    const double C1 = 6.5025, C2 = 58.5225;
/***************************** 初始化 **********************************/
    int d = CV_32F;
 
    Mat I1, I2;
i1.convertTo(I1, d); // 无法对 1 字节大小的值进行计算
i2.convertTo(I2, d);
 
    Mat I2_2 = I2.mul(I2); // I2^2
    Mat I1_2 = I1.mul(I1); // I1^2
    Mat I1_I2 = I1.mul(I2); // I1 * I2
 
/*************************** 初始化结束 **********************************/
 
    Mat mu1, mu2; // 预计算
    GaussianBlur(I1, mu1, Size(11, 11), 1.5);
    GaussianBlur(I2, mu2, Size(11, 11), 1.5);
 
    Mat mu1_2 = mu1.mul(mu1);
    Mat mu2_2 = mu2.mul(mu2);
    Mat mu1_mu2 = mu1.mul(mu2);
 
    Mat sigma1_2, sigma2_2, sigma12;
 
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
sigma1_2 -= mu1_2;
 
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
sigma2_2 -= mu2_2;
 
    GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
sigma12 -= mu1_mu2;
 
    Mat t1, t2, t3;
 
t1 = 2 * mu1_mu2 + C1;
t2 = 2 * sigma12 + C2;
t3 = t1.mul(t2); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
t1 = mu1_2 + mu2_2 + C1;
t2 = sigma1_2 + sigma2_2 + C2;
t1 = t1.mul(t2); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
    Mat ssim_map;
    divide(t3, t1, ssim_map); // ssim_map = t3./t1;
 
    Scalar mssim = mean( ssim_map ); // mssim = ssim map 的平均值
    return mssim;
}

Scalar getMSSIM_CUDA( const Mat& i1, const Mat& i2)
{
    const float C1 = 6.5025f, C2 = 58.5225f;
/***************************** 初始化 **********************************/
    cuda::GpuMat gI1, gI2, gs1, tmp1,tmp2;
 
gI1.upload(i1);
gI2.upload(i2);
 
gI1.convertTo(tmp1, CV_MAKE_TYPE(CV_32F, gI1.channels()));
gI2.convertTo(tmp2, CV_MAKE_TYPE(CV_32F, gI2.channels()));
 
vector<cuda::GpuMat> vI1, vI2;
cuda::split(tmp1, vI1);
cuda::split(tmp2, vI2);
    Scalar mssim;
 
    Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(vI2[0].type(), -1, Size(11, 11), 1.5);
 
    for( int i = 0; i < gI1.channels(); ++i )
    {
        cuda::GpuMat I2_2, I1_2, I1_I2;
 
cuda::multiply(vI2[i], vI2[i], I2_2); // I2^2
cuda::multiply(vI1[i], vI1[i], I1_2); // I1^2
cuda::multiply(vI1[i], vI2[i], I1_I2); // I1 * I2
 
/*************************** 初始化结束 **********************************/
        cuda::GpuMat mu1, mu2; // 预计算
gauss->apply(vI1[i], mu1);
gauss->apply(vI2[i], mu2);
 
        cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
cuda::multiply(mu1, mu1, mu1_2);
cuda::multiply(mu2, mu2, mu2_2);
cuda::multiply(mu1, mu2, mu1_mu2);
 
        cuda::GpuMat sigma1_2, sigma2_2, sigma12;
 
gauss->apply(I1_2, sigma1_2);
cuda::subtract(sigma1_2, mu1_2, sigma1_2); // sigma1_2 -= mu1_2;
 
gauss->apply(I2_2, sigma2_2);
cuda::subtract(sigma2_2, mu2_2, sigma2_2); // sigma2_2 -= mu2_2;
 
gauss->apply(I1_I2, sigma12);
cuda::subtract(sigma12, mu1_mu2, sigma12); // sigma12 -= mu1_mu2;
 
        cuda::GpuMat t1, t2, t3;
 
mu1_mu2.convertTo(t1, -1, 2, C1); // t1 = 2 * mu1_mu2 + C1;
sigma12.convertTo(t2, -1, 2, C2); // t2 = 2 * sigma12 + C2;
cuda::multiply(t1, t2, t3); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
cuda::addWeighted(mu1_2, 1.0, mu2_2, 1.0, C1, t1); // t1 = mu1_2 + mu2_2 + C1;
cuda::addWeighted(sigma1_2, 1.0, sigma2_2, 1.0, C2, t2); // t2 = sigma1_2 + sigma2_2 + C2;
cuda::multiply(t1, t2, t1); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
 
        cuda::GpuMat ssim_map;
cuda::divide(t3, t1, ssim_map); // ssim_map = t3./t1;
 
        Scalar s = cuda::sum(ssim_map);
mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols);
 
    }
    return mssim;
}

struct BufferMSSIM // 优化的 CUDA 版本
{ // CUDA 上的数据分配成本很高。使用缓冲区来解决：一次分配，稍后重用。
    cuda::GpuMat gI1, gI2, gs, t1,t2;
 
    cuda::GpuMat I1_2, I2_2, I1_I2;
vector<cuda::GpuMat> vI1, vI2;
 
    cuda::GpuMat mu1, mu2;
    cuda::GpuMat mu1_2, mu2_2, mu1_mu2;
 
    cuda::GpuMat sigma1_2, sigma2_2, sigma12;
    cuda::GpuMat t3;
 
    cuda::GpuMat ssim_map;
 
    cuda::GpuMat buf;
};

Scalar getMSSIM_CUDA_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b)
{
    const float C1 = 6.5025f, C2 = 58.5225f;
/***************************** 初始化 **********************************/
 
b.gI1.upload(i1);
b.gI2.upload(i2);
 
    cuda::Stream stream;
 
b.gI1.convertTo(b.t1, CV_32F, stream);
b.gI2.convertTo(b.t2, CV_32F, stream);
 
cuda::split(b.t1, b.vI1, stream);
cuda::split(b.t2, b.vI2, stream);
    Scalar mssim;
 
    Ptr<cuda::Filter> gauss = cuda::createGaussianFilter(b.vI1[0].type(), -1, Size(11, 11), 1.5);
 
    for( int i = 0; i < b.gI1.channels(); ++i )
    {
cuda::multiply(b.vI2[i], b.vI2[i], b.I2_2, 1, -1, stream); // I2^2
cuda::multiply(b.vI1[i], b.vI1[i], b.I1_2, 1, -1, stream); // I1^2
cuda::multiply(b.vI1[i], b.vI2[i], b.I1_I2, 1, -1, stream); // I1 * I2
 
gauss->apply(b.vI1[i], b.mu1, stream);
gauss->apply(b.vI2[i], b.mu2, stream);
 
cuda::multiply(b.mu1, b.mu1, b.mu1_2, 1, -1, stream);
cuda::multiply(b.mu2, b.mu2, b.mu2_2, 1, -1, stream);
cuda::multiply(b.mu1, b.mu2, b.mu1_mu2, 1, -1, stream);
 
gauss->apply(b.I1_2, b.sigma1_2, stream);
cuda::subtract(b.sigma1_2, b.mu1_2, b.sigma1_2, cuda::GpuMat(), -1, stream);
        //b.sigma1_2 -= b.mu1_2; - 这将导致额外的读取操作
 
gauss->apply(b.I2_2, b.sigma2_2, stream);
cuda::subtract(b.sigma2_2, b.mu2_2, b.sigma2_2, cuda::GpuMat(), -1, stream);
        //b.sigma2_2 -= b.mu2_2;
 
gauss->apply(b.I1_I2, b.sigma12, stream);
cuda::subtract(b.sigma12, b.mu1_mu2, b.sigma12, cuda::GpuMat(), -1, stream);
        //b.sigma12 -= b.mu1_mu2;
 
        //此处调用 operator*(Scalar, Mat) 也会导致额外的读取操作
cuda::multiply(b.mu1_mu2, 2, b.t1, 1, -1, stream); //b.t1 = 2 * b.mu1_mu2 + C1;
cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
cuda::multiply(b.sigma12, 2, b.t2, 1, -1, stream); //b.t2 = 2 * b.sigma12 + C2;
cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -12, stream);
 
cuda::multiply(b.t1, b.t2, b.t3, 1, -1, stream); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
 
cuda::add(b.mu1_2, b.mu2_2, b.t1, cuda::GpuMat(), -1, stream);
cuda::add(b.t1, C1, b.t1, cuda::GpuMat(), -1, stream);
 
cuda::add(b.sigma1_2, b.sigma2_2, b.t2, cuda::GpuMat(), -1, stream);
cuda::add(b.t2, C2, b.t2, cuda::GpuMat(), -1, stream);
 
 
cuda::multiply(b.t1, b.t2, b.t1, 1, -1, stream); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
cuda::divide(b.t3, b.t1, b.ssim_map, 1, -1, stream); // ssim_map = t3./t1;
 
stream.waitForCompletion();
 
        Scalar s = cuda::sum(b.ssim_map, b.buf);
mssim.val[i] = s.val[0] / (b.ssim_map.rows * b.ssim_map.cols);
 
    }
    return mssim;
}

如何实现？ - GPU

如上所述，我们为每个操作都有三种类型的函数。一种用于 CPU，两种用于 GPU。我为 GPU 创建两种函数的原因是说明，简单地将 CPU 代码移植到 GPU 实际上会使其变慢。如果你想获得性能提升，你需要记住一些规则，我稍后会详细介绍。

GPU 模块的开发使其尽可能地与其 CPU 对应物相似。这使得移植过程更加容易。在编写任何代码之前，你需要做的第一件事是将 GPU 模块链接到你的项目，并包含该模块的头文件。GPU 的所有函数和数据结构都位于*cv* 命名空间的*gpu* 子命名空间中。你可以通过*use namespace* 关键字将其添加到默认命名空间，或者到处明确地用 cv:: 标记它以避免混淆。我将采用后一种方式。

#include <opencv2/gpu.hpp> // GPU 结构和方法

GPU 是“图形处理单元”的缩写。它最初用于渲染图形场景。这些场景在某种程度上建立在大量数据的基础上。然而，这些数据并非都以顺序方式相互依赖，而是可以并行处理。因此，GPU 将包含多个较小的处理单元。这些不是最先进的处理器，一对一测试会输给 CPU。然而，它的优势在于数量。近年来，人们越来越倾向于利用 GPU 的大规模并行能力来渲染非图形场景。这催生了通用图形处理单元（GPGPU）计算。

GPU 有自己的内存。当你在 OpenCV 中使用*Mat* 对象从硬盘读取数据时，这些数据会存储在系统的内存中。CPU 直接（通过其缓存）对这些数据进行操作，但是 GPU 不能。它必须将计算所需的信息从系统内存传输到自己的内存。这通过上传过程完成，并且非常耗时。最后，结果必须下载回你的系统内存，以便你的 CPU 可以看到并使用它。不建议移植小型函数到 GPU，因为上传/下载时间将大于并行执行所获得的增益。

*Mat* 对象仅存储在系统内存（或 CPU 缓存）中。要将 OpenCV 矩阵放到 GPU 上，你需要使用其 GPU 对应物 cv::cuda::GpuMat。它的工作方式类似于 Mat，但有 2D 限制，并且其函数没有引用返回（不能将 GPU 引用与 CPU 引用混合）。要将 Mat 对象上传到 GPU，你需要在创建类实例后调用 upload 函数。要下载，你可以简单地赋值给 Mat 对象或使用 download 函数。

Mat I1; // 主内存项 - 例如使用 imread 将图像读入
gpu::GpuMat gI; // GPU 矩阵 - 目前为空
gI1.upload(I1); // 将数据从系统内存上传到 GPU 内存
 
I1 = gI1; // 下载，gI1.download(I1) 也能工作

一旦你的数据上传到 GPU 内存，你就可以调用 OpenCV 的 GPU 功能。大多数函数保留与 CPU 上的相同名称，不同之处在于它们只接受*GpuMat* 输入。

需要记住的另一件事是，并非所有通道数都能在 GPU 上实现高效算法。通常，我发现 GPU 图像的输入图像需要是单通道或四通道的，并且项类型为 char 或 float。不支持 double，抱歉。为某些函数传递其他类型的对象将导致异常抛出，并在错误输出中显示错误消息。文档在大多数地方详细说明了接受的输入类型。如果你有三通道图像作为输入，你可以做两件事：要么添加一个新通道（并使用 char 元素），要么分割图像并为每个通道调用函数。第一种方法并不推荐，因为它浪费内存。

对于某些函数，其中元素的位置（相邻元素）无关紧要，快速的解决方案是将其重塑为单通道图像。PSNR 实现就是这种情况，对于*absdiff* 方法，邻居的值并不重要。但是，对于*GaussianBlur*，这不是一个选项，因此需要使用 split 方法。有了这些知识，你就可以编写一个 GPU 可用的代码（像我的 GPU 代码一样）并运行它。你可能会惊讶地发现它可能比你的 CPU 实现慢。

优化

原因在于你忽略了内存分配和数据传输的成本。而在 GPU 上，这些成本非常高。另一种优化可能性是借助 cv::cuda::Stream 来引入异步 OpenCV GPU 调用。

GPU 上的内存分配是相当大的。因此，如果可能，应尽可能少地分配新内存。如果你创建一个打算多次调用的函数，最好只在第一次调用时分配该函数的所有局部参数。为此，你可以创建一个包含你将使用的所有局部变量的数据结构。例如，对于 PSNR，这些是：
struct BufferPSNR // 优化的 GPU 版本

{ // GPU 上的数据分配成本非常高。使用缓冲区来解决：一次分配，稍后重用。

gpu::GpuMat gI1, gI2, gs, t1,t2;

gpu::GpuMat buf;

};

然后，在主程序中创建此实例：
BufferPSNR bufferPSNR;

最后，每次调用函数时都将其传递给函数：
double getPSNR_GPU_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b)

现在，你可以像这样访问这些局部参数：*b.gI1*、*b.buf* 等等。GpuMat 仅在新矩阵大小与上一次调用不同时才会重新分配自身。
避免不必要的数据传输。任何小型数据传输在 GPU 上都会变得显著。因此，如果可能，请尽可能在原地进行所有计算（换句话说，不要创建新的内存对象 - 出于先前原因的解释）。例如，虽然算术运算可能更容易用单行公式表达，但它会更慢。在 SSIM 的一个点上，我需要计算：
b.t1 = 2 * b.mu1_mu2 + C1;

虽然上述调用会成功，但请注意其中存在隐藏的数据传输。在执行加法之前，它需要在某个地方存储乘法的结果。因此，它会在后台创建一个局部矩阵，然后加上*C1*值，最后将其赋值给*t1*。为了避免这种情况，我们使用 gpu 函数而不是算术运算符：
gpu::multiply(b.mu1_mu2, 2, b.t1); //b.t1 = 2 * b.mu1_mu2 + C1;

gpu::add(b.t1, C1, b.t1);
使用异步调用（cv::cuda::Stream）。默认情况下，每当你调用 GPU 函数时，它会等待调用完成，然后返回结果。但是，可以进行异步调用，这意味着它会调用操作执行，进行耗时的算法数据分配，然后立即返回。现在你可以调用另一个函数，如果你愿意的话。对于 MSSIM，这是一个小的优化点。在我们默认的实现中，我们将图像分割成通道，并为每个通道调用 GPU 函数。通过 stream 可以实现一定程度的并行化。通过使用 stream，我们可以在 GPU 已经执行某个方法时进行数据分配和上传操作。例如，我们需要上传两张图像。我们将它们一个接一个地入队，然后调用处理它们的函数。函数将等待上传完成，但是在此过程中，它会为下一个要执行的函数进行输出缓冲区分配。
gpu::Stream stream;

stream.enqueueConvert(b.gI1, b.t1, CV_32F); // 上传

gpu::split(b.t1, b.vI1, stream); // 方法（将 stream 作为最后一个参数传递）。

gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2

CV_32F
#define CV_32F
定义位于 interface.h:78

结果与结论

在一台配备低端 NVIDIA GT220M 的 Intel P8700 笔记本 CPU 上，性能数据如下：

PSNR CPU 时间（平均 10 次运行）：41.4122 毫秒。结果为：19.2506
PSNR GPU 时间（平均 10 次运行）：158.977 毫秒。结果为：19.2506
GPU 初始调用：31.3418 毫秒。结果为：19.2506
GPU 优化版本（/ 10 次运行）：24.8171 毫秒。结果为：19.2506
 
MSSIM CPU 时间（平均 10 次运行）：484.343 毫秒。结果为 B0.890964 G0.903845 R0.936934
MSSIM GPU 时间（平均 10 次运行）：745.105 毫秒。结果为 B0.89922 G0.909051 R0.968223
MSSIM GPU 初始调用：357.746 毫秒。结果为 B0.890964 G0.903845 R0.936934
MSSIM GPU 优化版本（/ 10 次运行）：203.091 毫秒。结果为 B0.890964 G0.903845 R0.936934

在这两种情况下，我们都实现了与 CPU 实现相比近 100% 的性能提升。这可能正是你的应用程序所需的那种改进。你可以在 YouTube 这里观看此演示的运行时实例。