【圖像配準】基于sift算法實現(xiàn)圖像配準matlab源碼
一、簡介
SIFT即尺度不變特征變換,是用于圖像處理領域的一種描述。這種描述具有尺度不變性,可在圖像中檢測出關鍵點,是一種局部特征描述子。
1 SIFT算法特點:
(1)具有較好的穩(wěn)定性和不變性,能夠適應旋轉、尺度縮放、亮度的變化,能在一定程度上不受視角變化、仿射變換、噪聲的干擾。
(2)區(qū)分性好,能夠在海量特征數(shù)據(jù)庫中進行快速準確的區(qū)分信息進行匹配
(3)多量性,就算只有單個物體,也能產(chǎn)生大量特征向量
(4)高速性,能夠快速的進行特征向量匹配
(5)可擴展性,能夠與其它形式的特征向量進行聯(lián)合
2 SIFT算法實質(zhì)
在不同的尺度空間上查找關鍵點,并計算出關鍵點的方向。

3 SIFT算法實現(xiàn)特征匹配主要有以下三個流程:
(1)提取關鍵點:關鍵點是一些十分突出的不會因光照、尺度、旋轉等因素而消失的點,比如角點、邊緣點、暗區(qū)域的亮點以及亮區(qū)域的暗點。此步驟是搜索所有尺度空間上的圖像位置。通過高斯微分函數(shù)來識別潛在的具有尺度和旋轉不變的興趣點。
(2)定位關鍵點并確定特征方向:在每個候選的位置上,通過一個擬合精細的模型來確定位置和尺度。關鍵點的選擇依據(jù)于它們的穩(wěn)定程度。然后基于圖像局部的梯度方向,分配給每個關鍵點位置一個或多個方向。所有后面的對圖像數(shù)據(jù)的操作都相對于關鍵點的方向、尺度和位置進行變換,從而提供對于這些變換的不變性。
(3)通過各關鍵點的特征向量,進行兩兩比較找出相互匹配的若干對特征點,建立景物間的對應關系。
4 尺度空間
(1)概念
尺度空間即試圖在圖像領域中模擬人眼觀察物體的概念與方法。例如:觀察一顆樹,關鍵在于我們想要觀察是樹葉子還是整棵樹:如果是一整棵樹(相當于大尺度情況下觀察),那么就應該去除圖像的細節(jié)部分。如果是樹葉(小尺度情況下觀察),那么就該觀察局部細節(jié)特征。
SIFT算法在構建尺度空間時候采取高斯核函數(shù)進行濾波,使原始圖像保存最多的細節(jié)特征,經(jīng)過高斯濾波后細節(jié)特征逐漸減少來模擬大尺度情況下的特征表示。
利用高斯核函數(shù)進行濾波的主要原因有兩個:
a 高斯核函數(shù)是唯一的尺度不變核函數(shù)。
b DoG核函數(shù)可以近似為LoG函數(shù),這樣可以使特征提取更加簡單。同時,David. Lowe作者在論文中提出將原始圖像進行2倍上采樣后濾波能夠保留更多的信息便于后續(xù)特征提取與匹配。其實尺度空間圖像生成就是當前圖像與不同尺度核參數(shù)σ進行卷積運算后產(chǎn)生的圖像。
(2)表示
L(x, y, σ) ,定義為原始圖像 I(x, y)與一個可變尺度的2維高斯函數(shù)G(x, y, σ) 卷積運算。

5 高斯金字塔的構建
(1)概念
尺度空間在實現(xiàn)時使用高斯金字塔表示,高斯金字塔的構建分為兩步:
a 對圖像做高斯平滑;
b 對圖像做降采樣。

圖像的金字塔模型是指將原始圖像不斷降階采樣,得到一系列大小不一的圖像,由大到小,從下到上構成的塔狀模型。原圖像為金子塔的第一層,每次降采樣所得到的新圖像為金字塔的一層(每層一張圖像),每個金字塔共n層。為了讓尺度體現(xiàn)其連續(xù)性,高斯金字塔在簡單降采樣的基礎上加上了高斯濾波。如上圖所示,將圖像金字塔每層的一張圖像使用不同參數(shù)做高斯模糊,Octave表示一幅圖像可產(chǎn)生的圖像組數(shù),Interval表示一組圖像包括的圖像層數(shù)。另外,降采樣時,高斯金字塔上一組圖像的初始圖像(底層圖像)是由前一組圖像的倒數(shù)第三張圖像隔點采樣得到的。
(2)表示
高斯圖像金字塔共o組、s層,則有

6 DOG空間極值檢測
(1)DOG函數(shù)

(2)DoG高斯差分金字塔
a 對應DOG算子,需構建DOG金字塔。
可以通過高斯差分圖像看出圖像上的像素值變化情況。(如果沒有變化,也就沒有特征。特征必須是變化盡可能多的點。)DOG圖像描繪的是目標的輪廓。

