opencv8-GPU之相似性计算
Opencv支持GPU计算，并且包含成一个gpu类用来方便调用，所以不需要去加上什么__global__什么的很方便，不过同时这个类还是有不足的，待opencv小组的更新和完善。
这里先介绍在之前的《opencv4-highgui之视频的输入和输出以及滚动条》未介绍的图像的相似性检测，当然这是cpu版本，然后接着在介绍对应的gpu版本。这里只介绍了PSNR和SSIM两种用来进行对比图像的方法
原理：
PSNR: 当我们想检查压缩视频带来的细微差异的时候，就需要构建一个能够逐帧比较差视频差异的系统（也就是两幅图像帧之间的对比）。最常用的比较算法是PSNR( Peak signal-to-noise ratio)。这是个使用“局部均值误差”来判断差异的最简单的方法，假设有这两幅图像：I1和I2，它们的行列数分别是i，j，有c个通道（两幅图片同样的size）:
均值误差就是两个矩阵相减，然后将所有的元素相加除以元素的个数即可。PSNR公式如下:
每个像素的每个通道的值占用一个字节，值域[0,255]。这里每个像素会有 个有效的最大值
注意当两幅图像的相同的话，MSE的值会变成0。这样会导致PSNR的公式会除以0而变得没有意义。所以我们需要单独的处理这样的特殊情况。此外由于像素的动态范围很广，在处理时会使用对数变换来缩小范围。
在考察压缩后的视频时，这个值大约在30到50之间，数字越大则表明压缩质量越好。如果图像差异很明显，就可能会得到15甚至更低的值。PSNR算法简单，检查的速度也很快。但是其呈现的差异值有时候和人的主观感受不成比例。所以有另外一种称作 结构相似性 的算法做出了这方面的改进。
SSIM:是一种相比之前的方法更加直接的将参考图片与失真图片进行结构上对比的算法。得到的对比值范围在【0，1】，值越大，越靠近1，那么失真度越小。下面就是几个简单的例子：
上图中：不同图之间的的MSE都等于210，但是肉眼发现还是有很大不同的，具体来说 a、为原始图像；b、对比度拉伸图像，MSSIM= 0.9168；c、移去均值的图像，MSSIM= 0.9900；d、JPEG压缩后的图像，MSSIM= 0.6949；e、平滑后的图像 MSSIM= 0.7052；f、增加了椒盐噪声的图片，MSSIM= 0.7784.
按照论文中说的，其实就是将三个部分的值进行对比：
a) $\mu_x=\frac{1}{N}\sum_{i=1}^Nx_i$;
$l({\bf x},{\bf y})$;
$l({\bf x},{\bf y})=\frac{2\mu_x\mu_y+C_1}{\mu_x^2+\mu_y^2+C_1}$
就是求得两个图像块它们的均值，采用函数$l({\bf x},{\bf y})$来对比之间的不同；其中的C1的取值为：$C_1=(K_1L)^2$，L是像素值的动态范围，这里为255，对于8-bit的灰度图来说，其中的K1<=1,所以按照下面的示例代码来说，$C1= （0.01×255）^2 = 6.5025$；同理下面的$C2= （0.03×255）^2 = 58.5225$.而且这个公式对相对亮度是敏感的，为了改变这种情况，将R设置成相对差，使得$\mu_y=(1+R)\mu_x$，从而改写公式得到：$l({\bf x},{\bf y})=\frac{2(1+R)}{1+(1+R)^2+(C/\mu_x^2)}$；
b) $\sigma_x=\left(\frac{1}{N-1}\sum_{i=1}^N(x_i-\mu_x)^2\right)^{\frac{1}{2}}$;
$c({\bf x},{\bf y})$;
$c({\bf x},{\bf y}) = \frac{2\sigma_x\sigma_y+C_2}{\sigma_x^2+\sigma_y^2+C_2}$
就是求得两个图像块它们之间的标准差，采用函数$c(x,y)$来对比之间的不同；
c) $\frac{{\bf x}-\mu_x}{\sigma_x}$;
$s({\bf x},{\bf y})$;
$s({\bf x},{\bf y})=\frac{\sigma_{xy}+C_3}{\sigma_x\sigma_y+C_3}$;
$\sigma_{xy}=\frac{1}{N-1}\sum_{i=1}^N(x_i-\mu_x)(y_i-\mu_y)$
就是求得两个图像块它们求得的去均值和标准方差之后的两个值之间的对比。
z) $S(x,y)=f(l(x,y),c(x,y),s(x,y))$;
$SSIM({\bf x},{\bf y})=[l(x,y)]^\alpha \cdot [c(x,y)]^\beta \cdot [s(x,y)]^\gamma$
为了简化，就把三个指数都设置为1，就得到了
$SSIM({\bf x},{\bf y})=\frac{(2\mu_x\mu_y+C_1)(2\sigma_{xy}+C_2)}{(\mu_x^2+\mu_y^2+C_1)(\mu_x^2+\mu_y^2+C_2)}$
采用一个另外的函数来对比上面 三个不同值之间的对比来达到SSIM计算的目的。最后的S(x，y)必须要满足三个条件：
a、对称性$s({\bf x},{\bf y})=s({\bf x},{\bf y})$；
b、边界性$S({\bf x},{\bf y})\leq 1$；
c、唯一最大值有且只有当$x=y$时$S({\bf x},{\bf y})=1$。当然上面的方法是原理，但是计算颇为麻烦而且:
i）图像统计特征通常具有高度的不稳定性，
ii）图像的失真具有空间不变性，
iii）只有局部区域才是被人类感兴趣的，所以通常采用局部的SSIM指数代替全局的方法，在所提到的论文中采用11×11大小的图像块进行移动来获取权值，${\bf w}=\{w_i|i=1,2,...,N\}$，使用1.5采样的标准差，并且归一化到$\sum_{i=1}^Nw_i=1$
$\mu_x=\sum_{i=1}^Nw_ix_i$;
$\sigma_x=\left(\sum_{i=1}^Nw_i(x_i-\mu_x)^2\right)^{\frac{1}{2}}$;
$\sigma_{xy}=\sum_{i=1}^Nw_i(x_i-\mu_x)(y_i-\mu_y)$
最后使用mean来计算整个图像的SSIM指数。
SSIM的原理来自《Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004.》
一、cpu版本的PSNR和SSIM算法实现
PSNR的代码：
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]; // 叠加所有通道 if( sse <= 1e-10) // 如果值太小就直接等于0 return 0; else { double mse =sse /(double)(I1.channels() * I1.total());//求得mse，为后面做准备 double psnr = 10.0*log10((255*255)/mse);//公式中的MAX就等于像素的最大值为255. return psnr; } }
差分绝对值的函数原型：void absdiff(InputArray src1, InputArray src2, OutputArray dst)
第一个参数：输入的数组或者一个标量；
第二个参数：输入的数组或者一个标量；
第三个参数：输出的数组与输入的数组有着一样的size和typoe。
该函数的计算如下：
a、当它们有向图的size和type的时候基于两个数组的差分的绝对值：
b、当是一个数组和一个来自Scalar的标量或者有着与src1通道数一样多的元素的时候的差分绝对值是：
c、当与b的情况相反：
这里 I 是一个多维度索引的数组元素，在多通道数组下啊，每个通道是分别计算的
Notes:dang数组有着depth为CV_32S的时候saturation不会被使用，在这种情况下溢出有可能得到一个负值。
矩阵转换的函数原型：convertTo(OutputArray m, int rtype, double alpha=1, double beta=0 ) const
参数列表：
第一个参数：输出矩阵，如果在该操作之前没有一个合适的size或者type，那么就重新分配；
第二个参数：合适的输出矩阵的type和depth，不过与输入有着同样的通道；如果rtype是负的，那么输出矩阵的type就和输入的一样；
第三个参数：可选的缩放因子；
第四个参数：可选的添加到缩放后的值。
该函数将原像素值转换到目标数据类型，saturate_cast<>用来防止溢出：
矩阵的乘法的函数原型：MatExpr Mat::mul(InputArray m, double scale=1) const
参数列表：
第一个参数：另一个与*this 有着一样type和size的数组，或者一个矩阵表达式。
第二个参数：可选的缩放因子。
该函数返回一个临时对象来逐元素计算两个数组的乘法或者除法，顺带着可选的缩放。注意到这不是一个对应着一个更简单“×”操作符的矩阵乘法；
例子：
Mat C = A.mul(5/B); // equivalent to divide(A, B, C, 5)
矩阵求和的函数：CV_EXPORTS_AS(sumElems) Scalar sum(InputArray src);
该函数暂时不在refman中，返回的是一个Scalar值，共有四个元素，每个对应每个通道上这个矩阵的所有元素的和。
计算矩阵中元素个数的函数原型：size_t Mat::total() const
覆盖函数返回数组元素的个数（如果该数组表示一张图像，那么就返回像素的个数，也就是一个通道的矩阵上的元素个数）。
SSIM的代码：为了计算方便，使用高速平滑来得到均值的方法。
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); // 不能在单字节像素上进行计算，范围不够。 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);//mu1^2 Mat mu2_2 = mu2.mul(mu2);//mu2^2 Mat mu1_mu2 = mu1.mul(mu2);//mu1*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)),SSI的分子 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))，SSIM的分母Matssim_map;divide(t3,t1,ssim_map);
// ssim_map = t3./t1;Scalarmssim=mean(ssim_map);
// mssim = ssim_map的平均值，returnmssim;}
这个操作会针对图像的每个通道返回一个相似度，取值范围应该在0到1之间，取值为1时代表完全符合。然而尽管SSIM能产生更优秀的数据，但是由于高斯模糊很花时间，所以在一个实时系统（每秒24帧）中，人们还是更多地采用PSNR算法。
矩阵的表达式：
上图是在3.0beta中refman中的相关部分，这里是矩阵的表达式，通过将不同的操作符进行重载实现的。
Notes:逗号隔开的初始化事和一些可能的其他操作需要额外的显式初始化Mat()或者Mat_<T>()来解决可能造成的混乱：
下面是一些例子：
// compute pseudo-inverse of A, equivalent to A.inv(DECOMP_SVD) SVD svd(A); Mat pinvA = svd.vt.t()*Mat::diag(1./svd.w)*svd.u.t(); // compute the new vector of parameters in the Levenberg-Marquardt algorithm x -= (A.t()*A + lambda*Mat::eye(A.cols,A.cols,A.type())).inv(DECOMP_CHOLESKY)*(A.t()*err); // sharpen image using "unsharp mask" algorithm Mat blurred; double sigma = 1, threshold = 5, amount = 1; GaussianBlur(img, blurred, Size(), sigma, sigma); Mat lowConstrastMask = abs(img - blurred) < threshold; Mat sharpened = img*(1+amount) + blurred*(-amount); img.copyTo(sharpened, lowContrastMask);
矩阵的除法函数原型：void divide(InputArray src1, InputArray src2, OutputArray dst, double scale=1, int dtype=-1)
void divide(double scale, InputArray src2, OutputArray dst, int dtype=-1)；
参数列表：
src1：第一个输入的数组；
src2：第二个与src1有着同样size和type的素组
scale：缩放因子
dst：与src2有着一样size和type的输出数组
dtype：输出数组的可选depty；如果为 -1，dst的depty为src2.depth(),不过在逐数组除法情况下，当src1.depth（） == src2.depth(),你只能传递一个-1 给它。
该函数divide通过一个数组除以另一个：
或者当没有src1的时候一个数组被一个向量除：
当src2（I）为0，dst（I）同样也会是0.多通道数组的不同通道是独立处理的；
notes：当输出数组的depth为CV_32S的时候，saturation操作不会被使用，所以结果有可能得到一个不正确的符号（即在溢出之后变成负的了）。
正是这个原因，最开始的源码里，我们用PSNR算法去计算每一帧图像，而仅当PSNR算法计算出的结果低于我们想要的标准的时候，用SSIM算法去验证：
psnrV = getPSNR(srcFrame,dstFrame); cin>>psnrTriggerValue; //我们输入的值，假设想要验证两幅图片的相似性是否在50或者30 if ( （psnrV < psnrTriggerValue） && psnrV) { mssimV = getMSSIM(frameReference,frameUnderTest); cout << " MSSIM: " << " R " << setiosflags(ios::fixed) << setprecision(2) << mssimV.val[2] * 100 << "%" << " G " << setiosflags(ios::fixed) << setprecision(2) << mssimV.val[1] * 100 << "%" << " B " << setiosflags(ios::fixed) << setprecision(2) << mssimV.val[0] * 100 << "%"; } cout << endl; ////////////////////////////////// Show Image ///////////////////////////////////////////// imshow( "windowName1", srcFrame); imshow( "windowName2", dstFrame); c = cvWaitKey(delay); if (c == 27) break; }
二、GPU版本的PSNR和SSIM算法的实现
首先声明在3.0版本之前的refman中比如2.4.10中，为gpu:: ；而在3.0版本的refman中统一换成了cuda::，也就是名字空间换的更加贴近N卡了，不过奇怪的是在3.0beta的代码中还是使用的gpu名字空间，这当然一方面是3.0beta我下载的时间上刚好还是使用的2.4.10版本的代码（内部发现dll都是2.4.10版本和caffe的一些代码）不过按照refman的趋势看，肯定是换成cuda了以后，所以这里需要注意。
上面已经学习了检测两幅图像相似度的两种方法：PSNR和SSIM。正如我们所看到的，执行这些算法需要相当长的计算时间，其中SSIM（结构相似度）的算法代价相当高昂。然而，OpenCV现在支持Nvidia的CUDA加速，如果你有一块支持CUDA的的Nvidia显卡。您可以将算法改为使用GPU计算从而大幅提高效率。不过之所以做两个GPU函数是为了说明一个常识：简单的移植你的CPU程序到GPU反而会让速度变慢。如果你想要一些性能上的改善，你将需牢记住一些规则。
所需要开发的GPU模块尽可能多的与CPU对应，这样才能方便移植。在你编写任何代码之间要做的第一件事是连GPU模块到你的项目中，包括模块的头文件。所有GPU的函数和数据结构都在以 cv 命名空间的子空间 gpu 内。你可以将其添加为默认 use
namespace 关键字, 也可以显示的通过cv::来避免混乱。
#include <opencv2/gpu/gpu.hpp> // GPU结构和方法
GPU代表**图形处理单元**。最开始是为渲染各种图形场景而建立，这些场景是基于大量的矢量数据建立的。由于矢量图形的特殊性，数据不需要以串行的方式一步一步执行的，而是并行的方式一次性渲染大量的数据。从GPU的结构上来说，不像CPU基于数个寄存器和高速指令集，GPU一般有数百个较小的处理单元。这些处理单元每一个都比CPU的核心慢很多很多。然而，它的力量在于它的数目众多，能够同时进行大量运算。在过去的几年中已经有越来越多的人常识用GPU来进行图形渲染以外的工作。这就产生了GPGPU（general-purpose computation on graphics processing units)的概念。
图像处理器有它自己的内存，一般称呼为显存。当你从硬盘驱动器读数据并产生一个 Mat 对象的时候，数据是放在普通内存当中的（由CPU掌管）CPU可以直接操作内存, 然而GPU不能这样，对于电脑来说GPU只是个外设，它只能操作它自己的显存，当计算时，需要先让CPU将用于计算的信息转移到GPU掌管的显存上。 这是通过一个上传过程完成，需要比访问内存多很多的时间。而且最终的计算结果将要回传到你的系统内存处理器才能和其他代码交互，由于传输的代价高昂，所以注定移植太小的函数到GPU上并不会提高效率。
Mat对象仅仅存储在内存或者CPU缓存中。为了得到一个GPU能直接访问的opencv 矩阵你必须使用GPU对象 GpuMat 。它的工作方式类似于2维 Mat，唯一的限制是你不能直接引用GPU函数。（因为它们本质上是完全不同的代码，不能混合引用）。要传输*Mat*对象到*GPU*上并创建GpuMat时需要调用上传函数，在回传时，可以使用简单的重载赋值操作符或者调用下载函数。
Mat I1; // 内存对象，可以用imread来创建 gpu::GpuMat gI; // GPU 矩阵 - 现在为空 gI1.upload(I1); //将内存数据上传到显存中 I1 = gI1; //回传, gI1.download(I1) 也可以
记住：一旦你的数据被传到GPU显存中，你就只能调用GPU函数来操作，大部分gpu函数名字与原来CPU名字保持一致，不同的是他们只接收 GpuMat 输入。可以在:在线文档 中找到一些说明，或者查阅OpenCV参考手册。
另外要记住的是：并非所有的图像类型都可以用于GPU算法中。很多时候，GPU函数只能接收单通道或者4通道的uchar或者float类型（CV_8UC1,CV_8UC4,CV_32FC1,CV_32FC4），GPU函数不支持双精度。如果你试图向这些函数传入非指定类型的数据时，这些函数会抛出异常或者输出错误信息。这些函数的文档里大都说明了接收的数据类型。如果函数只接收4通道或者单通道的图像而你的数据类型刚好是3通道的话，你能做的就是两件事：1.增加一个新的通道（使用char类型）或 2.将图像的三个通道切分为三个独立的图像，为每个图像调用该函数。不推荐使用第一个方法，因为有些浪费内存。
对于不需要在意元素的坐标位置的某些函数，快速解决办法就是将它直接当作一个单通道图像处理。这种方法适用于PSNR中需要调用的absdiff方法。然而，对于GaussianBlur就不能这么用，对于SSIM，就需要使用分离通道的方法。知道这些知识就已经可以让你的GPU开始运行代码。但是你也可能发现GPU“加速”后的代码仍然可能会低于你的CPU执行速度。
优化GPU：
因为大量的时间都耗费在传输数据到显存这样的步骤了，大量的计算时间被消耗在内存分配和传输上，GPU的高计算能力根本发挥不出来，我们可以利用gpu::Stream 来进行异步传输，从而节约一些时间。GPU上的显存分配代价是相当大的。因此如果可能的话，应该尽可能少的分配新的内存。如果你需要创建一个调用多次的函数，就应该在一开始分配全部的内存，而且仅仅分配一次。在第一次调用的时候可以创建一个结构体包含所有将使用局部变量。
这里先介绍PSNR的基本GPU版本，然后介绍优化后的GPU版本：
double getPSNR_GPU(const Mat& I1, const Mat& I2) { gpu::GpuMat gI1, gI2, gs, t1,t2; //在gpu上分配几个矩阵 gI1.upload(I1); gI2.upload(I2);//上传数据 gI1.convertTo(t1, CV_32F);//数据的类型转换 gI2.convertTo(t2, CV_32F); gpu::absdiff(t1.reshape(1), t2.reshape(1), gs); //在gpu上执行差分绝对值函数 gpu::multiply(gs, gs, gs);//在gpu上执行乘法 Scalar s = gpu::sum(gs);//在gpu上执行求和并下载结果 double sse = s.val[0] + s.val[1] + s.val[2];//在cpu上操作 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; } }
gpu版本的Mat类：class gpu::GpuMat
有着引用计数的GPU内存的基本存储类。它与cpu版本的Mat的接口匹配的时候是有着如下的限制：
a、没有任意维度的支持（只有2D）；
b、不会返回对它们数据的引用的函数（因为在GPU上的引用对于cpu来说不合法）
c、不支持表达式模板技术
注意到后者的限制会导致那些会引起内存分配的矩阵操作符被重载。GpuMat类对于gpu::PtrStepSz和gpu::PtrStep更加的方便，所以可以直接传给核kernel（cuda中的术语）。
notes：与Mat对比来说，大部分情况下GpuMat::isContinuous() == false。这是因为在gpu中行与行之间是分隔开的，不过单一行GpuMat总是一个连续的矩阵（也就是一行向量）。
class CV_EXPORTS GpuMat { public: //! 默认构造函数 GpuMat(); //! 以指定的size和type来构造GpuMat GpuMat(int rows, int cols, int type); GpuMat(Size size, int type); ..... //! 从主机内存中构建GpuMat (Blocking call) explicit GpuMat(InputArray arr);（这是3.0版本的） //! 返回轻量级的 PtrStepSz structure for passing //to nvcc-compiled code. Contains size, data ptr and step. template <class T> operator PtrStepSz<T>() const; template <class T> operator PtrStep<T>() const; //! pefroms upload data to GpuMat (Blocking call) void upload(InputArray arr); //! pefroms upload data to GpuMat (Non-Blocking call) void upload(InputArray arr, Stream& stream); //! pefroms download data from device to host memory (Blocking call) void download(OutputArray dst) const; //! pefroms download data from device to host memory (Non-Blocking call) void download(OutputArray dst, Stream& stream) const; };
Notes:不推荐使用static或者全局GpuMat变量分配去依赖它们的析构函数，因为这些变量和CUDA上下文的析构的顺序是未定义的。如果CUDA上下文之前已经被毁了，那么GPU内存的释放函数会返回错误。
上传到GPU的函数
优化后的GPU版本：
struct BufferPSNR // GPU的优化版本 { // 基于GPU的数据分配代价是高昂的。所以使用一个缓冲区来改善这点：分配一次，以后重用 gpu::GpuMat gI1, gI2, gs, t1,t2; gpu::GpuMat buf; }; BufferPSNR bufferPSNR;
最后每次调用函数时都使用这个结构:当你传入的参数为 b.gI1 ， b.buf （下面的优化后的函数）等，除非类型发生变化，GpuMat都将直接使用原来的内容而不重新分配。
double getPSNR_GPU_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); gpu::absdiff(b.t1.reshape(1), b.t2.reshape(1), b.gs); gpu::multiply(b.gs, b.gs, b.gs); double sse = gpu::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的基本GPU版本，然后介绍优化后的GPU版本：
Scalar getMSSIM_GPU( const Mat& i1, const Mat& i2) { const float C1 = 6.5025f, C2 = 58.5225f; /***************************** INITS **********************************/ gpu::GpuMat gI1, gI2, gs1, t1,t2; gI1.upload(i1); gI2.upload(i2); gI1.convertTo(t1, CV_MAKE_TYPE(CV_32F, gI1.channels())); gI2.convertTo(t2, CV_MAKE_TYPE(CV_32F, gI2.channels())); vector<gpu::GpuMat> vI1, vI2; gpu::split(t1, vI1); gpu::split(t2, vI2); Scalar mssim; for( int i = 0; i < gI1.channels(); ++i ) { gpu::GpuMat I2_2, I1_2, I1_I2; gpu::multiply(vI2[i], vI2[i], I2_2); // I2^2 gpu::multiply(vI1[i], vI1[i], I1_2); // I1^2 gpu::multiply(vI1[i], vI2[i], I1_I2); // I1 * I2 /*************************** END INITS **********************************/ gpu::GpuMat mu1, mu2; // PRELIMINARY COMPUTING gpu::GaussianBlur(vI1[i], mu1, Size(11, 11), 1.5); gpu::GaussianBlur(vI2[i], mu2, Size(11, 11), 1.5); gpu::GpuMat mu1_2, mu2_2, mu1_mu2; gpu::multiply(mu1, mu1, mu1_2); gpu::multiply(mu2, mu2, mu2_2); gpu::multiply(mu1, mu2, mu1_mu2); gpu::GpuMat sigma1_2, sigma2_2, sigma12; gpu::GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5); sigma1_2 -= mu1_2; gpu::GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5); sigma2_2 -= mu2_2; gpu::GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5); sigma12 -= mu1_mu2; ///////////////////////////////// FORMULA //////////////////////////////// gpu::GpuMat t1, t2, t3; t1 = 2 * mu1_mu2 + C1; t2 = 2 * sigma12 + C2; gpu::multiply(t1, t2, t3); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2)) t1 = mu1_2 + mu2_2 + C1; t2 = sigma1_2 + sigma2_2 + C2; gpu::multiply(t1, t2, t1); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2)) gpu::GpuMat ssim_map; gpu::divide(t3, t1, ssim_map); // ssim_map = t3./t1; Scalar s = gpu::sum(ssim_map); mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols); } return mssim; }
优化后的GPU版本：
struct BufferMSSIM // 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 I1_2, I2_2, I1_I2; vector<gpu::GpuMat> vI1, vI2; gpu::GpuMat mu1, mu2; gpu::GpuMat mu1_2, mu2_2, mu1_mu2; gpu::GpuMat sigma1_2, sigma2_2, sigma12; gpu::GpuMat t3; gpu::GpuMat ssim_map; gpu::GpuMat buf; }; Scalar getMSSIM_GPU_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b) { int cn = i1.channels(); const float C1 = 6.5025f, C2 = 58.5225f; /***************************** INITS **********************************/ b.gI1.upload(i1); b.gI2.upload(i2); gpu::Stream stream; stream.enqueueConvert(b.gI1, b.t1, CV_32F); stream.enqueueConvert(b.gI2, b.t2, CV_32F); gpu::split(b.t1, b.vI1, stream); gpu::split(b.t2, b.vI2, stream); Scalar mssim; for( int i = 0; i < b.gI1.channels(); ++i ) { gpu::multiply(b.vI2[i], b.vI2[i], b.I2_2, stream); // I2^2 gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream); // I1^2 gpu::multiply(b.vI1[i], b.vI2[i], b.I1_I2, stream); // I1 * I2 gpu::GaussianBlur(b.vI1[i], b.mu1, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream); gpu::GaussianBlur(b.vI2[i], b.mu2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream); gpu::multiply(b.mu1, b.mu1, b.mu1_2, stream); gpu::multiply(b.mu2, b.mu2, b.mu2_2, stream); gpu::multiply(b.mu1, b.mu2, b.mu1_mu2, stream); gpu::GaussianBlur(b.I1_2, b.sigma1_2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream); gpu::subtract(b.sigma1_2, b.mu1_2, b.sigma1_2, stream); //b.sigma1_2 -= b.mu1_2; - This would result in an extra data transfer operation gpu::GaussianBlur(b.I2_2, b.sigma2_2, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream); gpu::subtract(b.sigma2_2, b.mu2_2, b.sigma2_2, stream); //b.sigma2_2 -= b.mu2_2; gpu::GaussianBlur(b.I1_I2, b.sigma12, Size(11, 11), 1.5, 0, BORDER_DEFAULT, -1, stream); gpu::subtract(b.sigma12, b.mu1_mu2, b.sigma12, stream); //b.sigma12 -= b.mu1_mu2; //here too it would be an extra data transfer due to call of operator*(Scalar, Mat) gpu::multiply(b.mu1_mu2, 2, b.t1, stream); //b.t1 = 2 * b.mu1_mu2 + C1; gpu::add(b.t1, C1, b.t1, stream); gpu::multiply(b.sigma12, 2, b.t2, stream); //b.t2 = 2 * b.sigma12 + C2; gpu::add(b.t2, C2, b.t2, stream); gpu::multiply(b.t1, b.t2, b.t3, stream); // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2)) gpu::add(b.mu1_2, b.mu2_2, b.t1, stream); gpu::add(b.t1, C1, b.t1, stream); gpu::add(b.sigma1_2, b.sigma2_2, b.t2, stream); gpu::add(b.t2, C2, b.t2, stream); gpu::multiply(b.t1, b.t2, b.t1, stream); // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2)) gpu::divide(b.t3, b.t1, b.ssim_map, stream); // ssim_map = t3./t1; stream.waitForCompletion(); Scalar s = gpu::sum(b.ssim_map, b.buf); mssim.val[i] = s.val[0] / (b.ssim_map.rows * b.ssim_map.cols); } return mssim; }
