--- /dev/null
+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.
--- /dev/null
+
+# 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)
+
--- /dev/null
+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;
+
+
--- /dev/null
+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;
--- /dev/null
+function a = fif(a)
+
+a = fftshift(fftshift(ifft2(fftshift(fftshift( a, 1),2)),1),2);
--- /dev/null
+function a = ifi(a)
+a = ifftshift(ifftshift(fft2(ifftshift( ifftshift( a, 1),2)),1),2);
--- /dev/null
+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;
+
--- /dev/null
+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
--- /dev/null
+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);
+
--- /dev/null
+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;
--- /dev/null
+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' ]);
--- /dev/null
+% 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);
+
--- /dev/null
+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])) ] )
+
--- /dev/null
+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])) ] )
+
--- /dev/null
+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])) ] )
+
--- /dev/null
+function result = vec( x );\r
+result = x(:);
\ No newline at end of file
projects / dpg / commitdiff