calculate_phase_quality_3D.m

function phase_quality_array = calculate_phase_quality_3D(wrapped_phase_angle, kernel_radius)


% Calculate the horizontal and vertical wrapped phase differences
% phase_diff_rows = phase_diff_kernel_wrapped(wrapped_phase_angle, 1, 'sobel');
% phase_diff_cols = phase_diff_kernel_wrapped(wrapped_phase_angle, 2, 'sobel');

phase_diff_rows = wrapped_phase_difference(wrapped_phase_angle, 1);
phase_diff_cols = wrapped_phase_difference(wrapped_phase_angle, 2);
phase_diff_slice = wrapped_phase_difference(wrapped_phase_angle, 3);

% Measure the array size
% and subtract 2 from each dimensions due to the differencing 
[M, N, O] = size(wrapped_phase_angle);

% Array height and width
array_height = M;
array_width  = N;
array_slice = O;

% Min and maxes of rows and columns
col_min = 1 + kernel_radius;
col_max = array_width - kernel_radius;

row_min = 1 + kernel_radius;
row_max = array_height - kernel_radius;

slice_min = 1 + kernel_radius;
slice_max = array_slice - kernel_radius;

% Allocate the phase quality array
phase_quality_array = zeros(array_height, array_width, array_slice);

% Loop over the rows
tic

for m = row_min : row_max
	for n = col_min : col_max
        for o = slice_min : slice_max
		
            % Extract the sub-region of phase differences
            row_diffs = phase_diff_rows(m - kernel_radius : m + kernel_radius,...
             n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius);

            col_diffs = phase_diff_cols(m - kernel_radius : m + kernel_radius,...
             n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius);

            slice_diffs = phase_diff_slice(m - kernel_radius : m + kernel_radius,...
             n - kernel_radius : n + kernel_radius, o - kernel_radius : o + kernel_radius);
            % Standard deviations
            % The 1 indicates normalization by N, not by N-1
            row_diff_std = std(row_diffs(:), 1);
            col_diff_std = std(col_diffs(:), 1);
            slice_diff_std = std(slice_diffs(:), 1);

            % Populate the phase quality array
            phase_quality_array(m, n, o) = row_diff_std + col_diff_std + slice_diff_std;
        end
		
	end
end

toc
	
end