收入囊中 理解blob特征会利用OpenCV API提取blob特征自己实现blob特征检测 首先要了解，什么是blob特征，我们来看下面两幅图片。
文章图片
文章图片
直观上来看，blob特征就是一团，一坨东西，它并不一定是圆形的，总之它就是那么一团独立存在的特征。
葵花宝典 我们看待世界万物的特征跟物体的大小（scale）有很大的关系。假如当前相机镜头能清楚看到一个人，相机往后移动，那么就能清楚看到一栋建筑，再往后，就是城市，国家，地球，宇宙...... 因此，我们必须对不同的scale-space进行一系列的研究，才能挖掘出特征。如下图：这样，我们就达到了一个covariant的方法,如下图。
文章图片
什么是scale-space?如下图。通过取不同的系数sigma和不同的核size,我们得到不同程序高斯模糊的图片，组成了scale-space.
文章图片
scale-space有很多的创建方法。 在本文的实现中，我采用了log scale-space view source print ? 01.
vector<double> create_sigma(double start, double step, double end) 02.
{ 03.
vector<double> sigma;
04.
while(start <= end+1e-8){ 05.
double s = exp(start);
06.
sigma.push_back(s);
07.
start += step;
08.
} 09.
return sigma;
10.
} 如果我这样调用 vector
这里的3是一个比较好的取值，原因我会在下面分析（个人理解）.
文章图片
左图是原图，下面是利用上面的sigma建立出来的scale-space 
文章图片
下一步，为了得到特征，我们对scale-space作用拉普拉斯算子。拉普拉斯梯度在blob的中心会取得最大值，如下图，假设拉普拉斯窗口半径介于1－2之间，那么在下面第4张图片，拉普拉斯窗口覆盖了整个短信号，使得梯度很大，因此，只要用合适的窗口在合适的scale-space操作，就能得到一系列极大值。
文章图片
下面是经过拉普拉斯处理后的图片。
文章图片
看到上面的图片，最右下叫最明显，因为sigma比较大，最后一张图片能找到比较大的blob特征，那些特别泛白的地方就是。再看sigma = 9.97的图片，sigma = 20.09中白色的那块在sigma = 9.97中已经不白了，在sigma = 9.97中能找到的白色是特征比较适中，不大不小的。 经过我上面的解释，下一步我们要做什么你应该已经知道了。就是遍历我们的scale-space,找到每一个局部最大值点（泛白），存储下来备选为我们的blob特征，blob特征有自己的半径，我们取r = 1.5*sigma. 我想过在上面我们取的窗口大小半径是3倍的sigma,这里取1.5大概就是这个原因。［其他地方有更专业的解释，拉普拉斯是(
x
2+
y
2?2σ2)
e
?(
x
2+
y
2)/2σ2 最大响应在σ =
r
/ 根号2
］ 我们找到了备选的blob还得做删除工作，因为两个blob可能重叠的比较厉害，因此我们需要把重叠面积比较大的删去一个，做完这一步，我们就得到了最后的所有blob特征。
之前提到 我们一开始先进行高斯模糊，再做拉普拉斯。其实这两步可以一起做，叫做LOG算子 我在http://blog.csdn.net/abcd1992719g/article/details/25559527有稍微介绍
在matlab里，下面一句可以很简单得到ksize = 3, sigma = 1.0的LOG算子模版 view source print ? 1.
lap_gauss = fspecial('log', 3, 1.0) lap_gauss =
0.1004 -0.0234 0.1004
-0.0234 -0.3079 -0.0234
0.1004 -0.0234 0.1004
但是在OpenCV就没那么方便了，下面是我自己实现的一个 LOG算子模版，跟matlab的函数接近的比较好，看一下我的实现。
文章图片
我依据的公式，这也有多种不同的格式
view source print ? 01.
Mat create_LOG(int size, double sigma) 02.
{ 03.
Mat H(size, size, CV_64F);
04.
int cx = (size-1)/2;
05.
int cy = (size-1)/2;
06.
double sum = 0;
07.
for(int i = 0;
i < H.rows;
i++) 08.
{ 09.
for(int j = 0;
j < H.cols;
j++) 10.
{ 11.
int nx = i - cx;
12.
int ny = j - cy;
13.
double s = -1/(3.1415926 * sigma*sigma*sigma*sigma);
14.
s = s* (1 - (nx*nx + ny*ny)/(2*sigma*sigma));
15.
s = s*exp(-(nx*nx+ny*ny)/(2*sigma*sigma));
16.
sum += s;
17.
H.at<double>(i,j) = s;
18.
} 19.
} 20. 21.
double mean = sum/(size*size);
22. 23.
for(int i = 0;
i < H.rows;
i++) 24.
{ 25.
for(int j = 0;
j < H.cols;
j++) 26.
{ 27.
H.at<double>(i,j) -= mean;
28.
} 29.
} 30. 31.
return H;
32.
}
但是利用LOG算子时，要注意
文章图片
为了使得response一致，我们对做卷积的图片要先乘上sigma^2，这一步叫做scale normalize,我们简单看下做不做normalize的对比。 
文章图片
初识API OpenCV提供了Simpleblobdetector的接口
view source print ? 01.
class SimpleBlobDetector : public FeatureDetector 02.
{ 03.
public: 04.
struct Params 05.
{ 06.
Params();
07.
float thresholdStep;
08.
float minThreshold;
09.
float maxThreshold;
10.
size_t minRepeatability;
11.
float minDistBetweenBlobs;
12. 13.
bool filterByColor;
14.
uchar blobColor;
15. 16.
bool filterByArea;
17.
float minArea, maxArea;
18. 19.
bool filterByCircularity;
20.
float minCircularity, maxCircularity;
21. 22.
bool filterByInertia;
23.
float minInertiaRatio, maxInertiaRatio;
24. 25.
bool filterByConvexity;
26.
float minConvexity, maxConvexity;
27.
};
28. 29.
SimpleBlobDetector(const SimpleBlobDetector::Params ¶meters = SimpleBlobDetector::Params());
30. 31.
protected: 32.
... 33.
};
OpenCV实现的算法如下：
【opencv——simpleblob】 你可以设置提取特征的方法，有上面5个选项。默认提取黑色圆形的blob特征。下面是一个示例用法。
文章图片
) between minCircularity (inclusive) and maxCircularity(exclusive).
view source print ? 01.
#include "opencv2/core/core.hpp" 02.
#include "opencv2/features2d/features2d.hpp" 03.
#include "opencv2/highgui/highgui.hpp" 04.
#include "opencv2/nonfree/features2d.hpp" 05.
#include 06.
#include 07.
using namespace cv;
08.
using namespace std;
09. 10.
int main( int, char** argv ) 11.
{ 12.
Mat image = imread( argv[1] );
13.
vector 14.
SimpleBlobDetector::Params params;
15.
params.filterByArea = true;
16.
params.minArea = 10;
17.
params.maxArea = 10000;
18. 19.
SimpleBlobDetector blobDetector( params );
20.
blobDetector.create("SimpleBlob");
21.
blobDetector.detect( image, keypoints );
22.
drawKeypoints(image, keypoints, image, Scalar(255,0,0));
23. 24.
namedWindow("result", 1);
25.
imshow("result", image);
26.
waitKey();
27.
return 0;
28.
}
荷枪实弹 我们来实现自己的blob特征提取～ 这是效果图。
文章图片
下面是我们的第一个函数，已经见过了，生成log sigma 而不是 linear sigma view source print ? 01.
vector<double> create_sigma(double start, double step, double end) 02.
{ 03.
vector<double> sigma;
04.
while(start <= end+1e-8){ 05.
double s = exp(start);
06.
sigma.push_back(s);
07.
start += step;
08.
} 09.
return sigma;
10.
}
这是第二个函数我们也已经见过了，生成我们的LOG算子模版 view source print ? 01.
Mat create_LOG(int size, double sigma) 02.
{ 03.
Mat H(size, size, CV_64F);
04.
int cx = (size-1)/2;
05.
int cy = (size-1)/2;
06.
double sum = 0;
07.
for(int i = 0;
i < H.rows;
i++) 08.
{ 09.
for(int j = 0;
j < H.cols;
j++) 10.
{ 11.
int nx = i - cx;
12.
int ny = j - cy;
13.
double s = -1/(3.1415926 * sigma*sigma*sigma*sigma);
14.
s = s* (1 - (nx*nx + ny*ny)/(2*sigma*sigma));
15.
s = s*exp(-(nx*nx+ny*ny)/(2*sigma*sigma));
16.
sum += s;
17.
H.at<double>(i,j) = s;
18.
} 19.
} 20. 21.
double mean = sum/(size*size);
22. 23.
for(int i = 0;
i < H.rows;
i++) 24.
{ 25.
for(int j = 0;
j < H.cols;
j++) 26.
{ 27.
H.at<double>(i,j) -= mean;
28.
} 29.
} 30. 31.
return H;
32.
}
这是我们的blob结构体，x,y是坐标，sigma用于计算半径r=1.5*sigma,intensity用于当两个 blob重叠度高时取intensity大的（梯度大）
view source print ? 1.
struct blob 2.
{ 3.
int x;
4.
int y;
5.
double sigma;
6.
int intensity;
7.
};
第三个函数，生成LOG卷积后的scale_space view source print ? 01.
vectordouble> &sigma) 02.
{ 03.
vector 04. 05.
for(int i = 0;
i < sigma.size();
i++) 06.
{ 07.
int n = ceil(sigma[i]*3)*2+1;
08.
Mat LOG = create_LOG(n, sigma[i]);
09.
Mat dst;
10. 11.
filter2D(im_in.mul(sigma[i]*sigma[i]), dst, -1 , LOG, Point(-1,-1));
12.
scale_space.push_back(dst);
13.
} 14. 15.
return scale_space;
16.
}
第四个函数，检测所有可能的blob view source print ? 01.
vector detect_blobs(Mat &im_in, double thresh, vector<double> et, int padding = 10) 02.
{ 03.
vector blobs;
04.
Mat addborder,norm;
05.
copyMakeBorder(im_in, addborder, padding, padding, padding, padding, BORDER_CONSTANT, Scalar(255));
06.
normalize(addborder, norm, 1, 0, NORM_MINMAX, CV_64F);
07.
vector 08. 09.
for(int i = 0;
i < et.size();
i++) 10.
{ 11.
Mat grayscale,mx;
12.
normalize( all_ims[i], grayscale, 0, 255, NORM_MINMAX, CV_8U, Mat() );
13.
Mat structedElement(3, 3, CV_8U, Scalar(1));
14.
dilate(grayscale, mx, structedElement);
15.
grayscale = ( grayscale == mx) & (all_ims[i] > thresh);
//取大于threshold并且是周围最大 16.
for(int j = 0;
j < norm.rows;
j++) 17.
{ 18.
for(int k = 0;
k < norm.cols;
k++) 19.
{ 20.
if(grayscale.atchar>(j,k) > 0) 21.
{ 22.
struct blob b;
23.
b.x = j - padding;
24.
b.y = k - padding;
25.
b.sigma = et[i];
26.
b.intensity = all_ims[i].at<double>(j,k);
27.
blobs.push_back(b);
28.
} 29.
} 30.
} 31.
} 32. 33.
return blobs;
34.
}
第五个函数，删除重叠度大的blob view source print ? 01.
vector prune_blobs(vector blobs_in) 02.
{ 03.
vectortrue);
04. 05.
for(int i = 0;
i < blobs_in.size();
i++) 06.
{ 07.
for(int j = 0;
j < blobs_in.size();
j++) 08.
{ 09.
if ( (i==j) || (keep[i] == false) || (keep[j] == false) ) 10.
continue;
11. 12.
struct blob blobi = blobs_in[i];
13.
struct blob blobj = blobs_in[j];
14. 15.
int xi = blobi.x;
16.
int yi = blobi.y;
17.
double ri = blobi.sigma;
18.
int xj = blobj.x;
19.
int yj = blobj.y;
20.
double rj = blobj.sigma;
21. 22.
double d = sqrt((xj-xi)*(xj-xi) + (yj-yi)*(yj-yi));
23.
double rirj = ri+rj;
24.
double overlap_percent = rirj/d;
25.
if (overlap_percent > 2.0) 26.
{ 27.
if (blobi.intensity > blobj.intensity) 28.
{ 29.
keep[i] = true;
30.
keep[j] = false;
31.
} else { 32.
keep[i] = false;
33.
keep[j] = true;
34.
} 35.
} 36.
} 37.
} 38. 39. 40.
vector blobs_out;
41. 42.
for(int i = 0;
i < blobs_in.size();
i++) 43.
if (keep[i]) 44.
blobs_out.push_back(blobs_in[i]);
45. 46.
return blobs_out;
47.
}
完整代码 view source print ? 001.
#include "opencv2/highgui/highgui.hpp"002.
#include "opencv2/imgproc/imgproc.hpp"003.
#include 004.
#include 005.
#include 006.
using namespace cv;
007.
using namespace std;
008. 009.
const double pi = 3.1415926f;
010. 011.
struct blob 012.
{ 013.
int x;
014.
int y;
015.
double sigma;
016.
int intensity;
017.
};
018. 019.
Mat create_LOG(int size, double sigma) 020.
{ 021.
Mat H(size, size, CV_64F);
022.
int cx = (size-1)/2;
023.
int cy = (size-1)/2;
024.
double sum = 0;
025.
for(int i = 0;
i < H.rows;
i++) 026.
{ 027.
for(int j = 0;
j < H.cols;
j++) 028.
{ 029.
int nx = i - cx;
030.
int ny = j - cy;
031.
double s = -1/(3.1415926 * sigma*sigma*sigma*sigma);
032.
s = s* (1 - (nx*nx + ny*ny)/(2*sigma*sigma));
033.
s = s*exp(-(nx*nx+ny*ny)/(2*sigma*sigma));
034.
sum += s;
035.
H.at<double>(i,j) = s;
036.
} 037.
} 038. 039.
double mean = sum/(size*size);
040. 041.
for(int i = 0;
i < H.rows;
i++) 042.
{ 043.
for(int j = 0;
j < H.cols;
j++) 044.
{ 045.
H.at<double>(i,j) -= mean;
046.
} 047.
} 048. 049.
return H;
050.
} 051. 052.
vectordouble> &sigma) 053.
{ 054.
vector 055. 056.
for(int i = 0;
i < sigma.size();
i++) 057.
{ 058.
int n = ceil(sigma[i]*3)*2+1;
059.
Mat LOG = create_LOG(n, sigma[i]);
060.
Mat dst;
061. 062.
filter2D(im_in.mul(sigma[i]*sigma[i]), dst, -1 , LOG, Point(-1,-1));
063.
scale_space.push_back(dst);
064.
} 065. 066.
return scale_space;
067.
} 068. 069.
vector<double> create_sigma(double start, double step, double end) 070.
{ 071.
vector<double> sigma;
072.
while(start <= end+1e-8){ 073.
double s = exp(start);
074.
sigma.push_back(s);
075.
start += step;
076.
} 077.
return sigma;
078.
} 079. 080.
vector detect_blobs(Mat &im_in, double thresh, vector<double> et, int padding = 10) 081.
{ 082.
vector blobs;
083.
Mat addborder,norm;
084.
copyMakeBorder(im_in, addborder, padding, padding, padding, padding, BORDER_CONSTANT, Scalar(255));
085.
normalize(addborder, norm, 1, 0, NORM_MINMAX, CV_64F);
086.
vector 087. 088.
for(int i = 0;
i < et.size();
i++) 089.
{ 090.
Mat grayscale,mx;
091.
normalize( all_ims[i], grayscale, 0, 255, NORM_MINMAX, CV_8U, Mat() );
092.
Mat structedElement(3, 3, CV_8U, Scalar(1));
093.
dilate(grayscale, mx, structedElement);
094.
grayscale = ( grayscale == mx) & (all_ims[i] > thresh);
095.
for(int j = 0;
j < norm.rows;
j++) 096.
{ 097.
for(int k = 0;
k < norm.cols;
k++) 098.
{ 099.
if(grayscale.atchar>(j,k) > 0) 100.
{ 101.
struct blob b;
102.
b.x = j - padding;
103.
b.y = k - padding;
104.
b.sigma = et[i];
105.
b.intensity = all_ims[i].at<double>(j,k);
106.
blobs.push_back(b);
107.
} 108.
} 109.
} 110.
} 111. 112.
return blobs;
113.
} 114. 115.
vector prune_blobs(vector blobs_in) 116.
{ 117.
vectortrue);
118. 119.
for(int i = 0;
i < blobs_in.size();
i++) 120.
{ 121.
for(int j = 0;
j < blobs_in.size();
j++) 122.
{ 123.
if ( (i==j) || (keep[i] == false) || (keep[j] == false) ) 124.
continue;
125. 126.
struct blob blobi = blobs_in[i];
127.
struct blob blobj = blobs_in[j];
128. 129.
int xi = blobi.x;
130.
int yi = blobi.y;
131.
double ri = blobi.sigma;
132.
int xj = blobj.x;
133.
int yj = blobj.y;
134.
double rj = blobj.sigma;
135. 136.
double d = sqrt((xj-xi)*(xj-xi) + (yj-yi)*(yj-yi));
137.
double rirj = ri+rj;
138.
double overlap_percent = rirj/d;
139.
if (overlap_percent > 2.0) 140.
{ 141.
if (blobi.intensity > blobj.intensity) 142.
{ 143.
keep[i] = true;
144.
keep[j] = false;
145.
} else { 146.
keep[i] = false;
147.
keep[j] = true;
148.
} 149.
} 150.
} 151.
} 152. 153. 154.
vector blobs_out;
155. 156.
for(int i = 0;
i < blobs_in.size();
i++) 157.
if (keep[i]) 158.
blobs_out.push_back(blobs_in[i]);
159. 160.
return blobs_out;
161.
} 162. 163.
int main( int, char** argv ) 164.
{ 165.
Mat src,gray;
166. 167.
src = https://www.it610.com/article/imread( argv[1] );
168.
cvtColor( src, gray, CV_RGB2GRAY );
169.
vector<double> sigma = create_sigma(1.0, 0.2, 3.0);
170.
vector blobs = detect_blobs(gray, 0.2, sigma);
171.
vector blobs_trimmed = prune_blobs(blobs);
172. 173.
for(int i = 0;
i < blobs_trimmed.size();
i++) 174.
{ 175.
struct blob b = blobs_trimmed[i];
176.
circle( src, Point( b.y, b.x ), 1.5*b.sigma,Scalar(0), 2, 8, 0 );
177.
} 178. 179.
namedWindow("result", 1);
180.
imshow("result", src);
181.
waitKey();
182.
return 0;
183.
}
举一反三
文章图片
DoG是LOG算子的近似实现，速度更快 有兴趣的同学可以用DoG去实现，然后告诉我效果怎么样！！！
福利 matlab版本 http://www.cs.utah.edu/~jfishbau/advimproc/project1/code/
推荐阅读
- opencv|图像处理之椒盐噪声的添加与去除
- OpenCV|【OpenCV 完整例程】89. 带阻滤波器的传递函数
- OpenCV|【OpenCV 完整例程】90. 频率域陷波滤波器
- OpenCV|【OpenCV 完整例程】22. 图像添加非中文文字
- OpenCV|【OpenCV 完整例程】91. 高斯噪声、瑞利噪声、爱尔兰噪声
- 灰度世界算法（Gray World Algorithm）和White Patch Retinex算法
- 安卓开发|通过RenderScript 实现 NV21转Bitmap、两张Bitmap按照透明度混合的工具类
- 基于SLIC的细胞分割与检测
- 科学计算|超像素分割
- MFC|gdal 图像金字塔