本文系《数字图像处理原理与实践（MATLAB版）》一书之代码系列的Part9，辑录该书第431至第438页之代码，供有需要读者下载研究使用。至此全书代码发布已经接近尾声，希望这些源码能够对有需要的读者有所帮助。代码执行结果请参见原书配图，建议下载代码前阅读下文：
关于《数字图像处理原理与实践（MATLAB版）》一书代码发布的说明
http://blog.csdn.net/baimafujinji/article/details/40987807
首先给出的是原书P438所列之程序源码，详情请参见原书之相关内容。
===========================================================
I1=imread(‘box.png‘);
I2=imread(‘box_in_scene.png‘);
Options.upright=true;
Options.tresh=0.0001;
Ipts1=OpenSurf(I1,Options);
Ipts2=OpenSurf(I2,Options);
D1 = reshape([Ipts1.descriptor],64,[]);
D2 = reshape([Ipts2.descriptor],64,[]);
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
[err, ind]=sort(err);
cor1=cor1(ind);
cor2=cor2(ind);
I = zeros([max(size(I1,1),size(I2,1)),size(I1,2)+size(I1,2)]);
I(1:size(I1,1),1:size(I1,2))=I1;
I(:,size(I1,2)+1:size(I1,2)+size(I2,2))=I2;
figure, imshow(I, []); hold on;
for i=1:30,
c=rand(1,3);
plot([Ipts1(cor1(i)).x Ipts2(cor2(i)).x+size(I1,2)],...
[Ipts1(cor1(i)).y Ipts2(cor2(i)).y],‘-‘,‘Color‘,c)
plot([Ipts1(cor1(i)).x Ipts2(cor2(i)).x+size(I1,2)],...
[Ipts1(cor1(i)).y Ipts2(cor2(i)).y],‘o‘,‘Color‘,c)
end
由于上述代码调用的函数所占篇幅过长，于是以下部分之代码并未辑录于书中，在此详细列出，以供有需要的读者参阅。特别说明，下面的程序源码由荷兰特温特大学的Dirk-Jan Kroon博士提供，使用条款请参见相关说明。
================================================================================
function ipts=OpenSurf(img,Options)
% This function OPENSURF, is an implementation of SURF (Speeded Up Robust
% Features). SURF will detect landmark points in an image, and describe
% the points by a vector which is robust against (a little bit) rotation
% ,scaling and noise. It can be used in the same way as SIFT (Scale-invariant
% feature transform) which is patented. Thus to align (register) two
% or more images based on corresponding points, or make 3D reconstructions.
%
% This Matlab implementation of Surf is a direct translation of the
% OpenSurf C# code of Chris Evans, and gives exactly the same answer.
% Chris Evans wrote one of the best, well structured all inclusive SURF
% implementations. On his site you can find Evaluations of OpenSURF
% and the C# and C++ code. http://www.chrisevansdev.com/opensurf/
% Chris Evans gave me permisson to publish this code under the (Mathworks)
% BSD license.
%
% Ipts = OpenSurf(I, Options)
%
% inputs,
% I : The 2D input image color or greyscale
% (optional)
% Options : A struct with options (see below)
%
% outputs,
% Ipts : A structure with the information about all detected Landmark points
% Ipts.x , ipts.y : The landmark position
% Ipts.scale : The scale of the detected landmark
% Ipts.laplacian : The laplacian of the landmark neighborhood
% Ipts.orientation : Orientation in radians
% Ipts.descriptor : The descriptor for corresponding point matching
%
% options,
% Options.verbose : If set to true then useful information is
% displayed (default false)
% Options.upright : Boolean which determines if we want a non-rotation
% invariant result (default false)
% Options.extended : Add extra landmark point information to the
% descriptor (default false)
% Options.tresh : Hessian response treshold (default 0.0002)
% Options.octaves : Number of octaves to analyse(default 5)
% Options.init_sample : Initial sampling step in the image (default 2)
%
% Example 1, Basic Surf Point Detection
% % Load image
% I=imread(‘TestImages/test.png‘);
% % Set this option to true if you want to see more information
% Options.verbose=false;
% % Get the Key Points
% Ipts=OpenSurf(I,Options);
% % Draw points on the image
% PaintSURF(I, Ipts);
%
% Example 2, Corresponding points
% % See, example2.m
%
% Example 3, Affine registration
% % See, example3.m
%
% Function is written by D.Kroon University of Twente (July 2010)
% Add subfunctions to Matlab Search path
functionname=‘OpenSurf.m‘;
functiondir=which(functionname);
functiondir=functiondir(1:end-length(functionname));
addpath([functiondir ‘/SubFunctions‘])
% Process inputs
defaultoptions=struct(‘tresh‘,0.0002,‘octaves‘,5,‘init_sample‘,2,‘upright‘,false,‘extended‘,false,‘verbose‘,false);
if(~exist(‘Options‘,‘var‘)),
Options=defaultoptions;
else
tags = fieldnames(defaultoptions);
for i=1:length(tags)
if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end
end
if(length(tags)~=length(fieldnames(Options))),
warning(‘register_volumes:unknownoption‘,‘unknown options found‘);
end
end
% Create Integral Image
iimg=IntegralImage_IntegralImage(img);
% Extract the interest points
FastHessianData.thresh = Options.tresh;
FastHessianData.octaves = Options.octaves;
FastHessianData.init_sample = Options.init_sample;
FastHessianData.img = iimg;
ipts = FastHessian_getIpoints(FastHessianData,Options.verbose);
% Describe the interest points
if(~isempty(ipts))
ipts = SurfDescriptor_DecribeInterestPoints(ipts,Options.upright, Options.extended, iimg, Options.verbose);
end
================================================================================
function D=FastHessian_BuildDerivative(r,c,t,m,b)
% This function FastHessian_BuildDerivative will ..
%
% [D] = FastHessian_BuildDerivative( r,c,t,m,b )
%
% inputs,
% r :
% c :
% t :
% m :
% b :
%
% outputs,
% D :
%
% Function is written by D.Kroon University of Twente ()
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 rl=FastHessian_buildResponseLayer(rl,FastHessianData)
% This function FastHessian_buildResponseLayer will ..
%
% [rl] = FastHessian_buildResponseLayer( rl,FastHessianData )
%
% inputs,
% rl :
% FastHessianData :
%
% outputs,
% rl :
%
% Function is written by D.Kroon University of Twente ()
step = fix( rl.step); % step size for this filter
b = fix((rl.filter - 1) / 2 + 1); % border for this filter
l = fix(rl.filter / 3); % lobe for this filter (filter size / 3)
w = fix(rl.filter); % filter size
inverse_area = 1 / double(w * w); % normalisation factor
img=FastHessianData.img;
[ac,ar]=ndgrid(0:rl.width-1,0:rl.height-1);
ar=ar(:); ac=ac(:);
% get the image coordinates
r = int32(ar * step);
c = int32(ac * step);
% Compute response components
Dxx = IntegralImage_BoxIntegral(r - l + 1, c - b, 2 * l - 1, w,img) - IntegralImage_BoxIntegral(r - l + 1, c - fix(l / 2), 2 * l - 1, l, img) * 3;
Dyy = IntegralImage_BoxIntegral(r - b, c - l + 1, w, 2 * l - 1,img) - IntegralImage_BoxIntegral(r - fix(l / 2), c - l + 1, l, 2 * l - 1,img) * 3;
Dxy = + IntegralImage_BoxIntegral(r - l, c + 1, l, l,img) + IntegralImage_BoxIntegral(r + 1, c - l, l, l,img) ...
- IntegralImage_BoxIntegral(r - l, c - l, l, l,img) - IntegralImage_BoxIntegral(r + 1, c + 1, l, l,img);
% Normalise the filter responses with respect to their size
Dxx = Dxx*inverse_area;
Dyy = Dyy*inverse_area;
Dxy = Dxy*inverse_area;
% Get the determinant of hessian response & laplacian sign
rl.responses = (Dxx .* Dyy - 0.81 * Dxy .* Dxy);
rl.laplacian = (Dxx + Dyy) >= 0;
================================================================================
function an=IntegralImage_BoxIntegral(row, col, rows,cols,img)
% Get integer coordinates
row=fix(row);
col=fix(col);
rows=fix(rows);
cols=fix(cols);
% Get the corner coordinates of the box integral
r1 = min(row, size(img,1));
c1 = min(col, size(img,2));
r2 = min(row + rows, size(img,1));
c2 = min(col + cols, size(img,2));
% Get the values at the cornes of the box integral (fast 1D index look up)
sx=size(img,1);
A = img(max(r1+(c1-1)*sx,1));
B = img(max(r1+(c2-1)*sx,1));
C = img(max(r2+(c1-1)*sx,1));
D = img(max(r2+(c2-1)*sx,1));
% If coordinates are outside at the top or left, the value must be zero
A((r1<1)|(c1<1))=0;
B((r1<1)|(c2<1))=0;
C((r2<1)|(c1<1))=0;
D((r2<1)|(c2<1))=0;
% Minimum value of the integral is zero
an=max(0, A - B - C + D);
================================================================================
function ResponseLayerData=FastHessian_ResponseLayer(width, height, step, filter)
% This function FastHessian_ResponseLayer will ..
%
% [ResponseLayerData] = FastHessian_ResponseLayer( width,height,step,filter )
%
% inputs,
% width :
% height :
% step :
% filter :
%
% outputs,
% ResponseLayerData :
%
% Function is written by D.Kroon University of Twente (July 2010)
width=floor(width);
height=floor(height);
step=floor(step);
filter=floor(filter);
ResponseLayerData.width = width;
ResponseLayerData.height = height;
ResponseLayerData.step = step;
ResponseLayerData.filter = filter;
ResponseLayerData.responses = zeros(width * height,1);
ResponseLayerData.laplacian = zeros(width * height,1);
================================================================================
function an=FastHessian_isExtremum(r, c, t, m, b,FastHessianData)
% This function FastHessian_isExtremum will ..
%
% [an] = FastHessian_isExtremum( r,c,t,m,b,FastHessianData )
%
% inputs,
% r :
% c :
% t :
% m :
% b :
% FastHessianData :
%
% outputs,
% an :
%
% Function is written by D.Kroon University of Twente (July 2010)
% bounds check
layerBorder = fix((t.filter + 1) / (2 * t.step));
bound_check_fail=(r <= layerBorder | r >= t.height - layerBorder | c <= layerBorder | c >= t.width - layerBorder);
% check the candidate point in the middle layer is above thresh
candidate = FastHessian_getResponse(m,r,c,t);
treshold_fail=candidate < FastHessianData.thresh;
an=(~bound_check_fail)&(~treshold_fail);
for rr = -1:1
for cc = -1:1
% if any response in 3x3x3 is greater then the candidate is not a maximum
check1=FastHessian_getResponse(t,r + rr, c + cc, t) >= candidate;
check2=FastHessian_getResponse(m,r + rr, c + cc, t) >= candidate;
check3=FastHessian_getResponse(b,r + rr, c + cc, t) >= candidate;
check4=(rr ~= 0 || cc ~= 0);
an3 = ~(check1 | (check4 & check2) | check3);
an=an&an3;
end
end
function an=FastHessian_getResponse(a,row, column,b)
scale=fix(a.width/b.width);
% Clamp to boundary
% (The orignal C# code, doesn‘t contain this boundary clamp because if you
% process one coordinate at the time you already returned on the boundary check)
index=fix(scale*row) * a.width + fix(scale*column)+1;
index(index<1)=1; index(index>length(a.responses))=length(a.responses);
an=a.responses(index);
================================================================================
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 an=FastHessian_getResponse(a,row, column,b)
% This function FastHessian_getResponse will ..
%
% [an] = FastHessian_getResponse( a,row,column,b )
%
% inputs,
% a :
% row :
% column :
% b :
%
% outputs,
% an :
%
% Function is written by D.Kroon University of Twente ()
if(nargin<4)
scale=1;
else
scale=fix(a.width/b.width);
end
an=a.responses(fix(scale*row) * a.width + fix(scale*column)+1);
================================================================================
function an=FastHessian_getLaplacian(a,row, column,b)
% This function FastHessian_getLaplacian will ..
%
% [an] = FastHessian_getLaplacian( a,row,column,b )
%
% inputs,
% a :
% row :
% column :
% b :
%
% outputs,
% an :
%
% Function is written by D.Kroon University of Twente ()
if(nargin<4)
scale=1;
else
scale=fix(a.width/b.width);
end
an=a.laplacian(fix(scale*row) * a.width + fix(scale*column)+1);
================================================================================
function ipts=FastHessian_getIpoints(FastHessianData,verbose)
% filter index map
filter_map = [0,1,2,3;
1,3,4,5;
3,5,6,7;
5,7,8,9;
7,9,10,11]+1;
np=0; ipts=struct;
% Build the response map
responseMap=FastHessian_buildResponseMap(FastHessianData);
% Find the maxima acrros scale and space
for o = 1:FastHessianData.octaves
for i = 1:2
b = responseMap{filter_map(o,i)};
m = responseMap{filter_map(o,i+1)};
t = responseMap{filter_map(o,i+2)};
% loop over middle response layer at density of the most
% sparse layer (always top), to find maxima across scale and space
[c,r]=ndgrid(0:t.width-1,0:t.height-1);
r=r(:); c=c(:);
p=find(FastHessian_isExtremum(r, c, t, m, b,FastHessianData));
for j=1:length(p);
ind=p(j);
[ipts,np]=FastHessian_interpolateExtremum(r(ind), c(ind), t, m, b, ipts,np);
end
end
end
% Show laplacian and response maps with found interest-points
if(verbose)
% Show the response map
if(verbose)
fig_h=ceil(length(responseMap)/3);
h=figure; set(h,‘name‘,‘Laplacian‘);
for i=1:length(responseMap),
pic=reshape(responseMap{i}.laplacian,[responseMap{i}.width responseMap{i}.height]);
subplot(3,fig_h,i); imshow(pic,[]); hold on;
end
h=figure; set(h,‘name‘,‘Responses‘);
h_res=zeros(1,length(responseMap));
for i=1:length(responseMap),
pic=reshape(responseMap{i}.responses,[responseMap{i}.width responseMap{i}.height]);
h_res(i)=subplot(3,fig_h,i); imshow(pic,[]); hold on;
end
end
% Show the maximum points
disp([‘Number of interest points found ‘ num2str(np)]);
scales=zeros(1,length(responseMap));
scaley=zeros(1,length(responseMap));
scalex=zeros(1,length(responseMap));
for i=1:length(responseMap)
scales(i)=responseMap{i}.filter*(2/15);
scalex(i)=responseMap{i}.width/size(FastHessianData.img,2);
scaley(i)=responseMap{i}.height/size(FastHessianData.img,1);
end
for i=1:np
[t,ind]=min((scales-ipts(i).scale).^2);
plot(h_res(ind),ipts(i).y*scaley(ind)+1,ipts(i).x*scalex(ind)+1,‘o‘,‘color‘,rand(1,3));
end
end
================================================================================
function responseMap=FastHessian_buildResponseMap(FastHessianData)
% Calculate responses for the first 4 FastHessianData.octaves:
% Oct1: 9, 15, 21, 27
% Oct2: 15, 27, 39, 51
% Oct3: 27, 51, 75, 99
% Oct4: 51, 99, 147,195
% Oct5: 99, 195,291,387
% Deallocate memory and clear any existing response layers
responseMap=[]; j=0;
% Get image attributes
w = (size(FastHessianData.img,2) / FastHessianData.init_sample);
h = (size(FastHessianData.img,1)/ FastHessianData.init_sample);
s = (FastHessianData.init_sample);
% Calculate approximated determinant of hessian values
if (FastHessianData.octaves >= 1)
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w, h, s, 9);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w, h, s, 15);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w, h, s, 21);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w, h, s, 27);
end
if (FastHessianData.octaves >= 2)
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 2, h / 2, s * 2, 39);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 2, h / 2, s * 2, 51);
end
if (FastHessianData.octaves >= 3)
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 4, h / 4, s * 4, 75);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 4, h / 4, s * 4, 99);
end
if (FastHessianData.octaves >= 4)
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 8, h / 8, s * 8, 147);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 8, h / 8, s * 8, 195);
end
if (FastHessianData.octaves >= 5)
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 16, h / 16, s * 16, 291);
j=j+1; responseMap{j}=FastHessian_ResponseLayer(w / 16, h / 16, s * 16, 387);
end
% Extract responses from the image
for i=1:length(responseMap);
responseMap{i}=FastHessian_buildResponseLayer(responseMap{i},FastHessianData);
end
================================================================================
function an=IntegralImage_HaarX(row, column, size, img)
% This function IntegralImage_HaarX will ..
%
% [an] = IntegralImage_HaarX( row,column,size,img )
%
% inputs,
% row :
% column :
% size :
% img :
%
% outputs,
% an :
%
% Function is written by D.Kroon University of Twente (July 2010)
an= IntegralImage_BoxIntegral(row - size / 2, column, size, size / 2, img) - IntegralImage_BoxIntegral(row - size / 2, column - size / 2, size, size / 2, img);
================================================================================
function an=IntegralImage_HaarY(row, column, size, img)
% This function IntegralImage_HaarY will ..
%
% [an] = IntegralImage_HaarY( row,column,size,img )
%
% inputs,
% row :
% column :
% size :
% img : The integral image
%
% outputs,
% an : The haar response in y-direction
%
% Function is written by D.Kroon University of Twente (July 2010)
an= IntegralImage_BoxIntegral(row, column - size / 2, size / 2, size , img) - IntegralImage_BoxIntegral(row - size / 2, column - size / 2, size / 2, size , img);
================================================================================
function pic=IntegralImage_IntegralImage(I)
% This function IntegralImage_IntegralImage will ..
%
% J = IntegralImage_IntegralImage( I )
%
% inputs,
% I : An 2D image color or greyscale
%
% outputs,
% J : The integral image
%
% Function is written by D.Kroon University of Twente (July 2010)
% Convert Image to double
switch(class(I));
case ‘uint8‘
I=double(I)/255;
case ‘uint16‘
I=double(I)/65535;
case ‘int8‘
I=(double(I)+128)/255;
case ‘int16‘
I=(double(I)+32768)/65535;
otherwise
I=double(I);
end
% Convert Image to greyscale
if(size(I,3)==3)
cR = .2989; cG = .5870; cB = .1140;
I=I(:,:,1)*cR+I(:,:,2)*cG+I(:,:,3)*cB;
end
% Make the integral image
pic = cumsum(cumsum(I,1),2);
================================================================================
function PaintSURF(I, ipts)
% This function PaintSURF will display the image with the found Interest points
%
% [] = PaintSURF( img,ipts )
%
% inputs,
% img : Image 2D color or greyscale
% ipts : The interest points
%
% Function is written by D.Kroon University of Twente (July 2010)
% Convert Image to double
switch(class(I));
case ‘uint8‘
I=double(I)/255;
case ‘uint16‘
I=double(I)/65535;
case ‘int8‘
I=(double(I)+128)/255;
case ‘int16‘
I=(double(I)+32768)/65535;
otherwise
I=double(I);
end
figure, imshow(I), hold on;
if (isempty(fields(ipts))), return; end
for i=1:length(ipts)
ip=ipts(i);
S = 2 * fix(2.5 * ip.scale);
R = fix(S / 2);
pt = [(ip.x), (ip.y)];
ptR = [(R * cos(ip.orientation)), (R * sin(ip.orientation))];
if(ip.laplacian >0), myPen =[0 0 1]; else myPen =[1 0 0]; end
rectangle(‘Curvature‘, [1 1],‘Position‘, [pt(1)-R, pt(2)-R, S, S],‘EdgeColor‘,myPen);
plot([pt(1), pt(1)+ptR(1)]+1,[pt(2), pt(2)+ptR(2)]+1,‘g‘);
end
================================================================================
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));
================================================================================
function pic=IntegralImage_IntegralImage(I)
% This function IntegralImage_IntegralImage will ..
%
% J = IntegralImage_IntegralImage( I )
%
% inputs,
% I : An 2D image color or greyscale
%
% outputs,
% J : The integral image
%
% Function is written by D.Kroon University of Twente (July 2010)
% Convert Image to double
switch(class(I));
case ‘uint8‘
I=double(I)/255;
case ‘uint16‘
I=double(I)/65535;
case ‘int8‘
I=(double(I)+128)/255;
case ‘int16‘
I=(double(I)+32768)/65535;
otherwise
I=double(I);
end
% Convert Image to greyscale
if(size(I,3)==3)
cR = .2989; cG = .5870; cB = .1140;
I=I(:,:,1)*cR+I(:,:,2)*cG+I(:,:,3)*cB;
end
% Make the integral image
pic = cumsum(cumsum(I,1),2);
================================================================================
function ipts = SurfDescriptor_DecribeInterestPoints(ipts, upright, extended, img, verbose)
% This function SurfDescriptor_DecribeInterestPoints will ..
%
% [ipts] = SurfDescriptor_DecribeInterestPoints( ipts,upright,extended,img )
%
% inputs,
% ipts : Interest Points (x,y,scale)
% bUpright : If true not rotation invariant descriptor
% bExtended : If true make a 128 values descriptor
% img : Integral image
% verbose : If true show useful information
%
% outputs,
% ipts : Interest Points (x,y,orientation,descriptor)
%
% Function is written by D.Kroon University of Twente (July 2010)
if (isempty(fields(ipts))), return; end
if(verbose), h_ang=figure; drawnow, set(h_ang,‘name‘,‘Angles‘); else h_ang=[]; end
if(verbose), h_des=figure; drawnow, set(h_des,‘name‘,‘Aligned Descriptor XY‘); end
for i=1:length(ipts)
% Display only information about the first 40 points
if(i>40), verbose=false; end
ip=ipts(i);
% determine descriptor size
if (extended), ip.descriptorLength = 128; else ip.descriptorLength = 64; end
% Get the orientation
if(verbose), figure(h_ang), subplot(5,8,i), end
ip.orientation=SurfDescriptor_GetOrientation(ip,img,verbose);
% Extract SURF descriptor
if(verbose), figure(h_des), subplot(10,4,i), end
ip.descriptor=SurfDescriptor_GetDescriptor(ip, upright, extended, img,verbose);
ipts(i).orientation=ip.orientation;
ipts(i).descriptor=ip.descriptor;
end
if(~isempty(h_ang)), figure(h_ang), colormap(jet); end
================================================================================
function orientation=SurfDescriptor_GetOrientation(ip,img,verbose)
% This function SurfDescriptor_GetOrientation will ..
%
% [orientation] = SurfDescriptor_GetOrientation( ip,img )
%
% inputs,
% ip : InterestPoint data, (x,y,scale)
% img : Integral Image
% verbose : If true, show additional information
%
% outputs,
% orientation : Orientation of intereset point (radians)
%
% Function is written by D.Kroon University of Twente (July 2010)
gauss25 = [0.02350693969273 0.01849121369071 0.01239503121241 0.00708015417522 0.00344628101733 0.00142945847484 0.00050524879060;
0.02169964028389 0.01706954162243 0.01144205592615 0.00653580605408 0.00318131834134 0.00131955648461 0.00046640341759;
0.01706954162243 0.01342737701584 0.00900063997939 0.00514124713667 0.00250251364222 0.00103799989504 0.00036688592278;
0.01144205592615 0.00900063997939 0.00603330940534 0.00344628101733 0.00167748505986 0.00069579213743 0.00024593098864;
0.00653580605408 0.00514124713667 0.00344628101733 0.00196854695367 0.00095819467066 0.00039744277546 0.00014047800980;
0.00318131834134 0.00250251364222 0.00167748505986 0.00095819467066 0.00046640341759 0.00019345616757 0.00006837798818;
0.00131955648461 0.00103799989504 0.00069579213743 0.00039744277546 0.00019345616757 0.00008024231247 0.00002836202103];
gauss25=gauss25(:);
% Get rounded InterestPoint data
X = round(ip.x);
Y = round(ip.y);
S = round(ip.scale);
% calculate haar responses for points within radius of 6*scale
[j,i]=ndgrid(-6:6,-6:6);
j=j(:); i=i(:); check=(i.^2 + j.^2 < 36); j=j(check); i=i(check);
% Get gaussian filter (by mirroring gauss25)
id = [ 6, 5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5, 6 ];
gauss = gauss25(id(i + 6 + 1) + id(j + 6 + 1) *7+1);
resX = gauss .* IntegralImage_HaarX(Y + j * S, X + i * S, 4 * S, img);
resY = gauss .* IntegralImage_HaarY(Y + j * S, X + i * S, 4 * S, img);
Ang = mod(atan2(resY, resX),2*pi);
% loop slides pi/3 window around feature point
ang1 = 0:0.15:(2 * pi);
ang2 = mod(ang1+pi/3,2*pi);
% Repmat is used to check for all angles (x direction) and
% all responses (y direction) without for-loops.
cx=length(Ang); cy=length(ang1);
ang1=repmat(ang1,[cx 1]);
ang2=repmat(ang2,[cx 1]);
Ang =repmat(Ang,[1 cy]);
resX =repmat(resX,[1 cy]);
resY =repmat(resY,[1 cy]);
% determine whether the point is within the window
check1= (ang1 < ang2) & (ang1 < Ang) & (Ang < ang2);
check2= (ang2 < ang1) & ( ((Ang > 0) & (Ang < ang2)) | ((Ang > ang1) & (Ang < pi)) );
check=check1|check2;
sumX = sum(resX.*check,1);
sumY = sum(resY.*check,1);
% Find the most dominant direction
R=sumX.^2+ sumY.^2;
[t,ind]=max(R);
orientation = mod(atan2(sumY(ind), sumX(ind)),2*pi);
if(verbose)
pica=zeros(13,13);
pica(i+7+(j+6)*13)=Ang(:,1);
imshow(pica,[0 2*pi]);
end
《数字图像处理原理与实践（MATLAB版）》一书之代码Part9
原文地址：http://blog.csdn.net/baimafujinji/article/details/44815599
