1.原理
Difference of Gaussian(DOG)是高斯函數的差分。將兩幅圖像在不同參數下的高斯濾波結果相減，得到DoG圖。步驟：
-
處理一幅圖像在不同高斯參數下的DoG
用兩個不同的5x5高斯核對圖像進行卷積，然后再相減的操作。重復三次得到三個差分圖A,B,C。 -
根據DoG求角點
計算出的A,B,C三個DOG圖中求圖B中是極值的點。圖B的點在當前由A,B,C共27個點組成的block中是否為極大值或者極小值。若滿足此條件則認為是角點。
2.實現細節
2.1 差分得到DoG圖
Mat Process(Mat &A, double sig1, double sig2, Size Ksize)
{
Mat AF, out, out1, out2;
A.convertTo(AF, CV_32FC1);
GaussianBlur(AF, out1, Ksize, sig1, 0);
GaussianBlur(AF, out2, Ksize, sig2, 0);
subtract(out1, out2, out);
return out;
}
- converTo()
直接讀取灰度圖會得到CV_8UC1類型的mat，是單通道uchar型矩陣，因此高斯濾波后相減都是整型非負數據，影響后面進行比較取極值的步驟。所以需要將原本的數據類型轉化為CV_32FC1即單通道float型數據。再進行后續的操作 - subtract()
兩個mat類型相減直接用這個就好了。剛開始查到cvSub函數是老版本的，里面的指針亂七八糟的都挺復雜的，后在官方說明上找到這個c++版本的函數。現在的操作基本都是這個了。
2.2 求取當前block極值得到角點
Mat getExtrema(Mat &A, Mat &B, Mat &C, int thresh)
{
Mat block;
Mat extr = Mat::zeros(A.rows, A.cols, CV_32FC1);
double minv, maxv;
Point minLoc, maxLoc;
for (int i = 1; i < A.rows - 2; i++)
{
for (int j = 1; j < A.cols - 2; j++)
{
block.release();
block.push_back(A(Range(i - 1, i + 2), Range(j - 1, j + 2)));
block.push_back(B(Range(i - 1, i + 2), Range(j - 1, j + 2)));
block.push_back(C(Range(i - 1, i + 2), Range(j - 1, j + 2)));
minMaxLoc(block, &minv, &maxv, &minLoc, &maxLoc);
if ((maxLoc.x == 1 && maxLoc.y == 4) && maxv >(double)thresh)
extr.at<float>(i, j) = 1;
if ((minLoc.x == 1 && minLoc.y == 4) && minv < (double)(-thresh))
extr.at<float>(i, j) = -1;
}
}
return extr;
}
取出A，B，C中3x3的小塊存放在一個mat里，然后利用自帶的求取最大最小值函數minMaxLoc得到最大最小值和其位置。然后判斷該點\((i,j)\)是否為極值點。每個循環刷新釋放block。
2.3 求取當前block極值得到角點
void drawExtrema(Mat &imgB, Mat &extr)
將得到的極值點在原圖像中顯示。由於imread(img,0)得到的是單通道灰度圖 ，opencv提供的cirle()和line()函數在img上無法顯示，所以用imread(imgB)得到一個三通道的灰度圖，然后在這個imgB上畫出標記點。
2.4 討論
用樣例測試圖600x1200像素點的時候，運算速度很慢，逐步運行發現時間主要消耗在getExtrema()函數上。原因可能是：
- mat.push_back()截取矩陣然后合並得到新矩陣用這個操作比較費時間。如果改進的話可以直接用一個數組，將三個3x3矩陣的數據直接放進去。
- minMaxLoc()求取當前矩陣最大最小值和響應位置。這個函數涉及到排序可能也費時間。因為只用判斷當前\((i,j)\)點是否是極值，其實可以將這27個數據直接減去mat.at(i,j)根據是否都大於或者小於零來判斷極值。
然而我已經懶得改了……需要優化速度的話再動手。
3.測試結果
輸入圖像：
參考
關於opencv的函數，官方說明十分好，直接google之看一下就知道用法了。調試時將mat輸出到文件用到了參考鏈接3。
5.代碼
#include <iostream>
#include <fstream>
#include "opencv2/core/core.hpp"
#include "highgui.h"
#include "opencv2/imgproc/imgproc.hpp"
using namespace cv;
using namespace std;
Mat Process(Mat &A, double sig1, double sig2, Size Ksize)
{
Mat AF, out, out1, out2;
A.convertTo(AF, CV_32FC1);
GaussianBlur(AF, out1, Ksize, sig1, 0);
GaussianBlur(AF, out2, Ksize, sig2, 0);
subtract(out1, out2, out);
return out;
}
Mat getExtrema(Mat &A, Mat &B, Mat &C, int thresh)
{
Mat block;
Mat extr = Mat::zeros(A.rows, A.cols, CV_32FC1);
double minv, maxv;
Point minLoc, maxLoc;
for (int i = 1; i < A.rows - 2; i++)
{
for (int j = 1; j < A.cols - 2; j++)
{
block.release();
block.push_back(A(Range(i - 1, i + 2), Range(j - 1, j + 2)));
block.push_back(B(Range(i - 1, i + 2), Range(j - 1, j + 2)));
block.push_back(C(Range(i - 1, i + 2), Range(j - 1, j + 2)));
minMaxLoc(block, &minv, &maxv, &minLoc, &maxLoc);
if ((maxLoc.x == 1 && maxLoc.y == 4) && maxv >(double)thresh)
extr.at<float>(i, j) = 1;
if ((minLoc.x == 1 && minLoc.y == 4) && minv < (double)(-thresh))
extr.at<float>(i, j) = -1;
}
}
return extr;
}
void drawExtrema(Mat &img, Mat &extr)
{
Point center;
int r = 2;
for (int i = 0; i < img.rows; i++)
{
for (int j = 0; j < img.cols; j++)
{
if (extr.at<float>(i, j) == -1)
{
center = Point(j, i);
circle(img, center, r, Scalar(0, 255, 0));
}
if (extr.at<float>(i, j) == 1)
{
center = Point(j, i);
circle(img, center, r, Scalar(0, 0, 255));
}
}
}
namedWindow("DoG");
imshow("DoG", img);
waitKey(0);
}
int main(int argc, char** argv)
{
char file_name[] = "result1.png";
Mat img = imread(file_name, 0);
Mat A, B, C, a;
int threshold = 3;
A = Process(img, 0.3, 0.4, Size(5, 5));
B = Process(img, 0.6, 0.7, Size(5, 5));
C = Process(img, 0.7, 0.8, Size(5, 5));
a = getExtrema(A, B, C, threshold);
//const char* filename = "output.txt";
//writeMatToFile(B, filename);
Mat imgB = imread(file_name);
drawExtrema(imgB, a);
return 0;
}