【影像配準】基于matlab OpenSUFT影像配準【含Matlab原始碼 1232期】-有解無憂

一、SUFT配準簡介

0 引言
目前, 紅外檢測技術廣泛應用于電氣設備和電路板卡的故障檢測，影像配準技術是紅外影像處理中最關鍵的技術之一, 配準的結果直接影響到故障的檢測與定位，影像配準可分為基于灰度的影像配準和基于特征的影像配準，基于灰度的影像配準一般要求影像的相關性強, 而且計算量大, 很難達到實時性的需求;基于特征的影像配準計算量小、運算速度快, 且具有較強的魯棒性, 成為影像配準研究的主要方向，

常用的特征提取演算法有Harris, SUSAN, SIFT (Scale Invariant Features Transform) 和SUFT(Speeded-Up Robust Features) 等，SIFT算子最早由Lowe David G提出, 是建立在DoG (Difference of Gaussian) 尺度空間理論基礎上的一種演算法，該演算法采取鄰域方向性資訊聯合的思想, 從空間域和尺度域兩個方面對影像進行特征分析, 對檢測到的關鍵點用128維的特征向量表征, 具有尺度不變性和較強的魯棒性，由對比分析知, SIFT演算法性能優于Harris、SUSAN等角點演算法, 但SIFT算子比較耗時, 不能滿足實時性的要求，因此, Bay等人提出了一種基于快速魯棒特征的SUFT演算法, 它在特征點檢測的準確性、魯棒性以及實時性方面較其他演算法有很大優勢，

本文利用SUFT演算法進行特征點檢測, 采取粗匹配與精匹配結合的匹配策略選取特征點對, 設計了一種快速、有效、高精度的紅外影像配準演算法，

1 特征點提取
1.1 尺度空間特征點檢測
SUFT特征點檢測是基于Hessian矩陣進行的, 給定影像I (x, y) 中一點s=[x, y], 則在尺度σ的Hessian矩陣為
在這里插入圖片描述
式中, Lxx (x, σ) 為影像I (x, y) 和高斯函式G (x, y, σ) 在x方向上的二階導數在x的卷積，即

Lxy (x, σ) 、Lyy (x, σ) 與之類似，

高斯函式對尺度空間分析是最優的, 盒濾波器可近似為高斯函式二階導數, 由式 (2) 易知SUFT演算法可用盒濾波器近似Hessian矩陣，為了計算方便, 采用盒濾波器與輸入影像的卷積Dxx 、Dyy 、Dxy 代替Lxx 、Lyy 、Lxy ，把9×9的盒濾波器近似為σ=1.2的二階高斯導數Dxx、 Dxy與Lxx、Lxy之間關系為
在這里插入圖片描述
使用了盒濾波器和積分影像, 不必迭代地應用相同濾波器到前一個已濾層的輸出上, 而是應用不同大小的盒濾波器以相同的速度直接作用到原始影像上, 見圖1，

圖1 尺度空間的金字塔示意圖
因此, 對尺度空間的分析是通過增大濾波器的尺寸而不是迭代地降低影像的尺寸, 所需的計算時間是獨立于濾波器尺寸的, 從而大大降低了運算時間，

1.2 特征點描述子的生成
為保證特征點旋轉不變性, 要確定特征點主方向并建立坐標，以特征點為中心, 計算半徑為6s (s為特征點所在尺度值) 鄰域內的點在x、y方向的哈爾小波回應, 回應可表示為水平回應和垂直回應矢量和，按距離賦予回應值不同的高斯權重系數, 遠離特征點的回應貢獻小, 靠近特征點的回應貢獻大，通過計算π/3滑動方向視窗內所有回應的和來估計主方向, 視窗內的水平回應和垂直回應分別被求和，兩個回應的和產生一個區域方向向量, 定義視窗中最長的向量為特征點的主方向，

選定特征點主方向后, 以特征點為中心, 將坐標軸旋轉到主方向上，選取周圍邊長為20s×20s的正方形區域, 并將該區域劃分為4×4共16個子區域，在子區域內計算每個像素點x方向和y方向的哈爾小波回應, 記為dx、dy, 為了增強特征點描述子對幾何變換及區域誤差的魯棒性, 可對dx、dy進行高斯權重系數賦值，然后對每個子區域的dx、dy回應以及響
在這里插入圖片描述
2 影像配準方法及流程
利用SUFT演算法提取紅外影像中的特征點并生成特征點描述子, 采取歐氏距離最近鄰粗匹配和相似四邊形精匹配的方式提高配準精度, 然后使用8引數的平面透視變換模型描述匹配影像序列間的相對變換關系, 依據精匹配得到的匹配點對求解模型變換引數, 從而實作紅外影像的配準，演算法流程如圖2所示，
在這里插入圖片描述
圖2 本文演算法流程圖
演算法具體實作程序如下，

