阅读(2708) 赞(1)

GPU上的相似性检查（PNSR和SSIM）

2017-10-25 10:22:22 更新

目标

在OpenCV和相似度测量教程的视频输入中，我已经介绍了用于检查两个图像之间的相似性的PSNR和SSIM方法。正如你所看到的，执行过程需要相当长的一段时间，特别是在SSIM的情况下。然而，如果CPU的OpenCV实现的性能数字不能满足您的需求，并且您的系统中恰好有一个NVidia CUDA GPU设备，则所有这些都不会丢失。您可以尝试为视频卡端口或写入您的owm算法。

本教程将很好地掌握如何使用OpenCV的GPU模块来处理编码。作为先决条件，您应该已经知道如何处理core，highgui和imgproc模块。所以我们的主要目标是：

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

源代码

您也可以将源代码和视频文件放在samples/cpp/tutorial_code/gpu/gpu-basics-similarity/gpu-basics-similarityOpenCV源库的目录中或从这里下载。完整的源代码相当长（由于通过命令行参数和性能测量来控制应用程序）。因此，为了避免这些部分与这些部分混乱，这里只有功能本身。

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。这也是零和一之间的浮点数（更高更好），但是每个通道都有一个。因此，我们返回一个Scalar 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;
}

怎么做？ - 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。它的作用类似于仅具有2D限制的Mat，并且没有参考返回其功能（不能将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输入不同。您可以在文档中找到这些的完整列表：在这里或源代码附带的OpenCV参考手册。

要记住的另一件事是，不是所有的频道号码都可以在GPU上制作高效的算法。一般来说，我发现GPU图像的输入图像需要是一个或四个通道图像，而一个char或float类型用于项目大小。对GPU没有双重支持，对不起。传递某些功能的其他类型的对象将导致异常抛出，并在错误输出上显示错误消息。大多数地方的输入文件接受的类型的文档详细信息。如果您有三个通道图像作为输入，您可以做两件事情：添加新的通道（并使用char元素）或分割图像并调用每个图像的功能。第一个不是真的推荐，因为这浪费了浪费存储空间。

对于某些功能，元素（邻居项目）的位置无关紧要，快速解决方案是将其重新形成单个通道图像。这是PSNR实现的情况，对于absdiff方法，邻居的值不重要。然而，对于GaussianBlur，这不是一个选项，因此需要使用分割方法来进行SSIM。有了这个知识，你可以制作一个GPU可行的代码（像我的GPU一个），并运行它。您会惊讶地发现，它可能会比您的CPU执行速度慢。

优化

这样做的原因是您在窗口上扔出内存分配和数据传输的价格。在GPU上，这是很高的。优化的另一种可能性是在cv :: cuda :: Stream的帮助下引入异步OpenCV GPU调用。

GPU上的内存分配相当可观。因此，如果可能的话，尽可能多地分配新的内存。如果你创建一个函数你打算多次调用，最好在第一次调用期间为函数分配一次本地参数一次。为此，您将创建一个包含将使用的所有局部变量的数据结构。例如在PSNR的情况下，这些是：

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）;

使用asynchronous调用（cv :: cuda :: Stream）。默认情况下，只要您调用GPU功能，它将等待调用完成，然后返回结果。然而，可以进行asynchronous调用，这意味着它将调用操作执行，从而为算法提供昂贵的数据分配并立即返回。现在你可以调用另一个函数，如果你愿意的话。对于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

结果和结论

在配备低端NVidia GT220M的Intel P8700笔记本电脑CPU上，性能数字如下：

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％。这可能只是您的应用程序工作所需的改进

← 高级拼接API（Stitcher类）

使用cv :: cuda :: GpuMat推力（thrust） →

GPU上的相似性检查（PNSR和SSIM）

目标

源代码

怎么做？ - GPU

优化

结果和结论

推荐文章

推荐教程

技术文章