makeBrainMaskFromRecon.m

function [maskBrain, info] = makeBrainMaskFromRecon(x, thresholdScale)
%MAKEBRAINMASKFROMRECON Create a 3-D support mask from an MRI reconstruction.
%
% INPUTS
%   x              : complex or magnitude 3-D reconstruction
%   thresholdScale : multiplier applied to the Otsu threshold
%                    lower value -> larger mask
%                    higher value -> smaller mask
%
% OUTPUTS
%   maskBrain : logical 3-D mask with the same dimensions as x
%   info      : structure containing mask-generation diagnostics
%
% This creates a subject-specific signal-support mask. It is appropriate
% for signal normalization and can be used as a fixed evaluation ROI after
% visual verification. It is not a tissue-classification algorithm.

if nargin < 2
    thresholdScale = 0.75;
end

assert(ndims(x) == 3, ...
    'Input reconstruction must be a 3-D volume.');

%% 1. Magnitude and numerical cleanup

img = abs(single(x));
img(~isfinite(img)) = 0;

if ~any(img(:) > 0)
    error('The reconstruction contains no positive magnitude values.');
end

%% 2. Suppress isolated noise

imgSmooth = imgaussfilt3(img, 1);

%% 3. Robust intensity normalization

positiveValues = imgSmooth(imgSmooth > 0);
robustMaximum = prctile(positiveValues, 99.5);

if ~isfinite(robustMaximum) || robustMaximum <= 0
    error('Could not calculate a valid image-intensity scale.');
end

imgNormalized = imgSmooth / robustMaximum;
imgNormalized = min(max(imgNormalized, 0), 1);

%% 4. Automatic threshold

otsuLevel = graythresh(imgNormalized);

% Scaling below 1 makes the mask more inclusive, reducing the chance that
% lower-intensity brain regions are omitted.
threshold = thresholdScale * otsuLevel;

% Prevent an extremely low threshold from retaining most background noise.
threshold = max(threshold, 0.03);

maskInitial = imgNormalized > threshold;

%% 5. Three-dimensional morphological cleanup

% Join small gaps between neighbouring foreground regions.
maskInitial = imclose(maskInitial, strel('sphere', 2));

% Fill enclosed cavities.
maskInitial = imfill(maskInitial, 'holes');

% Remove small disconnected objects.
minimumObjectSize = max(100, round(0.0005 * numel(maskInitial)));

maskInitial = bwareaopen( ...
    maskInitial, minimumObjectSize, 26);

%% 6. Retain the largest connected component

components = bwconncomp(maskInitial, 26);

if components.NumObjects == 0
    error(['Mask generation produced no connected foreground object. ' ...
        'Try reducing thresholdScale.']);
end

componentSizes = cellfun(@numel, components.PixelIdxList);
[~, largestIndex] = max(componentSizes);

maskBrain = false(size(maskInitial));
maskBrain(components.PixelIdxList{largestIndex}) = true;

%% 7. Final cleanup

maskBrain = imclose(maskBrain, strel('sphere', 1));
maskBrain = imfill(maskBrain, 'holes');

%% 8. Diagnostics

info.thresholdScale = thresholdScale;
info.otsuLevel       = otsuLevel;
info.threshold       = threshold;
info.robustMaximum   = robustMaximum;
info.nVoxels         = nnz(maskBrain);
info.occupancy       = nnz(maskBrain) / numel(maskBrain);

end