-
Notifications
You must be signed in to change notification settings - Fork 0
/
register_rotation.m
122 lines (118 loc) · 3.99 KB
/
register_rotation.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
function [final_R,im2] = register_rotation( image1, image2, thresh, buff );
% S = register_rotation( image1, image2, thresh, buff );
%
% A. Schaum Registration Algorithm using Prewitt Operator.
% Uses Peleg's iterative estimation technique to treat large shift values
% ->Analytic Methods of Image Registration Displacement Estimation and
% Resampling
% See paper high resolution image reconstruction from a sequence
% of rotated and translated frames and its application to an infrared
% imaging system by russell hardie
% image1 Reference image
% image2 Unknown shift image
% thresh The threshold for the largest shift to accept, suggested=.1
% buff The border buffer to ignore due to border effects and
% shift effects, suggested=max expected shift
%
% final_R=[sx,sy,theta] Output shift and rotation vector
% sx horizontal shift
% sy vertical shift
% theta rotation
%
% Author: Jeff Jenkins
% Based on register_shift by Russell Hardie et. al
% June 1996, modified 12/12/2017
% size(image1)
% size(image2r)
%---------------------------------%
% Estimate the discrete gradients %
%---------------------------------%
image1 =double(image1);
image2 =double(image2);
xkernel=(1/6)*[1 0 -1
1 0 -1
1 0 -1];
ykernel=(1/6)*[1 1 1
0 0 0
-1 -1 -1];
gx=conv2(image1,xkernel,'same');
gy=conv2(image1,ykernel,'same');
%-----------------------------------------------%
% Cut out center because of later shift effects %
%-----------------------------------------------%
[fully,fullx]=size(gx); % old size
gx=gx(buff:fully-buff+1,buff:fullx-buff+1);
gy=gy(buff:fully-buff+1,buff:fullx-buff+1);
im1=image1(buff:fully-buff+1,buff:fullx-buff+1);
[ydim,xdim]=size(im1); % new smaller size
%-----------------------%
% Generate the M matrix %
%-----------------------%
ssx=buff:fully-buff+1; ssy=buff:fullx-buff+1;
% We can do this because we are performing rotation in pixel spacing of 1
T1 = 1; T2 = 1;
hb=floor(max(size(ssy(:)))/2);
ha=floor(max(size(ssx(:)))/2);
for b=1:max(size(ssy(:)))
for a=1:max(size(ssx(:)))
ghaty(a,b) = (b-hb)*T1*gy(a,b);
ghatx(a,b) = (a-ha)*T2*gx(a,b);
end
end
g_bar = ghaty - ghatx;
gbar_gx = sum(sum(g_bar.*gx));
gbar_gy = sum(sum(g_bar.*gy));
gx_2 = sum(sum(gx.^2));
gy_2 = sum(sum(gy.^2));
gbar_2 = sum(sum(g_bar.^2));
cross=sum(sum(gx.*gy));
M=[ gx_2, cross, gbar_gx;
cross, gy_2, gbar_gy;
gbar_gx, gbar_gy, gbar_2];
%---------------------------%
% Initialize loop constants %
%---------------------------%
count=0; % Number of times through loop
stop=0; % Loop stop flag
im2=(image2); % start with the ``warped'' image being the full second image
Rn=zeros(3,1); % set initial shift estimates to zero
Rold=Rn;
ykprev=zeros(size(im1));
%---------------------------%
% Begin iterative estimates %
%---------------------------%
while stop~=1
count=count+1;
% calculate the difference image over the interior region
yk=double(im2(buff:fully-buff+1,buff:fullx-buff+1)-im1);
figure(9),
imagesc(yk),colormap 'gray';
act_yk=abs(double(ykprev-yk));
% Generate the V matrix
V=[ sum(sum( yk.*gx ));
sum(sum( yk.*gy ));
sum(sum( yk.*g_bar ));];
local = M\V;
Rn = Rn + local;
[sy,sx]=size(im2);
[W,Z]=meshgrid( [1:sx],[1:sy] );
nearest_x = (Rn(1,1));
nearest_y = (Rn(2,1));
XI=[1:fullx];
YI=[1:fully]';
check_val = sqrt(sum((Rn-Rold).^2)) / sqrt(sum(Rold.^2));
% See if its time to stop or warp and continue
if count > 750 | (check_val <= thresh) | sum(act_yk(:))<eps %| (nnz(isnan(check_val)) >0)
stop=1;
sum(act_yk(:));
im2(find(isnan(im2))) = 0;
% im2 = interp2( W,Z,im2, XI,YI, 'bilinear' );
final_R = Rn;
else % sub-pixelly rotate and then shift image2 according to latest estimate
% Get the values at the new coordinates
a=imrotate(image2, -Rn(3,1)*1.5, 'bilinear', 'crop');
im2=imtranslate(a, [nearest_x, nearest_y]);
im2(find(isnan(im2))) = 0;
ykprev = yk;
end
end