b DOG局部極值檢測
特征點是由DOG空間的局部極值點組成的。為了尋找DoG函數(shù)的極值點,每一個像素點要和它所有的相鄰點比較,看其是否比它的圖像域和尺度域的相鄰點大或者小。特征點是由DOG空間的局部極值點組成的。為了尋找DoG函數(shù)的極值點,每一個像素點要和它所有的相鄰點比較,看其是否比它的圖像域和尺度域的相鄰點大或者小。如下圖,中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點比較,以確保在尺度空間和二維圖像空間都檢測到極值點。

b 去除邊緣效應
在邊緣梯度的方向上主曲率值比較大,而沿著邊緣方向則主曲率值較小。候選特征點的DoG函數(shù)D(x)的主曲率與2×2Hessian矩陣H的特征值成正比。


7 關鍵點方向分配
(1)通過尺度不變性求極值點,需要利用圖像的局部特征為給每一個關鍵點分配一個基準方向,使描述子對圖像旋轉具有不變性。對于在DOG金字塔中檢測出的關鍵點,采集其所在高斯金字塔圖像3σ鄰域窗口內(nèi)像素的梯度和方向分布特征。梯度的模值和方向如下:

(2)本算法采用梯度直方圖統(tǒng)計法,統(tǒng)計以關鍵點為原點,一定區(qū)域內(nèi)的圖像像素點確定關鍵點方向。在完成關鍵點的梯度計算后,使用直方圖統(tǒng)計鄰域內(nèi)像素的梯度和方向。梯度直方圖將0~360度的方向范圍分為36個柱,其中每柱10度。如下圖所示,直方圖的峰值方向代表了關鍵點的主方向,方向直方圖的峰值則代表了該特征點處鄰域梯度的方向,以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性,只保留峰值大于主方向峰值80%的方向作為該關鍵點的輔方向。

8 關鍵點描述
對于每一個關鍵點,都擁有位置、尺度以及方向三個信息。為每個關鍵點建立一個描述符,用一組向量將這個關鍵點描述出來,使其不隨各種變化而改變,比如光照變化、視角變化等等。這個描述子不但包括關鍵點,也包含關鍵點周圍對其有貢獻的像素點,并且描述符應該有較高的獨特性,以便于提高特征點正確匹配的概率。

Lowe實驗結果表明:描述子采用4×4×8=128維向量表征,綜合效果最優(yōu)(不變性與獨特性)。
9 關鍵點匹配
(1)分別對模板圖(參考圖,reference image)和實時圖(觀測圖,
observation image)建立關鍵點描述子集合。目標的識別是通過兩點集內(nèi)關鍵點描述子的比對來完成。具有128維的關鍵點描述子的相似性度量采用歐式距離。
(3)匹配可采取窮舉法完成,但所花費的時間太多。所以一般采用kd樹的數(shù)據(jù)結構來完成搜索。搜索的內(nèi)容是以目標圖像的關鍵點為基準,搜索與目標圖像的特征點最鄰近的原圖像特征點和次鄰近的原圖像特征點。
Kd樹如下如所示,是個平衡二叉樹

