sigprocexperts.com Git - dpg/blob - recongrappa.m

   1 function [F,I,mask] = recongrappa( varargin );
   2 % [F,I,mask ] = recongrappa(Wk, Sk, Y, [opt args]);
   3 %
   4 % opt args:
   5 %  'v', {0,1}              - verbosity off/on
   6 %  'kernel', string        - grappa kernel size, x by y. e.g. '2x5', '4x1'
   7 %  'N', (xy)-by-L matrix   - grappa reconstruction parameters
   8 %  'acs', array of indices - index locations to use for ACS lines
   9 %  'dk', delta-k           - distance between phase encode lines used by
  10 %                            the GRAPPA kernel.  default: 2
  11 %
  12 % This is an implementation of
  13 % "Generalized autocalibrating partially parallel acquisitions"
  14 %   by M. A. Griswold, P. M. Jakob, et. al..
  15 %  Mag. Reson. Med.   47(6):1202-1210.  Jun 2002
  16
  17 % Copyright 2004-2007  W. Scott Hoge  (wsh032580 at proton dot me)
  18 %
  19 % Licensed under the terms of the MIT License
  20 % (https://opensource.org/licenses/MIT)
  21
  22 %
  23 % uncombined images are generated for each coil in the array, by
  24 % filling in k-space lines.
  25 %
  26 %
  27 % GRAPPA: data from all coils used to patch data in one coil
  28 %
  29 %  coil 1     o * o * - * o *
  30 %  coil 2     o * o * - * o *
  31 %  coil 3     o * o * - * o *
  32 %  coil 4     o * o * x * o *
  33 %
  34 %  *: acquired line, o: unacquired line, -: auto-callibration line
  35 %  x: reconstructed line
  36 %
  37 %  S_j (k_y - m\Delta k_y ) = \sum_{l=1}^L \sum_{b=0}^{N_b-1}
  38 %                                  n(j,b,l,m) S_l(k_y - b A \Delta k_y)
  39 %
  40 % in a variable density sampling scheme, the center of k-space is fully
  41 % sampled and can be used to determine the linear fit parameters.
  42 %
  43 % A   - acceleration factor
  44 % N_b - number of blocks used in the reconstruction
  45 %
  46 % N_b = 1 ----> GRAPPA = VD-AUTO-SMASH
  47
  48 global prnt
  49
  50 Wk = varargin{1};
  51 s  = varargin{2};
  52 Y  = varargin{3};
  53
  54 acs = []; N = []; verbose=0; kernel = '2x3'; dk = 2;
  55 prnt = 0;   % set to '1' to output some debugging information
  56 Nbias = [];
  57 chi = 1e-8;
  58
  59 vararg = varargin(4:end);
  60 for cnt=1:2:length(vararg),
  61   if strcmp( vararg{cnt}, 'v' ),
  62     verbose = cell2mat(vararg(cnt+1));
  63   elseif strcmp( vararg{cnt}, 'N' ),
  64     N = vararg{cnt+1};
  65   elseif strcmp( vararg{cnt}, 'Nbias' ),
  66     Nbias = vararg{cnt+1};
  67   elseif strcmp( vararg{cnt}, 'kernel' ),
  68     kernel = vararg{cnt+1};
  69   elseif strcmp( vararg{cnt}, 'dk' ),
  70     dk = vararg{cnt+1};
  71   elseif strcmp( vararg{cnt}, 'chi' ),
  72     chi = vararg{cnt+1};
  73   elseif strcmp( vararg{cnt}, 'acs' ),
  74     acs = vararg{cnt+1};
  75   elseif (strcmp( vararg{cnt}, 'report' ) | strcmp( vararg{cnt}, 'r' )),
  76     prnt = vararg{cnt+1};
  77   end;
  78 end;
  79
  80 if isempty(acs),
  81   d = diff(Y);
  82   acs = Y(find(d == 1));
  83   % acs = acs(2:length(acs));
  84 end;
  85
  86 if length(size(Wk)) == 3,
  87   [Wm,Wn,Wc] = size(Wk);
  88 else,
  89   Wm = Wk(1); Wn = Wk(2); Wc = Wk(3);
  90 end;
  91
  92 %% zero-pad the k-space data
  93 Sk = zeros( Wm, Wn, Wc );
  94 Sk(Y,:,:) = s;
  95
  96 % set the kernel spacings based on 'kernel' string
  97 nky = str2double(kernel(1));               % should be even
  98 if (nky==1),
  99   indx0 = 1;
 100 else,
 101   indx0 = floor( [-(nky-1):2:(nky-1)] * dk/2 );
 102 end;
 103 % if (nky == '2'), indx0 = [-1 1];
 104 % elseif (nky == '6'), indx0 = [ -5 -3 -1 1 3 5];
 105 % else % (nky == '4'), indx0 = [ -3 -1 1 3];
 106 % end;
 107
 108 nkx = str2double(kernel(3));               % should be odd
 109 indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
 110
 111 rptupd( sprintf(['kernel size: ' num2str(length(indx0)) 'x' num2str(length(indx3)) ';']) );
 112 rptupd( sprintf([' indx0 = [' sprintf(' %d',indx0) ']; ']) );
 113 rptupd( sprintf([' indx3 = [' sprintf(' %d',indx3) '];\n']) );
 114
 115 stime = clock;
 116
 117 plist = zeros( length(acs), 1 );
 118
 119 for cnt2 = 1:length(acs);
 120   %% check if the GRAPPA pattern...
 121   indx = acs(cnt2) + indx0;
 122
 123   %% can be accomodated by the sampling pattern
 124   tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
 125   plist(cnt2) =  ( length(tmp) == length(indx) );
 126 end;
 127 %% use only those ACS lines that fit the GRAPPA kernel
 128 acs = vec( acs(logical(plist)) ).';
 129 rptupd([ ' ACS lines:' sprintf(' %d ',acs) '\n' ]);
 130
 131 rptupd('processing coil: ');
 132 for l=1:size(s,3),
 133   rptupd([' ' num2str(l)]);
 134   %% determine the k-space filling coefficients for each coil:
 135
 136   %%% setup the system matrix
 137   A2 = zeros( size(s,2), length(indx3)*length(indx0)*size(s,3) );
 138
 139   if size(N,2) < l,                     % if coefficients are not supplied
 140   %%% Griswolds way:
 141   %   tic;
 142   %   A2 = zeros( size(s,2), length(indx3)*length(indx)*size(s,3) );% 20*size(s,3) );
 143   %   for x0=1:size(s,2),
 144   %     indx2 = x0 + indx3;
 145   %     indx2( indx2 < 1 ) = indx2( indx2 < 1 ) + size(s,2);
 146   %     indx2( indx2 > size(s,2) ) = indx2( indx2 > size(s,2) ) - size(s,2);
 147   %
 148   %     A2(x0,:) = vec( permute(squeeze( Sk(indx,indx2,:)),[2 3 1 ]) ).';
 149   %   end;
 150   %   toc
 151
 152   %%% My way: (heh, 8x faster)
 153   % tic;
 154   A = zeros( length(acs)*size(s,2), length(indx3)*length(indx0)*size(s,3) );
 155   A2 = zeros( size(s,2), length(indx3)*length(indx0)*size(s,3) );
 156   if (l==1), rptupd(['\n size(A) = ' num2str(size(A)) '\n']), end;
 157
 158   for cnt2=1:length(acs),
 159     indx = acs(cnt2) + indx0;
 160     for cnt=1:length(indx3),
 161       y0 = mod( indx3(cnt) + [1:size(Sk,2)] , size(Sk,2) );
 162       y0( y0 == 0 ) = size(Sk,2);
 163
 164       A( size(s,2)*(cnt2-1) + (1:size(s,2)), ...
 165          cnt:length(indx3):size(A,2) ) = ...
 166           reshape( permute(squeeze( Sk(indx,y0,:) ),[ 2 3 1 ]), ...
 167                    size(Sk,2), size(Sk,3)*length(indx) );
 168     end;
 169   end;
 170   % toc
 171
 172   %%% setup the objective vector
 173   b = vec( Sk( acs, :, l ).' );
 174
 175   %%% find the coefficients
 176   % n = A \ b; % pinv( A, norm(A)*0.01 ) * b;
 177
 178   if (size(Nbias,1)==0), % isempty('Nbias'),
 179     % n = cgsolv( A'*A, A'*b, zeros(size(A,2),1), size(A,2) );
 180     [ n tt ] = bicg( A'*A + chi*eye(size(A,2)), A'*b );
 181   else,
 182     n = cgsolv( A'*A, A'*b, Nbias(:,l), size(A,2)/4 );
 183   end;
 184
 185   N(:,l) = n;                           % store them
 186   if verbose, pltcmplx( A*n, b ); keyboard; end;
 187   end;
 188
 189   % extract the coefficients for this coil
 190   n = N(:,l);
 191
 192   %% now go recon the missing lines of k-space:
 193   mss = zeros(Wm,1);            % set the missing lines index vector
 194   mss(Y) = 1;
 195
 196   for cnt0=1:length(mss),
 197     if (mss(cnt0) == 1);   continue; end;
 198
 199     indx = cnt0 + indx0;
 200     if ( (min(indx) < 1) | (max(indx) > size(Sk,1))); continue; end;
 201
 202     if ( sum(sum( Y(:)*ones(size(indx)) == ones(size(Y(:)))*indx )) < length(indx) ),
 203       continue;
 204     end;
 205
 206     %     for x0=1:size(s,2),
 207     %       indx2 = x0 ; % + [-2:2];
 208     %       indx2( indx2 < 1 ) = indx2( indx2 < 1 ) + size(s,2);
 209     %       indx2( indx2 > size(s,2) ) = indx2( indx2 > size(s,2) ) - size(s,2);
 210     %
 211     %     end;
 212     for cnt=1:length(indx3),
 213       y0 = mod( indx3(cnt) + [1:size(Sk,2)] , size(Sk,2) );
 214       y0( y0 == 0 ) = size(Sk,2);
 215
 216       A2( (1:size(s,2)), cnt:length(indx3):size(A2,2) ) = ...
 217           reshape( permute(squeeze( Sk(indx,y0,:) ),[ 2 3 1 ]), ...
 218                    size(Sk,2), size(Sk,3)*length(indx) );
 219     end;
 220     Sk( cnt0, :, l ) = A2 * n;
 221     % figure(1); imagesc( log10(abs(Sk(:,:,l)) + 1e-1) ); pause(0.5);
 222   end;
 223
 224 end;
 225 rptupd('\n');
 226 rptupd( sprintf( 'Elapsed time: %8.5g seconds.\n', etime(clock,stime) ) );
 227
 228 F = Sk;
 229 I = sqrt( sum(abs(ifft2(Sk)).^2,3) );
 230 mask = N;
 231
 232 clear global prnt
 233
 234 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 235 function rptupd( str )
 236 global prnt
 237
 238 if (prnt)
 239   fprintf( str );
 240 end;
 241
 242 return;