基于各向同性高斯核的斑点检测算法

本算法针对于具有高度重叠区域斑点的检测。采用迭代LoG处理图像削弱重叠区域，随后采用单边高斯核（USG核）继续处理，抑制重叠区域和鞍区，从而使斑点完全分离，随后采用阈值分割技术，提取斑点，最终实现斑点检测。 

1.迭代LoG核的演变

高斯拉普拉斯（Laplacian of Gaussian, LoG）算子，其中 是拉普拉斯算子，G表示为标准差为σ的二维高斯函数

                           

LoG算子的核实质上是测量了图像中的所有方向上的平均局部对比度，但是相对于有部分重叠的区域来说，重叠区域的斑点都存在有一个突起部分和一个重叠部分。与重叠区域相比较来说，突起区域在图像的内部以及周围存在较大的正局部对比。

LoG核对突起区域往往比对重叠区域产生的响应更加明显，并且这种差异可以通过迭代的方式进行累积。通过利用这一事实，可以通过对图像进行迭代高斯拉普拉斯滤波，保留孤立点的同时削弱重叠斑点的重叠区域。

首先需要建立一个二维高斯函数来模拟斑点，其中x0表示中心位置，p0∈[0,1]表示突出度，ω0表示斑点规模，b0∈[0,1]表示背景强度，为了方便计算，将x0设置为[0,0]T，因此构建的二维高斯函数为：

                                    

高斯斑点的空间拓展一般可以用半径r0来表示，通常是由尺度来决定的，因此可以近似的用r0 ≈ 3ω0来表示。

 

为了准确的实现高斯拉普拉斯滤波的迭代卷积功能，因此需要重新调整高斯拉普拉斯核，如下：                  

重新调整后的核与图像信号卷积后的响应更能精确的反应要检测的斑点的突出性，更加符合本文的算法要求。

将通过高斯函数模拟的斑点与重新调整过后的核进行卷积，得到的响应为：

      

         

通过设置x0=0得到点的中心响应为：

                                             

为了找到所需要的特征尺度，需要在L0(0,σ)中找到所产生的最大响应，因此需要通过计算L0(0,σ)对于σ的一阶导数并且将其设置为0以此来确定最大响应，因此所得到的的式子为：

                                     

 

通过将上述式子的偏导数设置为0，由此可以得出σ=ω0这就相当于，对于一个给定尺度为的ω0的高斯斑点，得到了其在LoG尺度空间中所产生的最大响应σ∗ = ω0通过σ∗ = ω0可以计算得出L0(0,σ)在尺度空间下所产生的的最大响应为L0(0;σ∗) = p0由此可以看出，将模拟出来的高斯斑点与调整之后的高斯核进行卷积，由此得到了最大响应p0，这就相当于是原始斑点中的突出，有上述结论可知，LoG核在σ∗ = ω0尺度上的响应为：

                             

由于本文所需要的斑点检测为亮斑点，所以对于上式中的结果只保留正值作为斑点的强度图。因此需要A1(x) = max (L0(x;σ∗), 0)。由于A1(x)是旋转对称的，所以它的半径r1是要从y轴上的L0(x；σ∗)的零点开始计算的。因此需要设置L0(x；σ∗)=0以及y=0。由此可以计算得出A1(x)的半径为：r1 = 2ω0由于之前已经得出高斯斑点的半径r0 ≈ 3ω0，可以将A1(x)近似为一个高斯斑点，由此可以得出A1(x)的近似尺度为：ω1 ≈ r1/3 = 2ω0 /3。

这就意味着在高斯拉普拉斯滤波的卷积过程中，高斯斑点的尺度在有规律的进行缩减，由此可以根据σ∗ = ω0和ω1 ≈ r1/3 = 2ω0 /3得出，在A1(x)与高斯拉普拉斯核进行卷积之后，所得到的特征尺度应该为σ∗∗ ≈ ω1 = 2σ∗/3。

更加通俗的来说，在第一次高斯拉普拉斯滤波中使用的特征尺度为σ∗，则第二次高斯拉普拉斯滤波进行卷积时，应该将其特征尺度设置为σ∗∗=2/3σ∗。所以在进行算法实现的时候，通过设置了一个循环结构，使得每次迭代卷积的时候，将上一次的特征尺度缩小为2/3，以便进行接下来的卷积操作。所以通过设置一个比例集，St=Kt-1×S1,其中t∈{1,2,……，T}表示t的迭代次数，T为总的迭代次数，K表示尺度衰退因子，由于上面的研究表明，K的值应该是2/3。

迭代LoG核代码展示：

