SIFT原理與原始碼分析：DoG尺度空間構造
《SIFT原理與原始碼分析》系列文章索引：http://blog.csdn.net/xiaowei_cqu/article/details/8069548
尺度空間理論
自然界中的物體隨著觀測尺度不同有不同的表現形態。例如我們形容建築物用“米”，觀測分子、原子等用“納米”。更形象的例子比如Google地圖，滑動滑鼠輪可以改變觀測地圖的尺度，看到的地圖繪製也不同；還有電影中的拉伸鏡頭等等……
尺度空間中各尺度影象的模糊程度逐漸變大，能夠模擬人在距離目標由近到遠時目標在視網膜上的形成過程。
尺度越大影象越模糊。
為什麼要討論尺度空間？
用機器視覺系統分析未知場景時，計算機並不預先知道影象中物體的尺度。我們需要同時考慮影象在多尺度下的描述，獲知感興趣物體的最佳尺度
影象的尺度空間表達就是影象在所有尺度下的描述。
尺度空間表達與金字塔多解析度表達
高斯模糊
高斯核是唯一可以產生多尺度空間的核（《Scale-space theory: A basic tool for analysing structures at different scales》）。一個影象的尺度空間L(x,y,σ) ,定義為原始影象I(x,y)與一個可變尺度的2維高斯函式G(x,y,σ)卷積運算。
二維空間高斯函式：
尺度空間：
尺度是自然客觀存在的，不是主觀創造的。高斯卷積只是表現尺度空間的一種形式。
二維空間高斯函式是等高線從中心成正太分佈的同心圓：
分佈不為零的點組成卷積陣與原始影象做變換，即每個畫素值是周圍相鄰畫素值的高斯平均。一個5*5的高斯模版如下所示：
高斯模版是圓對稱的，且卷積的結果使原始畫素值有最大的權重，距離中心越遠的相鄰畫素值權重也越小。
在實際應用中，在計算高斯函式的離散近似時，在大概3σ距離之外的畫素都可以看作不起作用，這些畫素的計算也就可以忽略。所以，通常程式只計算(6σ+1)*(6σ+1)
高斯模糊另一個很厲害的性質就是線性可分：使用二維矩陣變換的高斯模糊可以通過在水平和豎直方向各進行一維高斯矩陣變換相加得到。
O(N^2*m*n)次乘法就縮減成了O(N*m*n)+O(N*m*n)次乘法。（N為高斯核大小，m,n為二維影象高和寬）
其實高斯這一部分只需要簡單瞭解就可以了，在OpenCV也只需要一句程式碼：
GaussianBlur(dbl, dbl, Size(), sig_diff, sig_diff);我這裡詳寫了一下是因為這塊兒對分析演算法效率比較有用，而且高斯模糊的演算法真的很漂亮~
金字塔多解析度
金字塔是早期影象多尺度的表示形式。影象金字塔化一般包括兩個步驟：使用低通濾波器平滑影象；對平滑影象進行降取樣（通常是水平，豎直方向1/2），從而得到一系列尺寸縮小的影象。
上圖中（a）是對原始訊號進行低通濾波，（b）是降取樣得到的訊號。
而對於二維影象，一個傳統的金字塔中，每一層影象由上一層解析度的長、寬各一半，也就是四分之一的畫素組成：
多尺度和多解析度
尺度空間表達和金字塔多解析度表達之間最大的不同是：
- 尺度空間表達是由不同高斯核平滑卷積得到，在所有尺度上有相同的解析度；
- 而金字塔多解析度表達每層解析度減少固定比率。
所以，金字塔多解析度生成較快，且佔用儲存空間少；而多尺度表達隨著尺度引數的增加冗餘資訊也變多。
多尺度表達的優點在於影象的區域性特徵可以用簡單的形式在不同尺度上描述；而金字塔表達沒有理論基礎，難以分析影象區域性特徵。
DoG（Difference of Gaussian）
高斯拉普拉斯LoG金字塔
結合尺度空間表達和金字塔多解析度表達，就是在使用尺度空間時使用金字塔表示，也就是計算機視覺中最有名的拉普拉斯金子塔（《The Laplacian pyramid as a compact image code》）。
高斯拉普拉斯LoG（Laplace of Guassian）運算元就是對高斯函式進行拉普拉斯變換：
核心思想還是高斯，這個不多敘述。
高斯差分DoG金字塔
DoG（Difference of Gaussian）其實是對高斯拉普拉斯LoG的近似，也就是對
的近似。SIFT演算法建議，在某一尺度上的特徵檢測可以通過對兩個相鄰高斯尺度空間的影象相減，得到DoG的響應值影象D(x,y,σ)。然後仿照LoG方法，通過對響應值影象D(x,y,σ)進行區域性最大值搜尋，在空間位置和尺度空間定位區域性特徵點。其中：
k為相鄰兩個尺度空間倍數的常數。
上圖中（a）是DoG的三維圖，（b）是DoG與LoG的對比。
金字塔構建
構建高斯金字塔
為了得到DoG影象，先要構造高斯金字塔。我們回過頭來繼續說高斯金字塔~
高斯金字塔在多解析度金字塔簡單降取樣基礎上加了高斯濾波，也就是對金字塔每層影象用不同引數的σ做高斯模糊，使得每層金字塔有多張高斯模糊影象。金字塔每層多張影象合稱為一組（Octave），每組有多張（也叫層Interval）影象。另外，降取樣時，金字塔上邊一組影象的第一張影象（最底層的一張）是由前一組（金字塔下面一組）影象的倒數第三張隔點取樣得到。
以下是OpenCV中構建高斯金字塔的程式碼，我加了相應的註釋：
-
// 構建nOctaves組（每組nOctaves+3層）高斯金字塔 -
void SIFT::buildGaussianPyramid( const Mat& base, vector<Mat>& pyr, int nOctaves ) const -
{ -
vector<double> sig(nOctaveLayers + 3); -
pyr.resize(nOctaves*(nOctaveLayers + 3)); -
// precompute Gaussian sigmas using the following formula: -
// \sigma_{total}^2 = \sigma_{i}^2 + \sigma_{i-1}^2、 -
// 計算對影象做不同尺度高斯模糊的尺度因子 -
sig[0] = sigma; -
double k = pow( 2., 1. / nOctaveLayers ); -
for( int i = 1; i < nOctaveLayers + 3; i++ ) -
{ -
double sig_prev = pow(k, (double)(i-1))*sigma; -
double sig_total = sig_prev*k; -
sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev); -
} -
for( int o = 0; o < nOctaves; o++ ) -
{ -
// DoG金子塔需要nOctaveLayers+2層影象來檢測nOctaves層尺度 -
// 所以高斯金字塔需要nOctaveLayers+3層影象得到nOctaveLayers+2層DoG金字塔 -
for( int i = 0; i < nOctaveLayers + 3; i++ ) -
{ -
// dst為第o組（Octave）金字塔 -
Mat& dst = pyr[o*(nOctaveLayers + 3) + i]; -
// 第0組第0層為原始影象 -
if( o == 0 && i == 0 ) -
dst = base; -
// base of new octave is halved image from end of previous octave -
// 每一組第0副影象時上一組倒數第三幅影象隔點取樣得到 -
else if( i == 0 ) -
{ -
const Mat& src = pyr[(o-1)*(nOctaveLayers + 3) + nOctaveLayers]; -
resize(src, dst, Size(src.cols/2, src.rows/2), -
0, 0, INTER_NEAREST); -
} -
// 每一組第i副影象是由第i-1副影象進行sig[i]的高斯模糊得到 -
// 也就是本組影象在sig[i]的尺度空間下的影象 -
else -
{ -
const Mat& src = pyr[o*(nOctaveLayers + 3) + i-1]; -
GaussianBlur(src, dst, Size(), sig[i], sig[i]); -
} -
} -
} -
}
高斯金字塔的組數為：
程式碼10-17行是計算高斯模糊的係數σ，具體關係如下：
其中，σ為尺度空間座標，s為每組中層座標，σ0為初始尺度，S為每組層數（一般為3~5）。根據這個公式，我們可以得到金字塔組內各層尺度以及組間各影象尺度關係。
組內相鄰影象尺度關係：
相鄰組間尺度關係：
所以，相鄰兩組的同一層尺度為2倍的關係。
最終尺度序列總結為：
o為金字塔組數，n為每組金字塔層數。
構建DoG金字塔
構建高斯金字塔之後，就是用金字塔相鄰影象相減構造DoG金字塔。
下面為構造DoG的程式碼：
-
// 構建nOctaves組（每組nOctaves+2層）高斯差分金字塔 -
void SIFT::buildDoGPyramid( const vector<Mat>& gpyr, vector<Mat>& dogpyr ) const -
{ -
int nOctaves = (int)gpyr.size()/(nOctaveLayers + 3); -
dogpyr.resize( nOctaves*(nOctaveLayers + 2) ); -
for( int o = 0; o < nOctaves; o++ ) -
{ -
for( int i = 0; i < nOctaveLayers + 2; i++ ) -
{ -
// 第o組第i副影象為高斯金字塔中第o組第i+1和i組影象相減得到 -
const Mat& src1 = gpyr[o*(nOctaveLayers + 3) + i]; -
const Mat& src2 = gpyr[o*(nOctaveLayers + 3) + i + 1]; -
Mat& dst = dogpyr[o*(nOctaveLayers + 2) + i]; -
subtract(src2, src1, dst, noArray(), CV_16S); -
} -
} -
}
這個比較簡單，就是一個subtract()函式。
至此，SIFT第一步就完成了。參見《SIFT原理與原始碼分析》
（轉載請註明作者和出處：http://blog.csdn.net/xiaowei_cqu 未經允許請勿用於商業用途）
--------------------- 本文來自 xiaowei_cqu 的CSDN 部落格 ，全文地址請點選：https://blog.csdn.net/xiaowei_cqu/article/details/8067881?utm_source=copy