图像特征提取，斑点检测

斑点是数字图像的主要特征，是区域检测的一种特例，是许多特征生成、目标识别等方法的重要预处理环节。斑点通常和关键点(keypoint)，兴趣点(intrestpoint)以及特征点(featurepoint)表示同一个概念。

斑点通常指与周围有着颜色和灰度区别的区域。从远处看，一棵树是一个斑点，一块草地、一栋房子也可看成斑点。主要有LoG斑点，DoH斑点以及Gilles斑点。

LOG斑点检测

利用高斯拉普通拉斯（Laplace of Gaussian,LOG）算子检测图像斑点是一种十分常用的方法。

其实从更直观的角度去解释为什么LOG算子可以检测图像中的斑点是：

图像与某一个二维函数进行卷积运算实际就是求取图像与这一函数的相似性。

同理，图像与高斯拉普拉斯函数的卷积实际就是求取图像与高斯拉普拉斯函数的相似性。当图像中的斑点尺寸与高斯拉普拉斯函数的形状趋近一致时，图像的拉普拉斯响应达到最大。

从概率的角度解释为：假设原图像是一个与位置有关的随机变量X的密度函数，而LOG为随机变量Y的密度函数，则随机变量X+Y的密度分布函数即为两个函数的卷积形式。如果想让X+Y能取到最大值，则X与Y能保持步调一致最好，即X上升时，Y也上升，X最大时，Y也最大。

那么LOG算子是怎么被构想出来的呢？

事实上我们知道Laplace可以用来检测图像中的局部极值点，但是对噪声敏感，所以在我们对图像进行Laplace卷积之前，我们用一个高斯低通滤波对图像进行卷积，目标是去除图像中的噪声点。这一过程 可以描述为:
先对图像 f(x,y) 用方差为 $\sigma$ 的高斯核进行高斯滤波，去除图像中的噪点。

$$L(x,y;\sigma) = f(x,y) * G(x,y;\sigma)$$

然后对图像的拉普拉斯图像则为：

$$\Delta^2 = \frac{\partial^2L}{\partial x^2}+\frac{\partial^2L}{\partial y^2}$$

而实际上有下面等式：

$$\Delta^2[G(x,y)*f(x,y)] = \Delta^2[G(x,y)]*f(x,y)$$

所以，我们可以先求高斯核的拉普拉斯算子，再对图像进行卷积。

对于图像中二值化的圆形斑点，在尺度为$\sigma=r/\sqrt{2}$时，高斯拉普拉斯响应值达到最大。

同理，如果图像中的圆形斑点黑点反向，则在$\sigma=r/\sqrt{2}$时，高斯拉普拉斯响应值达到最小。将高斯拉普拉斯响应达到峰值时的尺度$\sigma$值，称为特征尺度。

如前所述，同时在空间和尺度上达到最大值（最小值）的点就是我们期望的斑点。对于计算二维图像I（x，y），计算图像在不同尺度下的离散拉普拉斯响应值，然后，检查位置空间中的每个点，如果该点的拉普拉斯响应值都大于或小于其他26个立方空间领域的值，那么，就是被检测的图像斑点。

LoG斑点的Matlab程序：

function [points]=LoG_Blob(img,num_blobs)
%功能：提取LoG斑点
%img——输入图像
%num——需要检测斑点数目
%point——检测出的斑点
img=double(img(:,:,1));
if nargin==1    %如果输入参数仅有一个（img）
num=120;    %则将检测斑点数设置为120
else
num=num_blobs;
end
%设定LoG参数
sigma_begin=2;
sigma_end=15;
sigma_step=1;
sigma_array=sigma_begin:sigma_step:sigma_end;
sigma_nb=numel(sigma_array);
%n = numel(A) returns the number of elements, n, in array A
%equivalent to prod(size(A)).
img_height=size(img,1);
img_width=size(img,2);
%计算尺度规范化高斯拉普拉斯算子
snlo=zeros(img_height,img_width,sigma_nb);
for i=1:sigma_nb
sigma=sigma_array(i);
snlo(:,:,i)=sigma*sigma*imfilter(img,fspecial('log',...
floor(6*sigma+1),sigma),'replicate');
end
%搜索局部极值
snlo_dil=imdilate(snlo,ones(3,3,3));
blob_candidate_index=find(snlo==snlo_dil);
blob_candidate_value=snlo(blob_candidate_index);
[temp,index]=sort(blob_candidate_value,'descend');
blob_index=blob_candidate_index(index(1:min(num,numel(index))));
[lig,col,sca]=ind2sub([img_height,img_width,sigma_nb],blob_index);
points=[lig,col,3*reshape(sigma_array(sca),[size(lig,1),1])];
end

function draw(img,pt,str)
%功能：在图像中绘制特征点
%img——输入图像
%pt——特征点坐标
%str——图上显示的名称
figure('Name',str);
imshow(img);
hold on;
axis off;
switch size(pt,2)
case 2
s=2;
for i=1:size(pt,1)
rectangle('Position',[pt(i,2)-s,pt(i,1)-s,2*s,2*s],'Curvature'...
,[0,0],'EdgeColor','b','LineWidth',2);
end
case 3
for i=1:size(pt,1)
rectangle('Position',[pt(i,2)-pt(i,3),pt(i,1)-pt(i,3),...
2*pt(i,3),2*pt(i,3)],'Curvature',[1,1],'EdgeColor',...
'w','LineWidth',2);
end
end
end

在matlab指令窗口输入：

img=imread('向日葵.jpg');

imshow(img);

pt=LoG_Blob(rgb2gray(img));

draw(img,pt,'LOG')

即可得到：