57 lines
2.5 KiB
Matlab
57 lines
2.5 KiB
Matlab
function [motion_uw,dyn_range_uw,conf_uw,stdv]=omme(single_enc_motion,dyn_ranges)
|
|
|
|
%% Fast OPTIMAL MULTIPLE ENCODING RECONSTRUCTION (OMME-fast).
|
|
% Also leads to best resuts when phase-contrast measurements are NOT
|
|
% i.i.d., e.g. in case that reference phase measurement is measured only
|
|
% once (as in 4 points 4D flow, for instance)
|
|
%
|
|
% Inputs:
|
|
% single_enc_motion [Nvoxels,Ndynranges]: Example: vel in flow MRI. dphi*encEff/pi in MRE
|
|
% dyn_ranges [1,Ndynranges]: relative dynamic ranges for each single_enc_motion data. Example: vencs in flow MRI, 1/encEff in MRE.
|
|
%
|
|
% Outputs:
|
|
% motion_uw [Nvoxels,1]: unwrapped motion
|
|
% dyn_range_uw [1,1]: new dynamic range after unwrapping
|
|
% conf_uw [Nvoxels,1]: confidence image (should give 1 when measurements do
|
|
% not have noise, and small when the noise is large)
|
|
% stdv [1,1]: theoretical confidence in the estimated unwrapped motion
|
|
|
|
%% INIT
|
|
|
|
% Check correctness of input
|
|
if length(dyn_ranges)~=size(single_enc_motion,2)
|
|
error('Different number of dynamic ranges and encoded phases')
|
|
end
|
|
|
|
if any(dyn_ranges<=0)
|
|
error('Some dynamic ranges are smaller or equal to zero. Dynamic ranges need to be positive!')
|
|
end
|
|
|
|
N=size(single_enc_motion,1); % number of voxels to be unwrapped
|
|
motion_uw = zeros(N,1); % init unwrapped phases
|
|
conf_uw = zeros(N,1); % init confidence
|
|
|
|
%% Optimal Multiple Encoding algorithm
|
|
% New dynamic range (resulting from the combination)
|
|
dyn_range_uw = double(lcm(sym(abs(dyn_ranges))));
|
|
|
|
% Sampling of u, according only to the smallest dyn_range
|
|
[min_dynrange , ind] = min(abs(dyn_ranges)) ;
|
|
nb_of_samples = ceil(dyn_range_uw/min_dynrange)+2 ; % This is for covering the whole range [-dyn_range_uw,dyn_range_uw]
|
|
% Phase-contrast motion for smallest dyn_range +- min_dyn_range*k candidates
|
|
range_u = (-nb_of_samples:nb_of_samples)*min_dynrange;
|
|
u = single_enc_motion(:,ind)*ones(1,length(range_u)) + ones(N,1)*range_u ;
|
|
|
|
for k=1:N % loop over all voxels
|
|
Jmulti = 0 ;
|
|
u_k = u(k,abs(u(k,:))<=dyn_range_uw); % Candidates only in the effective dynamic range
|
|
for i=1:length(dyn_ranges)
|
|
Jmulti = Jmulti - cos( pi*( single_enc_motion(k,i) - u_k )/dyn_ranges(i) );
|
|
end
|
|
[~,ind_k] = min(Jmulti);
|
|
motion_uw(k) = u_k(ind_k);
|
|
% Evaluate confidence of the estimation by computing second derivative of the cost function
|
|
conf_uw(k) = sum( cos( pi*( single_enc_motion(k,:) - motion_uw(k) )./dyn_ranges )./(dyn_ranges.^2) )/sum( 1./(dyn_ranges.^2) ) ;
|
|
end
|
|
stdv = min_dynrange; % Theoretical standard deviation of the unwrapped estimator (= to empirical for inf realizations)
|