Similarity check (PNSR and SSIM) on the GPU的简单翻译，原文地址
Global
在使用OpenCV和相似性度量的视频输入教程中，我已经介绍了PSNR和SSIM方法来检查两幅图像之间的相似性。如您所见，执行过程需要相当长的时间，特别是在SSIM的情况下。但是，如果CPU的OpenCV实现的性能数字不能让您满意，并且您的系统中碰巧有一个NVIDIA CUDA GPU设备，那么一切都不会丢失。您可以尝试为视频卡移植或编写owm算法。
本教程将很好地掌握如何使用OpenCV的GPU模块进行编码。作为先决条件，您应该已经知道如何处理core、highgui和imgproc模块。因此，我们的主要目标是：
-
与CPU相比有什么不同？
-
为PSNR和SSIM创建GPU代码
-
优化代码以获得最佳性能
源码
您还可以在OpenCV源代码库的samples/cpp/tutorial\u 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); // cannot make a square on 8 bits
s1 = s1.mul(s1); // |I1 - I2|^2
Scalar s = sum(s1); // sum elements per channel
double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels
if( sse <= 1e-10) // for small values return zero
return 0;
else
{
double mse =sse /(double)(I1.channels() * I1.total());
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];
if( sse <= 1e-10) // for small values return zero
return 0;
else
{
double mse =sse /(double)(gI1.channels() * I1.total());
double psnr = 10.0*log10((255*255)/mse);
return psnr;
}
}
struct BufferPSNR // Optimized CUDA versions
{ // Data allocations are very expensive on CUDA. Use a buffer to solve: allocate once reuse later.
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];
if( sse <= 1e-10) // for small values return zero
return 0;
else
{
double mse = sse /(double)(I1.channels() * I1.total());
double psnr = 10.0*log10((255*255)/mse);
return psnr;
}
}
SSIM返回图像的MSSIM。这也是一个介于0和1之间的浮点数（越高越好），但是每个通道有一个浮点数。因此，我们返回一个标量OpenCV数据结构：
Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
const double C1 = 6.5025, C2 = 58.5225;
/***************************** INITS **********************************/
int d = CV_32F;
Mat I1, I2;
i1.convertTo(I1, d); // cannot calculate on one byte large values
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
/*************************** END INITS **********************************/
Mat mu1, mu2; // PRELIMINARY COMPUTING
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 = average of ssim map
return mssim;
}
Scalar getMSSIM_CUDA( const Mat& i1, const Mat& i2)
{
const float C1 = 6.5025f, C2 = 58.5225f;
/***************************** INITS **********************************/
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
/*************************** END INITS **********************************/
cuda::GpuMat mu1, mu2; // PRELIMINARY COMPUTING
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 // Optimized CUDA versions
{ // Data allocations are very expensive on CUDA. Use a buffer to solve: allocate once reuse later.
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;
/***************************** INITS **********************************/
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; - This would result in an extra data transfer operation
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;
//here too it would be an extra data transfer due to call of 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;
}
How to do it? -The GPU
如上所述，每个操作有三种类型的函数。一个用于CPU，两个用于GPU。我为GPU做了两个示例的原因是，通常简单地将CPU移植到GPU实际上会使它变慢。如果你想获得一些性能提升，你需要记住一些规则，稍后我将详细介绍这些规则。
GPU模块的开发是为了尽可能地与CPU模块相似。这使得移植过程更容易。在编写任何代码之前，您需要做的第一件事是将GPU模块链接到您的项目，并包含模块的头文件。GPU的所有函数和数据结构都位于cv命名空间的GPU子命名空间中。您可以通过use namespace关键字将其添加到默认名称空间中，或者通过cv::显式地将其标记到任意位置以避免混淆。我晚点再做。
#include <opencv2/gpu.hpp> // GPU structures and methods
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对象的简单赋值或使用下载功能。
Mat I1; // Main memory item - read image into with imread for example
gpu::GpuMat gI; // GPU matrix - for now empty
gI1.upload(I1); // Upload a data from the system memory to the GPU memory
I1 = gI1; // Download, gI1.download(I1) will work too
一旦你的数据在GPU内存中，你可以调用OpenCV的GPU功能。大多数函数的名称与CPU上的相同，不同之处在于它们只接受GpuMat输入。
另一件事要记住的是，不是所有的频道号码，你可以作出高效的算法在GPU上。一般来说，我发现GPU图像的输入图像需要是一个或四个通道的图像，并且元素大小需要是char或float类型。GPU上没有double支持，抱歉。为某些函数传递其他类型的对象将导致异常抛出，并在错误输出上显示错误消息。文档在大多数地方详细说明了输入所接受的类型。如果您有三个通道图像作为输入，您可以做两件事：要么添加一个新通道（并使用char元素），要么分割图像并为每个图像调用函数。第一个并不推荐，因为这会浪费内存。
对于某些函数，元素（相邻项）的位置无关紧要，快速的解决方法是将其重塑为单通道图像。这是PSNR实现的情况，其中对于absdiff方法，邻居的值并不重要。但是，对于GaussianBlur，这不是一个选项，因此需要对SSIM使用分割方法。有了这些知识，你可以制作一个可行的GPU代码（比如我的GPU代码）并运行它。您会惊讶地发现，它可能会比您的CPU实现慢。
Optimization
这样做的原因是你把内存分配和数据传输的价格抛到了窗外。在GPU上，这非常高。优化的另一种可能性是在cv::cuda:Stream的帮助下引入异步OpenCV GPU调用。
- GPU上的内存分配相当可观。因此，如果可能的话，尽可能少地分配新内存。如果创建一个要多次调用的函数，最好只在第一次调用期间为该函数分配一次本地参数。为此，需要创建一个包含所有要使用的局部变量的数据结构。例如，在峰值信噪比的情况下，这些是：
struct BufferPSNR // Optimized GPU versions
{ // Data allocations are very expensive on GPU. Use a buffer to solve: allocate once reuse later.
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功能。流的并行化程度很小。通过使用流，我们可以在GPU已经执行给定方法的情况下进行数据分配、上传操作。例如，我们需要上传两张图片。我们把这些一个接一个地排队，并调用处理它的函数。这些函数将等待上载完成，但是在这种情况下，它会为下一个要执行的函数分配输出缓冲区。
gpu::Stream stream;
stream.enqueueConvert(b.gI1, b.t1, CV_32F); // Upload
gpu::split(b.t1, b.vI1, stream); // Methods (pass the stream as final parameter).
gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2
Result and conclusion
在Intel P8700笔记本电脑CPU与低端NVIDIA GT220M的搭配上，以下是性能数据：
Time of PSNR CPU (averaged for 10 runs): 41.4122 milliseconds. With result of: 19.2506
Time of PSNR GPU (averaged for 10 runs): 158.977 milliseconds. With result of: 19.2506
Initial call GPU optimized: 31.3418 milliseconds. With result of: 19.2506
Time of PSNR GPU OPTIMIZED ( / 10 runs): 24.8171 milliseconds. With result of: 19.2506
Time of MSSIM CPU (averaged for 10 runs): 484.343 milliseconds. With result of B0.890964 G0.903845 R0.936934
Time of MSSIM GPU (averaged for 10 runs): 745.105 milliseconds. With result of B0.89922 G0.909051 R0.968223
Time of MSSIM GPU Initial Call 357.746 milliseconds. With result of B0.890964 G0.903845 R0.936934
Time of MSSIM GPU OPTIMIZED ( / 10 runs): 203.091 milliseconds. With result of B0.890964 G0.903845 R0.936934
在这两种情况下，与CPU实现相比，我们的性能几乎提高了100%。它可能只是应用程序工作所需的改进。你可以在YouTube上看到一个运行时实例。
698
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
