From: WSHoge <>
Date: Sat, 28 Jan 2023 16:33:26 +0000 (-0500)
Subject: DPG demo version from 2019-08-20
X-Git-Url: https://sigprocexperts.com/demo_code/dpg/commitdiff_plain

DPG demo version from 2019-08-20
---

82df107ae0fe7e8430852935a91203370d249136
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..0b312ce
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023, W. Scott Hoge
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/ReadMe.md b/ReadMe.md
new file mode 100644
index 0000000..34f6fc5
--- /dev/null
+++ b/ReadMe.md
@@ -0,0 +1,49 @@
+
+# Dual-Polarity GRAPPA (DPG) demo
+
+This collection of Matlab code and data presents one implementation of the method described in
+"Dual-Polarity GRAPPA for simultaneous reconstruction and ghost correction of EPI data."
+_Magn Reson Med_. 76(1):32-44, Jul 2016.
+http://dx.doi.org/10.1002/mrm.25839
+
+
+## Data source
+
+EPI gradient-echo BOLD data
+
+- 3T Siemens Trio
+- 32 Channel Head coil
+- TR 3sec
+- TE 30msec
+- Echo Spacing 770usec
+- 2.0mm isotropic resolution
+
+- 96x60   matrix size
+- 192x120 (mm) FOV
+- 32 slices, 2mm thick
+
+
+collected on 2015-05-18, data file MID902
+
+## Demo Filelist:
+
+    data/dpg_acs1.nd     ACS data, with the first readout acquired using RO+
+    data/dpg_acs2.nd     ACS data, with the first readout acquired using RO-
+    data/dpg_C.nd        The PLACE ACS image, formed from a combination of dpg_acs1 and dpg_acs2
+    
+    data/dpg_k.nd        one frame of R=2 EPI data
+    data/dpg_k_pcref.nd  navigator data for dpg_k, to perform traditional Nyquist ghost correction
+    data/y.nd            NGC coefficients, derived from dpg_k_pcref and used with phzapply.m
+    data/y_acs2.nd       NGC coefficients, derived from A and B and used with phzapply.m
+    
+    dpg_cal.m       the DPG calibration code
+    dpg_recon.m     the DPG reconstruction code
+    phzshift.m	    calculates the coordinate shift between two k-space data sets using SIEPCM
+    phzapply.m      a script to apply standard Nyquist ghost correction coefficients
+    fif.m           a 2D FFT macro script
+    ifi.m           a 2D FFT macro script
+    readnd.m        script to read the .nd data files
+    
+    runthis.m       the main demo script
+    runthis2.m      an alternate demo script (w/o GESTE/PLACE acs data)
+
diff --git a/data/dpg_C.nd b/data/dpg_C.nd
new file mode 100644
index 0000000..8f964d7
Binary files /dev/null and b/data/dpg_C.nd differ
diff --git a/data/dpg_acs1.nd b/data/dpg_acs1.nd
new file mode 100644
index 0000000..4260119
Binary files /dev/null and b/data/dpg_acs1.nd differ
diff --git a/data/dpg_acs2.nd b/data/dpg_acs2.nd
new file mode 100644
index 0000000..53bf7ec
Binary files /dev/null and b/data/dpg_acs2.nd differ
diff --git a/data/dpg_k.nd b/data/dpg_k.nd
new file mode 100644
index 0000000..ef32bd1
Binary files /dev/null and b/data/dpg_k.nd differ
diff --git a/data/dpg_k_pcref.nd b/data/dpg_k_pcref.nd
new file mode 100644
index 0000000..2169b31
Binary files /dev/null and b/data/dpg_k_pcref.nd differ
diff --git a/data/y.nd b/data/y.nd
new file mode 100644
index 0000000..cc357ea
Binary files /dev/null and b/data/y.nd differ
diff --git a/data/y_acs2.nd b/data/y_acs2.nd
new file mode 100644
index 0000000..1ea5158
Binary files /dev/null and b/data/y_acs2.nd differ
diff --git a/dpg_cal.m b/dpg_cal.m
new file mode 100644
index 0000000..7d3597b
--- /dev/null
+++ b/dpg_cal.m
@@ -0,0 +1,479 @@
+function [F,I,N] = dpg_cal( varargin );
+
+% [F,I,N ] = dpg_cal( Wk, Sk, Y, [opt args]);
+% 
+% opt args:
+%  'v', {0,1}              - verbosity off/on
+%  'kernel', string        - grappa kernel size, x by y. e.g. '2x5', '4x1'
+%  'N', (xy)-by-num_patterns-L matrix   - grappa reconstruction parameters
+%  'acs', array of indices - index locations to use for ACS lines
+% 
+% This implementation is derived from 
+% "Generalized autocalibrating partially parallel acquisitions"
+%   by M. A. Griswold, P. M. Jakob, et. al..
+%  Mag. Reson. Med, 47(6):1202-1210, Jun 2002
+%
+
+%
+% Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+%
+% GRAPPA: data from all coils used to patch data in one coil
+%
+%  coil 1     o * o * a * o *
+%  coil 2     o * o * a * o *
+%  coil 3     o * o * a * o *
+%  coil 4     o * o * x * o *
+%
+%  *: acquired line, o: unacquired line, a: auto-callibration line
+%  x: reconstructed line
+%
+%  S_j (k_y - m\Delta k_y ) = \sum_{l=1}^L \sum_{b=0}^{N_b-1}
+%                                  n(j,b,l,m) S_l(k_y - b A \Delta k_y)
+%
+
+Wk  = varargin{1};
+kin = varargin{2};
+Y   = varargin{3};
+
+acs = []; N = []; verbose=0; kernel = '2x5'; dks = []; prnt = 0;
+eta = 1.0;
+chi = 0.0001;
+
+vararg = varargin(4:end);
+for cnt=1:2:length(vararg),
+  if strcmp( vararg{cnt}, 'v' ),
+    verbose = cell2mat(vararg(cnt+1));
+  elseif strcmp( vararg{cnt}, 'N' ),
+    N = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'dks' ),
+    dks = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'kernel' ),
+    kernel = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'acs' ),
+    acs = vararg{cnt+1};
+  elseif (strcmp( vararg{cnt}, 'report' ) || strcmp( vararg{cnt}, 'r' )),
+    prnt = vararg{cnt+1};
+  end;
+end;
+
+if length(size(Wk)) == 3, 
+  [Wm,Wn,ncoils] = size(Wk);
+else,
+  Wm = Wk(1); Wn = Wk(2); ncoils = Wk(3);
+end;
+
+
+%% zero-pad the k-space data
+Sk = zeros( Wm, Wn, ncoils );
+if iscell(kin)
+  dualks = 1;
+  s = kin{3};
+else
+  dualks = 0;
+  s = kin; clear kin;
+end;
+Sk(Y,:,:) = s;
+
+rptupd(prnt,['* ' mfilename ': kernel size = ' kernel '\n']);
+
+d = diff(Y(:));
+if ( isempty(acs) ),
+  prfl = max( sum(abs( Sk(:,:,:)),3).' );
+
+  if length( d == 1) == 1,
+    acs = find( prfl == max(prfl(Y(find(d==1)))) );
+  else
+    acs = Y([ min(find(d==1)); find( d == 1)+1 ]);
+  end;
+  rptupd(prnt,['* ' mfilename ': acs =']);
+  rptupd(prnt,sprintf(' %d', acs) );
+  rptupd(prnt,'\n');
+end;
+
+if ( isempty(dks) ),
+  %% determine which delta-k's to repair.
+  dks = [];
+  for cnt=max(d):-1:2,
+    if (sum( cnt == d ) > 0),
+      dks = [ cnt; dks ];
+    end;
+  end;
+end;
+
+% set the kernel spacings based on 'kernel' string
+nky = str2double(kernel(1));           
+nkx = str2double(kernel(3:end)); % should be odd
+
+indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
+
+pattern = {};
+
+if ~isempty(N),
+  for cnt=1:length(N)
+    pattern{cnt} = N(cnt).pattern;
+  end;
+end;
+
+if isempty(pattern)
+%% foreach dk spacing
+for cnt=1:length(dks),
+  
+  if ( dks(cnt) == 1 )
+    % for no acceleration, just specify the pattern
+    tmppat = repmat('*',[1 nky]);
+    for patcnt = 1:nky,
+      pattern{length(pattern)+1} = tmppat;
+      pattern{length(pattern)}(patcnt) = '<';
+    end;
+
+    for patcnt = 1:nky,
+      pattern{length(pattern)+1} = tmppat;
+      pattern{length(pattern)}(patcnt) = '>';
+    end;
+
+    continue;
+  end;
+
+  tmp = [ repmat([ '*' repmat('o',1,dks(cnt)-1) ],1,nky-1) '*' ];
+  z2 = strfind(tmp,'o');
+  if isempty(z2),                       
+    % pad with zeros, if the kernel is 1-by-N
+    tmp = [ 'o' tmp 'o' ];   
+    z2 = strfind(tmp,'o');
+  end;
+  [tmp2,indx] = sort( abs(z2 - z2( floor(median(1:length(z2))) )) );
+
+  olength = length(pattern);
+  if (dualks == 1)
+    % for dual-kernels, we want to genereate coefficients for the sampled
+    % lines as well.
+    npat = 2*length(tmp);
+  else
+    npat = length(z2);
+  end
+  for cnt2 = 1:npat, % (dks(cnt)-1), % 
+
+    if (dualks == 0 )
+      cnt3 = indx(cnt2);
+      pattern{length(pattern)+1} = tmp;
+      pattern{length(pattern)}(z2(cnt3)) = 'x';
+    else
+      cnt3 = mod( cnt2-1, length(tmp) ) + 1;
+      pattern{length(pattern)+1} = tmp;
+      if ( cnt2 <= npat/2 )
+        pattern{length(pattern)}( cnt3 ) = '<';
+      else
+        pattern{length(pattern)}( cnt3 ) = '>';
+      end;
+    end;
+  
+  end;
+
+end;
+end;
+
+
+% find the recon patterns we can use 
+plist = cell(length(pattern),1);
+
+for cnt = 1:length(pattern),
+  v = extrapolate_pattern(pattern{cnt},nky);     
+
+  % if cnt == 14, keyboard; end;
+
+  for cnt2 = 1:length(acs);
+    %% check if the GRAPPA pattern...
+    indx = acs(cnt2) + v;
+
+    %% can be accomodated by the sampling pattern 
+    tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
+    if ( length(tmp) == length(indx) ),
+      plist{cnt} =  [ plist{cnt} acs(cnt2) ];
+    end;
+  end;
+end;
+
+rptupd(prnt,['* ' mfilename ': sampling patterns in use = \n' ]);
+for cnt1=1:length(pattern);
+  rptupd(prnt,sprintf('  %s: %d acs lines =', pattern{cnt1}, length(plist{cnt1})) );
+  rptupd(prnt,sprintf(' %3d', plist{cnt1}) );
+  % for cnt2 = 1:length( plist{cnt1} ),
+  %   rptupd(prnt,('  %s', pattern{plist{cnt1}} );
+  %   % v = extrapolate_pattern(pattern{plist{cnt1}(cnt2)},nky);
+  %   % norm(vec(Sk( acs(cnt1)+v, :, : )));
+  % end;
+  rptupd(prnt,'\n');
+end;
+if (length(plist) == 0), 
+  error(['sorry, acquisition pattern is incompatible with available GRAPPA recon patterns']);
+end;
+
+st = clock; % tic;
+
+if verbose, rptupd(prnt,'\n'); end;
+
+mss = zeros( size(Sk,1), 1 );         % set the missing lines index vector 
+mss(Y) = 1;
+
+if isempty(N),
+  % N = zeros( length(indx3)*nky*size(s,3),length(pattern),size(s,3));
+
+  rptupd(prnt,'calc recon param\n');
+  for cnt = 1:length(pattern), 
+    v = extrapolate_pattern(pattern{cnt},nky);
+
+    %% and for each pattern ... 
+    %% A2 = [];     b2 = [];
+    
+    %%  A is data from the expanded input data matrix.
+    %%     size(A) = [ kxnum * nacs, ncoils * nky * nkx ]
+    %%
+    %% in the matlab (and c-code) structure is:
+    %%    for each element of indx3
+    %%                   coil 1   |    coil 2   |    coil 3   | ...
+    %%                -3 -1  1  3 | -3 -1  1  3 | -3 -1  1  3 | ...
+    %%                 o  o  o  o |  o  o  o  o |  o  o  o  o | ...
+    %%             ^
+    %%  ACS(1) xres|
+    %%             v
+    %%             ^
+    %%  ACS(2) xres|
+    %%             v
+    %%             ^
+    %%  ACS(3) xres|
+    %%             v
+    %% 
+    %%  the older version of matlab code sorted the columns by
+    %% 
+    %%     | c1 c2 c3 ... | c1 c2 c3 ... | c1 c2 c3 ...  | c1 c2 c3 ...
+    %%     | -3 -3 -3 ... | -1 -1 -1 ... |  1  1  1 ...  |  3  3  3 ...
+    %%     |  o  o  o ... |  o  o  o ... |  o  o  o ...  |  o  o  o ...
+    %% 
+    %%  the order is immaterial, as long at it is consistent...
+    %%  UPDATE (12 Jan 2011): the order does matter for GROG, so we have
+    %%                        now changed it here so that the c-code and matlab code match.
+    %%   
+    %%  The old way to switch between the two was:
+    %%    Accode(:, permute( reshape(1:size(Accode,2),YKS,Ncoils,XKS), [3 2 1]) ) == Amatlab
+    %%    Accode = Amatlab(:, permute( reshape(1:size(Accode,2),XKS,Ncoils,YKS), [3 2 1]))
+    A = zeros( length(plist{cnt})*size(s,2), length(indx3)*nky*size(s,3) );
+    
+    rptupd(prnt,sprintf('%3d  %s: ', cnt, pattern{cnt}));
+
+    % v = extrapolate_pattern(pattern{cnt});
+    N(length(N)+1).pattern = pattern{cnt};
+
+    n = zeros( nky*size(s,3)*length(indx3),size(s,3));
+
+    for cnt2 = 1:length(plist{cnt}),
+      %% and ACS line ...
+
+      %% if the sampling pattern accomodates the GRAPPA pattern...
+      indx  = plist{cnt}(cnt2) + v;         % feedforward indices
+      % indx2 = plist{cnt}(cnt2) + f;         % feedback indices
+
+      %% determine the k-space filling coefficients 
+      for cnt3=1:length(indx3),
+        y0 = mod( indx3(cnt3) + [1:size(Sk,2)] , size(Sk,2) );
+        y0( y0 == 0 ) = size(Sk,2);
+        
+        % generate a hybrid source matrix
+        if (dualks)
+          tmpSksub = zeros( size( Sk(indx,y0,:) ) );
+          if ( sum(pattern{cnt}=='<') ), 
+            tmpSksub( 1:2:length(indx), :, : ) = kin{1}( indx(1:2:end),y0,: );
+            tmpSksub( 2:2:length(indx), :, : ) = kin{2}( indx(2:2:end),y0,: );
+%            tmpSksub( (v < 0), :, : )  = kin{2}( indx( (v < 0) ),y0,: );
+%            tmpSksub( (v >= 0), :, : ) = kin{3}( indx( (v >= 0) ),y0,: );
+          else
+            tmpSksub( 1:2:length(indx), :, : ) = kin{2}( indx(1:2:end),y0,: );
+            tmpSksub( 2:2:length(indx), :, : ) = kin{1}( indx(2:2:end),y0,: );
+%            tmpSksub( (v >= 0), :, : ) = kin{2}( indx( (v >= 0) ),y0,: );
+%            tmpSksub( (v < 0), :, : )  = kin{3}( indx( (v < 0) ),y0,: );
+          end
+        else
+          tmpSksub = Sk(indx,y0,:);
+        end;
+        % A( size(s,2)*(cnt2-1) + (1:size(s,2)), cnt3:length(indx3):size(A,2) ) = ...
+        A( size(s,2)*(cnt2-1) + (1:size(s,2)), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
+            reshape( permute( tmpSksub, [ 2 1 3 ]), size(Sk,2), size(Sk,3)*length(indx) );
+
+      end;
+      %% for an N-by-1 kernel, these next lines are enough...
+      % A2 = [ A2; ...
+      %         reshape( permute(squeeze( Sk(indx,:,:) ),[ 2 3 1 ]), ...
+      %                  size(s,2), size(s,3)*length(indx) )
+      %       ];
+      % 
+      % b2 = [b2; squeeze( Sk( plist{cnt}(cnt2), :, l ) ).' ];
+    end;
+
+    % normalize, following to Algo 4.2 of M. Schneider, "GPGPU for Accelerated GRAPPA Autocalibration in
+    %    Magnetic Resonance Imaging," Masters Thesis, U of Erlangen, Apr 2008
+    % A_bak = A;
+    p_src = sum(abs(A),2).^(-eta/2);    % compute the power in each line of the linsys
+    A = tmult(A,diag(p_src),1);         % normalize by the power in each linsys line
+
+    % Apinv = pinv(A);
+    % AA = A'*A;
+    AA = [A ]'*[A ];
+    At = [A ]';
+    clear b2;
+    for l=1:size(Sk,3),
+      rptupd(prnt,'.');
+      % A is the same for each coil ...
+      b = vec( squeeze(Sk( plist{cnt}, :, l )).' ); % for dual-kernel, this defaults
+                                                    % to the target data
+
+      b = b.*p_src(:);                  % normalize by the power in each linsys line
+
+      b2(:,l) = b;
+      n(:,l) = A \ b;
+
+      if (verbose>2), pltcmplx( At'*n(:,l), b ); keyboard; end;
+    end;
+    b1 = b2;
+    N(length(N)).n = n;
+    
+    rptupd(prnt,'\n');
+  end;  
+end;
+
+rptupd(prnt,'recon missing lines\n');
+%% now go recon the missing lines of k-space:
+
+if (verbose); keyboard; end;
+
+kypts = 1:length(mss);
+kypts( 1:floor(length(mss)/2) ) = kypts( (floor(length(mss)/2)):-1:1 );
+
+linerpt = repmat(' ',1,length(mss));
+
+for kycnt=1:length(mss),
+  cnt2 = kypts(kycnt);
+  
+  if verbose, rptupd(prnt,sprintf(' line %3d, from lines ',cnt2)); end;
+
+  if ( find( Y == cnt2) ); 
+    % rptupd(prnt,'*'); 
+    linerpt(cnt2) = '*';
+    if verbose, rptupd(prnt,'\n'); end;
+    continue; 
+  end;
+
+  tmp = [];
+  %% determine which recon pattern we can use:
+  for cnt4 = 1:length(pattern),
+    v = extrapolate_pattern(pattern{cnt4},nky);
+
+    if ( (kycnt < length(mss)/2) && sum(find( pattern{cnt4} == '-')) ), continue, end;
+    if ( (kycnt > length(mss)/2) && sum(find( pattern{cnt4} == '+')) ), continue, end;
+
+    indx  = cnt2 + v;
+    % indx2 = cnt2 + f;
+    %     tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
+    %     if  ( ... % ( (min(indx) > 1) & (max(indx) < size(Sk,1)) ) & ...
+    %         ( length(tmp) == length(indx) ) ),
+    %       break;                          % jump out of the loop
+    %     end;
+    if ( indx(1) < 1) | (indx(end) > size(Sk,1)); continue; end;
+    
+    tmp = (mss( indx ) == ones(length(indx),1));
+    if ( sum(tmp) == length(tmp)  ),
+      break; 
+    end;
+
+  end;  
+  
+  if ( ( (min(indx) < 1) | (max(indx) > size(Sk,1)) ) | ...
+       ( sum(tmp) ~= length(indx) ) ),
+    % 
+    linerpt(cnt2) = 'o';
+    if verbose, rptupd(prnt,'o'); rptupd(prnt,'\n'); end;
+    continue; 
+  end;
+  
+  %% recon k-space line based on that pattern:
+  rsz = size(Sk,2);
+  A = zeros( rsz, length(indx3)*length(v)*size(s,3) );
+  B = []; % zeros( rsz, length(indx3)*length(f)*size(s,3) );
+  
+  for cnt3=1:length(indx3),
+    y0 = mod( indx3(cnt3) + [(size(Sk,2)-rsz)/2+(1:rsz)] , size(Sk,2) );
+    y0( y0 == 0 ) = size(Sk,2);
+    
+    A( (1:rsz), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
+        reshape( permute( Sk(indx,y0,:), [ 2 1 3 ]), ...
+                 rsz, size(Sk,3)*length(indx) );
+  end;
+
+  if (verbose), % ( sum(Sk(Y(mss(cnt2))+1,:,l)) ~= 0 ); 
+    rptupd(prnt,sprintf(' %3d',[ indx(:) ]));
+    rptupd(prnt,sprintf(' (pattern: %s)', pattern{cnt4}) );
+    rptupd(prnt,'\n');
+    if (verbose>1),
+      figure(1); imagesc( log10(abs(Sk(:,:,1)) + 1e-1) ); pause(0.5);
+      keyboard;
+    end;
+  end;
+  linerpt(cnt2) = 'x';
+
+  for cnt3=1:length(N),
+    if strcmp( N(cnt3).pattern, pattern{cnt4} ),
+      break;                            % jump out of the loop
+    end;
+  end;
+  n = N(cnt3).n;
+  
+  for l=1:size(Sk,3),
+    % for each coil ...
+
+    % if cnt2==109, keyboard; end;
+    % rptupd(prnt,(' line %3d, pattern %3d\n', cnt2, plist(cnt3) );
+    Sk( cnt2, :, l ) = [A B] * n(:,l);
+  end;
+
+  if (verbose>1), keyboard; end;
+end;
+rptupd(prnt,linerpt);
+rptupd(prnt,'  \n');
+
+% Sk(1:2:size(Sk,1),:,:) = -Sk(1:2:size(Sk,1),:,:);
+% Sk(:,1:2:size(Sk,2),:) = -Sk(:,1:2:size(Sk,2),:);
+
+F = Sk;
+I = sqrt( sum(abs(ifft2(Sk)).^2,3) ); 
+I = I./(prod(size(I)));
+mask = 0;
+
+rpttime = etime(clock,st);
+rptupd(prnt,sprintf(' GRAPPA recon time: %g (avg per kernel: %g)\n',  ... 
+        rpttime, rpttime/length(pattern)) ); 
+
+return;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function rptupd( prnt, str )
+
+if (prnt)
+  fprintf( str );
+end;
+
+return;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function v = extrapolate_pattern(p,YKS)
+ 
+r = find( (p == 'x') | (p == '<') | (p == '>') );
+v = find( (p == '*') | (p == '<') | (p == '>') ) - r;
+
+if length(v)>YKS,
+  v = v( find( v~= 0 ) );
+end;
+return;
+
+
diff --git a/dpg_recon.m b/dpg_recon.m
new file mode 100644
index 0000000..2f3757a
--- /dev/null
+++ b/dpg_recon.m
@@ -0,0 +1,147 @@
+function [Sk,Sk2] = dpg_recon( Fin, N, dk, nky );
+% reconstruct measured EPI data using Dual-Polarity GRAPPA kernels
+
+%
+% Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+Sk  = zeros(size(Fin{1}));
+if nargout > 1,
+  Sk2 = zeros(size(Fin{1}));
+end;
+sz = size(Sk);
+
+
+if nargin < 4,
+  nky = 2;
+end;
+if isstr(nky),
+  kernel = nky;
+  nky = str2num(kernel(1));
+end;
+
+if ( ~isfield( N, 'n') & isfield( N, 'coef' ) )
+  for cnt=1:length(N),
+    N(cnt).n = N(cnt).coef;
+  end;
+end;
+
+nkx = size(N(1).n,1) / size(Fin{1},3) / nky;
+ncoils = size(Sk,3);
+
+% fprintf(' ky:%d kx:%d nc:%d\n',nky,nkx,ncoils);
+indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
+
+verbose=1;
+prnt = 0;
+
+% for kycnt=2:2*dk:(size(Sk,1)-2*dk),
+for kycnt=2:2*dk:(size(Sk,1)),
+
+  % fprintf('---%d ---\n',kycnt);
+  for patcnt = 1:length(N); % [ 2:length(N)/2  (length(N)/2+(2:length(N)/2)) ];
+    
+    % fprintf(' %d',patcnt)
+    v = extrapolate_pattern(N(patcnt).pattern,nky);
+
+    if (patcnt <= length(N)/2)
+      indx  = kycnt + (patcnt-1) + v;
+    else
+      indx  = kycnt + (patcnt-2) + v;
+    end;        
+
+    % fprintf(': '); fprintf(' %d', indx); fprintf('\n');
+    if sum(indx>sz(1));
+      indx = mod( indx, sz(1) );
+      indx( indx==0 ) = sz(1);
+    end;
+    % fprintf('# '); fprintf(' %d', indx); fprintf('\n');
+
+    %% recon k-space line based on that pattern:
+    A = zeros( size(Sk,2), length(indx3)*length(v)*size(Sk,3) );
+  
+    for cnt3=1:length(indx3),
+      y0 = mod( indx3(cnt3) + [1:size(Sk,2)] , size(Sk,2) );
+      y0( y0 == 0 ) = size(Sk,2);
+
+      tmpSk = zeros([ length(indx) length(y0) size(Sk,3) ]);
+  
+      if (patcnt <= length(N)/2)
+        tmpSk(1:2:end,:,:) = Fin{1}( indx(1:2:end), y0, :);
+        tmpSk(2:2:end,:,:) = Fin{2}( indx(2:2:end), y0, :);
+      else
+        tmpSk(1:2:end,:,:) = Fin{2}( indx(1:2:end), y0, :);
+        tmpSk(2:2:end,:,:) = Fin{1}( indx(2:2:end), y0, :);
+      end;
+
+      A( (1:size(Sk,2)), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
+          reshape( permute( tmpSk, [ 2 1 3 ]), ...
+                   size(Sk,2), size(Sk,3)*length(indx) );
+
+    end;
+
+    if (verbose), % ( sum(Sk(Y(mss(cnt2))+1,:,l)) ~= 0 ); 
+      tmp = kycnt+patcnt; 
+      if (patcnt > length(N)/2); tmp = tmp-1; end;
+      rptupd(prnt,sprintf(' %3d',[ tmp; indx(:) ]));
+      rptupd(prnt,sprintf(' (pattern: %s)', N(patcnt).pattern) );
+      rptupd(prnt,'\n');
+      if (verbose>1),
+        figure(1); imagesc( log10(abs(Sk(:,:,1)) + 1e-1) ); pause(0.5);
+        keyboard;
+      end;
+    end;
+    % linerpt(cnt2) = 'x';
+
+    n = N(patcnt).n;
+    
+    for l=1:size(Sk,3),
+      % for each coil ...
+
+      % if cnt2==109, keyboard; end;
+      % rptupd(prnt,(' line %3d, pattern %3d\n', cnt2, plist(cnt3) );
+
+      if (nargout == 1),
+        tmp = kycnt+(patcnt);
+        if (patcnt > length(N)/2), tmp = tmp-1; end;
+        Sk( tmp, :, l ) = [A] * n(:,l);
+      else
+        if (patcnt > length(N)/2 ), 
+          Sk2( kycnt+(patcnt-1), :, l ) = [A] * n(:,l);
+        else
+          Sk( kycnt+(patcnt), :, l ) = [A] * n(:,l);
+        end;
+      end;
+    end;
+
+    % if (verbose>1), keyboard; end;
+  end;
+
+end;
+
+Sk = Sk(1:sz(1),:,:);
+
+%% need to fix this above, but for now...
+Sk = Sk([2:sz(1) 1],:,:);               % rotate the output back 1 kspace line
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function rptupd( prnt, str )
+
+if (prnt)
+  fprintf( str );
+end;
+
+return;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function v = extrapolate_pattern(p,YKS)
+ 
+r = find( (p == 'x') | (p == '<' ) | (p == '>') );
+v = find( (p == '*') | (p == '<' ) | (p == '>') ) - r;
+
+if length(v)>YKS,
+  v = v( find( v~= 0 ) );
+end;
+return;
diff --git a/fif.m b/fif.m
new file mode 100644
index 0000000..565b7d8
--- /dev/null
+++ b/fif.m
@@ -0,0 +1,3 @@
+function a = fif(a)
+
+a = fftshift(fftshift(ifft2(fftshift(fftshift( a, 1),2)),1),2);
diff --git a/ifi.m b/ifi.m
new file mode 100644
index 0000000..fc55937
--- /dev/null
+++ b/ifi.m
@@ -0,0 +1,2 @@
+function a = ifi(a)
+a = ifftshift(ifftshift(fft2(ifftshift( ifftshift( a, 1),2)),1),2);
diff --git a/phzapply.m b/phzapply.m
new file mode 100644
index 0000000..0a7e397
--- /dev/null
+++ b/phzapply.m
@@ -0,0 +1,50 @@
+function [a_out] = phzapply(a,y);
+% [a_out] = phzapply(a,y);
+%
+% apply the EPI phase shift parameters calcuated by phzshift or
+% comp_local_pc to the original k-space data
+%
+
+%
+% Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+
+[m n c]=size(a);
+
+[pm pc ] = size(y);
+
+A1 = a(1:2:end,:,1:c);
+A1a = fftshift(ifft(fftshift(A1,[ 2]),[],2),[ 2]);
+
+A2 = a(2:2:end,:,1:c);
+A2a = fftshift(ifft(fftshift(A2,[ 2]),[],2),[ 2]);
+
+A1b = zeros(size(A1a));
+A2b = zeros(size(A2a));
+for cnt=1:c,
+  A1b(:,:,cnt) = A1a(:,:,cnt) * diag( exp( +j*y(1,cnt)/2*((1:size(A2a,2))-size(A2a,2)/2)));
+  A2b(:,:,cnt) = A2a(:,:,cnt) * diag( exp( -j*y(1,cnt)/2*((1:size(A2a,2))-size(A2a,2)/2)));
+end;
+
+A1c = ifftshift(fft(ifftshift(A1b,[ 2]),[],2),[ 2]);
+A2c = ifftshift(fft(ifftshift(A2b,[ 2]),[],2),[ 2]);
+
+%% keyboard;
+
+for cnt=1:c,
+  if ~isreal( y(2,cnt) ),
+    phz = angle(y(2,cnt));
+  else,
+    phz = y(2,cnt);
+  end;
+  A1c(:,:,cnt) = A1c(:,:,cnt) .* exp(-j*phz/2); % angle(y(2,cnt)));
+  A2c(:,:,cnt) = A2c(:,:,cnt) .* exp(j*phz/2); % angle(y(2,cnt)));
+end;
+
+a_out = zeros(size(a));
+a_out(1:2:end,:,:) = A1c;
+a_out(2:2:end,:,:) = A2c;
+
diff --git a/phzshift.m b/phzshift.m
new file mode 100644
index 0000000..fcdce8a
--- /dev/null
+++ b/phzshift.m
@@ -0,0 +1,282 @@
+function varargout = phzshift(a,b,mode,verbose,phz,thrsh);
+% [c_out,y_out] = phzshift( a, b, [mode, verbose, phz]);
+%
+% calculates the coordinate shift between two k-space data sets, a and b.
+% Using the Phase Correlation Method, each k-space set is transformed
+% to the image domain, the phase difference is calculated and then
+% applied to the second set.  After this correction, the two data sets
+% are added together and output.
+%
+% inputs can be any number of dimensions, but 'a' and 'b' must be the
+% same size
+%
+% c_out is the combination of (a+b)/2, after phase correction.
+% 
+% y_out contains the estimated phase shift between 'a' and 'b' for
+% each frame.
+%
+% Advanced (optional) options:
+%    mode - a cell variable containing flags to various settings.
+%           'nofft'    : The input data is already in the spatial domain
+%           'nocombo'  : Returns  b with the phase correction applied.
+%           'nofilter' : use the full data to estimate the phase shift, 
+%                        instead of the default setting which uses a
+%                        low-passed filtered version.
+%           'center'   : center the low-pass filter.
+% verbose - set to '1' to show debugging information
+%     phz~ - if nonempty, use this matrix instead of estimating the phase
+%           shift.  (i.e. y_out from a previous run)
+%   thrsh - the mask threshold for the svd calculation. range: [0:1]
+
+%
+% Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+
+% 2013/nov/07: jonathan polimeni <jonp at nmr dot mgh dot harvard dot edu>
+% fixed non-integer index, fixed missing conditional for "nofilter" option
+
+% 2014/sep/25: jonathan polimeni <jonp at nmr dot mgh dot harvard dot edu>
+% replaced phase difference fitting with Local Phase Correction (LPC) method
+% based on the methods described in:
+%
+%   Feiweier T. Magnetic Resonance Method and Apparatus to Determine Phase
+%   Correction Parameters. U.S. Patent No. 8,497,681, Issued 7/30/2013.
+
+  
+[m,n,c]=size(a);
+c_out = zeros(m,n,c);
+if nargout > 1,
+  y_out = zeros(2,c);
+end;
+
+if nargin < 3,
+  mode = {''};
+end;
+
+if nargin < 4,
+  verbose = 0;
+end;
+
+if nargin < 5,
+  phz = [];
+end;
+if ~isempty(phz),
+  phz = reshape(phz,[2 c]);
+end;
+
+if nargin < 6,
+  thrsh = 0.2;
+end;
+
+
+if ( sum(strcmp(mode,'nofilter')) || sum(strcmp(mode,'nofft')) ),
+  f1 = 1; f2 = 1;
+else
+  % define a low-pass envelope to apply in kspace
+  
+  alph = 20;
+  
+  if (alph > m), alph2 = m; else, alph2 = alph; end;
+  f1 = [ zeros(floor((m-alph2)/2),1); gausswin(alph2); zeros(ceil((m-alph2)/2),1)]; 
+
+  if (alph > n), alph2 = n; else, alph2 = alph; end;
+  f2 = [ zeros(floor((n-alph2)/2),1); gausswin(alph2); zeros(ceil((n-alph2)/2),1)]; 
+
+  [tmp,maxx] = max(max(sum(abs(a(:,:,:)),3) ));
+  [tmp,maxy] = max(max(sum(abs(a(:,:,:)),3)'));
+
+  %%% correct the window coordinates to set ontop of DC (which may be off center)
+  [ans,j2]=max(f1);
+  [ans,j3]=max(f2);
+  
+  if sum( strcmp(mode,'center') )
+    c1 = 0; c2 = 0;
+  else,
+    c1 = j2 - maxy;
+    c2 = j3 - maxx;
+  end;
+  % fprintf('c1: %d c2: %d\n',c1,c2);
+  
+  if c1>0;
+    f1 = f1( [ c1:length(f1) 1:(c1-1)]);
+  elseif c1<0
+    f1 = f1( [ length(f1)+(c1:0) 1:(length(f1)+c1-1)]);
+  end;
+  
+  if c2>0;
+    f2 = f2( [ c2:length(f2) 1:(c2-1)]);
+  elseif c2<0
+    f2 = f2( [ length(f2)+(c2:0) 1:(length(f2)+c2-1)]);
+  end;
+
+end;
+
+
+% convert k-space to image-space
+if sum(strcmp(mode,'nofft'));
+  A = a;
+  B = b;
+elseif sum(strcmp(mode,'nofilter')),
+  A = fftshifthift(ifft2(fftshifthift( reshape(a,[m n c]) )));
+  B = fftshifthift(ifft2(fftshifthift( reshape(b,[m n c]) )));
+else,
+  A = fftshifthift(ifft2(fftshifthift( repmat( (f1*f2'),[1 1 c]).*reshape(a,[m n c]) )));
+  B = fftshifthift(ifft2(fftshifthift( repmat( (f1*f2'),[1 1 c]).*reshape(b,[m n c]) )));
+end;
+
+mask = ( abs(A) > thrsh*max( abs(A(:)) ) );
+
+A1 = reshape( permute(A,[2 1 3 4]),[n m*c]).'; 
+B1 = reshape( permute(B,[2 1 3 4]),[n m*c]).';
+M1 = reshape( permute(mask,[2 1 3 4]),[n m*c]).';
+
+if (1),
+  if isempty(phz)
+    % calculate the normalized cross-correlation
+    C = zeros(size(A1));
+    % [JRP] why does next definition of "tt" below use a threshold of 0.1*max?
+    tt = find( abs(A1(:)) > 0.01*max(abs(A1(:))) );
+    C(tt) = A1(tt).*conj(B1(tt))./abs(A1(tt).*conj(A1(tt)));
+  
+    % extract out the dominant x and y shift vectors
+    [u d v]=svds(M1.*C,1);
+    % v = rsvd( mask.*C, 1 );
+  
+
+    m2 = sum(M1,1).';
+    t2 = find(m2~=0);
+
+    % [JRP] previously the vector "v" was unwrapped after extracting a
+    % subset defined by vector "t2", but since the subset of samples may be
+    % non-contiguous the unwrapping could add erroroneous phase
+    % jumps---better to unwrap the vector first then extract the subset
+    % afterwards
+    unv = unwrap(angle(v));
+    
+    % calculate the shift from the slope of the phase.
+    % here, we are only concerned with shifts along rows
+    if (verbose),
+      fprintf(' %d ',t2);
+      fprintf('\n');
+      Atmp = [ ones(length(t2),1) t2(:)-size(B,2) ];
+      fprintf('A''*A  =');
+      fprintf(' %d ', Atmp'*Atmp );
+      fprintf('\n');
+      fprintf('A''*v2 =');
+      fprintf(' %f ', Atmp'*unv(t2) );
+      fprintf('\n');
+    end;
+    if (0),
+      % [JRP] "robustfit" performs better than backslash when outliers are
+      % present and are not properly masked
+      %y2 = ( [ ones(length(t2),1) t2(:)-size(B,2) ] ) \ unv(t2);
+      y2 = robustfit(t2(:)-size(B,2), unv(t2) );
+    else,
+      % [JRP] based on "Local Phase Correction" of Feiweier
+      dphz = angle(sum(v(2:end) .* conj(v(1:end-1))));
+      y2 = [1; dphz];
+    end;
+    if (verbose)
+      figure; plot( t2-size(B,2), unv(t2), 'bx-', t2-size(B,2), [ ones(length(t2),1) t2(:) ]*y2,'k-.' );
+
+      fprintf(' x = %f  %f \n', y2 );
+      keyboard;
+    end;
+  else,
+    y2(2) = phz(1);
+  end; % if exist('phz')
+  
+  % correct the phase of the second data set to match the first
+  % (i.e. shift the k-space grids to match)
+  if sum(strcmp(mode,'nofft'));
+    A = a;% (:,:,cnt);
+    B = b;% (:,:,cnt);
+  else,
+    A = fftshift(ifft2(fftshift( a,1:2 )), 1:2 );
+    B = fftshift(ifft2(fftshift( b,1:2 )), 1:2 );
+  end;
+  for cnt=1:c,
+    if sum(strcmp(mode,'nocombo'))
+      B(:,:,cnt) = B(:,:,cnt) * diag( exp( -j*y2(2)*( (1:size(B,2))-size(B,2)/2 ) ) );
+    else 
+      A(:,:,cnt) = A(:,:,cnt) * diag( exp( +j*y2(2)/2*( (1:size(A,2))-size(A,2)/2 ) ) );
+      B(:,:,cnt) = B(:,:,cnt) * diag( exp( -j*y2(2)/2*( (1:size(B,2))-size(B,2)/2 ) ) );
+   end;
+  end;
+%  m1 = sum(mask,2).';
+%  t1 = find(m1~=0);
+%  % calculate the shift from the slope of the phase.
+%  % here, we are only concerned with shifts along rows
+%  y1 = ( [ ones(length(t1),1) t1(:) ] ) \ unv(t1);
+%  
+%  % correct the phase of the second data set to match the first
+%  % (i.e. shift the k-space grids to match)
+%  B = diag( exp( j*y1(2)*((1:size(B,1))-1) ) ) * B;
+
+  % calculate the scalar phase difference between the sets ...
+  if isempty(phz)
+    if (0),
+    P = zeros(size(A));
+    % [JRP] why does previous definition of "tt" above use a threshold of 0.01*max?
+    tt = find(A(:)>0.1*max(A(:)));
+    P(tt) = A(tt).*conj(B(tt))./abs(A(tt).*conj(A(tt)));
+    indx1 = round((size(P,1) - 10)/2) + (1:10); if ( sum(indx1<=0) > 0 ) indx1 = 1:size(P,1); end;
+    x = mean(vec(P(  indx1, (end-10)/2 +(1:10), : ) ));
+    else,
+      % [JRP] based on "Local Phase Correction" of Feiweier
+      x = sum(vec(A.*conj(B)));
+    end;
+  else,
+    x = phz(2);
+  end;
+
+  % and apply a scalar correction as well.
+  B = B.*exp(j*angle(x));
+
+  if nargout > 1,
+    y_out(:,cnt) = [ y2(2); exp(j*angle(x)) ];
+    if (verbose)
+      fprintf(' p = %f  %f + j %f\n',y_out(1,cnt), ...
+      real(y_out(2,cnt)),imag(y_out(2,cnt)) );
+    end;
+  end;
+
+  % add the two sets together, and convert back to k-space for output.
+  if sum(strcmp(mode,'nocombo'))
+    c_out = B;
+  else,
+    c_out = (A+B)/2;
+  end;
+  if sum(strcmp(mode,'nofft')),
+  else,
+    c_out = ifftshift(fft2(ifftshift(c_out,1:2)),1:2);
+  end;
+
+end;
+
+
+c_out = reshape( c_out, size(a) );
+sz = size(c_out);
+if ( (nargout > 1) && length(sz)>2),
+  y_out = reshape( y_out, [ 2 sz(3:end) ] );
+end;
+
+% return output
+if nargout > 0
+  varargout{1}= c_out;
+end
+
+if nargout > 1
+  varargout{2}= y_out;
+end
+ 
+if nargout > 2
+  varargout{3}= y2;
+end
+
+if nargout > 3
+  varargout{4}= x;
+end
diff --git a/readnd.m b/readnd.m
new file mode 100755
index 0000000..4730333
--- /dev/null
+++ b/readnd.m
@@ -0,0 +1,81 @@
+function [A,meta] = readnd(fname);
+%
+% READND - read a binary N-dimensional data file
+%
+% The function reads N-dimensional data saved by either the C or Matlab
+% form of savend.
+%
+%    See Also SAVEND
+%
+
+%
+% Copyright 2004-2008  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+
+fid = fopen(fname,'rb');
+if fid <0, error(['Couldn''t open ' fname]); end;
+
+fseek(fid,0,'eof');
+maxfln = ftell(fid);
+fseek(fid,0,'bof');
+
+fln = fread(fid,1,'integer*4');
+hdr = fscanf(fid,'%c',4);          
+if ~strcmp(hdr,'nddf')
+  warning(['This does not appear to be a .nd file.  The file type field is ' ...
+           'not ''nddf'' ']);
+  return;
+end;
+ 
+szln = fread(fid,1,'integer*4');    
+  sz = fread(fid,szln,'integer*4');
+ hdr = fscanf(fid,'%c',4);
+
+if ( strcmp(hdr(4),'c') ),
+  cmplx = 2;
+else
+  cmplx = 1;
+end;
+   
+if strcmp( hdr(1:3), 'dbl' ),
+  data = fread(fid,cmplx*prod(sz),'double');
+elseif strcmp( hdr(1:3), 'flt' ),
+  data = fread(fid,cmplx*prod(sz),'float');
+elseif strcmp( hdr(1:3), 'int' ),
+  data = fread(fid,cmplx*prod(sz),'int16');
+elseif strcmp( hdr(1:3), 'cha' ),
+  data = fread(fid,cmplx*prod(sz),'uint8');
+end;
+fprintf('format: %s\n',hdr);
+
+if ( strcmp(hdr(4),'c') ),
+  data = data(1:2:length(data)) + j * data(2:2:length(data));
+end;
+if ( prod(size(data)) ~= prod(sz) ), 
+  fprintf('declared: %10d\n',prod(sz));
+  fprintf('    read: %10d\n',length(data));
+  warning('Error: size of read data doesn''t match declared size');
+else
+  if length(sz)==1, sz = size(data)'; end;
+  A = reshape( data, sz' );
+end;
+
+if ( ftell(fid) ~= maxfln ), 
+
+  metaln = fread(fid,1,'integer*4');
+  hdr = fscanf(fid,'%c',4);  
+  
+  if strcmp( hdr, 'meta' ),
+    meta = fread(fid,metaln,'char');
+  else,
+    warning('Warning: number of bytes read from file is suspicious'); 
+    disp([ ftell(fid) fln ]);
+    keyboard;
+  end;
+end;
+
+fclose(fid); 
+
diff --git a/recongrappa.m b/recongrappa.m
new file mode 100755
index 0000000..d4e77d1
--- /dev/null
+++ b/recongrappa.m
@@ -0,0 +1,242 @@
+function [F,I,mask] = recongrappa( varargin );
+% [F,I,mask ] = recongrappa(Wk, Sk, Y, [opt args]);
+% 
+% opt args:
+%  'v', {0,1}              - verbosity off/on
+%  'kernel', string        - grappa kernel size, x by y. e.g. '2x5', '4x1'
+%  'N', (xy)-by-L matrix   - grappa reconstruction parameters
+%  'acs', array of indices - index locations to use for ACS lines
+%  'dk', delta-k           - distance between phase encode lines used by
+%                            the GRAPPA kernel.  default: 2
+% 
+% This is an implementation of     
+% "Generalized autocalibrating partially parallel acquisitions"
+%   by M. A. Griswold, P. M. Jakob, et. al..
+%  Mag. Reson. Med.   47(6):1202-1210.  Jun 2002
+
+% Copyright 2004-2007  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+%
+% uncombined images are generated for each coil in the array, by
+% filling in k-space lines.
+%
+%
+% GRAPPA: data from all coils used to patch data in one coil
+%
+%  coil 1     o * o * - * o *
+%  coil 2     o * o * - * o *
+%  coil 3     o * o * - * o *
+%  coil 4     o * o * x * o *
+%
+%  *: acquired line, o: unacquired line, -: auto-callibration line
+%  x: reconstructed line
+%
+%  S_j (k_y - m\Delta k_y ) = \sum_{l=1}^L \sum_{b=0}^{N_b-1}
+%                                  n(j,b,l,m) S_l(k_y - b A \Delta k_y)
+%
+% in a variable density sampling scheme, the center of k-space is fully
+% sampled and can be used to determine the linear fit parameters.
+%
+% A   - acceleration factor
+% N_b - number of blocks used in the reconstruction
+% 
+% N_b = 1 ----> GRAPPA = VD-AUTO-SMASH
+
+global prnt
+
+Wk = varargin{1};
+s  = varargin{2};
+Y  = varargin{3};
+
+acs = []; N = []; verbose=0; kernel = '2x3'; dk = 2;
+prnt = 0;   % set to '1' to output some debugging information
+Nbias = [];
+chi = 1e-8;
+
+vararg = varargin(4:end);
+for cnt=1:2:length(vararg),
+  if strcmp( vararg{cnt}, 'v' ),
+    verbose = cell2mat(vararg(cnt+1));
+  elseif strcmp( vararg{cnt}, 'N' ),
+    N = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'Nbias' ),
+    Nbias = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'kernel' ),
+    kernel = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'dk' ),
+    dk = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'chi' ),
+    chi = vararg{cnt+1};
+  elseif strcmp( vararg{cnt}, 'acs' ),
+    acs = vararg{cnt+1};
+  elseif (strcmp( vararg{cnt}, 'report' ) | strcmp( vararg{cnt}, 'r' )),
+    prnt = vararg{cnt+1};
+  end;
+end;
+
+if isempty(acs), 
+  d = diff(Y);
+  acs = Y(find(d == 1));
+  % acs = acs(2:length(acs));
+end;
+
+if length(size(Wk)) == 3, 
+  [Wm,Wn,Wc] = size(Wk);
+else,
+  Wm = Wk(1); Wn = Wk(2); Wc = Wk(3);
+end;
+
+%% zero-pad the k-space data
+Sk = zeros( Wm, Wn, Wc );
+Sk(Y,:,:) = s;
+
+% set the kernel spacings based on 'kernel' string
+nky = str2double(kernel(1));               % should be even
+if (nky==1),
+  indx0 = 1;
+else,
+  indx0 = floor( [-(nky-1):2:(nky-1)] * dk/2 );
+end;
+% if (nky == '2'), indx0 = [-1 1];
+% elseif (nky == '6'), indx0 = [ -5 -3 -1 1 3 5];
+% else % (nky == '4'), indx0 = [ -3 -1 1 3];
+% end;
+
+nkx = str2double(kernel(3));               % should be odd
+indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
+
+rptupd( sprintf(['kernel size: ' num2str(length(indx0)) 'x' num2str(length(indx3)) ';']) );
+rptupd( sprintf([' indx0 = [' sprintf(' %d',indx0) ']; ']) );
+rptupd( sprintf([' indx3 = [' sprintf(' %d',indx3) '];\n']) );
+
+stime = clock;
+
+plist = zeros( length(acs), 1 );
+
+for cnt2 = 1:length(acs);
+  %% check if the GRAPPA pattern...
+  indx = acs(cnt2) + indx0;
+  
+  %% can be accomodated by the sampling pattern 
+  tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
+  plist(cnt2) =  ( length(tmp) == length(indx) );
+end;
+%% use only those ACS lines that fit the GRAPPA kernel
+acs = vec( acs(logical(plist)) ).';
+rptupd([ ' ACS lines:' sprintf(' %d ',acs) '\n' ]);
+
+rptupd('processing coil: ');
+for l=1:size(s,3),
+  rptupd([' ' num2str(l)]);
+  %% determine the k-space filling coefficients for each coil:
+
+  %%% setup the system matrix
+  A2 = zeros( size(s,2), length(indx3)*length(indx0)*size(s,3) );
+
+  if size(N,2) < l,                     % if coefficients are not supplied
+  %%% Griswolds way:
+  %   tic;
+  %   A2 = zeros( size(s,2), length(indx3)*length(indx)*size(s,3) );% 20*size(s,3) );
+  %   for x0=1:size(s,2),
+  %     indx2 = x0 + indx3;
+  %     indx2( indx2 < 1 ) = indx2( indx2 < 1 ) + size(s,2);
+  %     indx2( indx2 > size(s,2) ) = indx2( indx2 > size(s,2) ) - size(s,2);
+  %     
+  %     A2(x0,:) = vec( permute(squeeze( Sk(indx,indx2,:)),[2 3 1 ]) ).';
+  %   end;
+  %   toc
+  
+  %%% My way: (heh, 8x faster)
+  % tic;
+  A = zeros( length(acs)*size(s,2), length(indx3)*length(indx0)*size(s,3) );
+  A2 = zeros( size(s,2), length(indx3)*length(indx0)*size(s,3) );
+  if (l==1), rptupd(['\n size(A) = ' num2str(size(A)) '\n']), end;
+
+  for cnt2=1:length(acs),
+    indx = acs(cnt2) + indx0;
+    for cnt=1:length(indx3),
+      y0 = mod( indx3(cnt) + [1:size(Sk,2)] , size(Sk,2) );
+      y0( y0 == 0 ) = size(Sk,2);
+      
+      A( size(s,2)*(cnt2-1) + (1:size(s,2)), ...
+         cnt:length(indx3):size(A,2) ) = ...
+          reshape( permute(squeeze( Sk(indx,y0,:) ),[ 2 3 1 ]), ...
+                   size(Sk,2), size(Sk,3)*length(indx) );
+    end;
+  end;
+  % toc
+
+  %%% setup the objective vector
+  b = vec( Sk( acs, :, l ).' );
+  
+  %%% find the coefficients
+  % n = A \ b; % pinv( A, norm(A)*0.01 ) * b;  
+  
+  if (size(Nbias,1)==0), % isempty('Nbias'),
+    % n = cgsolv( A'*A, A'*b, zeros(size(A,2),1), size(A,2) );
+    [ n tt ] = bicg( A'*A + chi*eye(size(A,2)), A'*b );
+  else,
+    n = cgsolv( A'*A, A'*b, Nbias(:,l), size(A,2)/4 );
+  end;
+  
+  N(:,l) = n;                           % store them
+  if verbose, pltcmplx( A*n, b ); keyboard; end;
+  end;
+
+  % extract the coefficients for this coil
+  n = N(:,l);
+
+  %% now go recon the missing lines of k-space:
+  mss = zeros(Wm,1);            % set the missing lines index vector 
+  mss(Y) = 1;
+
+  for cnt0=1:length(mss),
+    if (mss(cnt0) == 1);   continue; end;
+
+    indx = cnt0 + indx0;
+    if ( (min(indx) < 1) | (max(indx) > size(Sk,1))); continue; end;
+
+    if ( sum(sum( Y(:)*ones(size(indx)) == ones(size(Y(:)))*indx )) < length(indx) ),
+      continue; 
+    end;
+
+    %     for x0=1:size(s,2),
+    %       indx2 = x0 ; % + [-2:2];
+    %       indx2( indx2 < 1 ) = indx2( indx2 < 1 ) + size(s,2);
+    %       indx2( indx2 > size(s,2) ) = indx2( indx2 > size(s,2) ) - size(s,2);
+    %       
+    %     end;
+    for cnt=1:length(indx3),
+      y0 = mod( indx3(cnt) + [1:size(Sk,2)] , size(Sk,2) );
+      y0( y0 == 0 ) = size(Sk,2);
+        
+      A2( (1:size(s,2)), cnt:length(indx3):size(A2,2) ) = ...
+          reshape( permute(squeeze( Sk(indx,y0,:) ),[ 2 3 1 ]), ...
+                   size(Sk,2), size(Sk,3)*length(indx) );
+    end;
+    Sk( cnt0, :, l ) = A2 * n;
+    % figure(1); imagesc( log10(abs(Sk(:,:,l)) + 1e-1) ); pause(0.5);
+  end;
+
+end;
+rptupd('\n');
+rptupd( sprintf( 'Elapsed time: %8.5g seconds.\n', etime(clock,stime) ) );
+
+F = Sk;
+I = sqrt( sum(abs(ifft2(Sk)).^2,3) ); 
+mask = N;
+
+clear global prnt
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function rptupd( str )
+global prnt
+
+if (prnt)
+  fprintf( str );
+end;
+
+return;
diff --git a/runthis.m b/runthis.m
new file mode 100644
index 0000000..776504b
--- /dev/null
+++ b/runthis.m
@@ -0,0 +1,59 @@
+addpath('data');
+addpath('tensor');
+
+% load the calibration data
+acs1 = readnd('dpg_acs1.nd');               % RO+ lead data
+acs2 = readnd('dpg_acs2.nd');               % RO- lead data
+
+% interleave to form the RO+ image
+A = zeros(size(acs1));
+A(:,1:2:end,:) = acs1(:,1:2:end,:);
+A(:,2:2:end,:) = acs2(:,2:2:end,:);
+A = permute(A,[2 1 3]);
+
+% interleave to form the RO- image
+B = zeros(size(acs1));
+B(:,1:2:end,:) = acs2(:,1:2:end,:);
+B(:,2:2:end,:) = acs1(:,2:2:end,:);
+B = permute(B,[2 1 3]);
+
+% combine coherently, to form the PLACE image
+% the lines below calls a MEX function from the Fast Imaging Library 
+% C = ifi(pos_neg_add( fif(A),fif(B) ))/2;
+C = ifi(phzshift(  fif(A), fif(B), {'nofft'}));
+
+% if the FIL function is not available, you can replace the above line(s) with
+% the call below, to read in the result from a file.
+% C = readnd('dpg_C.nd');
+
+% calibrate the DPG coefficients
+kin{1} = A; 
+kin{2} = B; 
+kin{3} = C;
+[~,~,N]=dpg_cal( size(kin{3}), kin, 1:size(kin{3},1), 'kernel','2x5','dks',2);
+
+% read in the raw EPI data (R=2)
+k = readnd('dpg_k.nd');
+
+% reconstruct the data using DPG
+Fin{1} = zeros([60 96 32]); Fin{1}(2:4:end,:,:) = permute(k(:,1:4:end,:),[2 1 3]);
+Fin{2} = zeros([60 96 32]); Fin{2}(4:4:end,:,:) = permute(k(:,3:4:end,:),[2 1 3]);
+Fdpg = dpg_recon( Fin, N, 2, '2x5');
+
+% compare with an LPC+GRAPPA recon:
+[~,~,Ng]=recongrappa(size(kin{3}),kin{3},vec(1:30),'kernel','2x5','dks',[2;4]);
+
+% compute the standard Nyquist ghost correction coefficients:
+s = readnd('dpg_k_pcref.nd');
+S0 = fif(mean(s(:,[1 3],:),2));
+S1 = fif(mean(s(:,[ 2 ],:),2));
+for cnt=1:size(S0,3); [~,y(:,cnt)]=phzshift( S1(:,:,cnt).', S0(:,:,cnt).', {'nofft'}); end;
+% y = readnd('y.nd');                     % derived from dpg_k_pcref.nd
+
+k_gc = phzapply( permute(k(:,1:2:end,:),[2 1 3]), y);
+Flpc = recongrappa([60 96 32],k_gc,vec(1:2:60),'kernel','2x5','dks',2,'N',Ng);
+
+imagesc(sqrt(sum(abs([ flipdim(fif(Flpc),1); flipdim(fif(Fdpg),1) ]).^2,3)));
+axis('image');
+colormap(gray(256))
+title([ 'top: 3-nav-line correction. bottom: DPG correction' ]);
diff --git a/runthis2.m b/runthis2.m
new file mode 100644
index 0000000..ce9f8ba
--- /dev/null
+++ b/runthis2.m
@@ -0,0 +1,32 @@
+% an example of DPG w/o PLACE calibration data
+addpath('data');
+addpath('tensor');
+
+acs2 = readnd('dpg_acs2.nd');           % RO- lead data
+y2 = readnd('y_acs2.nd');               % ghost correction coefficients
+
+% ghost-corrected ACS data
+c2_gc = permute( phzapply( permute(acs2,[2 1 3]), y2),[2 1 3]);
+
+% generate the GRAPPA coefficients 
+[~,~,Ng2]=grappa(size(c2_gc),c2_gc,vec(1:30),'2x5',2);
+
+% synthesize the RO+ and RO- source data from the ACS data using GRAPPA
+acs2n = grappa(size(c2_gc),acs2(:,1:2:end,:),vec(1:2:30),'2x5',2,Ng2);
+acs2p = grappa(size(c2_gc),acs2(:,2:2:end,:),vec(2:2:30),'2x5',2,Ng2);
+acs2c = ifi(pos_neg_add( fif(acs2n),fif(acs2p) ))/2;
+
+% calibrate DPG
+kin{1} = permute(acs2p,[2 1 3]);        % RO+ in the first cell
+kin{2} = permute(acs2n,[2 1 3]);        % RO- in the second cell
+kin{3} = permute(acs2c,[2 1 3]);
+[~,~,N2]=dpg_cal( size(kin{3}), kin, 1:size(kin{3},1), 'kernel','2x5','dks',2);
+
+
+% reconstruct the accelerated data using DPG
+k = readnd('dpg_k.nd');
+
+Fin{1} = zeros([60 96 32]); Fin{1}(2:4:end,:,:) = permute(k(:,1:4:end,:),[2 1 3]);
+Fin{2} = zeros([60 96 32]); Fin{2}(4:4:end,:,:) = permute(k(:,3:4:end,:),[2 1 3]);
+Fdpg2 = dpg_recon( Fin, N2, 2, 2);
+
diff --git a/tensor/tmult.m b/tensor/tmult.m
new file mode 100755
index 0000000..da2be49
--- /dev/null
+++ b/tensor/tmult.m
@@ -0,0 +1,68 @@
+function result = tmult(A,U,d)
+%
+% TMULT - tensor multiply (S = A x_i U) of A by U along dimension i
+%
+% usage: S = tmult( A, U, i )
+
+%
+% Copyright 2008  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+sz = size(A);
+tmp = U * local_unfold(A,d);
+
+sz(d) = size(U,1);
+N = length(sz);
+
+result = permute( reshape( tmp, sz( [ d:-1:1 N:-1:(d+1) ] ) ), ...
+                  [ d:-1:1 N:-1:(d+1) ] );
+
+%%%% the following was used with CFW's permute version of the unfolding
+%
+% result = permute(  ...
+%     reshape( tmp, sz([ d:N 1:(d-1) ]) ), ...
+%     [  (N-d+2):N 1:(N-d+1) ] );
+
+function result = local_unfold(A,d)
+%
+% UNFOLD - reduces the dimension of a tensor by "unfolding" along one direction
+%          (performs unfolding as described by deLathauwer in SIMAX 21(4):1253.
+%
+% Example: given a 3D tensor, A, 
+%
+% unfold(A,2) =>        I1
+%                     +----+ +----+ +----+
+%                  I2 |    | |    | |    |
+%                     |    | |    | |    |
+%                     +----+ +----+ +----+
+%                        \ ____|____ /
+%                              I3
+%
+% for a 3D tensor, 
+%
+% unfold(A,1) =  [ A(:,1,:)  A(:,2,:)  ... A(:,n2,:)  ],
+% unfold(A,2) = ( A(:,:,1)' A(:,:,2)' ... A(:,:,n3)' ),
+% unfold(A,3) = ( A(1,:,:)' A(2,:,:)' ... A(n1,:,:)' )
+%
+% To fold the matrix back (i.e. undo the unfold), one *must* know the
+% size of the final matrix and use the same permutation as the unfolding.  
+% An example:
+%    sz = size(A);
+%    A == permute( reshape( unfold(A,d) , sz([ d:-1:1 N:-1:(d+1) ]) ), ...
+%                  [  d:-1:1 N:-1:(d+1) ] );
+%
+
+
+sz = size(A);
+N = length(sz);
+if (d > length(sz)), return; end;
+
+result = reshape( permute( A, [ d:-1:1 N:-1:(d+1) ]), ...
+                  [ sz(d) prod(sz([1:(d-1) (d+1):N])) ] );
+
+% unfold(A,1) = reshape( permute( A, [ 1 3 2 ]), [sz(1) prod(sz(3:2)) ] )
+% unfold(A,2) = reshape( permute( A, [ 2 1 3 ]), [sz(2) prod(sz([1 3])) ] )
+% unfold(A,3) = reshape( permute( A, [ 3 2 1 ]), [sz(3) prod(sz([1 2])) ] )
+
diff --git a/tensor/trefold.m b/tensor/trefold.m
new file mode 100644
index 0000000..9c36d43
--- /dev/null
+++ b/tensor/trefold.m
@@ -0,0 +1,50 @@
+function result = trefold(tmp,sz,d)
+%
+% TREFOLD - tensor re-fold
+%           expands the dimension of a tensor by "refolding" along one direction
+%          (un-does the unfolding as described by deLathauwer in SIMAX 21(4):1253)
+%
+% Usage: A = trefold(unfA,sz,d)
+%
+% Example: given a 3D tensor, A, 
+%
+% unfold(A,2) =>        I1
+%                     +----+ +----+ +----+
+%                  I2 |    | |    | |    |
+%                     |    | |    | |    |
+%                     +----+ +----+ +----+
+%                        \ ____|____ /
+%                              I3
+%
+% for a 3D tensor, 
+%
+% unfold(A,1) =  [ A(:,1,:)  A(:,2,:)  ... A(:,n2,:)  ],
+% unfold(A,2) = ( A(:,:,1)' A(:,:,2)' ... A(:,:,n3)' ),
+% unfold(A,3) = ( A(1,:,:)' A(2,:,:)' ... A(n1,:,:)' )
+%
+% To fold the matrix back (i.e. undo the unfold), one *must* know the
+% size of the final matrix and use the same permutation as the unfolding.  
+% An example:
+%    sz = size(A);
+%    A == permute( reshape( unfold(A,d) , sz([ d:-1:1 N:-1:(d+1) ]) ), ...
+%                  [  d:-1:1 N:-1:(d+1) ] );
+%
+
+%
+% Copyright 2008  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+% sz = size(A);
+N = length(sz);
+if (d > length(sz)), result = A; return; end;
+
+result = permute( reshape( tmp, sz( [ d:-1:1 N:-1:(d+1) ] ) ), ...
+                  [ d:-1:1 N:-1:(d+1) ] );
+
+
+% unfold(A,1) = reshape( permute( A, [ 1 3 2 ]), [sz(1) prod(sz(3:2)) ] )
+% unfold(A,2) = reshape( permute( A, [ 2 1 3 ]), [sz(2) prod(sz([1 3])) ] )
+% unfold(A,3) = reshape( permute( A, [ 3 2 1 ]), [sz(3) prod(sz([1 2])) ] )
+
diff --git a/tensor/tunfold.m b/tensor/tunfold.m
new file mode 100644
index 0000000..dc41a91
--- /dev/null
+++ b/tensor/tunfold.m
@@ -0,0 +1,45 @@
+function result = tunfold(A,d)
+%
+% TUNFOLD - tensor unfold - reduces the dimension of a tensor by "unfolding" along one direction
+%          (performs unfolding as described by deLathauwer in SIMAX 21(4):1253.)
+%
+% Example: given a 3D tensor, A, 
+%
+% unfold(A,2) =>        I1
+%                     +----+ +----+ +----+
+%                  I2 |    | |    | |    |
+%                     |    | |    | |    |
+%                     +----+ +----+ +----+
+%                        \ ____|____ /
+%                              I3
+%
+% for a 3D tensor, 
+%
+% unfold(A,1) =  [ A(:,1,:)  A(:,2,:)  ... A(:,n2,:)  ],
+% unfold(A,2) = ( A(:,:,1)' A(:,:,2)' ... A(:,:,n3)' ),
+% unfold(A,3) = ( A(1,:,:)' A(2,:,:)' ... A(n1,:,:)' )
+%
+% To fold the matrix back (i.e. undo the unfold), one *must* know the
+% size of the final matrix and use the same permutation as the unfolding.  
+% An example:
+%    sz = size(A);
+%    A == permute( reshape( unfold(A,d) , sz([ d:-1:1 N:-1:(d+1) ]) ), ...
+%                  [  d:-1:1 N:-1:(d+1) ] );
+%
+
+% Copyright 2008  W. Scott Hoge  (wsh032580 at proton dot me)
+%
+% Licensed under the terms of the MIT License 
+% (https://opensource.org/licenses/MIT)
+
+sz = size(A);
+N = length(sz);
+if (d > length(sz)), result = A; return; end;
+
+result = reshape( permute( A, [ d:-1:1 N:-1:(d+1) ]), ...
+                  [ sz(d) prod(sz([1:(d-1) (d+1):N])) ] );
+
+% unfold(A,1) = reshape( permute( A, [ 1 3 2 ]), [sz(1) prod(sz(3:2)) ] )
+% unfold(A,2) = reshape( permute( A, [ 2 1 3 ]), [sz(2) prod(sz([1 3])) ] )
+% unfold(A,3) = reshape( permute( A, [ 3 2 1 ]), [sz(3) prod(sz([1 2])) ] )
+
diff --git a/vec.m b/vec.m
new file mode 100755
index 0000000..3c748e5
--- /dev/null
+++ b/vec.m
@@ -0,0 +1,2 @@
+function result = vec( x );
+result = x(:);
\ No newline at end of file