SIE-PCM demo from 2003

diff --git a/ReadMe.md b/ReadMe.md
new file mode 100644 (file)
index 0000000..5396c61
--- /dev/null
+++ b/ReadMe.md
@@ -0,0 +1,99 @@
+
+# Subspace Identification Extension to the Phase Correlation Method (SIE-PCM) demo
+
+This short bit of  Matlab code was the original implementation of the method described in
+``A subspace identification extension to the phase correlation method,''
+IEEE Trans. Medical Imaging, 22(2):277-280, Feb. 2003.
+http://dx.doi.org/10.1109/TMI.2002.808359
+
+Images from the original 5 test frame data are included as well.
+
+# File List
+
+    run_example.m - start here
+    siepcm.m      - a function to estimate the row and column translation between two images
+
+The images are from 5 different MRI acquisition of a grapefruit, with
+the grapefruit at different locations within the field-of-view.  In
+the table below, the coordinates of one corner of the field of view
+are given.
+
+     Image File | Source file |  dim1 loc    | dim2 loc
+    ------------|-------------|--------------|------------
+     img1.png   | P64512.7    | y = 8.30 cm  | z= 2.00 cm
+     img2.png   | P00000.7    | y = 8.15 cm  | z= 2.25 cm
+     img3.png   | P00512.7    | y = 8.00 cm  | z= 2.50 cm
+     img4.png   | P01024.7    | y = 7.85 cm  | z= 2.27 cm
+     img5.png   | P01536.7    | y = 7.85 cm  | z= 2.75 cm
+
+To determine the expected translational shift between the images, use the
+difference between the coordinates divided by the distance per pixel.
+
+                     dy = (y2 - y1) / n  where n=0.0625
+
+The original data was acquired on a GE Signa Lx 1.5T MR Scanner using
+the following parameters:
+ TR: 500  TE: 15   protocol: FSE  coil: HEAD
+
+----------------------------------------------------------------------
+
+FOV length : Image Size :: 160mm : 256 pixels; => 0.0625 mm/pixel
+
+ compare images 1 and 2.  S: -2.5 mm, A: 1.5 mm
+       expect z =  -4.0000,  y =  +2.4000
+
+ compare images 1 and 3.  S: -5 mm, A: 3 mm
+       expect z =  -8.0000,  y =  +4.8000
+
+ compare images 1 and 4.  S: -2.7 mm, A: 4.5 mm
+       expect z =  -4.3200,  y =  +7.2000
+
+ compare images 1 and 5.  S: -7.5 mm, A: 4.5 mm
+       expect z = -12.0000,  y =  +7.2000
+
+ compare images 2 and 3.  S: -2.5 mm, A: 1.5 mm
+       expect z =  -4.0000,  y =  +2.4000
+
+ compare images 2 and 4.  S: -0.2 mm, A: 3 mm
+       expect z =  -0.3200,  y =  +4.8000
+
+ compare images 2 and 5.  S: -5 mm, A: 3 mm
+       expect z =  -8.0000,  y =  +4.8000
+
+ compare images 3 and 4.  S: 2.3 mm, A: 1.5 mm
+       expect z =  +3.6800,  y =  +2.4000
+
+ compare images 3 and 5.  S: -2.5 mm, A: 1.5 mm
+       expect z =  -4.0000,  y =  +2.4000
+
+ compare images 4 and 5.  S: -4.8 mm, A: 0 mm
+       expect z =  -7.6800,  y =  +0.0000
+
+----------------------------------------------------------------------
+
+Copyright 2003, W. Scott Hoge (wsh032580 at proton dot me)
+
+This images and text in this zip file are licensed under the Creative
+Commons Attribution License. To view a copy of this license, visit
+http://creativecommons.org/licenses/by/1.0/ or send a letter to Creative
+Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
+
+You are free:
+
+    * to copy, distribute, display, and perform the work
+    * to make derivative works
+    * to make commercial use of the work
+
+Under the following conditions:
+by Attribution. You must give the original author credit.
+
+    * For any reuse or distribution, you must make clear to others the
+      license terms of this work.
+    * Any of these conditions can be waived if you get permission from the
+      author. 
+
+Your fair use and other rights are in no way affected by the above.
+
+This is a human-readable summary of the Legal Code 
+(the full license: http://creativecommons.org/licenses/by/1.0/legalcode)
+

diff --git a/img1.png b/img1.png
new file mode 100755 (executable)
index 0000000..21fac8c
Binary files /dev/null and b/img1.png differ

diff --git a/img2.png b/img2.png
new file mode 100755 (executable)
index 0000000..d0bbe54
Binary files /dev/null and b/img2.png differ

diff --git a/img3.png b/img3.png
new file mode 100755 (executable)
index 0000000..3df03e2
Binary files /dev/null and b/img3.png differ

diff --git a/img4.png b/img4.png
new file mode 100755 (executable)
index 0000000..d49292e
Binary files /dev/null and b/img4.png differ

diff --git a/img5.png b/img5.png
new file mode 100755 (executable)
index 0000000..bcd9dfa
Binary files /dev/null and b/img5.png differ

diff --git a/run_example.m b/run_example.m
new file mode 100755 (executable)
index 0000000..d70ebc8
--- /dev/null
+++ b/run_example.m
@@ -0,0 +1,11 @@
+% load images
+a = imread('img1.png');
+b = imread('img2.png');
+
+% transform image data to Fourier domain
+A = fftshift(fft2(double(a)));
+B = fftshift(fft2(double(b)));
+
+% call pcm estimation code
+[x,y]=siepcm(A,B,1);                    % 3rd arg == 1: display debug images
+

diff --git a/siepcm.m b/siepcm.m
new file mode 100755 (executable)
index 0000000..7a34228
--- /dev/null
+++ b/siepcm.m
@@ -0,0 +1,102 @@
+function [x,y,verbose] = siepcm(a,b,verbose,threshold,radius)
+% SIEPCM - estimation of sub-pixel shifts between images using 
+%          SVD of phase correlation 
+%
+% Inputs: 
+%        a, b - input image data, in Fourier domain with DC centered
+%     verbose - if non-zero, display debug plots and images
+%    thrshold - a percentage of maximum value, used to set initial mask
+%      radius - the width of data to use for LMS fit of phase slopes (<=0.25)
+
+if nargin < 3, verbose = 0;     end;    % plot intermediate images
+if nargin < 4, threshold = 0.0015; end;    % a percentage of max to set mask
+if nargin < 5, radius = 0.25;   end;    % size of post-mask radius (<= 0.25)
+
+
+[m,n] = size(a);
+
+s = abs( a .* conj(b) ); s( s == 0 ) = 1;
+c = (a .* conj(b)) ./ s; % abs( a .* conj(a) + 1e-10 );
+
+s = ones(size(c));
+s( abs(c) ~= 0 ) = abs(c( abs(c) ~= 0 ));
+
+
+%%%% apply pre-masking here.
+
+if exist('radius'), r=radius(1); else, r = 0.3;  end;
+
+mask = zeros(size(c));
+t = max(abs(a(:))).*threshold;
+
+mask( abs(a) > t ) = 1;
+mask = medfilt2( mask, [10 10] );       % 2D median filter, from image
+                                        % processing toolbox
+if (sum(mask(:)) == 0);
+  error(sprintf([' mask in empty.  Need to lower the threshold (currently=' ...
+  ' %f)'], threshold)); 
+end;
+                                        
+[u d v]=svds(mask.*c,1);
+% [u d v]=svds( c,1);
+
+r2 = min(r,0.5);
+% solve an LMS fit to the angles of the rank one matrix components...
+if n > 3,
+  n2 = length(v);
+  %% apply post-masking:
+  t_n = ceil( (0.5-r2)*length(v) ):floor( (0.5+r2)*length(v) );
+
+  tmp = unwrap( angle(v(t_n)) ); A = [ (t_n-1)' ones(length(t_n),1) ];
+  sys_y = A \ tmp;
+  y = sys_y(1) * n / (2*pi);
+  fullphz1 = unwrap(angle(v));
+  fullphz1 = fullphz1 - fullphz1(floor(n2/2)) + ([floor(n2/2)-1 1]*sys_y);
+  pltcmd1 = '1:length(v), fullphz1, t_n, [ unwrap( angle(v(t_n)) ) ( sys_y(1) * (t_n-1)'' + sys_y(2) ) ]';
+else,
+  pltcmd1 = ' 0,0 '; sys_y = [ 0 0]; y = 0;
+end;
+
+if m > 3,
+  m2 = length(u);
+
+  %% apply post-masking:
+  t_m = ceil( (0.5-r2)*length(u) ):floor( (0.5+r2)*length(u) );
+  
+  tmp = unwrap( angle(u(t_m)) ); A = [ (t_m-1)' ones(length(t_m),1) ];
+  sys_x = A \ tmp;
+  x = sys_x(1) * m / (2*pi);
+  fullphz2 = unwrap(angle(u));
+  fullphz2 = fullphz2 - fullphz2(floor(m2/2)) + ([floor(m2/2)-1 1]*sys_x);
+  pltcmd2 = '1:length(u), fullphz2, t_m, [ unwrap( angle(u(t_m)) ) ( sys_x(1) * (t_m-1)'' + sys_x(2) ) ]';
+else,
+  pltcmd2 = ' 0, 0 '; sys_x = [ 0 0]; x = 0;
+end;
+
+% if verbose mode is on, show the estimation quality images
+if verbose,  
+  figure(1);
+  disp([ 'plot( ' pltcmd1 ', ''o-'', ' pltcmd2 ', ''*-'');' ]);
+  eval([ 'plot( ' pltcmd1 ', ''o-'', ' pltcmd2 ', ''*-'');' ]);
+  title(' plot of least squares linear fit ');
+  
+  
+  Q = exp( j*2*pi* x * [0:m-1]'/m ) * exp( -j*2*pi* y * [0:n-1]/n );
+  q = u*v'; 
+  
+  % phase align each of these beasts so that the abs( c - Q ) looks good.
+  mc = floor(m/2); nc = floor(n/2);
+  Q = Q * exp(angle(Q(mc,nc)'*c(mc,nc))*j);
+  q = q * exp(angle(q(mc,nc)'*c(mc,nc))*j);
+  figure(2);
+  imagesc( angle([ c mask.*c; u*v' Q ]) ); colormap('hsv');
+  title('[ Q | mask \circ Q ; u v^H | est of (u v^H) ]');
+  
+  disp([ '    x_0: ' num2str(x) '       y_0: ' num2str(y) ]);
+  keyboard;
+end;
+
+% Copyright 2002-2005 W. Scott Hoge (wsh032580 at proton dot me) 
+%
+% Licensed under the terms of the MIT License
+% (https://opensource.org/licenses/MIT)
\ No newline at end of file