20557人阅读
计算机视觉（42）
图形图像（59）
OpenCV（60）
转载请注明出处：，
对计算图像相似度的方法，本文做了如下总结，主要有三种办法：
1.PSNR峰值信噪比
PSNR（Peak Signal to Noise Ratio），一种全参考的图像质量评价指标。
PSNR是最普遍和使用最为广泛的一种图像客观评价指标，然而它是基于对应像素点间的误差，即基于误差敏感的图像质量评价。由于并未考虑到人眼的视觉特性（人眼对空间频率较低的对比差异敏感度较高，人眼对亮度对比差异的敏感度较色度高，人眼对一个区域的感知结果会受到其周围邻近区域的影响等），因而经常出现评价结果与人的主观感觉不一致的情况。
SSIM（structural similarity）结构相似性，也是一种全参考的图像质量评价指标，它分别从亮度、对比度、结构三方面度量图像相似性。
SSIM取值范围[0,1]，值越大，表示图像失真越小.
在实际应用中，可以利用滑动窗将图像分块，令分块总数为N，考虑到窗口形状对分块的影响，采用高斯加权计算每一窗口的均值、方差以及协方差，然后计算对应块的结构相似度SSIM，最后将平均值作为两图像的结构相似性度量，即平均结构相似性MSSIM：
[1] 峰值信噪比-维基百科
[2] 王宇庆，刘维亚，王勇. 一种基于局部方差和结构相似度的图像质量评价方法[J]. 光电子激光，2008。
官方文档的说明，不过是GPU版本的，我们可以修改不用gpu不然还得重新编译
#include "stdafx.h"
#include &iostream&
#include &sstream&
#include &opencv2/core/core.hpp&
#include &opencv2/imgproc/imgproc.hpp&
#include &opencv2/highgui/highgui.hpp&
#include &opencv2/gpu/gpu.hpp&
using namespace std;
using namespace
double getPSNR(const Mat& I1, const Mat& I2);
Scalar getMSSIM( const Mat& I1, const Mat& I2);
double getPSNR_GPU(const Mat& I1, const Mat& I2);
Scalar getMSSIM_GPU( const Mat& I1, const Mat& I2);
struct BufferPSNR
gpu::GpuMat gI1, gI2, gs, t1,t2;
double getPSNR_GPU_optimized(const Mat& I1, const Mat& I2, BufferPSNR& b);
struct BufferMSSIM
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_
Scalar getMSSIM_GPU_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b);
void help()
&& "\n--------------------------------------------------------------------------" && endl
&& "This program shows how to port your CPU code to GPU or write that from scratch." && endl
&& "You can see the performance improvement for the similarity check methods (PSNR and SSIM)."
&& "Usage:"
&& "./gpu-basics-similarity referenceImage comparedImage numberOfTimesToRunTest(like 10)." && endl
&& "--------------------------------------------------------------------------"
int main(int argc, char *argv[])
Mat I1 = imread("swan1.jpg",1);
Mat I2 = imread("swan2.jpg",1);
if (!I1.data || !I2.data)
cout && "Couldn't read the image";
BufferPSNR bufferPSNR;
BufferMSSIM bufferMSSIM;
int TIMES;
stringstream sstr("500");
sstr && TIMES;
double time,
time = (double)getTickCount();
for (int i = 0; i & TIMES; ++i)
result = getPSNR(I1,I2);
time = 1000*((double)getTickCount() - time)/getTickFrequency();
time /= TIMES;
cout && "Time of PSNR CPU (averaged for " && TIMES && " runs): " && time && " milliseconds."
&& " With result of: " &&
getchar();
double getPSNR(const Mat& I1, const Mat& I2)
absdiff(I1, I2, s1);
s1.convertTo(s1, CV_32F);
s1 = s1.mul(s1);
Scalar s = sum(s1);
double sse = s.val[0] + s.val[1] + s.val[2];
if( sse &= 1e-10)
mse =sse /(double)(I1.channels() * I1.total());
double psnr = 10.0*log10((255*255)/mse);
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)
double mse = sse /(double)(I1.channels() * I1.total());
double psnr = 10.0*log10((255*255)/mse);
double getPSNR_GPU(const Mat& I1, const Mat& I2)
gpu::GpuMat gI1, gI2, gs, t1,t2;
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::multiply(gs, gs, gs);
Scalar s = gpu::sum(gs);
double sse = s.val[0] + s.val[1] + s.val[2];
if( sse &= 1e-10)
mse =sse /(double)(gI1.channels() * I1.total());
double psnr = 10.0*log10((255*255)/mse);
Scalar getMSSIM( const Mat& i1, const Mat& i2)
const double C1 = 6.5025, C2 = 58.5225;
Mat I1, I2;
i1.convertTo(I1, d);
i2.convertTo(I2, d);
= I2.mul(I2);
= I1.mul(I1);
= I1.mul(I2);
Mat mu1, mu2;
GaussianBlur(I1, mu1, Size(11, 11), 1.5);
GaussianBlur(I2, mu2, Size(11, 11), 1.5);
mu1.mul(mu1);
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);
t1 = mu1_2 + mu2_2 + C1;
t2 = sigma1_2 + sigma2_2 + C2;
t1 = t1.mul(t2);
divide(t3, t1, ssim_map);
Scalar mssim = mean( ssim_map );
Scalar getMSSIM_GPU( const Mat& i1, const Mat& i2)
const float C1 = 6.5025f, C2 = 58.5225f;
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);
for( int i = 0; i & gI1.channels(); ++i )
gpu::GpuMat I2_2, I1_2, I1_I2;
gpu::multiply(vI2[i], vI2[i], I2_2);
gpu::multiply(vI1[i], vI1[i], I1_2);
gpu::multiply(vI1[i], vI2[i], I1_I2);
gpu::GpuMat mu1, mu2;
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);
gpu::subtract(sigma1_2,mu1_2,sigma1_2);
gpu::GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
gpu::GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
(Mat)sigma12 =(Mat)sigma12 - (Mat)mu1_mu2;
gpu::GpuMat t1, t2, t3;
gpu::GpuMat ssim_
gpu::divide(t3, t1, ssim_map);
Scalar s = gpu::sum(ssim_map);
mssim.val[i] = s.val[0] / (ssim_map.rows * ssim_map.cols);
Scalar getMSSIM_GPU_optimized( const Mat& i1, const Mat& i2, BufferMSSIM& b)
int cn = i1.channels();
const float C1 = 6.5025f, C2 = 58.5225f;
b.gI1.upload(i1);
b.gI2.upload(i2);
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);
for( int i = 0; i & b.gI1.channels(); ++i )
gpu::multiply(b.vI2[i], b.vI2[i], b.I2_2, stream);
gpu::multiply(b.vI1[i], b.vI1[i], b.I1_2, stream);
gpu::multiply(b.vI1[i], b.vI2[i], b.I1_I2, 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::multiply(b.mu1_mu2, 2, b.t1, stream);
gpu::multiply(b.sigma12, 2, b.t2, stream);
gpu::multiply(b.t1, b.t2, b.t3, stream);
gpu::multiply(b.t1, b.t2, b.t1, stream);
gpu::divide(b.t3, b.t1, b.ssim_map, stream);
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);
两幅一样的图片，对比结果：
2.感知哈希算法
(perceptual hash algorithm）
感知哈希算法(perceptual hash algorithm)，它的作用是对每张图像生成一个“指纹”(fingerprint)字符串，然后比较不同图像的指纹。结果越接近，就说明图像越相似。
实现步骤：
缩小尺寸：将图像缩小到8*8的尺寸，总共64个像素。这一步的作用是去除图像的细节，只保留结构/明暗等基本信息，摒弃不同尺寸/比例带来的图像差异；
简化色彩：将缩小后的图像，转为64级灰度，即所有像素点总共只有64种颜色；
计算平均值：计算所有64个像素的灰度平均值；
比较像素的灰度：将每个像素的灰度，与平均值进行比较，大于或等于平均值记为1，小于平均值记为0；
计算哈希值：将上一步的比较结果，组合在一起，就构成了一个64位的整数，这就是这张图像的指纹。组合的次序并不重要，只要保证所有图像都采用同样次序就行了；
得到指纹以后，就可以对比不同的图像，看看64位中有多少位是不一样的。在理论上，这等同于”汉明距离”(Hamming distance,在信息论中，两个等长字符串之间的汉明距离是两个字符串对应位置的不同字符的个数)。如果不相同的数据位数不超过5，就说明两张图像很相似；如果大于10，就说明这是两张不同的图像。
以上内容摘自：
#include "stdafx.h"
#include &iostream&
#include &opencv2/core/core.hpp&
#include &opencv2/highgui/highgui.hpp&
#include &opencv2/imgproc/imgproc.hpp&
#pragma comment(lib,"opencv_core2410d.lib")
#pragma comment(lib,"opencv_highgui2410d.lib")
#pragma comment(lib,"opencv_imgproc2410d.lib")
using namespace std;
int _tmain(int argc, _TCHAR* argv[])
string strSrcImageName = "swan.jpg";
cv::Mat matSrc, matSrc1, matSrc2;
matSrc = cv::imread(strSrcImageName, CV_LOAD_IMAGE_COLOR);
CV_Assert(matSrc.channels() == 3);
cv::resize(matSrc, matSrc1, cv::Size(357, 419), 0, 0, cv::INTER_NEAREST);
cv::resize(matSrc, matSrc2, cv::Size(2177, 3233), 0, 0, cv::INTER_LANCZOS4);
cv::Mat matDst1, matDst2;
cv::resize(matSrc1, matDst1, cv::Size(8, 8), 0, 0, cv::INTER_CUBIC);
cv::resize(matSrc2, matDst2, cv::Size(8, 8), 0, 0, cv::INTER_CUBIC);
cv::cvtColor(matDst1, matDst1, CV_BGR2GRAY);
cv::cvtColor(matDst2, matDst2, CV_BGR2GRAY);
int iAvg1 = 0, iAvg2 = 0;
int arr1[64], arr2[64];
for (int i = 0; i & 8; i++)
uchar* data1 = matDst1.ptr&uchar&(i);
uchar* data2 = matDst2.ptr&uchar&(i);
int tmp = i * 8;
for (int j = 0; j & 8; j++)
int tmp1 = tmp +
arr1[tmp1] = data1[j] / 4 * 4;
arr2[tmp1] = data2[j] / 4 * 4;
iAvg1 += arr1[tmp1];
iAvg2 += arr2[tmp1];
iAvg1 /= 64;
iAvg2 /= 64;
for (int i = 0; i & 64; i++)
arr1[i] = (arr1[i] &= iAvg1) ? 1 : 0;
arr2[i] = (arr2[i] &= iAvg2) ? 1 : 0;
int iDiffNum = 0;
for (int i = 0; i & 64; i++)
if (arr1[i] != arr2[i])
cout&&"iDiffNum = "&&iDiffNum&&
if (iDiffNum &= 5)
cout&&"two images are very similar!"&&
else if (iDiffNum & 10)
cout&&"they are two different images!"&&
cout&&"two image are somewhat similar!"&&
getchar();
一幅图片自己对比：
3.计算特征点
OpenCV的feature2d module中提供了从局部图像特征（Local image feature）的检测、特征向量（feature vector）的提取，到特征匹配的实现。其中的局部图像特征包括了常用的几种局部图像特征检测与描述算子，如FAST、SURF、SIFT、以及ORB。对于高维特征向量之间的匹配，OpenCV主要有两种方式：1）BruteForce穷举法；2）FLANN近似K近邻算法（包含了多种高维特征向量匹配的算法，例如随机森林等）。
feature2d module：
OpenCV FLANN:
#ifndef _FEATURE_H_
#define _FEATURE_H_
#include &iostream&
#include &vector&
#include &string&
#include &opencv2/opencv.hpp&
#include &opencv2/core/core.hpp&
#include &opencv2/highgui/highgui.hpp&
#include &opencv2/imgproc/imgproc.hpp&
#include &opencv2/nonfree/nonfree.hpp&
#include &opencv2/nonfree/features2d.hpp&
using namespace
using namespace std;
class Feature
Feature();
~Feature();
Feature(const string& detectType, const string& extractType, const string& matchType);
void detectKeypoints(const Mat& image, vector&KeyPoint&& keypoints);
void extractDescriptors(const Mat& image, vector&KeyPoint&& keypoints, Mat& descriptor);
void bestMatch(const Mat& queryDescriptor, Mat& trainDescriptor, vector&DMatch&& matches);
void knnMatch(const Mat& queryDescriptor, Mat& trainDescriptor, vector&vector&DMatch&&& matches, int k);
void saveKeypoints(const Mat& image, const vector&KeyPoint&& keypoints, const string& saveFileName = "");
void saveMatches(const Mat& queryImage,
const vector&KeyPoint&& queryKeypoints,
const Mat& trainImage,
const vector&KeyPoint&& trainKeypoints,
const vector&DMatch&& matches,
const string& saveFileName = "");
Ptr&FeatureDetector& m_
Ptr&DescriptorExtractor& m_
Ptr&DescriptorMatcher& m_
string m_detectT
string m_extractT
string m_matchT
#include "stdafx.h"
#include "localfeature.h"
Feature::Feature()
m_detectType = "SIFT";
m_extractType = "SIFT";
m_matchType = "BruteForce";
Feature::~Feature()
Feature::Feature(const string& detectType, const string& extractType, const string& matchType)
assert(!detectType.empty());
assert(!extractType.empty());
assert(!matchType.empty());
m_detectType = detectT
m_extractType = extractT
m_matchType = matchT
void Feature::detectKeypoints(const Mat& image, std::vector&KeyPoint&& keypoints)
assert(image.type() == CV_8UC1);
assert(!m_detectType.empty());
keypoints.clear();
initModule_nonfree();
m_detector = FeatureDetector::create(m_detectType);
m_detector-&detect(image, keypoints);
void Feature::extractDescriptors(const Mat& image, std::vector&KeyPoint&& keypoints, Mat& descriptor)
assert(image.type() == CV_8UC1);
assert(!m_extractType.empty());
initModule_nonfree();
m_extractor = DescriptorExtractor::create(m_extractType);
m_extractor-&compute(image, keypoints, descriptor);
void Feature::bestMatch(const Mat& queryDescriptor, Mat& trainDescriptor, std::vector&DMatch&& matches)
assert(!queryDescriptor.empty());
assert(!trainDescriptor.empty());
assert(!m_matchType.empty());
matches.clear();
m_matcher = DescriptorMatcher::create(m_matchType);
m_matcher-&add(std::vector&Mat&(1, trainDescriptor));
m_matcher-&train();
m_matcher-&match(queryDescriptor, matches);
void Feature::knnMatch(const Mat& queryDescriptor, Mat& trainDescriptor, std::vector&std::vector&DMatch&&& matches, int k)
assert(k & 0);
assert(!queryDescriptor.empty());
assert(!trainDescriptor.empty());
assert(!m_matchType.empty());
matches.clear();
m_matcher = DescriptorMatcher::create(m_matchType);
m_matcher-&add(std::vector&Mat&(1, trainDescriptor));
m_matcher-&train();
m_matcher-&knnMatch(queryDescriptor, matches, k);
void Feature::saveKeypoints(const Mat& image, const vector&KeyPoint&& keypoints, const string& saveFileName)
assert(!saveFileName.empty());
cv::drawKeypoints(image, keypoints, outImage, Scalar(255,255,0), DrawMatchesFlags::DRAW_RICH_KEYPOINTS );
string saveKeypointsImgName = saveFileName + "_" + m_detectType + ".jpg";
imwrite(saveKeypointsImgName, outImage);
void Feature::saveMatches(const Mat& queryImage,
const vector&KeyPoint&& queryKeypoints,
const Mat& trainImage,
const vector&KeyPoint&& trainKeypoints,
const vector&DMatch&& matches,
const string& saveFileName)
assert(!saveFileName.empty());
cv::drawMatches(queryImage, queryKeypoints, trainImage, trainKeypoints, matches, outImage,
Scalar(255, 0, 0), Scalar(0, 255, 255), vector&char&(),
DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
string saveMatchImgName = saveFileName + "_" + m_detectType + "_" + m_extractType + "_" + m_matchType + ".jpg";
imwrite(saveMatchImgName, outImage);
#include "stdafx.h"
#include &iostream&
#include &opencv2/core/core.hpp&
#include &opencv2/highgui/highgui.hpp&
#include &opencv2/imgproc/imgproc.hpp&
#include &opencv2/nonfree/nonfree.hpp&
#include &opencv2/nonfree/features2d.hpp&
#include "localfeature.h"
#pragma comment(lib,"opencv_core2410d.lib")
#pragma comment(lib,"opencv_highgui2410d.lib")
#pragma comment(lib,"opencv_imgproc2410d.lib")
#pragma comment(lib,"opencv_nonfree2410d.lib")
#pragma comment(lib,"opencv_features2d2410d.lib")
using namespace std;
int main(int argc, char** argv)
string detectorType = "SIFT";
string extractorType = "SIFT";
string matchType = "BruteForce";
string queryImagePath = "swan.jpg";
string trainImagePath = "swan.jpg";
Mat queryImage = imread(queryImagePath, CV_LOAD_IMAGE_GRAYSCALE);
if (queryImage.empty())
cout&&"read failed"&&
return -1;
Mat trainImage = imread(trainImagePath, CV_LOAD_IMAGE_GRAYSCALE);
if (trainImage.empty())
cout&&"read failed"&&
return -1;
Feature feature(detectorType, extractorType, matchType);
vector&KeyPoint& queryKeypoints, trainK
feature.detectKeypoints(queryImage, queryKeypoints);
feature.detectKeypoints(trainImage, trainKeypoints);
Mat queryDescriptor, trainD
feature.extractDescriptors(queryImage, queryKeypoints, queryDescriptor);
feature.extractDescriptors(trainImage, trainKeypoints, trainDescriptor);
vector&DMatch&
feature.bestMatch(queryDescriptor, trainDescriptor, matches);
vector&vector&DMatch&&
feature.knnMatch(queryDescriptor, trainDescriptor, knnmatches, 2);
feature.saveMatches(queryImage, queryKeypoints, trainImage, trainKeypoints, matches, "../");
两幅同样图片结果：
转载请注明出处：，
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
公众号：&&&&&老王和他的IT界朋友们
欢迎投稿：&&
我们想用一段音乐，几张图片，
些许文字绘制的IT圈心路历程，
为攻城狮，为程序员带来更多的人文关怀。
访问：467265次
积分：7421
积分：7421
排名：第1881名
原创：150篇
转载：47篇
译文：61篇
评论：370条
文章：35篇
阅读：27013
阅读：24890
阅读：21849
阅读：12122
文章：45篇
阅读：175781
(5)(4)(7)(8)(4)(4)(4)(4)(4)(4)(5)(4)(7)(13)(14)(17)(12)(5)(4)(22)(14)(12)(5)(15)(1)(3)(5)(1)(3)(1)(3)(3)(3)(1)(1)(2)(1)(5)(4)(5)(1)(1)(9)解析opencv中Box Filter的实现并提出进一步加速的方案（源码共享）。 - 推酷
解析opencv中Box Filter的实现并提出进一步加速的方案（源码共享）。
说明：本文所有算法的涉及到的优化均指在PC上进行的，对于其他构架是否合适未知，请自行试验。
Box Filter，最经典的一种领域操作，在无数的场合中都有着广泛的应用，作为一个很基础的函数，其性能的好坏也直接影响着其他相关函数的性能，最典型莫如现在很好的EPF滤波器：GuideFilter。因此其优化的档次和程度是非常重要的，网络上有很多相关的代码和博客对该算法进行讲解和优化，提出了不少O(1)算法，但所谓的0(1)算法也有优劣之分，0(1)只是表示执行时间和某个参数无关，但本身的耗时还是有区别。比较流行的就是累积直方图法，其实这个是非常不可取的，因为稍大的图就会导致溢出，这里我们解析下 opencv的相关代码，看看他是如何进行优化的。
首先找到opencv的代码的位置，其在\opencv\sources\modules\imgproc\src\smooth.cpp中。
Box Filter 是一种行列可分离的滤波，因此，opencv也是这样处理的，先进行行方向的滤波，得到中间结果，然后再对中间结果进行列方向的处理，得到最终的结果。
opencv 行方向处理的相关代码如下：
template&typename T, typename ST&
struct RowSum :
public BaseRowFilter
RowSum( int _ksize, int _anchor ) :
BaseRowFilter()
anchor = _
virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
const T* S = (const T*)
ST* D = (ST*)
int i = 0, k, ksz_cn = ksize*
width = (width - 1)*
for( k = 0; k & k++, S++, D++ )
for( i = 0; i & ksz_ i += cn )
s += S[i];
for( i = 0; i & i += cn )
s += S[i + ksz_cn] - S[i];
这个代码考虑了多个通道以及多种数据类型的情况，为了分析方便我们重写下单通道时的核心代码。
for(Z = 0, Value = 0; Z & S Z++)
Value += RowData[Z];
LinePD[0] = V
for(X = 1; X & W X ++)
Value += RowData[X + Size - 1] - RowData[X - 1];
LinePD[X] = V
上述代码中RowData是指对某一行像素进行扩展后的像素值，其宽度为Width + Radius + Radius, Radius是值滤波的半径， 而代码中Size = 2 * Radius + 1，LinePD即所谓的中间结果，考虑数据类型能表达的范围，必须使用 int类型。
对于每行第一个点很简单，直接用for计算出行方向的指定半径内的累加值。而之后所有点，用前一个累计值加上新加入的点，然后去除掉移出的哪一个点，得到新的累加值。这个方法在很多文献中都有提及，并没有什么新鲜之处，这里不多说。
按照正常的思维，在列方向的处理应该和行方向完全相同，只是处理的方向的不同、处理的数据源的不同以及最后需要对计算结果多一个归一化的过程而已。事实也是如此，有很多算法都是这样做的，但是由于CPU构架的一些原因（主要是cache miss的增加），同样的算法沿列方向处理总是会比沿行方向慢一个档次，解决方式有很多，例如先对中间结果进行转置，然后再按照行方向的规则进行处理，处理完后在将数据转置回去。转置过程是有非常高效的处理方式的，借助于SSE以及Cache优化，能实现比原始两重for循环快4倍以上的效果。还有一种方式就正如opencv中列处理过程所示，这正是下面将要重点描述的。
在opencv的代码中，并没有直接沿列方向处理，而是继续沿着行的方向一行一行处理，我先贴出我自己翻译的有关纯C的代码进行解说：
for (Y = 0; Y & Size - 1; Y++)
注意没有最后一项哦
int *LinePS = (int *)(Sum-&Data + ColPos[Y] * Sum-&WidthStep);
for(X = 0; X & W X++)
ColSum[X] += LinePS[X];
for (Y = 0; Y & H Y++)
unsigned char* LinePD
= Dest-&Data + Y * Dest-&WidthS
int *AddPos
= (int*)(Sum-&Data + ColPos[Y + Size - 1] * Sum-&WidthStep);
int *SubPos
= (int*)(Sum-&Data + ColPos[Y] * Sum-&WidthStep);
for(X = 0; X & W X++)
Value = ColSum[X] + AddPos[X];
LinePD[X] = Value * S
ColSum[X] = Value - SubPos[X];
上述代码中定义了一个ColSum用于保存每行某个位置处在列方向上指定半径内各中间元素的累加值，对于第一行，做特殊处理，其他行的每个元素的处理方式和BaseRowFilter里的处理方式很像，只是在书写上有所区别，特别注意的是对第一行的累加时，最后一个元素并没有计算在内，这个处理技巧在下面的X循环里有体现，请大家仔细体味下。
上述代码这样做的好处是，能有效的减少CPU的cache miss，但是总的计算量是没有改变的，因此能有效的提高速度。
针对PC，在opencv内部，其在列方向上采用了SSE进行进一步的优化，我们贴出其处理uchar类型的代码：
代码比较多，我稍微精简下（注意，精简后的是不可运行的，只是更清晰的表达了思路）。
virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
bool haveScale = scale != 1;
double _scale =
if( width != (int)sum.size() )
sum.resize(width);
sumCount = 0;
SUM = ∑[0];
if( sumCount == 0 )
memset((void*)SUM, 0, width*sizeof(int));
for( ; sumCount & ksize - 1; sumCount++, src++ )
const int* Sp = (const int*)src[0];
for( ; i &= width-4; i+=4 )
__m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
__m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
_mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
for( ; i & i++ )
SUM[i] += Sp[i];
src += ksize-1;
for( ; count--; src++ )
const int* Sp = (const int*)src[0];
const int* Sm = (const int*)src[1-ksize];
uchar* D = (uchar*)
const __m128 scale4 = _mm_set1_ps((float)_scale);
for( ; i &= width-8; i+=8 )
__m128i _sm
= _mm_loadu_si128((const __m128i*)(Sm+i));
__m128i _sm1
= _mm_loadu_si128((const __m128i*)(Sm+i+4));
__m128i _s0
= _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
_mm_loadu_si128((const __m128i*)(Sp+i)));
__m128i _s01
= _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
_mm_loadu_si128((const __m128i*)(Sp+i+4)));
__m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
__m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));
_s0T = _mm_packs_epi32(_s0T, _s0T1);
_mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));
_mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
_mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
for( ; i & i++ )
int s0 = SUM[i] + Sp[i];
D[i] = saturate_cast&uchar&(s0*_scale);
SUM[i] = s0 - Sm[i];
在行方向上，ColSum每个元素的更新相互之间是没有任何关系的，因此，借助于SSE可以实现一次性的四个元素的更新，并且上述代码还将第一行的特殊计算也用SSE实现了，虽然这部分计算量是非常小的。
在具体的SSE细节上，对于uchar类型，虽然中间结果是用int类型表达的，但是由于SSE没有整形的除法指令（是否用？），因此上面借用了浮点数的乘法SSE指令实现，当然就多了_mm_cvtepi32_ps 以及&_mm_cvtps_epi32这样的类型转换的函数。如果有整形除法，那就能更好了。
SSE的实现，无非也就是用_mm_loadu_si128加载数据，用_mm_add_epi32，&_mm_mul_ps之类的函数进行基本的运算，用_mm_storeu_si128来保存数据，并没有什么特别复杂的部分，注意到考虑到数据的普遍性，不一定都是16字节对齐的，因此在加载和保存是需要使用u方面的函数，其实现在的_mm_loadu_si128和_mm_load_si128在处理速度上已经没有特别明显的区别了。
注意到在每个SSE代码后面，总还有部分C代码，这是因为我们实际数据宽度并不一定是4的整数倍，未被SSE处理的部分必须使用C实现，这其实对读者来说是非常好的事情，因为我们能从这部分C代码中搞明白上述的SSE代码是干什么的。
以上就是opencv的Box Filter实现的关键代码，如果读者去看具体细节，opencv还有针对手机上的Neon优化，其实这个和SSE的意思基本是一样的。
那是否还有改进的空间呢，从我的认知中，在对垂直方向的处理上，应该没有了，但是水平方向呢， SSE有没有用武之地，经过我的思考我认为还是有的。
在行方向的计算中，这个for循环是主要的计算。
for(X = 1; X & W X ++)
Value += RowData[X + Size - 1] - RowData[X - 1];
LinePD[X] = V
本人认为虽然这里的操作是前后依赖的，全局无法并行化，但是观察这一行代码：
Value += RowData[X + Size - 1] - RowData[X - 1];
其中RowData[X + Size - 1] - RowData[X - 1]; 并不是前后依赖的，是可以并行的，因此如果提前快速的计算出所有的差值，那么提速的空间还比较大，即需要提前计算出下面的函数：
for(X = 0; X & (Width - 1); X++)
Diff[X] = AddPos[X] - SubPos[X];
这里用SSE实现则非常简单了：
unsigned char *AddPos = RowData + Size * C
unsigned char *SubPos = RowD
注意这个赋值在下面的循环外部，这可以避免当Width&8时第二个for循环循环变量未初始化
__m128i Zero = _mm_setzero_si128();
for (; X &= (Width - 1) * Channel - 8; X += 8)
__m128i Add = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(AddPos + X)), Zero);
__m128i Sub = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(SubPos + X)), Zero);
_mm_store_si128((__m128i *)(Diff + X + 0), _mm_sub_epi32(_mm_unpacklo_epi16(Add, Zero), _mm_unpacklo_epi16(Sub, Zero)));
由于采用了_aligned_malloc函数分配内存，可是使用_mm_store_si128
_mm_store_si128((__m128i *)(Diff + X + 4), _mm_sub_epi32(_mm_unpackhi_epi16(Add, Zero), _mm_unpackhi_epi16(Sub, Zero)));
for(; X & (Width - 1) * C X++)
Diff[X] = AddPos[X] - SubPos[X];
和列方向的SSE处理代码不同的是，这里加载的是uchar数据，因此加载的函数就有所不同，处理的方式也有区别，上面几个SSE函数各位查查MSDN就能理解其意思，还是很有味道的。
经过这样的处理，经过测试发现，速度能够比opencv的版本在提高30%，也是额外的惊喜。
再有的优化可能提速有限，比如把剩下的一些for循环内部分成四路等等。
在我的I3笔记本中，使用上述算法对的三通道彩色图像进行半径为5的测试，进行100次的计算纯算法部分的耗时为800ms，如果是纯C版本大概为1800ms；当半径为200时，SSE版本约为950ms， C版本约1900ms；当半径为400时，SSE版本约为1000ms， C版本约2100ms；可见，半径增大，耗时稍有增加，这主要是由于算法中有部分初始化的计算和半径有关，但是这些都是不重要的。
在不使用多线程（虽然本算法非常适合多线程计算），不使用GPU，只使用单线程\CPU进行计算的情况下， 个人觉得目前我这个代码是很难有质的超越的， 从某个方面讲，代码中的用于计算的时间占用的时间比从内存等待数据的时间可能还要短，而似乎也没有更好的算法能进一步减少内存数据的访问量。
本人在这里邀请各位高手对目前我优化的这个代码进行进一步的优化，希望高人不要谦虚。
源代码下载地址（VS2010编写）：Box Filter
运行界面：
本文的代码是针对常用的图像数据进行的优化处理，在很多场合下，需要对int或者float类型进行处理，比如GuideFilter，如果读者理解了本文解读的代码的原理，更改代码以便适应他们则是一件非常简单的事情，如果您不会，那么也请你不要给我留言，或者我可以有偿提供，因为从本质上讲我喜欢提供渔，而不是鱼，渔具已经送给你了，你却找我要鱼，那只能..............。
****************************作者： laviewpbt & 时间：
& &联系QQ: & 转载请保留本行信息**********************
已发表评论数()
请填写推刊名
描述不能大于100个字符!
权限设置： 公开
仅自己可见
正文不准确
标题不准确
排版有问题
主题不准确
没有分页内容
图片无法显示
视频无法显示
与原文不一致}