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)