Harris Corner Detector生成的特征对于缩放不是不变的。对于特征跟踪，我们需要对仿射变换不变的特征。拉普拉斯斑点检测器是生成对缩放不变的特征的基本方法之一。

拉普拉斯斑点检测器的想法是将图像与多个尺度的“斑点滤波器”进行卷积，并在得到的尺度空间中寻找滤波器响应的极值。

斑点滤波器（Blob Filter）：

这个滤波器是由沿x轴和y轴的double derivating Gaussian filter产生并相加。以上也称为高斯拉普拉斯算子。

当图像中有一个角时，LOG滤波器的输出最大。

我们想通过将图像与不同的sigma卷积来找到blob的特征尺度。但随着我们增加sigma，blob滤波器对图像的响应会降低，因此我们可能无法找到更好的尺度（scale）。为了减少这种影响，我们将LOG滤波器σ平方相乘。这个过程称为尺度归一化。

拉普拉斯档案对半径为r的binary圆的最大响应是σ= 1.414 * r

以上是blob滤波器的一些基础知识。整个过程归结为两个步骤

1. 在几个尺度上用尺度归一化拉普拉斯变换卷积图像（不同尺度意味着不同的西格玛）
2. 在尺度空间中找到平方拉普拉斯响应的最大值

拉普拉斯Blob检测的Python代码

该代码有四个主要步骤：

1. 生成LOG滤波器
2. 用高斯滤波器对图像进行卷积
3. 找到最大峰值
4. 绘制blob。

生成LOG滤波器的Python代码

import cv2
from pylab import *
import numpy as np
import matplotlib.pyplot as plt
img = cv2.imread("sunflowers.jpg",0) #gray scale conversion
from scipy import ndimage
from scipy.ndimage import filters
from scipy import spatial
k = 1.414
sigma = 1.0
img = img/255.0 #image normalization
def LoG(sigma):
 #window size 
 n = np.ceil(sigma*6)
 y,x = np.ogrid[-n//2:n//2+1,-n//2:n//2+1]
 y_filter = np.exp(-(y*y/(2.*sigma*sigma)))
 x_filter = np.exp(-(x*x/(2.*sigma*sigma)))
 final_filter = (-(2*sigma**2) + (x*x + y*y) ) * (x_filter*y_filter) * (1/(2*np.pi*sigma**4))
 return final_filter

LoG函数将sigma作为输入，并生成大小为6 * sigma（最佳选项）的滤波器。

用高斯滤波器对图像进行卷积

图像与不同的LOG滤波器进行卷积，Python代码如下：

def LoG_convolve(img):
 log_images = [] #to store responses
 for i in range(0,9):
 y = np.power(k,i) 
 sigma_1 = sigma*y #sigma 
 filter_log = LoG(sigma_1) #filter generation
 image = cv2.filter2D(img,-1,filter_log) # convolving image
 image = np.pad(image,((1,1),(1,1)),'constant') #padding 
 image = np.square(image) # squaring the response
 log_images.append(image)
 log_image_np = np.array([i for i in log_images]) # storing the #in numpy array
 return log_image_np
log_image_np = LoG_convolve(img)
#print(log_image_np.shape)

要选择不同的sigma，我将sigma与k ^ i相乘。log_image_np数组的维度为z * image_height * image_width，其中z表示与sigma相乘的第k次方。

在上面的函数中，我只生成了9个滤波器，图中的9个滤波器是

找到最大峰值

def detect_blob(log_image_np):
 co_ordinates = [] #to store co ordinates
 (h,w) = img.shape
 for i in range(1,h):
 for j in range(1,w):
 slice_img = log_image_np[:,i-1:i+2,j-1:j+2] #9*3*3 slice
 result = np.amax(slice_img) #finding maximum
 if result >= 0.03: #threshold
 z,x,y = np.unravel_index(slice_img.argmax(),slice_img.shape)
 co_ordinates.append((i+x-1,j+y-1,k**z*sigma)) #finding co-rdinates
 return co_ordinates
co_ordinates = list(set(detect_blob(log_image_np)))

在每个像素处，在所有尺度上考虑3x3邻域，其给出9x3x3矩阵。矩阵的最大峰值由矩阵的最大值给出。

对于每个像素，都会有一个最大值。但并非所有像素都有助于blob。所以我们需要通过阈值来消除它们。

如果最大值大于阈值，则考虑该点并存储坐标。

（z，x，y）是切片图像的坐标，因此我们需要将它们更改为原始图像坐标。在更改它们之后，我们将它们存储在co_ordinates列表中并绘制它们。Python代码如下：

fig, ax = plt.subplots()
nh,nw = img.shape
count = 0
ax.imshow(img, interpolation='nearest',cmap="gray")
for blob in co_ordinates:
 y,x,r = blob
 c = plt.Circle((x, y), r*1.414, color='red', linewidth=1.5, fill=False)
 ax.add_patch(c)
ax.plot() 
plt.show()

上面的代码在存储列表的每个点绘制blob。

结果

在上面的图像中，我们可以发现许多blob都是重叠的。为了减少重叠的斑点，我们将找到斑点之间的重叠区域。如果该区域大于阈值，则将丢弃具有较小半径的斑点。python实现如下：

#provided in scipy doucumentaion
def blob_overlap(blob1, blob2):
 n_dim = len(blob1) - 1
 root_ndim = sqrt(n_dim)
 #print(n_dim)
 
 # radius of two blobs
 r1 = blob1[-1] * root_ndim
 r2 = blob2[-1] * root_ndim
 
 d = sqrt(np.sum((blob1[:-1] - blob2[:-1])**2))
 
 #no overlap between two blobs
 if d > r1 + r2:
 return 0
 # one blob is inside the other, the smaller blob must die
 elif d <= abs(r1 - r2):
 return 1
 else:
 #computing the area of overlap between blobs
 ratio1 = (d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d * r1)
 ratio1 = np.clip(ratio1, -1, 1)
 acos1 = math.acos(ratio1)
ratio2 = (d ** 2 + r2 ** 2 - r1 ** 2) / (2 * d * r2)
 ratio2 = np.clip(ratio2, -1, 1)
 acos2 = math.acos(ratio2)
a = -d + r2 + r1
 b = d - r2 + r1
 c = d + r2 - r1
 d = d + r2 + r1
area = (r1 ** 2 * acos1 + r2 ** 2 * acos2 -0.5 * sqrt(abs(a * b * c * d)))
 return area/(math.pi * (min(r1, r2) ** 2))
 
def redundancy(blobs_array, overlap):
 sigma = blobs_array[:, -1].max()
 distance = 2 * sigma * sqrt(blobs_array.shape[1] - 1)
 tree = spatial.cKDTree(blobs_array[:, :-1])
 pairs = np.array(list(tree.query_pairs(distance)))
 if len(pairs) == 0:
 return blobs_array
 else:
 for (i, j) in pairs:
 blob1, blob2 = blobs_array[i], blobs_array[j]
 if blob_overlap(blob1, blob2) > overlap:
 if blob1[-1] > blob2[-1]:
 blob2[-1] = 0
 else:
 blob1[-1] = 0
return np.array([b for b in blobs_array if b[-1] > 0])
co_ordinates = np.array(co_ordinates)
co_ordinates = redundancy(co_ordinates,0.5)

删除重叠blob后输出为

我们可以看到大多数斑点已被消除。