利用SUFT演算法分別檢測標準影像F與待配準影像F′的特征點, 形成64維的特征點描述子，
最近鄰匹配，以歐氏距離作為兩個特征點描述子的相似性度量進行粗匹配, 算式為

式中:Xik表示影像F中第i個特征點對應特征向量的第k個元素;Xjk表示影像F′中第j個特征點對應特征向量的第k個元素;n為特征向量的維數，計算每個特征點對應特征向量的歐氏距離, 按照從小到大的順序排列形成距離集合，設定閾值T1, 當特征點最小歐氏距離與次小歐氏距離的比值小于T1時, 認為這兩個特征點匹配，T1越小, 匹配點對數目越少, 但更加穩定，
特征點精匹配，利用景物幾何結構間的相似性, 在粗匹配點對中尋找相似四邊形進行精匹配, 從而減少誤匹配的概率，選取粗匹配點對中的一對匹配點作為四邊形的一個固定頂點, 在粗匹配點對中隨機選取三對匹配點組成兩對四邊形，由相似四邊形性質可知, 若兩四邊形相似, 則對應四條邊和兩條對角線互成比例, 即滿足

依據四邊形相似的性質構造歸一化的均方誤差運算式e, 設定閾值T2 , 當e小于T2時認為兩四邊形相似，T2越小, 匹配精度越高, 由于SUFT提取特征點時存在誤差, 故閾值不能設定太小, 通常設定比計算精度高出2～3個數量級, 本文T2設為0.05,

由影像F和F′中固定頂點組成相似四邊形的數量判斷是否為匹配點對, 剔除粗匹配點對中不能組成相似形的點對, 從而實作精匹配，
求解變換模型引數，由于紅外熱像儀拍攝位置的不固定性, 采用更加符合實際情況的8引數平面透視變換模型, 則影像F和F′對應像素點的關系為:X′=HX, 表示為

① 當n<4時, 適度增大閾值T1、T2, 重復步驟2) 、3) 以獲得更多的精匹配點對;
② 當n=4時, 依據式 (9) 求解矩陣H;
③ 當n>4時, 屬于過約束的情況, 利用最小二乘法求解H, 尋找一個最佳解使得每對匹配點的平面透視模型均方誤差最小, 達到多個配準點擬合最優引數解的目的，
利用矩陣H對紅外影像進行模型變換, 并通過線性插值得到變換后的影像, 從而實作影像配準，

二、部分源代碼

% Example 3, Affine registration
% Load images
I1=im2double(imread('TestImages/lena1.png'));
I2=im2double(imread('TestImages/lena2.png'));

% Get the Key Points
Options.upright=true;
Options.tresh=0.0001;
Ipts1=OpenSurf(I1,Options);
Ipts2=OpenSurf(I2,Options);

% Put the landmark descriptors in a matrix
D1 = reshape([Ipts1.descriptor],64,[]);
D2 = reshape([Ipts2.descriptor],64,[]);

% Find the best matches
err=zeros(1,length(Ipts1));
cor1=1:length(Ipts1);
cor2=zeros(1,length(Ipts1));
for i=1:length(Ipts1),
    distance=sum((D2-repmat(D1(:,i),[1 length(Ipts2)])).^2,1);
    [err(i),cor2(i)]=min(distance);
end

% Sort matches on vector distance
[err, ind]=sort(err);
cor1=cor1(ind);
cor2=cor2(ind);

% Make vectors with the coordinates of the best matches
Pos1=[[Ipts1(cor1).y]',[Ipts1(cor1).x]'];
Pos2=[[Ipts2(cor2).y]',[Ipts2(cor2).x]'];
Pos1=Pos1(1:30,:);
Pos2=Pos2(1:30,:);

% Show both images
I = zeros([size(I1,1) size(I1,2)*2 size(I1,3)]);
I(:,1:size(I1,2),:)=I1; I(:,size(I1,2)+1:size(I1,2)+size(I2,2),:)=I2;
figure, imshow(I); hold on;

% Show the best matches
plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','-');
plot([Pos1(:,2) Pos2(:,2)+size(I1,2)]',[Pos1(:,1) Pos2(:,1)]','o');

% Calculate affine matrix
Pos1(:,3)=1; Pos2(:,3)=1;
M=Pos1'/Pos2';

% Add subfunctions to Matlab Search path
functionname='OpenSurf.m';
functiondir=which(functionname);
functiondir=functiondir(1:end-length(functionname));
addpath([functiondir '/WarpFunctions'])

% Warp the image
I1_warped=affine_warp(I1,M,'bicubic');

