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

   1 function [Sk,Sk2] = dpg_recon( Fin, N, dk, nky );
   2 % reconstruct measured EPI data using Dual-Polarity GRAPPA kernels
   3
   4 %
   5 % Copyright 2004-2016  W. Scott Hoge  (wsh032580 at proton dot me)
   6 %
   7 % Licensed under the terms of the MIT License
   8 % (https://opensource.org/licenses/MIT)
   9
  10 Sk  = zeros(size(Fin{1}));
  11 if nargout > 1,
  12   Sk2 = zeros(size(Fin{1}));
  13 end;
  14 sz = size(Sk);
  15
  16
  17 if nargin < 4,
  18   nky = 2;
  19 end;
  20 if isstr(nky),
  21   kernel = nky;
  22   nky = str2num(kernel(1));
  23 end;
  24
  25 if ( ~isfield( N, 'n') & isfield( N, 'coef' ) )
  26   for cnt=1:length(N),
  27     N(cnt).n = N(cnt).coef;
  28   end;
  29 end;
  30
  31 nkx = size(N(1).n,1) / size(Fin{1},3) / nky;
  32 ncoils = size(Sk,3);
  33
  34 % fprintf(' ky:%d kx:%d nc:%d\n',nky,nkx,ncoils);
  35 indx3 = ceil((0-nkx)/2) -1 + (1:nkx); % [-2:2];
  36
  37 verbose=1;
  38 prnt = 0;
  39
  40 % for kycnt=2:2*dk:(size(Sk,1)-2*dk),
  41 for kycnt=2:2*dk:(size(Sk,1)),
  42
  43   % fprintf('---%d ---\n',kycnt);
  44   for patcnt = 1:length(N); % [ 2:length(N)/2  (length(N)/2+(2:length(N)/2)) ];
  45
  46     % fprintf(' %d',patcnt)
  47     v = extrapolate_pattern(N(patcnt).pattern,nky);
  48
  49     if (patcnt <= length(N)/2)
  50       indx  = kycnt + (patcnt-1) + v;
  51     else
  52       indx  = kycnt + (patcnt-2) + v;
  53     end;
  54
  55     % fprintf(': '); fprintf(' %d', indx); fprintf('\n');
  56     if sum(indx>sz(1));
  57       indx = mod( indx, sz(1) );
  58       indx( indx==0 ) = sz(1);
  59     end;
  60     % fprintf('# '); fprintf(' %d', indx); fprintf('\n');
  61
  62     %% recon k-space line based on that pattern:
  63     A = zeros( size(Sk,2), length(indx3)*length(v)*size(Sk,3) );
  64
  65     for cnt3=1:length(indx3),
  66       y0 = mod( indx3(cnt3) + [1:size(Sk,2)] , size(Sk,2) );
  67       y0( y0 == 0 ) = size(Sk,2);
  68
  69       tmpSk = zeros([ length(indx) length(y0) size(Sk,3) ]);
  70
  71       if (patcnt <= length(N)/2)
  72         tmpSk(1:2:end,:,:) = Fin{1}( indx(1:2:end), y0, :);
  73         tmpSk(2:2:end,:,:) = Fin{2}( indx(2:2:end), y0, :);
  74       else
  75         tmpSk(1:2:end,:,:) = Fin{2}( indx(1:2:end), y0, :);
  76         tmpSk(2:2:end,:,:) = Fin{1}( indx(2:2:end), y0, :);
  77       end;
  78
  79       A( (1:size(Sk,2)), (cnt3-1)*nky*ncoils + (1:nky*ncoils) ) = ...
  80           reshape( permute( tmpSk, [ 2 1 3 ]), ...
  81                    size(Sk,2), size(Sk,3)*length(indx) );
  82
  83     end;
  84
  85     if (verbose), % ( sum(Sk(Y(mss(cnt2))+1,:,l)) ~= 0 );
  86       tmp = kycnt+patcnt;
  87       if (patcnt > length(N)/2); tmp = tmp-1; end;
  88       rptupd(prnt,sprintf(' %3d',[ tmp; indx(:) ]));
  89       rptupd(prnt,sprintf(' (pattern: %s)', N(patcnt).pattern) );
  90       rptupd(prnt,'\n');
  91       if (verbose>1),
  92         figure(1); imagesc( log10(abs(Sk(:,:,1)) + 1e-1) ); pause(0.5);
  93         keyboard;
  94       end;
  95     end;
  96     % linerpt(cnt2) = 'x';
  97
  98     n = N(patcnt).n;
  99
 100     for l=1:size(Sk,3),
 101       % for each coil ...
 102
 103       % if cnt2==109, keyboard; end;
 104       % rptupd(prnt,(' line %3d, pattern %3d\n', cnt2, plist(cnt3) );
 105
 106       if (nargout == 1),
 107         tmp = kycnt+(patcnt);
 108         if (patcnt > length(N)/2), tmp = tmp-1; end;
 109         Sk( tmp, :, l ) = [A] * n(:,l);
 110       else
 111         if (patcnt > length(N)/2 ),
 112           Sk2( kycnt+(patcnt-1), :, l ) = [A] * n(:,l);
 113         else
 114           Sk( kycnt+(patcnt), :, l ) = [A] * n(:,l);
 115         end;
 116       end;
 117     end;
 118
 119     % if (verbose>1), keyboard; end;
 120   end;
 121
 122 end;
 123
 124 Sk = Sk(1:sz(1),:,:);
 125
 126 %% need to fix this above, but for now...
 127 Sk = Sk([2:sz(1) 1],:,:);               % rotate the output back 1 kspace line
 128
 129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 130 function rptupd( prnt, str )
 131
 132 if (prnt)
 133   fprintf( str );
 134 end;
 135
 136 return;
 137
 138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 139 function v = extrapolate_pattern(p,YKS)
 140
 141 r = find( (p == 'x') | (p == '<' ) | (p == '>') );
 142 v = find( (p == '*') | (p == '<' ) | (p == '>') ) - r;
 143
 144 if length(v)>YKS,
 145   v = v( find( v~= 0 ) );
 146 end;
 147 return;