Code covered by the BSD License

# 2D - 2D Projective Homography (3x3) Estimation

### SasiKanth (view profile)

19 Sep 2010 (Updated )

This function estimates 2D-2D projective homography between two images.

homography( i1, i2)
```function [ h wim1 ] = homography( i1, i2)

% This function estimates 2D-2D plane projective homography between two
% perspective images using Direct Linear Transformation, RANSAC
% and Levenberg Marquardt optimisation.

% The format for calling upon the function is as follows:
%
% [h wim] = homography(im1, im2);
%
% where
%
% im1 -> 1st Image
% im2 -> 2nd Image
% h -> Returned homography matrix
% wim -> Warped version of im1 w.r.t. im2

% Code written by B S SasiKanth (bsasikanth@gmail.com) of the Indian
% Institute of Technology (IIT), Guwahati.

% Switch off warnings for non-singular cases of RANSAC to be encountered

warning off all;

% Finding Correspondences

a = matching(i1,i2);

% Computing Homography

h = autohomest2d(a);

% Computing Left Warped Image

wim1 = (imTrans(i1,(h)));
wim1 = imresize(wim1,[size(i1,1) size(i1,2)],'bicubic');

end

%%

function [ hmatrixauto1 ] = autohomest2d( correspondence_input )
%AUTOHOMEST2D This function computes the automatic homography between two
%perspective images using RANSAC and Levenberg Marquardt optimization.

corr_req = 4;
[corr_num, tempx, tempx] = size(correspondence_input);
max_no_inliers = 0; p = 0.99; iter_no = 0;

while(1)

iter_no = iter_no + 1;
inlier_detail = zeros(1,corr_num);

% Random number generation

stream = RandStream('mt19937ar','seed',sum(100*clock));
randomvector = rand(stream,1,corr_req);

% Choosing the correspondences

selected_correspondences = zeros(corr_req,2,2);
selected_indices = ceil(corr_num*randomvector);

for i=1:corr_req

selected_correspondences(i,:,:) = correspondence_input(selected_indices(i),:,:);

end

% Computing the temporary H matrix

hmatrixtemp = homest2d(selected_correspondences);
hmatrixtempi = inv(hmatrixtemp);

% Calculating the symmetric cost and checking for inliers

no_inliers = 0;

for i = 1:corr_num

temp1 = hmatrixtemp*[correspondence_input(i,1,1); correspondence_input(i,1,2); 1];
temp1 = temp1/temp1(3);

temp2 = hmatrixtempi*[correspondence_input(i,2,1); correspondence_input(i,2,2); 1]; %#ok<MINV>
temp2 = temp2/temp2(3);

error1 = sum(abs([correspondence_input(i,2,1); correspondence_input(i,2,2); 1] - temp1));
error2 = sum(abs([correspondence_input(i,1,1); correspondence_input(i,1,2); 1] - temp2));

error = error1 + error2;

if(error < 20)
no_inliers = no_inliers + 1;
inlier_detail(i) = 1;
end

end

if(no_inliers > max_no_inliers)
max_no_inliers = no_inliers;
max_inlier_detail = inlier_detail;
hmatrixauto = hmatrixtemp;
end

% Dynamic updatation of the no of iterations, N

w = max_no_inliers/corr_num;
N = abs((log(1-p))/(log(1-w^4)));

if(N < iter_no)
break;
end

end

% Isolating inlier correspondences

inlier_corr = zeros(max_no_inliers,2,2);
k = 1;
for i=1:corr_num

if(max_inlier_detail(i) == 1)
inlier_corr(k,:,:) = correspondence_input(i,:,:);
k = k+1;
end

end

% LevMar Optimization

hvector(1:3) = hmatrixauto(1,:);
hvector(4:6) = hmatrixauto(2,:);
hvector(7:9) = hmatrixauto(3,:);

ydata = zeros(1,max_no_inliers);

hoptim = lsqcurvefit(@myfun,hvector,inlier_corr,ydata);

% Reassigning the H matrix

hmatrixauto1 = [hoptim(1), hoptim(2), hoptim(3); hoptim(4), hoptim(5), hoptim(6); hoptim(7), hoptim(8), hoptim(9)];
hmatrixauto1 = hmatrixauto1./hmatrixauto1(3,3);

% clc;
%
% disp('The computation exited in ');
% disp(iter_no);
% disp(' iterations');

end

%%
function [ hmatrix ] = homest2d( correspondence_input )
%HOMEST2D This function computes the 2 dimensional homography for two
%different perspective images which have been inputted.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FORMAT OF THE CORRESPONDENCE MATRIX IS AS FOLLOWS:

% correspondence_input[i][j][k]

% i = correspondence_index

% j = 1 for image1
% j = 2 for image2

% k = 1 for xcoordinate
% k = 2 for ycoordinate

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[corr_num, tempx, tempx] = size(correspondence_input);

if(corr_num < 4)
error('Invalid no. of correspondences for homography estimation');
end

% NORMALIZE THE CORRESPONDENCES

% Find the centroid of the correspondences

xcentroid1 = 0;
ycentroid1 = 0;
xcentroid2 = 0;
ycentroid2 = 0;

for i = 1:corr_num

xcentroid1 = xcentroid1 + correspondence_input(i,1,1);
ycentroid1 = ycentroid1 + correspondence_input(i,1,2);
xcentroid2 = xcentroid2 + correspondence_input(i,2,1);
ycentroid2 = ycentroid2 + correspondence_input(i,2,2);

end

xcentroid1 = xcentroid1/corr_num;
ycentroid1 = ycentroid1/corr_num;
xcentroid2 = xcentroid2/corr_num;
ycentroid2 = ycentroid2/corr_num;

% Generating the Translation Matrices

tr1 = [1, 0, -xcentroid1; 0, 1, -ycentroid1; 0, 0, 1];
tr2 = [1, 0, -xcentroid2; 0, 1, -ycentroid2; 0, 0, 1];

% Generating Translated Correspondences

correspondence_input_tr = zeros(corr_num,2,2);

for i = 1:corr_num

temp_res = tr1*[correspondence_input(i,1,1); correspondence_input(i,1,2); 1];
correspondence_input_tr(i,1,1) = temp_res(1);
correspondence_input_tr(i,1,2) = temp_res(2);

temp_res = tr2*[correspondence_input(i,2,1); correspondence_input(i,2,2); 1];
correspondence_input_tr(i,2,1) = temp_res(1);
correspondence_input_tr(i,2,2) = temp_res(2);

end

% Computing the average RMS distance of translated correspondences

rms1 = 0;
rms2 = 0;

for i = 1:corr_num

rms1 = rms1 + sqrt(correspondence_input_tr(i,1,1)^2 + correspondence_input_tr(i,1,2)^2);
rms2 = rms2 + sqrt(correspondence_input_tr(i,2,1)^2 + correspondence_input_tr(i,2,2)^2);

end

rms1 = rms1/corr_num;
rms2 = rms2/corr_num;

% Generating the scaling matrices

scale_factor1 = sqrt(2)/rms1;
scale_factor2 = sqrt(2)/rms2;

sc1 = [scale_factor1, 0, 0; 0, scale_factor1, 0; 0, 0, 1];
sc2 = [scale_factor2, 0, 0; 0, scale_factor2, 0; 0, 0, 1];

% Normalizing correspondences to average RMS distance sqrt(2)

correspondence_input_norm = zeros(corr_num,2,2);

for i = 1:corr_num

temp_res = sc1*[correspondence_input_tr(i,1,1); correspondence_input_tr(i,1,2); 1];
correspondence_input_norm(i,1,1) = temp_res(1);
correspondence_input_norm(i,1,2) = temp_res(2);

temp_res = sc2*[correspondence_input_tr(i,2,1); correspondence_input_tr(i,2,2); 1];
correspondence_input_norm(i,2,1) = temp_res(1);
correspondence_input_norm(i,2,2) = temp_res(2);

end

% Generating Transformation matrices and required inverse for
% denormalization

T1 = sc1*tr1;
T2 = sc2*tr2;

% T2i = inv(T2);

% Generation of the A matrix for computation of the hmatrix

A = zeros(2*corr_num, 9);

for i = 1:corr_num

%A(2*i-1,:) = [0, 0, 0, -correspondence_input_norm(i,1,1), -correspondence_input_norm(i,1,2), -1, correspondence_input_norm(i,2,2)*correspondence_input_norm(i,1,1), correspondence_input_norm(i,2,2)*correspondence_input_norm(i,1,2), correspondence_input_norm(i,2,2)];
%A(2*i,:) = [correspondence_input_norm(i,1,1), correspondence_input_norm(i,1,2), 1, 0, 0, 0, -correspondence_input_norm(i,2,1)*correspondence_input_norm(i,1,1), -correspondence_input_norm(i,2,1)*correspondence_input_norm(i,1,2), -correspondence_input_norm(i,2,1)];
A(2*i-1,:) = [correspondence_input_norm(i,1,1), correspondence_input_norm(i,1,2), 1, 0, 0, 0, -correspondence_input_norm(i,2,1)*correspondence_input_norm(i,1,1), -correspondence_input_norm(i,2,1)*correspondence_input_norm(i,1,2), -correspondence_input_norm(i,2,1)];
A(2*i,:) = [0, 0, 0, correspondence_input_norm(i,1,1), correspondence_input_norm(i,1,2), 1, -correspondence_input_norm(i,2,2)*correspondence_input_norm(i,1,1), -correspondence_input_norm(i,2,2)*correspondence_input_norm(i,1,2), -correspondence_input_norm(i,2,2)];

end

% Solving the linear equation Ax = 0 using the SVD

if ((sum(sum(isnan(A))) + (sum(sum(isinf(A))) ~= 0)))
hvector = zeros(1,9);
else
[tempx, tempx, V] = svd(A, 0);
hvector = V(:,9);
end

hmatrix_norm = [hvector(1), hvector(2), hvector(3); hvector(4), hvector(5), hvector(6); hvector(7), hvector(8), hvector(9)];

% Denormalizing the hmatrix_norm

hmatrix_temp = T2\hmatrix_norm;
hmatrix = hmatrix_temp*T1;

end

%%

function [ ydata ] = myfun( x, xdata )
%MYFUN LevMar cost calculator

ydata = zeros(1,length(xdata));

hmatrixtemp = [x(1), x(2), x(3); x(4), x(5), x(6); x(7), x(8), x(9)];

for i = 1:length(xdata)

temp1 = hmatrixtemp*[xdata(i,1,1); xdata(i,1,2); 1];
temp1 = temp1/temp1(3);

temp2 = hmatrixtemp\[xdata(i,2,1); xdata(i,2,2); 1];
temp2 = temp2/temp2(3);

error1 = sum(abs([xdata(i,2,1); xdata(i,2,2); 1] - temp1));
error2 = sum(abs([xdata(i,1,1); xdata(i,1,2); 1] - temp2));

ydata(i) = error1 + error2;

end

end

%%

function [ a ] = matching( I1, I2 )
%MATCHING SURF matching

% Get the Key Points
Options.upright=true;
Options.tresh=0.000001;
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
[tempx, ind]=sort(err);
cor1=cor1(ind);
cor2=cor2(ind);

a = zeros(length(Ipts1),2,2);

for i=1:length(Ipts1)
a(i,1,1) = Ipts1(cor1(i)).x;
a(i,2,1) = Ipts2(cor2(i)).x;
a(i,1,2) = Ipts1(cor1(i)).y;
a(i,2,2) = Ipts2(cor2(i)).y;
end

end

%%

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)
%
% 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
% % 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));

% 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

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];
end

%%

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;

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

end

%%

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
end

%%

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);

end

%%

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);
end

%%

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];
end

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;

end
end

%%

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);
end
end

%%

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);
end

%%

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);

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);
end

%%

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);

end

%%

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);
end

%%

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
end

%%

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
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));

end
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

end

%%

% IMTRANS - Homogeneous transformation of an image.
%
% Applies a geometric transform to an image
%
%  [newim, newT] = imTrans(im, T, region, sze);
%
%  Arguments:
%        im     - The image to be transformed.
%        T      - The 3x3 homogeneous transformation matrix.
%        region - An optional 4 element vector specifying
%                 [minrow maxrow mincol maxcol] to transform.
%                 This defaults to the whole image if you omit it
%                 or specify it as an empty array [].
%        sze    - An optional desired size of the transformed image
%                 (this is the maximum No of rows or columns).
%                 This defaults to the maximum of the rows and columns
%                 of the original image.
%
%  Returns:
%        newim  - The transformed image.
%        newT   - The transformation matrix that relates transformed image
%                 coordinates to the reference coordinates for use in a
%                 function such as DIGIPLANE.
%
%  The region argument is used when one is inverting a perspective
%  transformation of a plane and the vanishing line of the plane lies within
%  the image.  Attempts to transform any part of the vanishing line will
%  position you at infinity.  Accordingly one should specify a region that
%  excludes any part of the vanishing line.
%
%  The sze parameter is optionally used to control the size of the
%  output image.  When inverting a perpective or affine transformation
%  the scale parameter is unknown/arbitrary, and without specifying
%  it explicitly the transformed image can end up being very small
%  or very large.
%
%  Problems: If your transformed image ends up as being two small bits of
%  image separated by a large black area then the chances are that you have
%  included the vanishing line of the plane within the specified region to
%  transform.  If your image degenerates to a very thin triangular shape
%  part of your region is probably very close to the vanishing line of the
%  plane.

% Copyright (c) 2000-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% April 2000 - original version.
% July 2001  - transformation of region boundaries corrected.

function [newim, newT] = imTrans(im, T, region, sze);

if isa(im,'uint8')
im = double(im);  % Make sure image is double
end

% Set up default region and transformed image size values
if ndims(im) == 3
[rows cols depth] = size(im);
else
[rows cols] = size(im);
depth = 1;
end

if nargin == 2
region = [1 rows 1 cols];
sze = max([rows cols]);
elseif nargin == 3
sze = max([rows cols]);
end

if isempty(region)
region = [1 rows 1 cols];
end

threeD = (ndims(im)==3);  % A colour image
if threeD    % Transform red, green, blue components separately
im = im/255;
[r, newT] = transformImage(im(:,:,1), T, region, sze);
[g, newT] = transformImage(im(:,:,2), T, region, sze);
[b, newT] = transformImage(im(:,:,3), T, region, sze);

newim = repmat(uint8(0),[size(r),3]);
newim(:,:,1) = uint8(round(r*255));
newim(:,:,2) = uint8(round(g*255));
newim(:,:,3) = uint8(round(b*255));

else                % Assume the image is greyscale
[newim, newT] = transformImage(im, T, region, sze);
end
end

%------------------------------------------------------------

% The internal function that does all the work

function [newim, newT] = transformImage(im, T, region, sze);

[rows, cols] = size(im);

if 0
% Determine default parameters if needed
if nargin == 2
region = [1 rows 1 cols];
sze = max(rows,cols);
elseif nargin == 3
sze = max(rows,cols);
elseif nargin ~= 4
error('Incorrect arguments to imtrans');
end
end
% Cut the image down to the specified region
%if nargin == 3 | nargin == 4
im = im(region(1):region(2), region(3):region(4));
[rows, cols] = size(im);
%end

% Find where corners go - this sets the bounds on the final image
B = bounds(T,region);
nrows = B(2) - B(1);
ncols = B(4) - B(3);

% Determine any rescaling needed
s = sze/max(nrows,ncols);

S = [s 0 0        % Scaling matrix
0 s 0
0 0 1];

T = S*T;
Tinv = inv(T);

% Recalculate the bounds of the new (scaled) image to be generated
B = bounds(T,region);
nrows = B(2) - B(1);
ncols = B(4) - B(3);

% Construct a transformation matrix that relates transformed image
% coordinates to the reference coordinates for use in a function such as
% DIGIPLANE.  This transformation is just an inverse of a scaling and
% origin shift.
newT=inv(S - [0 0 B(3); 0 0 B(1); 0 0 0]);

% Set things up for the image transformation.
newim = zeros(nrows,ncols);
[xi,yi] = meshgrid(1:ncols,1:nrows);    % All possible xy coords in the image.

% Transform these xy coords to determine where to interpolate values
% from. Note we have to work relative to x=B(3) and y=B(1).
sxy = homoTrans(Tinv, [xi(:)'+B(3) ; yi(:)'+B(1) ; ones(1,ncols*nrows)]);
xi = reshape(sxy(1,:),nrows,ncols);
yi = reshape(sxy(2,:),nrows,ncols);

[x,y] = meshgrid(1:cols,1:rows);
x = x+region(3)-1; % Offset x and y relative to region origin.
y = y+region(1)-1;
newim = interp2(x,y,double(im),xi,yi); % Interpolate values from source image.

% Plot bounding region
%P = [region(3) region(4) region(4) region(3)
%     region(1) region(1) region(2) region(2)
%      1    1    1    1   ];
%B = round(homoTrans(T,P));
%Bx = B(1,:);
%By = B(2,:);
%Bx = Bx-min(Bx); Bx(5)=Bx(1);
%By = By-min(By); By(5)=By(1);
%show(newim,2), axis xy
%line(Bx,By,'Color',[1 0 0],'LineWidth',2);
% end plot bounding region

end

%---------------------------------------------------------------------
%
% Internal function to find where the corners of a region, R
% defined by [minrow maxrow mincol maxcol] are transformed to
% by transform T and returns the bounds, B in the form
% [minrow maxrow mincol maxcol]

function B = bounds(T, R)

P = [R(3) R(4) R(4) R(3)      % homogeneous coords of region corners
R(1) R(1) R(2) R(2)
1    1    1    1   ];

PT = round(homoTrans(T,P));

B = [min(PT(2,:)) max(PT(2,:)) min(PT(1,:)) max(PT(1,:))];
%      minrow          maxrow      mincol       maxcol

end

%%

% HOMOTRANS - homogeneous transformation of points
%
% Function to perform a transformation on homogeneous points/lines
% The resulting points are normalised to have a homogeneous scale of 1
%
% Usage:
%           t = homoTrans(P,v);
%
% Arguments:
%           P  - 3 x 3 or 4 x 4 transformation matrix
%           v  - 3 x n or 4 x n matrix of points/lines

%  Peter Kovesi
%  School of Computer Science & Software Engineering
%  The University of Western Australia
%  pk @ csse uwa edu au
%  http://www.csse.uwa.edu.au/~pk
%
%  April 2000
%  September 2007

function t = homoTrans(P,v)

[dim,npts] = size(v);

if ~all(size(P)==dim)
error('Transformation matrix and point dimensions do not match');
end

t = P*v;  % Transform

for r = 1:dim-1     %  Now normalise
t(r,:) = t(r,:)./t(end,:);
end

t(end,:) = ones(1,npts);

end

%%

function [ BB ] = calcbb( H, im )
%CALCBB Compute bounding box

[imwidth imheight tempx] = size(im);

% Bounding box calculation

xmax=0;
ymax=0;
xmin=0;
ymin=0;

% Transform Point(1,1)
oldPoint = [1; 1; 1];
newPoint= H* oldPoint;
newPoint= newPoint/newPoint(3);

if (newPoint(1) < xmin)
xmin = newPoint(1);
end
if (newPoint(2) < ymin)
ymin = newPoint(2);
end
if (newPoint(1) > xmax)
xmax = newPoint(1);
end
if (newPoint(2) > ymax)
ymax = newPoint(2);
end

% Transform Point(1,imheight)
oldPoint = [1; imheight; 1];
newPoint= H* oldPoint;
newPoint= newPoint/newPoint(3);

if (newPoint(1) < xmin)
xmin = newPoint(1);
end
if (newPoint(2) < ymin)
ymin = newPoint(2);
end
if (newPoint(1) > xmax)
xmax = newPoint(1);
end
if (newPoint(2) > ymax)
ymax = newPoint(2);
end

% Transform Point(imwidth,1)
oldPoint = [imwidth; 1; 1];
newPoint= H* oldPoint;
newPoint= newPoint/newPoint(3);

if (newPoint(1) < xmin)
xmin = newPoint(1);
end
if (newPoint(2) < ymin)
ymin = newPoint(2);
end
if (newPoint(1) > xmax)
xmax = newPoint(1);
end
if (newPoint(2) > ymax)
ymax = newPoint(2);
end

% Transform Point(imwidth,imheight)
oldPoint = [imwidth; imheight; 1];
newPoint= H* oldPoint;
newPoint= newPoint/newPoint(3);

if (newPoint(1) < xmin)
xmin = newPoint(1);
end
if (newPoint(2) < ymin)
ymin = newPoint(2);
end
if (newPoint(1) > xmax)
xmax = newPoint(1);
end
if (newPoint(2) > ymax)
ymax = newPoint(2);
end

BB = [xmin, xmax; ymin, ymax];

% End of BB calculation

end

```