% Show the result
figure,
subplot(1,3,1), imshow(I1);title('Figure 1');
subplot(1,3,2), imshow(I2);title('Figure 2');
subplot(1,3,3), imshow(I1_warped);title('Warped Figure 1');
function [ipts, np]=FastHessian_interpolateExtremum(r, c,  t,  m,  b,  ipts, np)
% This function FastHessian_interpolateExtremum will ..
%
% [ipts,np] = FastHessian_interpolateExtremum( r,c,t,m,b,ipts,np )
%  
%  inputs,
%    r : 
%    c : 
%    t : 
%    m : 
%    b : 
%    ipts : 
%    np : 
%  
%  outputs,
%    ipts : 
%    np : 
%  
% Function is written by D.Kroon University of Twente (July 2010)
D = FastHessian_BuildDerivative(r, c, t, m, b);
H = FastHessian_BuildHessian(r, c, t, m, b);

%get the offsets from the interpolation
Of = - H\D;
O=[ Of(1, 1), Of(2, 1), Of(3, 1) ];

%get the step distance between filters
filterStep = fix((m.filter - b.filter));

%If point is sufficiently close to the actual extremum
if (abs(O(1)) < 0.5 && abs(O(2)) < 0.5 && abs(O(3)) < 0.5)
    np=np+1;
    ipts(np).x = double(((c + O(1))) * t.step);
    ipts(np).y = double(((r + O(2))) * t.step);
    ipts(np).scale = double(((2/15) * (m.filter + O(3) * filterStep)));
    ipts(np).laplacian = fix(FastHessian_getLaplacian(m,r,c,t));
end
  
function D=FastHessian_BuildDerivative(r,c,t,m,b)
dx = (FastHessian_getResponse(m,r, c + 1, t) - FastHessian_getResponse(m,r, c - 1, t)) / 2;
dy = (FastHessian_getResponse(m,r + 1, c, t) - FastHessian_getResponse(m,r - 1, c, t)) / 2;
ds = (FastHessian_getResponse(t,r, c) - FastHessian_getResponse(b,r, c, t)) / 2;
D = [dx;dy;ds];

function H=FastHessian_BuildHessian(r, c, t, m, b)
v = FastHessian_getResponse(m, r, c, t);
dxx = FastHessian_getResponse(m,r, c + 1, t) + FastHessian_getResponse(m,r, c - 1, t) - 2 * v;
dyy = FastHessian_getResponse(m,r + 1, c, t) + FastHessian_getResponse(m,r - 1, c, t) - 2 * v;
dss = FastHessian_getResponse(t,r, c) + FastHessian_getResponse(b,r, c, t) - 2 * v;
dxy = (FastHessian_getResponse(m,r + 1, c + 1, t) - FastHessian_getResponse(m,r + 1, c - 1, t) - FastHessian_getResponse(m,r - 1, c + 1, t) + FastHessian_getResponse(m,r - 1, c - 1, t)) / 4;
dxs = (FastHessian_getResponse(t,r, c + 1) - FastHessian_getResponse(t,r, c - 1) - FastHessian_getResponse(b,r, c + 1, t) + FastHessian_getResponse(b,r, c - 1, t)) / 4;
dys = (FastHessian_getResponse(t,r + 1, c) - FastHessian_getResponse(t,r - 1, c) - FastHessian_getResponse(b,r + 1, c, t) + FastHessian_getResponse(b,r - 1, c, t)) / 4;

H = zeros(3,3);
H(1, 1) = dxx;
H(1, 2) = dxy;
H(1, 3) = dxs;
H(2, 1) = dxy;
H(2, 2) = dyy;
H(2, 3) = dys;
H(3, 1) = dxs;
H(3, 2) = dys;
H(3, 3) = dss;

function descriptor=SurfDescriptor_GetDescriptor(ip, bUpright, bExtended, img, verbose)
% This function SurfDescriptor_GetDescriptor will ..
%
% [descriptor] = SurfDescriptor_GetDescriptor( ip,bUpright,bExtended,img )
%  
%  inputs,
%    ip : Interest Point (x,y,scale, orientation)
%    bUpright : If true not rotation invariant descriptor
%    bExtended :  If true make a 128 values descriptor
%    img : Integral image
%    verbose : If true show additional information
%  
%  outputs,
%    descriptor : Descriptor of interest point length 64 or 128 (extended)  
%  
% Function is written by D.Kroon University of Twente (July 2010)

% Get rounded InterestPoint data
X = round(ip.x);
Y = round(ip.y);
S = round(ip.scale);

if (bUpright)
    co = 1;
    si = 0;
else
    co = cos(ip.orientation);
    si = sin(ip.orientation);
end

% Basis coordinates of samples, if coordinate 0,0, and scale 1
[lb,kb]=ndgrid(-4:4,-4:4); lb=lb(:); kb=kb(:);

