Skip to content

Commit

Permalink
bug fix for R2*
Browse files Browse the repository at this point in the history
  • Loading branch information
Kwok-shing Chan committed Oct 11, 2022
1 parent 4d3b589 commit 759e89f
Show file tree
Hide file tree
Showing 2 changed files with 95 additions and 0 deletions.
54 changes: 54 additions & 0 deletions utils/ARLO.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
%% function u = ARLO(y,x)
%
% Description: Fast monoexponential fitting by auto-regression
% Samples have to be evenly sampled (i.e. even spacing)
% At least 3 samples are needed
% Assuming time series in the last dimension
% ref: Pei et al. MRM 73:843-850(2015)
% e.g. for function y=exp(-x/u), u can be estimated by ARLO
%
% Kwok-shing Chan @ DCCN
% [email protected]
% Date created: 8 October, 2016
% Date last modified:
%
function u = ARLO(y,x)
% ensure y is real
% y=abs(y);

% get dimension of y
ndim = ndims(y);
if ndim==2 && size(y,2)==1
y = permute(y,[2 1]);
end
matrixSize = size(y);
N = matrixSize(end);

% reshape y s.t. new y = [all y, time]
y = reshape(y,[numel(y)/matrixSize(end)],N);

% get the spacing
dx = x(2)-x(1);

% get sum of signal for i and delta i
Si = zeros(size(y,1),N-2);
deltai = zeros(size(y,1),N-2);
for k=1:N-2
[Si(:,k), deltai(:,k)]= Simpson(y(:,k:k+2),dx);
end

% analytical solution for minimiser to obatain u
a = sum(Si.^2,2);
b = sum(Si.*deltai,2);
u = (a + (dx/3)*b)./((dx/3)*a + b);

% reshape u based on input dimension
u = reshape(u,[matrixSize(1:end-1)]);

end

% Quadratic approximation of Simpson rule's in 4th order accuracy when J=2
function [Si, deltai] = Simpson(y,dx)
Si = (dx/3) * (y(:,1) + 4*y(:,2) + y(:,3));
deltai = y(:,1)-y(:,3);
end
41 changes: 41 additions & 0 deletions utils/computeFiter.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
%% function fiter = computeFiter(s,shat,NUM_MAGN)
% s - measured signal
% shat - simulated signal
% NUM_GAGN - no. of phase corrupted echoes:
% NUM_MAGN=0 : complex fitting
% NUM_MAGN=length(s) : magnitude fitting
% NUM_MAGN (0,length(s)) : mixed fitting
%
% Description: Compute the fitter for lsqnonlin
%
% Kwok-shing Chan @ DCCN
% [email protected]
% Date created:
% Date last modified:
%
function fiter = computeFiter(s,shat,NUM_MAGN)
if NUM_MAGN == length(s) % Magnitude fitting
shat1 = abs(shat);
s1 = abs(s);
fiter = shat1(:) - s1(:);
elseif NUM_MAGN == 0 % Complex fitting
fiter2 = shat(:) - s(:);
fiter2 = [real(fiter2); imag(fiter2)];
% fiter2 = [real(fiter2), imag(fiter2)];
fiter = fiter2;
% fiter = abs(fiter2);
else
% Compute mixed fitting fit error
shat1 = abs(shat(1:NUM_MAGN));
s1 = abs(s(1:NUM_MAGN));
shat2 = shat(NUM_MAGN+1:end);
s2 = s(NUM_MAGN+1:end);

fiter1 = shat1(:) - s1(:);
fiter2 = shat2(:) - s2(:);
fiter2 = [real(fiter2); imag(fiter2)];

fiter = [fiter1;fiter2];
end
fiter = double(fiter);
end

0 comments on commit 759e89f

Please sign in to comment.