calculate_phase_mask_ellipse_3D.m

function [PHASE_MASK] = calculate_phase_mask_ellipse_3D(PHASE_QUALITY, KERNEL_RADIUS)
	% This function computes a quality-based mask of the phase angle plane of a cross correlation.
		
    % Rename the kernel radius
    rad = KERNEL_RADIUS;
    
    % Extract the central portion of the phase quality
    phase_quality_interior = PHASE_QUALITY(rad + 1 : end - rad - 1, rad + 1 : end - rad - 1, rad + 1 : end - rad - 1);
    
    phase_quality_sub = (phase_quality_interior - ...
        min(phase_quality_interior(:)));
    
    % Rescale the phase quality so that its range is [0, 1]
    phase_quality_scaled = phase_quality_sub ./ max(phase_quality_sub(:));
    
    % Measure size
	[region_height, region_width, region_slice] = size(phase_quality_scaled);
    
	% Make coordinates (Cartesian)
	[x, y, z] = meshgrid(1 : region_width, 1 : region_height, 1 : region_slice);
	
	% Centroid
	xc = region_width  / 2;
	yc = region_height / 2;
    zc = region_slice / 2;
	
	% Make coordinates (polar),
	% with origin (r = 0) at geometric centroid of the region
	[~, r] = cart2pol(x - xc, y - yc, z - zc);
	
	% Threshold the phase quality map using histogram equalization
% 	phase_quality_bw = (histeq(phase_quality_scaled, 2));
    phase_quality_bw = im2bw(phase_quality_scaled, 0.5);
	
% 	Set border pixels to 1 to prevent connecting regions via the border
	phase_quality_bw(1, :, :) = 1;
	phase_quality_bw(end, :, :) = 1;
	phase_quality_bw(:, 1, :) = 1;
	phase_quality_bw(:, end, :) = 1;
    phase_quality_bw(:, :, 1) = 1;
    phase_quality_bw(:, :, end) = 1;
	
	% Find properties of all the connected regions in the thresholded image. 
	phase_quality_region_props = regionprops(~phase_quality_bw, 'Centroid', 'PixelIdxList');
	
	% Count the number of regions
	num_regions = length(phase_quality_region_props);
	
	% Allocate a vector that will contain the median values
	% of the radial coordinates of each region.
	region_weighted_centroid_radial = zeros(num_regions, 1);
	
	% Allocate a vector to hold the centroids of the detected regions
	region_centroid_radial = zeros(num_regions, 1);
	
	% Allocate a vector to hold the mean values of the radial coordinates
	% of the pixels comprising the detected regions.
	region_radius_median = zeros(num_regions, 1);
	
	% Measure the median radial coordinate of each region
	for k = 1 : num_regions
		
		% Find the radial coordinate of the region's centroid.
		% whose value is equal to the k'th iteration of this loop
		region_centroid_subs = phase_quality_region_props(k).Centroid;
		
		% Radial coordinates of the pixels in the region
		region_radius_median(k) = median(r(phase_quality_region_props(k).PixelIdxList));
		
		% Convert to radial
		region_centroid_radial(k) = sqrt((region_centroid_subs(1) - xc)^2 + (region_centroid_subs(2) - yc)^2 + (region_centroid_subs(3) - zc)^2);
		
		% Centroid coordinates weighted by the median radial coordinate of the pixels in the region
		region_weighted_centroid_radial(k) = region_radius_median(k) * region_centroid_radial(k);
	
	end
		
	% Find the minimum
	[~, region_weighted_centroid_idx] = min(region_weighted_centroid_radial);
	
	% Mask indices
	mask_idx = phase_quality_region_props(region_weighted_centroid_idx).PixelIdxList;
	
	% Allocate a mask matrix
	phase_mask_holder = zeros(region_height, region_width, region_slice);
	
	% Apply the mask
	phase_mask_holder(mask_idx) = 1;
	
	% Zero the border pixels just in cae
	phase_mask_holder(:, 1, :) = 0 ;
	phase_mask_holder(:, end, :) = 0;
	phase_mask_holder(1, :, :) = 0;
	phase_mask_holder(end, :, :) = 0;
    phase_mask_holder(:, :, 1) = 0;
	phase_mask_holder(:, :, end) = 0;

	% Fill the holes
	phase_mask_holder = imfill(phase_mask_holder);
    
    % Get the region properties again
    phase_mask_region_props = regionprops(phase_mask_holder,...
        'PixelIdxList', 'MajorAxisLength', 'MinorAxisLength', ...
        'Orientation');
    
    % Read the measured major axis of the ellipse
    % that best fits the binary quality mask
    ax_maj = phase_mask_region_props.MajorAxisLength/2;
    
    % Same for the minor axis
    ax_min = phase_mask_region_props.MinorAxisLength/2;
    
    % Measure the orientation angle of
    % the best fit ellipse and convert it to radians.
    ax_angle = 1 * deg2rad(phase_mask_region_props.Orientation);
    
    % Rotate the coordinates for the elliptical Gaussian
    x2 = (x - xc) * cos(ax_angle) - (y - yc) * sin(ax_angle);
    y2 = (x - xc) * sin(ax_angle) + (y - yc) * cos(ax_angle);
    
    % Calculate the Gaussian function
    elliptical_mask = zeros(size(phase_mask_holder));
    elliptical_mask((x2.^2 / ax_maj^2 + y2.^2 / ax_min^2) < 1) = 1;
    
    % Insert the cropped mask into the phase
    % quality matrix, so that it's the same
    % size as the original phase matrix.
    PHASE_MASK = zeros(size(PHASE_QUALITY));
    PHASE_MASK(rad + 1 : end - rad - 1,...
        rad + 1 : end - rad - 1) = elliptical_mask;
  
end