Yiwei_EFC.m

function efc = efc_masked(img, mask)
% EFC_MASKED  Entropy-Focus Criterion for 2D or 3D image
%
%   efc = efc_masked(img)
%   efc = efc_masked(img, mask)
%
% Inputs
%   img  : 2D or 3D image
%   mask : logical/binary mask, same size as img
%
% Output
%   efc  : normalized entropy-focus criterion
%          lower values indicate less ghosting / blurring

% -------------------------------------------------------------------------
% if no mask is given, the mask becomes the whole image with all voxels
% included
if nargin < 2 || isempty(mask)
    mask = true(size(img));
else
    mask = logical(mask);
end

xNorm = qc_normalize01(img, mask);
% use normalized magnitude inside mask only
x = xNorm(mask);             
% remove NaNs
x = x(~isnan(x));            

N = numel(x);
if N == 0
    efc = NaN;
    return;
end
%this is the energy of the masked intensities. It's what each voxel gets
%divided by, so the entropy is scale-independent.
bmax = sqrt(sum(x.^2));
if bmax == 0
    efc = 0;
    return;
end
%Energy-normalized intensities. Each p is a voxel's fraction of the total
%energy so the sum of p² = 1. It represents how concentrated the signal is
%in this voxel.
p = x / bmax;
small = 1e-16;

% same formula as MRIQC-style implementation
%Peaky, concentrated images give a less negative sum. 
% Spread-out, blurred/ghosted images give a more negative sum.
num = sum(p .* log(p + small));             %entropy sum
%What a maximally spread uniform image of N voxels would give and dividing
%by it makes EFC comparable across imaages of different sizes.
efc_max = N * (1/sqrt(N)) * log(1/sqrt(N)); %maximum-entropy normalizer
%entropy normalized by the max-entropy
%Lower value = sharper; higher = more blur/ghosting
efc = num / efc_max;                        %output

end

function x = qc_normalize01(img, mask)
% QC_NORMALIZE01  Convert magnitude image data to [0, 1] for QC metrics.
%
%   x = qc_normalize01(img)
%   x = qc_normalize01(img, mask)
%
% The min/max range is estimated inside mask when a non-empty mask is given.
% Non-finite output values are set to zero.
%double precision array
x = double(abs(img));
if nargin < 2 || isempty(mask)
    validMask = isfinite(x);
else
    validMask = logical(mask) & isfinite(x);
end

if ~any(validMask(:))
    x = zeros(size(x));
    return;
end
%smallest intensityinside a valid region
minValue = min(x(validMask));
%subtracts the offset from every voxel, shifting the whole distribution
%down before scaling.
x = x - minValue;
%largest intensity inside the valid region after the min is subtracted
maxValue = max(x(validMask));
if maxValue > 0
    x = x ./ maxValue;
else
    x(:) = 0;
end
x(~isfinite(x)) = 0;
x = min(max(x, 0), 1);
end

noMocoFile = 'I:\yannic.delisle\scripts\sub-10\baseline_nomoco\x_nomoco_1bin_mathilda_manual_finalNu240.mat'; 
mocoFile = 'I:\yannic.delisle\scripts\sub-10\sliding_L16.0_S3.2\x_moco_mathilda_smoothedMotion_from_x_cs_tres3.2s_Nu120x120x120_delta2_sliding.mat';
% stevaMocoFile = 'I:\yannic.delisle\scripts\sub-06\sliding_L3.5_S0.7\x_refined_moco_steva_delta0p1_from_x_cs_tres0.7s_Nu120x120x120_delta2_sliding.mat';
maskPath = 'I:\yannic.delisle\scripts\sub-10\sliding_L16.0_S3.2\prep_teva_tres3.2s_Nu120x120x120.mat';

S = load(noMocoFile);              
x = S.x_nomoco;  
T = load(mocoFile, 'x_moco');
moco = T.x_moco;
m = load(maskPath, 'maskBrain');
mask = m.maskBrain;
% U = load(stevaMocoFile, 'x_refined_moco');
% steva = U.x_refined_moco;

if iscell(x)
    x = x{1}; 
end 
%resize so the mask is the same size as the image
if ~isequal(size(mask), size(moco))
    mask = imresize3(single(mask), size(x), 'nearest') > 0.5;
end

format long
efc_masked(x, mask)