·主要思路:
斑点检测主要使检测出图像中比它周围像素灰度值大或比周围灰度值小的区域。
一般有两种方法:
1.基于求导的微分方法,也叫微分检测器(LOG算子)
2.基于局部极值的分水岭算法(SimpleBlobDetector斑点检测算子)
·LOG斑点检测
LOG算子–高斯拉普拉斯算子(Laplace of Gaussian,LOG)
·对于二维高斯函数:
exp():
·它的拉普拉斯变换为:
·规范化后:
如上图所示,规范化后算子在二维图像上使一个圆对称函数,可以用该算子去检测图像中的斑点。同时也可以通过改变σ的值,检测不同尺寸的斑点。
·LOG原理解析:
为什么LOG算子可以检测图像中斑点?
图像与某一个二维函数进行卷积运算实际就是求取图像与这一函数的相似性。
同理,图像与高斯拉普拉斯的卷积实际就是求取图像与高斯拉普拉斯函数的相似性。
当图像中斑点尺寸与函数形状趋向一致时,响应最大。
假设原图像是一个与位置有关的随机变量X的密度函数,LOG则是随机变量Y的密度函数。
则随机变量X+Y的密度分布函数时两个函数的卷积形式。
如果想让X+Y能取得最大值,则X和Y的步调越接近越好,即X上升,Y也上升。
Laplace可以用来检测图像中的局部极值点,但是对噪声敏感。所以在卷积之前,可以用高斯滤波对图像进行卷积。去除图像中的噪点:
先对图像f(x,y)用方差为σ的高斯核做高斯滤波:
然后图像的拉普拉斯图像即为:
·LOG算子实现:
根据以上公式,可以用openCV实现LOG算子的模板
//创建LOG算子模板
Mat create_LOG(int size, double sigma) {
Mat H(size, size, CV_64F);
//中心点位置
int cx = (size - 1) / 2;
int cy = (size - 1) / 2;
double sum = 0;
for (int i = 0; i < H.rows; i++)
{
for (int j = 0; j < H.cols; j++)
{
//LOG公式拆解
int nx = i - cx;//获得对应位置的x值
int ny = j - cy;//获得对应位置的y值
double s = -1 / (3.1415926 * sigma * sigma * sigma * sigma);
s = s * (1 - (nx * nx + ny * ny) / (2 * sigma * sigma));
s = s * exp(-(nx * nx + ny * ny) / (2 * sigma * sigma));
sum += s;
//将对应点值更改
H.at<double>(i, j) = s;
}
}
double mean = sum / (size * size);
for (int i = 0; i < H.rows; i++)
{
for (int j = 0; j < H.cols; j++) {
//每个点最终结果都要减去平均值
H.at<double>(i, j) -= mean;
}
}
return H;
}
·多尺度检测
当σ 尺度一定是,只能检测对应半径的斑点。具体检测的斑点半径,可以通过规范化的二维拉普拉斯算子求导。
规范化的高斯拉普拉斯函数为:
求上式的极值点等价于求下式:
得:
因此可知,在尺度σ=r/√2时,响应值达到最大。而将相应达到峰值时的尺度σ值,称为特征尺度。
在多尺度的情况下,同时在空间核尺度上达到最值便是半袋按。即对应点的拉普拉斯响应值都大于或小于其他26个立方空间领域(9+8+9)的值。
·openCV中的斑点检测
openCV中检测Blobs的类为SimpleBlobDetector,其定义如下:
class SimpleBlobDetector : public FeatureDetector
{
public:
struct Params
{
Params();
float thresholdStep;
float minThreshold;
float maxThreshold;
size_t minRepeatability;
float minDistBetweenBlobs;
bool filterByColor;
uchar blobColor;
bool filterByArea;
float minArea, maxArea;
bool filterByCircularity;
float minCircularity, maxCircularity;
bool filterByInertia;
float minInertiaRatio, maxInertiaRatio;
bool filterByConvexity;
float minConvexity, maxConvexity;
};
SimpleBlobDetector(const SimpleBlobDetector::Params ¶meters = SimpleBlobDetector::Params());
protected:
...
};
openCV实现的算法如下:
1.对[minThreshold,maxThreshold]区间,以thresholdStep为步长,做多次二值化。
同时也可以限定检测对应阈值区间内的斑点
2.对每张二值图片,使用findContours()提取连通域并计算每一个连通域的中心。
3.提取2得到的点的中心,全部放在一起。一些接近的点被归为一个group,对应一个bolb特征
(多少算接近由minDistBetweenBlobs控制)
4.从3得到的点中,估计最后的blob特征和相应半径,并以key points返回
·By color:通过将斑点中心的二值图像强度与斑点颜色进行比较
如果不同,则会过滤掉斑点。可以使用blobColor=0提取深色斑点,255提取浅色斑点
·By area:将提取的点限定在[minArea-maxArea]
·By ratio of the minimum inertia to maximum inertia:提取的点在对应的惯性比区间内
·By convexity:凸性区间
实例:
int main(int argc,char* argv[])
{
Mat srcImage = imread("F:\\opencv\\OpenCVImage\\FeatureDetectSrc1.jpg");
Mat srcGrayImage;
if (srcImage.channels() == 3)
{
cvtColor(srcImage,srcGrayImage,CV_RGB2GRAY);
}
else
{
srcImage.copyTo(srcGrayImage);
}
vector<KeyPoint>detectKeyPoint;
Mat keyPointImage1,keyPointImage2;
SimpleBlobDetector::Params params;
params.minArea = 10;
params.maxArea = 1000;
params.filterByCircularity = true;
Ptr<SimpleBlobDetector> sbd = SimpleBlobDetector::create(params);
sbd->detect(srcGrayImage,detectKeyPoint);
drawKeypoints(srcImage,detectKeyPoint,keyPointImage1,Scalar(0,0,255),DrawMatchesFlags::DRAW_RICH_KEYPOINTS);
drawKeypoints(srcImage,detectKeyPoint,keyPointImage2,Scalar(0,0,255),DrawMatchesFlags::DEFAULT);
imshow("src image",srcImage);
imshow("keyPoint image1",keyPointImage1);
imshow("keyPoint image2",keyPointImage2);
imwrite("F:\\opencv\\OpenCVImage\\FeatureDetectSrc1SimpleBlobDetectorKeyPointImageDefault.jpg",keyPointImage2);
waitKey(0);
return 0;
}
·手动实现blob特征提取:
#include<opencv2/opencv.hpp>
#include<vector>
#include<iostream>
#include<stdlib.h>
using namespace std;
using namespace cv;
Mat image4;
//设置sigma模板,定义检测范围大小
vector<double> create_sigma(double start, double step, double end) {
vector<double> sigma;
while (start<=end+1e-8)
{
double s = exp(start);
sigma.push_back(s);//将每步数据导入
start += step;//每步指数递增
}
return sigma;
}
//创建LOG算子模板
Mat create_LOG(int size, double sigma) {
Mat H(size, size, CV_64F);
//中心点位置
int cx = (size - 1) / 2;
int cy = (size - 1) / 2;
double sum = 0;
for (int i = 0; i < H.rows; i++)
{
for (int j = 0; j < H.cols; j++)
{
//LOG公式拆解
int nx = i - cx;//获得对应位置的x值
int ny = j - cy;//获得对应位置的y值
double s = -1 / (3.1415926 * sigma * sigma * sigma * sigma);
s = s * (1 - (nx * nx + ny * ny) / (2 * sigma * sigma));
s = s * exp(-(nx * nx + ny * ny) / (2 * sigma * sigma));
sum += s;
//将对应点值更改
H.at<double>(i, j) = s;
}
}
double mean = sum / (size * size);
for (int i = 0; i < H.rows; i++)
{
for (int j = 0; j < H.cols; j++) {
//每个点最终结果都要减去平均值
H.at<double>(i, j) -= mean;
}
}
return H;
}
//blob结构体
struct blob
{
int x;
int y;
double sigma;//用于计算半径 r = 1.5*sigma;
int intensity;//用于当两个blob重叠度高时,取intensity(梯度)更大的
};
//生成LOG卷积后的scale_space
//通过取不同的系数sigma和不同的核size,会得到不同程度高斯模糊的图片,这些图片组成了scale-space
vector<Mat> create_scale_space(Mat& im_in, vector<double>& sigma) {
vector<Mat> scale_space;
for (int i = 0; i < sigma.size(); i++)
{
int n = ceil(sigma[i] * 3) * 2 + 1;//向上取整
Mat LOG = create_LOG(n, sigma[i]);
Mat dst;
//使用高斯拉普拉斯算子进行滤波运算
filter2D(im_in.mul(sigma[i] * sigma[i]), dst, -1, LOG, Point(-1, -1));//mul()--两个相同大小矩阵中元素对应相乘
//https://blog.csdn.net/dcrmg/article/details/52404580
scale_space.push_back(dst);
}
return scale_space;
}
//检测所有可能的blob
vector<struct blob> detect_blobs(Mat& im_in, double thresh, vector<double> et, int padding = 10) {
vector<struct blob> blobs;
Mat addborder, norm;
copyMakeBorder(im_in, addborder, padding, padding, padding, padding, BORDER_CONSTANT, Scalar(255));
//扩充图像边界的函数,让LOG算子可以处理边界
normalize(addborder, norm, 1, 0, NORM_MINMAX, CV_64F);//标准化,防止数据溢出
//生成空间
vector<Mat> all_ims = create_scale_space(norm, et);
//开始检测
for (int i = 0; i < et.size(); i++)
{
Mat grayscale, mx;
//灰度化
normalize(all_ims[i], grayscale, 0, 255, NORM_MINMAX, CV_8U, Mat());
//形态学操作
Mat structedElement(3, 3, CV_8U, Scalar(1));
dilate(grayscale, mx, structedElement);
grayscale = (grayscale == mx) & (all_ims[i] > thresh);//取大于threshold并且使周围最大的点
for (int j = 0; j < norm.rows; j++)
{
for (int k = 0; k < norm.cols; k++)
{
//符合条件即为极值,保存进队列中
if (grayscale.at<unsigned char>(j,k)>0)
{
struct blob b;
b.x = j - padding;
b.y = k - padding;
b.sigma = et[i];
b.intensity = all_ims[i].at<double>(j, k);
blobs.push_back(b);
}
}
}
}
return blobs;
}
//删除重叠度最大的blob
vector<struct blob> prune_blobs(vector<struct blob> blobs_in) {
vector<bool> keep(blobs_in.size(), true);
for (int i = 0; i < blobs_in.size(); i++)
{
for (int j = 0; j < blobs_in.size(); j++)
{
if ((i == j) || (keep[i] == false) || (keep[j] == false)) {
continue;
}
struct blob blobi = blobs_in[i];
struct blob blobj = blobs_in[j];
int xi = blobi.x;
int yi = blobi.y;
double ri = blobi.sigma;
int xj = blobj.x;
int yj = blobj.y;
double rj = blobj.sigma;
double d = sqrt((xj - xi) * (xj - xi) + (yj - yi) * (yj - yi));
double rirj = ri + rj;
double overlap_percent = rirj / d;
//如果超出一定范围,则判断为重合度高,社区一个
if (overlap_percent>2.0)
{
if (blobi.intensity>blobj.intensity)
{
keep[i] = true;
keep[j] = false;
}
else {
keep[i] = false;
keep[j] = true;
}
}
}
}
vector<struct blob> blobs_out;
for (int i = 0; i < blobs_in.size(); i++)
{
if (keep[i])
{
blobs_out.push_back(blobs_in[i]);
}
}
return blobs_out;
}
int main(int argc, char** argv) {
image4 = imread("C:/Users/18929/Desktop/博客项目/项目图片/sf.jpg");
Mat gray_image4;
cvtColor(image4, gray_image4, COLOR_BGR2GRAY);
vector<double> sigma = create_sigma(1.0, 0.2, 3.0);
vector<struct blob> blobs = detect_blobs(gray_image4, 0.2, sigma);
vector<struct blob> blobs_trimmed = prune_blobs(blobs);
for (int i = 0; i < blobs_trimmed.size(); i++)
{
struct blob b = blobs_trimmed[i];
circle(image4, Point(b.y, b.x), 1.5 * b.sigma, Scalar(0), 2, 8, 0);//虽然不知道为什么,但是这里Point的x,y得调换
}
imshow("result_img", image4);
waitKey(0);
return 0;
}