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

   1 function [F,I,N] = dpg_cal( varargin );
   2
   3 % [F,I,N ] = dpg_cal( Wk, Sk, Y, [opt args]);
   4 %
   5 % opt args:
   6 %  'v', {0,1}              - verbosity off/on
   7 %  'kernel', string        - grappa kernel size, x by y. e.g. '2x5', '4x1'
   8 %  'N', (xy)-by-num_patterns-L matrix   - grappa reconstruction parameters
   9 %  'acs', array of indices - index locations to use for ACS lines
  10 %
  11 % This implementation is derived from
  12 % "Generalized autocalibrating partially parallel acquisitions"
  13 %   by M. A. Griswold, P. M. Jakob, et. al..
  14 %  Mag. Reson. Med, 47(6):1202-1210, Jun 2002
  15 %
  16
  17 %
  18 % Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
  19 %
  20 % Licensed under the terms of the MIT License
  21 % (https://opensource.org/licenses/MIT)
  22
  23 %
  24 % GRAPPA: data from all coils used to patch data in one coil
  25 %
  26 %  coil 1     o * o * a * o *
  27 %  coil 2     o * o * a * o *
  28 %  coil 3     o * o * a * o *
  29 %  coil 4     o * o * x * o *
  30 %
  31 %  *: acquired line, o: unacquired line, a: auto-callibration line
  32 %  x: reconstructed line
  33 %
  34 %  S_j (k_y - m\Delta k_y ) = \sum_{l=1}^L \sum_{b=0}^{N_b-1}
  35 %                                  n(j,b,l,m) S_l(k_y - b A \Delta k_y)
  36 %
  37
  38 Wk  = varargin{1};
  39 kin = varargin{2};
  40 Y   = varargin{3};
  41
  42 acs = []; N = []; verbose=0; kernel = '2x5'; dks = []; prnt = 0;
  43 eta = 1.0;
  44 chi = 0.0001;
  45
  46 vararg = varargin(4:end);
  47 for cnt=1:2:length(vararg),
  48   if strcmp( vararg{cnt}, 'v' ),
  49     verbose = cell2mat(vararg(cnt+1));
  50   elseif strcmp( vararg{cnt}, 'N' ),
  51     N = vararg{cnt+1};
  52   elseif strcmp( vararg{cnt}, 'dks' ),
  53     dks = vararg{cnt+1};
  54   elseif strcmp( vararg{cnt}, 'kernel' ),
  55     kernel = vararg{cnt+1};
  56   elseif strcmp( vararg{cnt}, 'acs' ),
  57     acs = vararg{cnt+1};
  58   elseif (strcmp( vararg{cnt}, 'report' ) || strcmp( vararg{cnt}, 'r' )),
  59     prnt = vararg{cnt+1};
  60   end;
  61 end;
  62
  63 if length(size(Wk)) == 3,
  64   [Wm,Wn,ncoils] = size(Wk);
  65 else,
  66   Wm = Wk(1); Wn = Wk(2); ncoils = Wk(3);
  67 end;
  68
  69
  70 %% zero-pad the k-space data
  71 Sk = zeros( Wm, Wn, ncoils );
  72 if iscell(kin)
  73   dualks = 1;
  74   s = kin{3};
  75 else
  76   dualks = 0;
  77   s = kin; clear kin;
  78 end;
  79 Sk(Y,:,:) = s;
  80
  81 rptupd(prnt,['* ' mfilename ': kernel size = ' kernel '\n']);
  82
  83 d = diff(Y(:));
  84 if ( isempty(acs) ),
  85   prfl = max( sum(abs( Sk(:,:,:)),3).' );
  86
  87   if length( d == 1) == 1,
  88     acs = find( prfl == max(prfl(Y(find(d==1)))) );
  89   else
  90     acs = Y([ min(find(d==1)); find( d == 1)+1 ]);
  91   end;
  92   rptupd(prnt,['* ' mfilename ': acs =']);
  93   rptupd(prnt,sprintf(' %d', acs) );
  94   rptupd(prnt,'\n');
  95 end;
  96
  97 if ( isempty(dks) ),
  98   %% determine which delta-k's to repair.
  99   dks = [];
 100   for cnt=max(d):-1:2,
 101     if (sum( cnt == d ) > 0),
 102       dks = [ cnt; dks ];
 103     end;
 104   end;
 105 end;
 106
 107 % set the kernel spacings based on 'kernel' string
 108 nky = str2double(kernel(1));
 109 nkx = str2double(kernel(3:end)); % should be odd
 110
 111 indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
 112
 113 pattern = {};
 114
 115 if ~isempty(N),
 116   for cnt=1:length(N)
 117     pattern{cnt} = N(cnt).pattern;
 118   end;
 119 end;
 120
 121 if isempty(pattern)
 122 %% foreach dk spacing
 123 for cnt=1:length(dks),
 124
 125   if ( dks(cnt) == 1 )
 126     % for no acceleration, just specify the pattern
 127     tmppat = repmat('*',[1 nky]);
 128     for patcnt = 1:nky,
 129       pattern{length(pattern)+1} = tmppat;
 130       pattern{length(pattern)}(patcnt) = '<';
 131     end;
 132
 133     for patcnt = 1:nky,
 134       pattern{length(pattern)+1} = tmppat;
 135       pattern{length(pattern)}(patcnt) = '>';
 136     end;
 137
 138     continue;
 139   end;
 140
 141   tmp = [ repmat([ '*' repmat('o',1,dks(cnt)-1) ],1,nky-1) '*' ];
 142   z2 = strfind(tmp,'o');
 143   if isempty(z2),
 144     % pad with zeros, if the kernel is 1-by-N
 145     tmp = [ 'o' tmp 'o' ];
 146     z2 = strfind(tmp,'o');
 147   end;
 148   [tmp2,indx] = sort( abs(z2 - z2( floor(median(1:length(z2))) )) );
 149
 150   olength = length(pattern);
 151   if (dualks == 1)
 152     % for dual-kernels, we want to genereate coefficients for the sampled
 153     % lines as well.
 154     npat = 2*length(tmp);
 155   else
 156     npat = length(z2);
 157   end
 158   for cnt2 = 1:npat, % (dks(cnt)-1), %
 159
 160     if (dualks == 0 )
 161       cnt3 = indx(cnt2);
 162       pattern{length(pattern)+1} = tmp;
 163       pattern{length(pattern)}(z2(cnt3)) = 'x';
 164     else
 165       cnt3 = mod( cnt2-1, length(tmp) ) + 1;
 166       pattern{length(pattern)+1} = tmp;
 167       if ( cnt2 <= npat/2 )
 168         pattern{length(pattern)}( cnt3 ) = '<';
 169       else
 170         pattern{length(pattern)}( cnt3 ) = '>';
 171       end;
 172     end;
 173
 174   end;
 175
 176 end;
 177 end;
 178
 179
 180 % find the recon patterns we can use
 181 plist = cell(length(pattern),1);
 182
 183 for cnt = 1:length(pattern),
 184   v = extrapolate_pattern(pattern{cnt},nky);
 185
 186   % if cnt == 14, keyboard; end;
 187
 188   for cnt2 = 1:length(acs);
 189     %% check if the GRAPPA pattern...
 190     indx = acs(cnt2) + v;
 191
 192     %% can be accomodated by the sampling pattern
 193     tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
 194     if ( length(tmp) == length(indx) ),
 195       plist{cnt} =  [ plist{cnt} acs(cnt2) ];
 196     end;
 197   end;
 198 end;
 199
 200 rptupd(prnt,['* ' mfilename ': sampling patterns in use = \n' ]);
 201 for cnt1=1:length(pattern);
 202   rptupd(prnt,sprintf('  %s: %d acs lines =', pattern{cnt1}, length(plist{cnt1})) );
 203   rptupd(prnt,sprintf(' %3d', plist{cnt1}) );
 204   % for cnt2 = 1:length( plist{cnt1} ),
 205   %   rptupd(prnt,('  %s', pattern{plist{cnt1}} );
 206   %   % v = extrapolate_pattern(pattern{plist{cnt1}(cnt2)},nky);
 207   %   % norm(vec(Sk( acs(cnt1)+v, :, : )));
 208   % end;
 209   rptupd(prnt,'\n');
 210 end;
 211 if (length(plist) == 0),
 212   error(['sorry, acquisition pattern is incompatible with available GRAPPA recon patterns']);
 213 end;
 214
 215 st = clock; % tic;
 216
 217 if verbose, rptupd(prnt,'\n'); end;
 218
 219 mss = zeros( size(Sk,1), 1 );         % set the missing lines index vector
 220 mss(Y) = 1;
 221
 222 if isempty(N),
 223   % N = zeros( length(indx3)*nky*size(s,3),length(pattern),size(s,3));
 224
 225   rptupd(prnt,'calc recon param\n');
 226   for cnt = 1:length(pattern),
 227     v = extrapolate_pattern(pattern{cnt},nky);
 228
 229     %% and for each pattern ...
 230     %% A2 = [];     b2 = [];
 231
 232     %%  A is data from the expanded input data matrix.
 233     %%     size(A) = [ kxnum * nacs, ncoils * nky * nkx ]
 234     %%
 235     %% in the matlab (and c-code) structure is:
 236     %%    for each element of indx3
 237     %%                   coil 1   |    coil 2   |    coil 3   | ...
 238     %%                -3 -1  1  3 | -3 -1  1  3 | -3 -1  1  3 | ...
 239     %%                 o  o  o  o |  o  o  o  o |  o  o  o  o | ...
 240     %%             ^
 241     %%  ACS(1) xres|
 242     %%             v
 243     %%             ^
 244     %%  ACS(2) xres|
 245     %%             v
 246     %%             ^
 247     %%  ACS(3) xres|
 248     %%             v
 249     %%
 250     %%  the older version of matlab code sorted the columns by
 251     %%
 252     %%     | c1 c2 c3 ... | c1 c2 c3 ... | c1 c2 c3 ...  | c1 c2 c3 ...
 253     %%     | -3 -3 -3 ... | -1 -1 -1 ... |  1  1  1 ...  |  3  3  3 ...
 254     %%     |  o  o  o ... |  o  o  o ... |  o  o  o ...  |  o  o  o ...
 255     %%
 256     %%  the order is immaterial, as long at it is consistent...
 257     %%  UPDATE (12 Jan 2011): the order does matter for GROG, so we have
 258     %%                        now changed it here so that the c-code and matlab code match.
 259     %%
 260     %%  The old way to switch between the two was:
 261     %%    Accode(:, permute( reshape(1:size(Accode,2),YKS,Ncoils,XKS), [3 2 1]) ) == Amatlab
 262     %%    Accode = Amatlab(:, permute( reshape(1:size(Accode,2),XKS,Ncoils,YKS), [3 2 1]))
 263     A = zeros( length(plist{cnt})*size(s,2), length(indx3)*nky*size(s,3) );
 264
 265     rptupd(prnt,sprintf('%3d  %s: ', cnt, pattern{cnt}));
 266
 267     % v = extrapolate_pattern(pattern{cnt});
 268     N(length(N)+1).pattern = pattern{cnt};
 269
 270     n = zeros( nky*size(s,3)*length(indx3),size(s,3));
 271
 272     for cnt2 = 1:length(plist{cnt}),
 273       %% and ACS line ...
 274
 275       %% if the sampling pattern accomodates the GRAPPA pattern...
 276       indx  = plist{cnt}(cnt2) + v;         % feedforward indices
 277       % indx2 = plist{cnt}(cnt2) + f;         % feedback indices
 278
 279       %% determine the k-space filling coefficients
 280       for cnt3=1:length(indx3),
 281         y0 = mod( indx3(cnt3) + [1:size(Sk,2)] , size(Sk,2) );
 282         y0( y0 == 0 ) = size(Sk,2);
 283
 284         % generate a hybrid source matrix
 285         if (dualks)
 286           tmpSksub = zeros( size( Sk(indx,y0,:) ) );
 287           if ( sum(pattern{cnt}=='<') ),
 288             tmpSksub( 1:2:length(indx), :, : ) = kin{1}( indx(1:2:end),y0,: );
 289             tmpSksub( 2:2:length(indx), :, : ) = kin{2}( indx(2:2:end),y0,: );
 290 %            tmpSksub( (v < 0), :, : )  = kin{2}( indx( (v < 0) ),y0,: );
 291 %            tmpSksub( (v >= 0), :, : ) = kin{3}( indx( (v >= 0) ),y0,: );
 292           else
 293             tmpSksub( 1:2:length(indx), :, : ) = kin{2}( indx(1:2:end),y0,: );
 294             tmpSksub( 2:2:length(indx), :, : ) = kin{1}( indx(2:2:end),y0,: );
 295 %            tmpSksub( (v >= 0), :, : ) = kin{2}( indx( (v >= 0) ),y0,: );
 296 %            tmpSksub( (v < 0), :, : )  = kin{3}( indx( (v < 0) ),y0,: );
 297           end
 298         else
 299           tmpSksub = Sk(indx,y0,:);
 300         end;
 301         % A( size(s,2)*(cnt2-1) + (1:size(s,2)), cnt3:length(indx3):size(A,2) ) = ...
 302         A( size(s,2)*(cnt2-1) + (1:size(s,2)), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
 303             reshape( permute( tmpSksub, [ 2 1 3 ]), size(Sk,2), size(Sk,3)*length(indx) );
 304
 305       end;
 306       %% for an N-by-1 kernel, these next lines are enough...
 307       % A2 = [ A2; ...
 308       %         reshape( permute(squeeze( Sk(indx,:,:) ),[ 2 3 1 ]), ...
 309       %                  size(s,2), size(s,3)*length(indx) )
 310       %       ];
 311       %
 312       % b2 = [b2; squeeze( Sk( plist{cnt}(cnt2), :, l ) ).' ];
 313     end;
 314
 315     % normalize, following to Algo 4.2 of M. Schneider, "GPGPU for Accelerated GRAPPA Autocalibration in
 316     %    Magnetic Resonance Imaging," Masters Thesis, U of Erlangen, Apr 2008
 317     % A_bak = A;
 318     p_src = sum(abs(A),2).^(-eta/2);    % compute the power in each line of the linsys
 319     A = tmult(A,diag(p_src),1);         % normalize by the power in each linsys line
 320
 321     % Apinv = pinv(A);
 322     % AA = A'*A;
 323     AA = [A ]'*[A ];
 324     At = [A ]';
 325     clear b2;
 326     for l=1:size(Sk,3),
 327       rptupd(prnt,'.');
 328       % A is the same for each coil ...
 329       b = vec( squeeze(Sk( plist{cnt}, :, l )).' ); % for dual-kernel, this defaults
 330                                                     % to the target data
 331
 332       b = b.*p_src(:);                  % normalize by the power in each linsys line
 333
 334       b2(:,l) = b;
 335       n(:,l) = A \ b;
 336
 337       if (verbose>2), pltcmplx( At'*n(:,l), b ); keyboard; end;
 338     end;
 339     b1 = b2;
 340     N(length(N)).n = n;
 341
 342     rptupd(prnt,'\n');
 343   end;
 344 end;
 345
 346 rptupd(prnt,'recon missing lines\n');
 347 %% now go recon the missing lines of k-space:
 348
 349 if (verbose); keyboard; end;
 350
 351 kypts = 1:length(mss);
 352 kypts( 1:floor(length(mss)/2) ) = kypts( (floor(length(mss)/2)):-1:1 );
 353
 354 linerpt = repmat(' ',1,length(mss));
 355
 356 for kycnt=1:length(mss),
 357   cnt2 = kypts(kycnt);
 358
 359   if verbose, rptupd(prnt,sprintf(' line %3d, from lines ',cnt2)); end;
 360
 361   if ( find( Y == cnt2) );
 362     % rptupd(prnt,'*');
 363     linerpt(cnt2) = '*';
 364     if verbose, rptupd(prnt,'\n'); end;
 365     continue;
 366   end;
 367
 368   tmp = [];
 369   %% determine which recon pattern we can use:
 370   for cnt4 = 1:length(pattern),
 371     v = extrapolate_pattern(pattern{cnt4},nky);
 372
 373     if ( (kycnt < length(mss)/2) && sum(find( pattern{cnt4} == '-')) ), continue, end;
 374     if ( (kycnt > length(mss)/2) && sum(find( pattern{cnt4} == '+')) ), continue, end;
 375
 376     indx  = cnt2 + v;
 377     % indx2 = cnt2 + f;
 378     %     tmp = find( Y(:)*ones(1,length(indx)) == ones(length(Y),1)*indx );
 379     %     if  ( ... % ( (min(indx) > 1) & (max(indx) < size(Sk,1)) ) & ...
 380     %         ( length(tmp) == length(indx) ) ),
 381     %       break;                          % jump out of the loop
 382     %     end;
 383     if ( indx(1) < 1) | (indx(end) > size(Sk,1)); continue; end;
 384
 385     tmp = (mss( indx ) == ones(length(indx),1));
 386     if ( sum(tmp) == length(tmp)  ),
 387       break;
 388     end;
 389
 390   end;
 391
 392   if ( ( (min(indx) < 1) | (max(indx) > size(Sk,1)) ) | ...
 393        ( sum(tmp) ~= length(indx) ) ),
 394     %
 395     linerpt(cnt2) = 'o';
 396     if verbose, rptupd(prnt,'o'); rptupd(prnt,'\n'); end;
 397     continue;
 398   end;
 399
 400   %% recon k-space line based on that pattern:
 401   rsz = size(Sk,2);
 402   A = zeros( rsz, length(indx3)*length(v)*size(s,3) );
 403   B = []; % zeros( rsz, length(indx3)*length(f)*size(s,3) );
 404
 405   for cnt3=1:length(indx3),
 406     y0 = mod( indx3(cnt3) + [(size(Sk,2)-rsz)/2+(1:rsz)] , size(Sk,2) );
 407     y0( y0 == 0 ) = size(Sk,2);
 408
 409     A( (1:rsz), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
 410         reshape( permute( Sk(indx,y0,:), [ 2 1 3 ]), ...
 411                  rsz, size(Sk,3)*length(indx) );
 412   end;
 413
 414   if (verbose), % ( sum(Sk(Y(mss(cnt2))+1,:,l)) ~= 0 );
 415     rptupd(prnt,sprintf(' %3d',[ indx(:) ]));
 416     rptupd(prnt,sprintf(' (pattern: %s)', pattern{cnt4}) );
 417     rptupd(prnt,'\n');
 418     if (verbose>1),
 419       figure(1); imagesc( log10(abs(Sk(:,:,1)) + 1e-1) ); pause(0.5);
 420       keyboard;
 421     end;
 422   end;
 423   linerpt(cnt2) = 'x';
 424
 425   for cnt3=1:length(N),
 426     if strcmp( N(cnt3).pattern, pattern{cnt4} ),
 427       break;                            % jump out of the loop
 428     end;
 429   end;
 430   n = N(cnt3).n;
 431
 432   for l=1:size(Sk,3),
 433     % for each coil ...
 434
 435     % if cnt2==109, keyboard; end;
 436     % rptupd(prnt,(' line %3d, pattern %3d\n', cnt2, plist(cnt3) );
 437     Sk( cnt2, :, l ) = [A B] * n(:,l);
 438   end;
 439
 440   if (verbose>1), keyboard; end;
 441 end;
 442 rptupd(prnt,linerpt);
 443 rptupd(prnt,'  \n');
 444
 445 % Sk(1:2:size(Sk,1),:,:) = -Sk(1:2:size(Sk,1),:,:);
 446 % Sk(:,1:2:size(Sk,2),:) = -Sk(:,1:2:size(Sk,2),:);
 447
 448 F = Sk;
 449 I = sqrt( sum(abs(ifft2(Sk)).^2,3) );
 450 I = I./(prod(size(I)));
 451 mask = 0;
 452
 453 rpttime = etime(clock,st);
 454 rptupd(prnt,sprintf(' GRAPPA recon time: %g (avg per kernel: %g)\n',  ...
 455         rpttime, rpttime/length(pattern)) );
 456
 457 return;
 458
 459 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 460 function rptupd( prnt, str )
 461
 462 if (prnt)
 463   fprintf( str );
 464 end;
 465
 466 return;
 467
 468 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 469 function v = extrapolate_pattern(p,YKS)
 470
 471 r = find( (p == 'x') | (p == '<') | (p == '>') );
 472 v = find( (p == '*') | (p == '<') | (p == '>') ) - r;
 473
 474 if length(v)>YKS,
 475   v = v( find( v~= 0 ) );
 476 end;
 477 return;
 478
 479