blob特征检测原理与实现

最新推荐文章于 2024-09-26 16:21:42 发布

南山牧笛

最新推荐文章于 2024-09-26 16:21:42 发布

阅读量1.2w

点赞数 4

分类专栏：计算机视觉

计算机视觉专栏收录该内容

50 篇文章 0 订阅

订阅专栏

计算机视觉讨论群162501053 转载请注明：http://blog.csdn.net/abcd1992719g/article/details/27071273
收入囊中理解blob特征会利用OpenCV API提取blob特征自己实现blob特征检测首先要了解，什么是blob特征，我们来看下面两幅图片。

直观上来看，blob特征就是一团，一坨东西，它并不一定是圆形的，总之它就是那么一团独立存在的特征。

葵花宝典我们看待世界万物的特征跟物体的大小（scale）有很大的关系。假如当前相机镜头能清楚看到一个人，相机往后移动，那么就能清楚看到一栋建筑，再往后，就是城市，国家，地球，宇宙...... 因此，我们必须对不同的scale-space进行一系列的研究，才能挖掘出特征。如下图：这样，我们就达到了一个covariant的方法,如下图。

什么是scale-space?如下图。通过取不同的系数sigma和不同的核size,我们得到不同程序高斯模糊的图片，组成了scale-space.

scale-space有很多的创建方法。在本文的实现中，我采用了log scale-space

 
   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<double> sigma = create_sigma(1.0, 0.1, 3.0); 就生成了21个sigma 高斯窗口大小 ksize = ceil(sigma[i]*3)*2+1;
这里的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算子模版

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的函数接近的比较好，看一下我的实现。

我依据的公式，这也有多种不同的格式

 
   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算子时，要注意

The response of a derivative of Gaussianfilter to a perfect step edge decreases asσincreases
To keep response the same (scale-invariant),must multiply Gaussian derivative byσ
Laplacian is the second Gaussian derivative,so it must be multiplied byσ2

为了使得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 &parameters = SimpleBlobDetector::Params());

30.

31. protected:

32. ...

33. };

OpenCV实现的算法如下：
1. 对[minThreshold,maxThreshold)区间，以thresholdStep为间隔，做多次二值化。
2. 对每张二值图片，使用findContours()提取连通域并计算每一个连通域的中心。
3. 根据2得到的中心，全部放在一起。一些很接近的点［由theminDistBetweenBlobs控制多少才算接近］被归为一个group,对应一个bolb特征..
4. 从3得到的那些点,估计最后的blob特征和相应半径，并以key points返回。
  - By color. This filter compares the intensity of a binary image at the center of a blob to blobColor. If they differ, the blob is filtered out. UseblobColor = 0 to extract dark blobs and blobColor = 255 to extract light blobs.
  - By area. Extracted blobs have an area between minArea (inclusive) and maxArea (exclusive).
  - By circularity. Extracted blobs have circularity ( $rac{4*pi*Area}{perimeter * perimeter}$ ) between minCircularity (inclusive) and maxCircularity(exclusive).
  - By ratio of the minimum inertia to maximum inertia. Extracted blobs have this ratio between minInertiaRatio (inclusive) andmaxInertiaRatio (exclusive).
  - By convexity. Extracted blobs have convexity (area / area of blob convex hull) between minConvexity (inclusive) and maxConvexity(exclusive).
  你可以设置提取特征的方法，有上面5个选项。默认提取黑色圆形的blob特征。下面是一个示例用法。
  
  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 <iostream>
  
  06. #include <vector>
  
  07. using namespace cv;
  
  08. using namespace std;
  
  09.
  
  10. int main( int, char** argv )
  
  11. {
  
  12. Mat image = imread( argv[1] );
  
  13. vector<KeyPoint> keypoints;
  
  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. vector<Mat> create_scale_space(Mat &im_in, vector<double> &sigma)
  
  02. {
  
  03. vector<Mat> scale_space;
  
  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<struct blob> detect_blobs(Mat &im_in, double thresh, vector<double> et, int padding = 10)
  
  02. {
  
  03. vector<struct blob> 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<Mat> all_ims = create_scale_space(norm, et);
  
  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.at<unsigned char>(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<struct blob> prune_blobs(vector<struct blob> blobs_in)
  
  02. {
  
  03. vector<bool>keep(blobs_in.size(), true);
  
  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<struct blob> 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 <cmath>
  
  004. #include <iostream>
  
  005. #include <vector>
  
  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. vector<Mat> create_scale_space(Mat &im_in, vector<double> &sigma)
  
  053. {
  
  054. vector<Mat> scale_space;
  
  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<struct blob> detect_blobs(Mat &im_in, double thresh, vector<double> et, int padding = 10)
  
  081. {
  
  082. vector<struct blob> 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<Mat> all_ims = create_scale_space(norm, et);
  
  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.at<unsigned char>(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<struct blob> prune_blobs(vector<struct blob> blobs_in)
  
  116. {
  
  117. vector<bool>keep(blobs_in.size(), true);
  
  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<struct blob> 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 = imread( argv[1] );
  
  168. cvtColor( src, gray, CV_RGB2GRAY );
  
  169. vector<double> sigma = create_sigma(1.0, 0.2, 3.0);
  
  170. vector<struct blob> blobs = detect_blobs(gray, 0.2, sigma);
  
  171. vector<struct blob> 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. }