sigprocexperts.com Git - sie_pcm/blob - siepcm.m

   1 function [x,y,verbose] = siepcm(a,b,verbose,threshold,radius)
   2 % SIEPCM - estimation of sub-pixel shifts between images using
   3 %          SVD of phase correlation
   4 %
   5 % Inputs:
   6 %        a, b - input image data, in Fourier domain with DC centered
   7 %     verbose - if non-zero, display debug plots and images
   8 %    thrshold - a percentage of maximum value, used to set initial mask
   9 %      radius - the width of data to use for LMS fit of phase slopes (<=0.25)
  10
  11 if nargin < 3, verbose = 0;     end;    % plot intermediate images
  12 if nargin < 4, threshold = 0.0015; end;    % a percentage of max to set mask
  13 if nargin < 5, radius = 0.25;   end;    % size of post-mask radius (<= 0.25)
  14
  15
  16 [m,n] = size(a);
  17
  18 s = abs( a .* conj(b) ); s( s == 0 ) = 1;
  19 c = (a .* conj(b)) ./ s; % abs( a .* conj(a) + 1e-10 );
  20
  21 s = ones(size(c));
  22 s( abs(c) ~= 0 ) = abs(c( abs(c) ~= 0 ));
  23
  24
  25 %%%% apply pre-masking here.
  26
  27 if exist('radius'), r=radius(1); else, r = 0.3;  end;
  28
  29 mask = zeros(size(c));
  30 t = max(abs(a(:))).*threshold;
  31
  32 mask( abs(a) > t ) = 1;
  33 mask = medfilt2( mask, [10 10] );       % 2D median filter, from image
  34                                         % processing toolbox
  35 if (sum(mask(:)) == 0);
  36   error(sprintf([' mask in empty.  Need to lower the threshold (currently=' ...
  37   ' %f)'], threshold));
  38 end;
  39
  40 [u d v]=svds(mask.*c,1);
  41 % [u d v]=svds( c,1);
  42
  43 r2 = min(r,0.5);
  44 % solve an LMS fit to the angles of the rank one matrix components...
  45 if n > 3,
  46   n2 = length(v);
  47   %% apply post-masking:
  48   t_n = ceil( (0.5-r2)*length(v) ):floor( (0.5+r2)*length(v) );
  49
  50   tmp = unwrap( angle(v(t_n)) ); A = [ (t_n-1)' ones(length(t_n),1) ];
  51   sys_y = A \ tmp;
  52   y = sys_y(1) * n / (2*pi);
  53   fullphz1 = unwrap(angle(v));
  54   fullphz1 = fullphz1 - fullphz1(floor(n2/2)) + ([floor(n2/2)-1 1]*sys_y);
  55   pltcmd1 = '1:length(v), fullphz1, t_n, [ unwrap( angle(v(t_n)) ) ( sys_y(1) * (t_n-1)'' + sys_y(2) ) ]';
  56 else,
  57   pltcmd1 = ' 0,0 '; sys_y = [ 0 0]; y = 0;
  58 end;
  59
  60 if m > 3,
  61   m2 = length(u);
  62
  63   %% apply post-masking:
  64   t_m = ceil( (0.5-r2)*length(u) ):floor( (0.5+r2)*length(u) );
  65
  66   tmp = unwrap( angle(u(t_m)) ); A = [ (t_m-1)' ones(length(t_m),1) ];
  67   sys_x = A \ tmp;
  68   x = sys_x(1) * m / (2*pi);
  69   fullphz2 = unwrap(angle(u));
  70   fullphz2 = fullphz2 - fullphz2(floor(m2/2)) + ([floor(m2/2)-1 1]*sys_x);
  71   pltcmd2 = '1:length(u), fullphz2, t_m, [ unwrap( angle(u(t_m)) ) ( sys_x(1) * (t_m-1)'' + sys_x(2) ) ]';
  72 else,
  73   pltcmd2 = ' 0, 0 '; sys_x = [ 0 0]; x = 0;
  74 end;
  75
  76 % if verbose mode is on, show the estimation quality images
  77 if verbose,
  78   figure(1);
  79   disp([ 'plot( ' pltcmd1 ', ''o-'', ' pltcmd2 ', ''*-'');' ]);
  80   eval([ 'plot( ' pltcmd1 ', ''o-'', ' pltcmd2 ', ''*-'');' ]);
  81   title(' plot of least squares linear fit ');
  82
  83
  84   Q = exp( j*2*pi* x * [0:m-1]'/m ) * exp( -j*2*pi* y * [0:n-1]/n );
  85   q = u*v';
  86
  87   % phase align each of these beasts so that the abs( c - Q ) looks good.
  88   mc = floor(m/2); nc = floor(n/2);
  89   Q = Q * exp(angle(Q(mc,nc)'*c(mc,nc))*j);
  90   q = q * exp(angle(q(mc,nc)'*c(mc,nc))*j);
  91   figure(2);
  92   imagesc( angle([ c mask.*c; u*v' Q ]) ); colormap('hsv');
  93   title('[ Q | mask \circ Q ; u v^H | est of (u v^H) ]');
  94
  95   disp([ '    x_0: ' num2str(x) '       y_0: ' num2str(y) ]);
  96   keyboard;
  97 end;
  98
  99 % Copyright 2002-2005 W. Scott Hoge (wsh032580 at proton dot me)
 100 %
 101 % Licensed under the terms of the MIT License
 102 % (https://opensource.org/licenses/MIT)