function im_log = log_filt(I,sigma) %构建高斯拉普拉斯核
radius = round(3*sigma); 
G1 = zeros(2*radius+1,2*radius+1);
for x = -radius:radius
    for y = -radius:radius
        G1(x+radius+1,y+radius+1) =exp(-(x^2+y^2)/(2*sigma^2))/(-(pi*sigma^4)/(x^2+y^2-2*sigma^2));
    end
end
K = imfilter(double(I),G1,'replicate')/255;  

im_log = max(K,0);

2.单边高斯核演变

单边高斯核核已经被用于去除斑点噪声的操作，在这里需要的是进一步利用它能够分离部分重叠斑点的能力。为了进一步方便计算，在此将之前等式中的二阶高斯核进行了重新的调整，调整后的结果为：

                             

使用这样的方式，就能构建出一个能够精确对斑点的突起部分作出响应的核。

通过旋转θ角度的方式可以得到一个新的归一化高斯核，由此表示为：

由此，可以将归一化的二阶高斯核在空间上被分为3个部分，即中心部分KC,以及两个对称的横向部分KL和KR，每一个都可以表示为一个单独的高斯核：

                

       

有上述式子可以得出结论，KL(x;σ;θ)和KR(x;σ;θ)具有相同的符号，而KC(x;σ;θ)具有相反的符号。

单边高斯核代码展示 

function U = log_usg(theta,sigma) %构建单边高斯核
    radius = round(3*sigma);
    K = zeros(2*radius+1,2*radius+1);
    Kc = K; Kr = K;
    R = [cos(theta),sin(theta);(-sin(theta)),(cos(theta))];
    for x = -radius:radius
        for y = -radius:radius
            X = [x;y];
            K(x+radius+1,y+radius+1) =(-2)*exp(-((R*X)'*(R*X))/(2*sigma^2))./(-(pi*sigma^4))*(((x*cos(theta)+y*sin(theta))^2-sigma^2));
            if abs(x*cos(theta)+y*sin(theta))<= sigma
                Kc(x+radius+1,y+radius+1)= K(x+radius+1,y+radius+1);
            end
            if (x*cos(theta)+y*sin(theta)) >= sigma
                Kr(x+radius+1,y+radius+1) = K(x+radius+1,y+radius+1);
            end
        end
    end
    U = Kc+2*Kr;

完整代码：

clear
clc
X = imread('图片存储地址'); %%调取图片

set(0,'defaultfigurecolor','w')
image(X);
[m,n,z] = size(X);%%得到图片像素大小
G = rgb2gray(X);
[m,n] = size(G);%%得到图片像素大小
iters = 3;
result_t = zeros(m,n,iters);
for sigma0 = 13:0.5:13.5
   for t = 1:3 %%迭代次数
        sigma = sigma0*(2/3)^(t-1);                                          
        result_t(:,:,t) = log_filt(G,sigma);  %调用高斯拉普拉斯核函数
        G = result_t(:,:,t);
        figure(t+1)
        imagesc(result_t(:,:,t)); 
        title(['第' num2str(t) '迭代高斯拉普拉斯核'])
   end
end
    for k = 1:4
        theta = 90*(k-1);           %对theta的取值限定为0°，45°,90°,270°
        sigma = sigma0*(2/3)^(3); 
        result_u_t(:,:,t) = log_usg(theta,sigma);  %调用单边高斯核函数
        E = result_u_t(:,:,t);
        P = conv2(G,E,'same')/255; 
        P(:,:,k) = P;
        figure
        imagesc(result_u_t(:,:,t)); 
        title( 'usg核')
        figure
        imagesc(P(:,:,k));  
        title( 'usg核卷积')
    end
L= min(P,[],3);
figure
imagesc(L)
%对L进行归一化处理
Lmax = max(max(L));    
Lmin = min(min(L));
Y = (L-Lmin)/(Lmax-Lmin);
%对归一化后的图像进行阈值分割
Q = im2double(Y);
c = graythresh(Q);
BW = imbinarize(Q,c);
figure
imshow(BW)
%反转二值图像
bw = 1 - BW ;
figure
imagesc(bw)
%对图像非连接区域标记
q = bwlabel(bw,4);
figure
imagesc(q)
%得到每个区域质心的坐标
STATS = regionprops(q,'Centroid');
centroids = cat(1,STATS.Centroid);
%标记坐标
image(X);
hold on;
plot(centroids(:,1),centroids(:,2), 'r.','Markersize',20);
hold off;

算法结果展示：

原图

 

经过迭代LoG处理后结果图

单边高斯核处理结果图

 

 

最终结果图

注意事项：本算法是基础的处理重叠斑点分离，并不够智能，因此每次设置新图片时需要手动更改sigma0处的数值大小，使其达到最优数值。