--- /dev/null
+
+# 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)
+
--- /dev/null
+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
projects / sie_pcm / commitdiff