10 總結
SIFT特征具有穩(wěn)定性和不變性,在圖像處理和計算機視覺領域有著很重要的作用,其本身也是非常復雜的,由于接觸SIFT不是很久,對其中的相關知識了解還很不足,經(jīng)多方查閱參考,寫得此文,內(nèi)容還不夠詳盡,望多多見諒。以下是SIFT算法的粗略總結。
(1)DoG尺度空間的極值檢測。
(2)刪除不穩(wěn)定的極值點。
(3)確定特征點的主方向
(4)生成特征點的描述子進行關鍵點匹配。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
%% image path
file_image='F:\class_file\圖像配準\圖像配準';
%% read images
[filename,pathname]=uigetfile({'*.*','All Files(*.*)'},'Reference image',...
? ? ? ? ? ? ? ? ? ? ? ? ?file_image);
image_1=imread(strcat(pathname,filename));
[filename,pathname]=uigetfile({'*.*','All Files(*.*)'},'Image to be registered',...
? ? ? ? ? ? ? ? ? ? ? ? ?file_image);
image_2=imread(strcat(pathname,filename));
figure;
subplot(1,2,1);
imshow(image_1);
title('Reference image');
subplot(1,2,2);
imshow(image_2);
title('Image to be registered');
%% make file for save images
if (exist('save_image','dir')==0)%如果文件夾不存在
? ?mkdir('save_image');
end
t1=clock;%Start time
%% Convert input image format
[~,~,num1]=size(image_1);
[~,~,num2]=size(image_2);
if(num1==3)
? ?image_11=rgb2gray(image_1);
else
? ?image_11=image_1;
end
if(num2==3)
? ?image_22=rgb2gray(image_2);
else
? ?image_22=image_2;
end
%Converted to floating point data
image_11=im2double(image_11);
image_22=im2double(image_22); ?
%% Define the constants used
sigma=1.6;%最底層高斯金字塔的尺度
dog_center_layer=3;%定義了DOG金字塔每組中間層數(shù),默認是3
contrast_threshold_1=0.03;%Contrast threshold
contrast_threshold_2=0.03;%Contrast threshold
edge_threshold=10;%Edge threshold
is_double_size=false;%expand image or not
change_form='affine';%change mode,'perspective','affine','similarity'
is_sift_or_log='GLOH-like';%Type of descriptor,it can be 'GLOH-like','SIFT'
%% The number of groups in Gauss Pyramid
nOctaves_1=num_octaves(image_11,is_double_size);
nOctaves_2=num_octaves(image_22,is_double_size);
%% Pyramid first layer image
image_11=create_initial_image(image_11,is_double_size,sigma);
image_22=create_initial_image(image_22,is_double_size,sigma);
%% ?Gauss Pyramid of Reference image
tic;
[gaussian_pyramid_1,gaussian_gradient_1,gaussian_angle_1]=...
build_gaussian_pyramid(image_11,nOctaves_1,dog_center_layer,sigma); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
disp(['參考圖像創(chuàng)建Gauss Pyramid花費時間是:',num2str(toc),'s']);
%% DOG Pyramid of Reference image
tic;
dog_pyramid_1=build_dog_pyramid(gaussian_pyramid_1,nOctaves_1,dog_center_layer);
disp(['參考圖像創(chuàng)建DOG Pyramid花費時間是:',num2str(toc),'s']);
%% display the Gauss Pyramid,DOG Pyramid,gradient of Reference image
display_product_image(gaussian_pyramid_1,dog_pyramid_1,gaussian_gradient_1,...
? ? ? ?gaussian_angle_1,nOctaves_1,dog_center_layer,'Reference image'); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
clear gaussian_pyramid_1;
%% Reference image DOG Pyramid extreme point detection
tic;
[key_point_array_1]=find_scale_space_extream...
(dog_pyramid_1,nOctaves_1,dog_center_layer,contrast_threshold_1,sigma,...
edge_threshold,gaussian_gradient_1,gaussian_angle_1);
disp(['參考圖像關鍵點定位花費時間是:',num2str(toc),'s']);
clear dog_pyramid_1;
%% descriptor generation of the reference image
tic;
[descriptors_1,locs_1]=calc_descriptors(gaussian_gradient_1,gaussian_angle_1,...
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?key_point_array_1,is_double_size,is_sift_or_log);
disp(['參考圖像描述符生成花費時間是:',num2str(toc),'s']);
clear gaussian_gradient_1;
clear gaussian_angle_1;
%% Gauss Pyramid of the image to be registered
tic;
[gaussian_pyramid_2,gaussian_gradient_2,gaussian_angle_2]=...
build_gaussian_pyramid(image_22,nOctaves_2,dog_center_layer,sigma); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
disp(['待配準圖像創(chuàng)建Gauss Pyramid花費時間是:',num2str(toc),'s']);
%% DOG of the image to be registered
tic;
dog_pyramid_2=build_dog_pyramid(gaussian_pyramid_2,nOctaves_2,dog_center_layer);
disp(['待配準圖像創(chuàng)建DOG Pyramid花費時間是:',num2str(toc),'s']);
display_product_image(gaussian_pyramid_2,dog_pyramid_2,gaussian_gradient_2,...
? ? ? ?gaussian_angle_2,nOctaves_2,dog_center_layer,'Image to be registered'); ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
clear gaussian_pyramid_2;
%% Image to be registered DOG Pyramid extreme point detection
tic;
[key_point_array_2]=find_scale_space_extream...
(dog_pyramid_2,nOctaves_2,dog_center_layer,contrast_threshold_2,sigma,....
edge_threshold,gaussian_gradient_2,gaussian_angle_2);
disp(['待配準圖像關鍵點定位花費時間是:',num2str(toc),'s']);
clear dog_pyramid_2;
%% descriptor generation of the Image to be registered
tic;
[descriptors_2,locs_2]=calc_descriptors(gaussian_gradient_2,gaussian_angle_2,...
? ? ? ? ? ? ? ? ? ? ? key_point_array_2,is_double_size,is_sift_or_log);
disp(['待配準圖像描述符生成花費時間是:',num2str(toc),'s']);
clear gaussian_gradient_2;
clear gaussian_angle_2;
%% match
tic;
[solution,rmse,cor1,cor2]=...
? ?match(image_2, image_1,descriptors_2,locs_2,descriptors_1,locs_1,change_form);
disp(['特征點匹配花費時間是:',num2str(toc),'s']);



?

?