1. 什么是斑點
斑點通常是指與周圍有着顏色和灰度差別的區域。在實際地圖中，往往存在着大量這樣的斑點，如一顆樹是一個斑點，一塊草地是一個斑點，一棟房子也可以是一個斑點。由於斑點代表的是一個區域，相比單純的角點，它的穩定性要好，抗噪聲能力要強，所以它在圖像配准上扮演了很重要的角色。
同時有時圖像中的斑點也是我們關心的區域，比如在醫學與生物領域，我們需要從一些X光照片或細胞顯微照片中提取一些具有特殊意義的斑點的位置或數量。
比如下圖中天空的飛機、向日葵的花盤、X線斷層圖像中的兩個斑點。
在視覺領域，斑點檢測的主要思路都是檢測出圖像中比它周圍像素灰度值大或比周圍灰度值小的區域。一般有兩種方法來實現這一目標：
- 基於求導的微分方法，這類的方法稱為微分檢測器；
- 基於局部極值的分水嶺算法。
這里我們重點介紹第一種方法，主要檢測LOG斑點。而OpenCV中SimpleBlobDetector斑點檢測算子就實現了第二種方法，我們這里也會介紹它的接口使用方法。
2. LOG斑點檢測
2.1 基本原理
利用高斯拉普通拉斯（Laplace of Gaussian,LOG）算子檢測圖像斑點是一種十分常用的方法，對於二維高斯函數：
$$G(x,y;\sigma) = \frac{1}{2\pi\sigma^2}exp(-\frac{x^2+y^2}{2\sigma^2})$$
它的拉普拉斯變換為:
$$\nabla^2g = \frac{\partial^2g}{\partial x^2}+\frac{\partial^2g}{\partial y^2}$$
規范化的高斯拉普變換為：
$$\nabla^2_{norm}=\sigma^2\nabla^2g=\sigma^2(\frac{\partial^2g}{\partial x^2}+\frac{\partial^2g}{\partial y^2}) = -\frac{1}{2\pi\sigma^2}[1-\frac{x^2+y^2}{\sigma^2}]\cdot exp(-\frac{x^2+y^2}{2\sigma^2})$$
規范化算法子在二維圖像上顯示是一個圓對稱函數，如下圖所示。我們可以用這個算子來檢測圖像中的斑點，並且可以通過改變$\sigma$的值，可以檢測不同尺寸的二維斑點。
2.2 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)$$
然后對圖像的拉普拉斯圖像則為：
$$\nabla^2 = \frac{\partial^2L}{\partial x^2}+\frac{\partial^2L}{\partial y^2}$$
而實際上有下面等式：
$$\nabla^2[G(x,y)*f(x,y)] = \nabla^2[G(x,y)]*f(x,y)$$
所以，我們可以先求高斯核的拉普拉斯算子，再對圖像進行卷積。也就是一開始描述的步驟。
2.3 LOG算子的實現
Mat Feat::getHOGKernel(Size& ksize, double sigma) { Mat kernel(ksize, CV_64F); Point centPoint = Point((ksize.width -1)/2, ((ksize.height -1)/2)); // first calculate Gaussian for (int i=0; i < kernel.rows; i++) { double* pData = kernel.ptr<double>(i); for (int j = 0; j < kernel.cols; j++) { double param = -((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x)) / (2*sigma*sigma); pData[j] = exp(param); } } double maxValue; minMaxLoc(kernel, NULL, &maxValue); for (int i=0; i < kernel.rows; i++) { double* pData = kernel.ptr<double>(i); for (int j = 0; j < kernel.cols; j++) { if (pData[j] < EPS* maxValue) { pData[j] = 0; } } } double sumKernel = sum(kernel)[0]; if (sumKernel != 0) { kernel = kernel / sumKernel; } // now calculate Laplacian for (int i=0; i < kernel.rows; i++) { double* pData = kernel.ptr<double>(i); for (int j = 0; j < kernel.cols; j++) { double addition = ((i - centPoint.y) * (i - centPoint.y) + (j - centPoint.x) * (j - centPoint.x) - 2*sigma*sigma)/(sigma*sigma*sigma*sigma); pData[j] *= addition; } } // make the filter sum to zero sumKernel = sum(kernel)[0]; kernel -= (sumKernel/(ksize.width * ksize.height)); return kernel; }
2.4 多尺度檢測
我們注意到當$\sigma$尺度一定時，只能檢測對應半徑的斑點，那么檢測的是多大半徑的斑點呢，我們可以通過對規范化的二維拉普拉斯高斯算子求導：
規范化的高斯拉普拉斯函數為：
$$\nabla^2_{norm} = -\frac{1}{2\pi\sigma^2}[1-\frac{x^2+y^2}{\sigma^2}]\cdot exp(-\frac{x^2+y^2}{2\sigma^2})$$
求$\nabla^2_{norm}$的極點值等價於求取下式：
$$\frac{\partial(\nabla^2_{norm})}{\partial\sigma} = 0$$
得到：
$$(x^2+y^2-2\sigma^2)\cdot exp(-\frac{(x^2+y^2)}{2\sigma^2})$$
$$r^2-2\sigma^2=0$$
對於圖像中的斑點，在尺度$\sigma=r/\sqrt{2}$時，高斯拉普拉斯響應值達到最大。同理，如果圖像中的圓形斑點黑白反向，那么，它的高斯拉普拉斯響應值在$\sigma=r/\sqrt{2}$時達到最小。將高斯拉普拉斯響應達到峰值時的尺度$\sigma$值，稱為特征尺度。
那么在多尺度的情況下，同時在空間和尺度上達到最大值（或最小值）的點就是我們所期望的斑點。對於二維圖像$I(x,y)$，計算圖像在不同尺度下的離散拉普拉斯響應值，然后檢查位置空間中的每個點；如果該點的拉普拉斯響應值都大小於或小於其他26個立方空間領域(9+8+9)的值，那么該點就是被檢測到的圖像斑點。
3 OpenCV進行斑點檢測
opencv中檢測Blobs的類為SimpleBlobDetector，這個類在opencv中的定義如下：
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: ... };
算法的大致步驟如下：
- 對[minThreshold,maxThreshold)區間，以thresholdStep為間隔，做多次二值化。
- 對每張二值圖片，使用findContours()提取連通域並計算每一個連通域的中心。
- 根據2得到的中心，全部放在一起。一些很接近的點［由theminDistBetweenBlobs控制多少才算接近］被歸為一個group,對應一個bolb特征..
- 從3得到的那些點,估計最后的blob特征和相應半徑，並以key points返回。
同時該支持提取特征的方法，一共有5個選項，這里就不多加描述了，默認是提取黑色圓形的Blob特征。下面是一個示例
int main(int argc, char** argv) { Mat image = imread(argv[1]); vector<KeyPoint> keyPoints; SimpleBlobDetector::Params params; SimpleBlobDetector blobDetect(params); blobDetect.create("SimpleBlob"); blobDetect.detect(image, keyPoints); cout << keyPoints.size() << endl; drawKeypoints(image, keyPoints, image, Scalar(255,0,0)); namedWindow("blobs"); imshow("blobs", image); waitKey(); return 0; }
總體來說，OpenCV的斑點檢測效果還算不錯，但是在有些圖像的效果上明顯不如LOG算子檢測的檢測效果。
4. 擴展閱讀
一個與LOG濾波核近似的是高斯差分DOG濾波核，它的定義為：
$$D(x,y,\sigma) = (G(x,y,k\sigma) – G(x,y,\sigma))*I(x,y) = L(x,y,k\sigma)-L(x,y,\sigma)$$
其中$k$為兩個相鄰尺度間的比例因子。
DOG可以看作為LOG的一個近似，但是它比LOG的效率更高。
前面介紹的微分算子在近圓的斑點檢測方面效果很好，但是這些檢測算子被限定於只能檢測圓形斑點，而且不能估計斑點的方向，因為LOG算子等都是中心對稱的。如果我們定義一種二維高斯核的變形，記它在X方向與Y方向上具有不同的方差，則這種算子可以用來檢測帶有方向的斑點。
$$G(x,y)=\mathcal{A}\cdot exp(-[(ax^2+2bxy+cy^2)])$$
$$a = \frac{cos^2\theta}{2\sigma^2_x}+\frac{sin^2\theta}{2\sigma^2_y}, b=-\frac{sin2\theta}{2\sigma^2_x}+\frac{sin2\theta}{4\sigma^2_y},c = \frac{sin^2\theta}{2\sigma^2_x}+\frac{cos^2\theta}{2\sigma_y^2}$$
其中$\mathcal{A}$是規一性因子。
5. 參考資料
1. 《現代數字圖像 -- 處理技術提高與應用案例詳解》
2. 《圖像局部不變性特征與描述》
3. Lindeberg, T. Feature Detection with Automatic Scale Selection
4. Hui Kong. A Generalized Laplacian Of Gaussian Filter for Blob Detection and Its Applications.