%Calculate descriptor for this interest point
[jl,il]=ndgrid(0:3,0:3); il=il(:)'; jl=jl(:)';

ix = (il*5-8);
jx = (jl*5-8);

% 2D matrices instead of double for-loops, il, jl
cx=length(lb); cy=length(ix);
lb=repmat(lb,[1 cy]); lb=lb(:);
kb=repmat(kb,[1 cy]); kb=kb(:);
ix=repmat(ix,[cx 1]); ix=ix(:);
jx=repmat(jx,[cx 1]); jx=jx(:);

% Coordinates of samples (not rotated)
l=lb+jx; k=kb+ix;

%Get coords of sample point on the rotated axis
sample_x = round(X + (-l * S * si + k * S * co)); 
sample_y = round(Y + (l * S * co + k * S * si));
                
%Get the gaussian weighted x and y responses
xs = round(X + (-(jx+1) * S * si + (ix+1) * S * co));
ys = round(Y + ((jx+1) * S * co + (ix+1) * S * si));

gauss_s1 = SurfDescriptor_Gaussian(xs - sample_x, ys - sample_y, 2.5 * S);
rx = IntegralImage_HaarX(sample_y, sample_x, 2 * S,img);
ry = IntegralImage_HaarY(sample_y, sample_x, 2 * S,img);
                
%Get the gaussian weighted x and y responses on the aligned axis
rrx = gauss_s1 .* (-rx * si + ry * co);  rrx=reshape(rrx,cx,cy);
rry = gauss_s1 .* ( rx * co + ry * si);  rry=reshape(rry,cx,cy);
        
% Get the gaussian scaling
cx = -0.5 + il + 1; cy = -0.5 + jl + 1;
gauss_s2 = SurfDescriptor_Gaussian(cx - 2, cy - 2, 1.5);

if (bExtended)
    % split x responses for different signs of y
    check=rry >= 0; rrx_p=rrx.*check;  rrx_n=rrx.*(~check);
    
    dx = sum(rrx_p); mdx = sum(abs(rrx_p),1);
    dx_yn = sum(rrx_n); mdx_yn = sum(abs(rrx_n),1);

    % split y responses for different signs of x
    check=(rrx >= 0); rry_p=rry.*check; rry_n=rry.*(~check);
    dy = sum(rry_p,1); 
    mdy = sum(abs(rry_p),1);
    dy_xn = sum(rry_n,1); 
    mdy_xn =  sum(abs(rry_n),1);
else
    dx = sum(rrx,1);
    dy = sum(rry,1);
    mdx = sum(abs(rrx),1);
    mdy = sum(abs(rry),1);
    dx_yn = 0; mdx_yn = 0; 
    dy_xn = 0; mdy_xn = 0;
end

if (bExtended)
    descriptor=[dx;dy;mdx;mdy;dx_yn;dy_xn;mdx_yn;mdy_xn].* repmat(gauss_s2,[8 1]);
else
    descriptor=[dx;dy;mdx;mdy].* repmat(gauss_s2,[4 1]);
end
  
len = sum((dx.^2 + dy.^2 + mdx.^2 + mdy.^2 + dx_yn + dy_xn + mdx_yn + mdy_xn) .* gauss_s2.^2);

%Convert to Unit Vector
descriptor= descriptor(:) / sqrt(len);
if(verbose)
    for i=1:size(rrx,2)
        p1=reshape(rrx(:,i),[9,9]);
        p2=reshape(rry(:,i),[9,9]);
        p=[p1;ones(1,9)*0.02;p2];
        if(i==1)
            pic=p;
        else
            pic=[pic ones(19,1)*0.02 p];
        end
    end
    imshow(pic,[]);
end

function an= SurfDescriptor_Gaussian(x, y, sig)
an = 1 / (2 * pi * sig^2) .* exp(-(x.^2 + y.^2) / (2 * sig^2));

三、運行結果

在這里插入圖片描述

四、matlab版本及參考文獻

1 matlab版本
2014a

2 參考文獻
[1] 蔡利梅.MATLAB影像處理——理論、演算法與實體分析[M].清華大學出版社，2020.
[2]楊丹,趙海濱,龍哲.MATLAB影像處理實體詳解[M].清華大學出版社，2013.
[3]周品.MATLAB影像處理與圖形用戶界面設計[M].清華大學出版社，2013.
[4]劉成龍.精通MATLAB影像處理[M].清華大學出版社，2015.
[5]魏新,馬麗華,李云霞,徐志燕,李大為.改進配準測度的SUFT紅外影像快速配準演算法[J].電光與控制. 2012,19(11)

轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/357222.html

標籤：其他

上一篇：MATLAB用數學形態學法提取影像邊界

下一篇：Caught FileNotFoundError in DataLoader worker process 0