www.gusucode.com > superresolution_v_2源码程序 > superresolution_v_2源码程序/superresolution_v_2.0_超分辨率图像处理_matlab源码_POCS/superresolution_v_2.0_超分辨率图像处理_matlab源码_POCS/superresolution_v_2.0/application/lucchese.m
function [delta_est, phi_est] = lucchese(s,M)
% LUCCHESE - estimate shift and rotation parameters using Lucchese and Cortelazzo algorithm
% [delta_est, phi_est] = lucchese(s,M)
% horizontal and vertical shifts DELTA_EST and rotations PHI_EST are
% estimated from the input images S (S{1},etc.). For the shift
% estimation, the Fourier transform images are interpolated by a factor
% M to increase precision
%% -----------------------------------------------------------------------
% SUPERRESOLUTION - Graphical User Interface for Super-Resolution Imaging
% Copyright (C) 2005-2007 Laboratory of Audiovisual Communications (LCAV),
% Ecole Polytechnique Federale de Lausanne (EPFL),
% CH-1015 Lausanne, Switzerland
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version. This software is distributed in the hope that
% it will be useful, but without any warranty; without even the implied
% warranty of merchantability or fitness for a particular purpose.
% See the GNU General Public License for more details
% (enclosed in the file GPL).
%
% Latest modifications: January 12, 2006, by Patrick Vandewalle
nr=length(s);
S = size(s{1});
if (nargin==1)
M = 4; % magnification factor to increase precision of shift estimation
end
% ROTATION ESTIMATION
z = zeros(2*S);
for i=1:nr
s2{i} = imresize(s{i},0.5,'bicubic');
end
% init stage 1 - coarse estimation of phi
Rmax = 80;
interpolation = 4;
histsize = 2;
IMREF = fftshift(fft2(s2{1},2*S(1),2*S(2)));
IMREF_SHIFT = fft2(s{1});
for i=2:nr
% stage 1 - coarse estimation of phi
s2{i}=s2{i}(:,end:-1:1);
IM = fftshift(fft2(s2{i},2*S(1),2*S(2)));
Delta = abs(IMREF.^2)/abs(IMREF(S(1)+1,S(2)+1)^2) ...
- abs(IM.^2)/abs(IM(S(1)+1,S(2)+1)^2); % S is even
Delta_part = Delta(S(1)-Rmax:S(1)+Rmax,S(2)-Rmax:S(2)+Rmax);
Delta_part = imresize(Delta_part,interpolation,'bilinear');
ds = size(Delta_part);
locus = contourc(Delta_part,[0 0]);
ind = find(locus(1,:)>0); % exclude the 'index' columns with 0 and the length of the contour
locus = locus(:,ind)-[326.5;326.5]*ones(1,length(ind));
[th,r] = cart2pol(locus(1,:),locus(2,:));
th2a = th(find(r<Rmax & th>-pi/4 & th<pi/4)); % only keep points for which r<Rmax
th2b = th(find(r<Rmax & th>pi/4 & th<3*pi/4)); % only keep points for which r<Rmax
H1 = hist(th2a*180/pi,[-45:1/2/histsize:45]);
H2 = hist(th2b*180/pi,[45:1/2/histsize:135]);
H = H1+H2;
H = filter(gausswin(11),1,H);
H = H(6:end); % need to cut off negative part of gaussian from filtering
% gaussian was not centered
[m,ind] = max(H);
phi_est1(i) = (ind-histsize*90-1)/histsize;
%figure; plot(H);
% stage 2 - refinement of stage 1
gamma = 4.5; % 1/2 aperture in degrees
histsize_ = 100;
th2a_ = th(find(r<Rmax & th>(phi_est1(i)/2-gamma)*pi/180 ...
& th<(phi_est1(i)/2+gamma)*pi/180)); % only keep points for which r<Rmax
th2b_ = th(find(r<Rmax & th>(phi_est1(i)/2-gamma)*pi/180+pi/2 ...
& th<(phi_est1(i)/2+gamma)*pi/180+pi/2)); % only keep points for which r<Rmax
H1_ = hist(th2a_*180/pi,[-gamma:1/2/histsize_:gamma]);
H2_ = hist(th2b_*180/pi,[-gamma:1/2/histsize_:gamma]);
H_ = H1_+H2_;
H_ = filter(gausswin(11),1,H_);
Hnorm = H_/sum(H_);
theta = [-gamma:1/2/histsize_:gamma];
phi_est2(i) = phi_est1(i)+sum(theta.*Hnorm);
sigma2 = sum(theta.*theta.*Hnorm);
% stage 3
omega = find(r<Rmax & ...
((th>(phi_est2(i)/2-sigma2)*pi/180 & th<(phi_est2(i)/2+sigma2)*pi/180) ...
| (th>(phi_est2(i)/2-sigma2)*pi/180 + pi & th<(phi_est2(i)/2+sigma2)*pi/180 + pi) ...
| (th>(phi_est2(i)/2-sigma2)*pi/180 - pi & th<(phi_est2(i)/2+sigma2)*pi/180 - pi)));
omega_ = find(r<Rmax & ...
((th>(phi_est2(i)/2-sigma2)*pi/180+pi/2 & th<(phi_est1(i)/2+sigma2)*pi/180+pi/2) ...
| (th>(phi_est2(i)/2-sigma2)*pi/180+3*pi/2 & th<(phi_est1(i)/2+sigma2)*pi/180+3*pi/2) ...
| (th>(phi_est2(i)/2-sigma2)*pi/180-pi/2 & th<(phi_est1(i)/2+sigma2)*pi/180-pi/2)));
kx = locus(1,omega);
ky = locus(2,omega);
kx_ = locus(1,omega_);
ky_ = locus(2,omega_);
%figure; plot(locus(1,:),locus(2,:),'.'); hold; plot(kx,ky,'ro'); plot(kx_,ky_,'rx');title('s3')
phi_est3(i) = atan(2*(sum(kx.*ky)/length(kx) - sum(kx_.*ky_)/length(kx_))/...
(sum(kx.*kx-ky.*ky)/length(kx) - sum(kx_.*kx_-ky_.*ky_)/length(kx_)))*180/pi;
phi_est = phi_est3;
% rotation compensation, required to estimate shifts
s2{i} = imrotate(s{i},-phi_est(i),'bilinear','crop');
% SHIFT ESTIMATION
IM = fft2(s2{i});
psi = IM./IMREF_SHIFT;
PSI = fft2(psi,M*S(1),M*S(2));
[m,ind] = max(PSI);
[mm,iind] = max(m);
delta_est(i,1) = (ind(iind)-1)/M;
delta_est(i,2) = (iind-1)/M;
if delta_est(i,1)>S(1)/2
delta_est(i,1) = delta_est(i,1)-S(1);
end
if delta_est(i,2)>S(2)/2
delta_est(i,2) = delta_est(i,2)-S(2);
end
end