function [freq_good, varargout] = PSF_classification(PSF, varargin)
% 
% PSF_CLASSIFICATION Determine frequencies at which PSF quality is 
% sufficient for use in image reconstruction.
%
%   [freq_good] = PSF_classification(PSF) returns the frequency indices of 
%   PSF that has CNR and/or SNR values above user threshold(s). The default
%   mode evaluates complex magnitude CNR with a threshold of 2 (linear 
%   ratio; 3 on dB scale).
%
%   INPUTS:         PSF:            PSF dataset, normalized (x,y,f,z)    
%
%   OUTPUTS:        freq_good:      Structure returning indices of good frequencies for a particular calculation. Field will be blank if calculation was not evaluated by the algorithm.
%                                           ie. if phase CNR was chosen to be evaluated, use freq_good.CNR_phase_i to access the good frequencies
%
%   ------------------------------------------------------------------
%
%   PSF_CLASSIFICATION uses the cluster ROI method by default. The user can 
%   use the following optional inputs, written in a 'key', value pairing 
%   during the function call:
%
%   OPT. INPUTS:    'CNR_thresholds', [a b c d]   Select CNR calculations to execute by assigning linear thresholds.
%                                                       a = magnitude, b = phase, c = average real/imag, d = complex magnitude 
%                                                       ie. to calculate magnitude CNR only, feed in 'CNR_thresholds', [a 0 0 0] where a is the integer threshold for magnitude CNR     
%                   'SNR_thresholds', [e f g h]   Select SNR calculations to execute by assigning linear thresholds.
%                                                       e = magnitude, f = phase, g = average real/imag, h = complex magnitude 
%                   'maxima_ROI', 1:              Select maxima method to estimate ROI; default is cluster method.
%                   'ROI_fixed_size', [x,y]:      Size of fixed ROI [x,y] in mm. Default is 10x10 (scattering probe size).
%                   'peak_signal', 1:             Select peak signal method of evaluating CNR/SNR; default is mean method.
%                   'plot_ROI_summation', 1:      Plot automated ROI against frequency summation.
%                   'cluster_threshold', i:       Select threshold cutoff for determining cluster ROI. Default is 3dB. 
%                                                       i = threshold in dB ie. i = 3.5dB, taking signal above 55% of baseline mininum as potential cluster points.
%                   'plot_CNR_SNR', 1:            Plot CNR and/or SNR against frequency index (default - see 'ftable').
%                   'dB_scale', 0:                Plot CNR and/or SNR in linear scale. Default is dB scale (10*log10(ratio)).
%                   'ftable', ftable:             Plot CNR and/or SNR against sampled frequencies. 
%                                                       ftable = vector of sampled frequencies (in GHz) ie. 3.0:0.1:8.0 
%
%   ------------------------------------------------------------------
%
%   OPT. OUPUTS:    ROI_mask:           Cell matrix containing logical ROI mask [f,z] 
%                   exclusion_zone:     Structure containing fields ...
%                                           1) 'ez_mask' - cell matrix containing logical exclusion zone mask [f,z] 
%                                           2) 'center_coords' - matrix containing x,y coordinates of center coordinates of circle
%                                           3) 'radius' - matrix containing radius of each circle [f,z]
%                   PSF_inverted:       Initial magnitude PSF dataset with inversion performed  
%
%   ------------------------------------------------------------------
%
%   Example 1
%   ---------
%   Evaluate phase SNR (threshold = 3dB) and complex SNR (threshold = 2dB) with ROI cluster method. Plot results on dB scale.
%       [freq_good] = PSF_classification(PSF, 'plot_ROI_summation', 1, 'plot_CNR_SNR', 1, 'SNR_thresholds', [0 2 0 1.5849], 'dB_scale', 1);
%       * Note: will also evaluate complex magnitude CNR (threshold = 3dB)
%               because user did not override
%
%   Example 2
%   ---------
%   Evaluate real/imag CNR (linear threshold = 5) with ROI maxima method and specific ROI size. Retrieve optional outputs.
%       [freq_good, ROI_mask, exclusion_zone, PSF_inversion] = PSF_classification(PSF, 'maxima_ROI', 1, 'probe_size', [20,10], 'CNR_thresholds', [0 0 5 0]);
%
%
%   ------------------------------------------------------------------
%
%   Written by Jessica Trac, McMaster University, EMVi Lab
%         Adapted from Justin McCombe's SNR GUI code
%         Directory must also have phase_unwrap.m
%
%   Created: May 3 2018
%   Last Revision:  Auust 24 2018
%
%   ------------------------------------------------------------------
%

%% (0) SETTING OPTIONAL PARAMETERS 
% ---------------------------------------------------------------------- 
% From user input, where an optional parameter can be set using a 'key',
% value pairing (ie. to evaluate magnitude SNR, pass in 'SNR_mag', x where
% x is the desired user-threshold for magnitude SNR classification).
% ---------------------------------------------------------------------- 
 
[nx, ny, nf, nz] = size(PSF);

fprintf('Classifying PSFs... \n');

var_input = struct('CNR_thresholds', [0 0 0 2], 'SNR_thresholds', [0 0 0 0], 'maxima_ROI', 0, 'ROI_fixed_size', [10,10],'plot_ROI_summation', 0, 'cluster_threshold', 3, 'peak_signal', 0, 'plot_CNR_SNR', 0, 'dB_scale', 1, 'ftable', 0); % Default values for each optional parameter
var_input_names = fieldnames(var_input); % Optional parameter names

for pair = reshape(varargin, 2, []) % 1x2 cell of parameter name, value
    var_name = pair{1}; % Parameter name, as specified by user in function call
    var_value = pair{2}; % Parameter value, as specified by user in function call
    if any(strcmp(var_name, var_input_names))
        var_input.(var_name) = var_value; % Override default value
    else
        error('%s is not a recognized parameter name - please re-enter', var_name); 
    end
end

rho_mag_CNR = var_input.CNR_thresholds(1); % Assign user-defined thresholds; 0 if not chosen for execution, so 0 cannot be used 
rho_phase_CNR = var_input.CNR_thresholds(2);
rho_real_imag_CNR = var_input.CNR_thresholds(3);  
rho_complex_mag_CNR = var_input.CNR_thresholds(4); % Used as default calculation mode if no user-input; linear threshold of 2 - user must override with input of 0, if not desired

rho_mag_SNR = var_input.SNR_thresholds(1); % Assign user-defined thresholds; 0 if not chosen for execution, so 0 cannot be used
rho_phase_SNR = var_input.SNR_thresholds(2);
rho_real_imag_SNR = var_input.SNR_thresholds(3);  
rho_complex_mag_SNR = var_input.SNR_thresholds(4); 

cluster_threshold = var_input.cluster_threshold; 

if var_input.ftable == 0 
    ftable = linspace(1,nf,nf); % Plot frequencies as general indices (ie. 1,2,3...etc)
    f_plot_title = 'Frequency Index'; 
else
    ftable = var_input.ftable; % Plot frequencies as actual sampling frequencies
    f_plot_title = 'Frequency (GHz)';
end

flag_maxima = 0; % Flag to track algorithm switch from cluster ROI to maxima ROI due to noisy dataset
center_xy = [round(nx/2), round(ny/2)]; % Center-point of image in x,y (x = columns, y = rows)

%% (1) INVERSION
% ---------------------------------------------------------------------- 
% Inversion is performed if mean(ROI) < mean(background). Here, mean(ROI) 
% is evaluated as the center 10% of the image as an initial estimate.
% ---------------------------------------------------------------------- 

padding = ceil(nx*0.05); % Take 5% of image range to set initial PSF estimate location to center 10% of image (ie. center point is [50,50] for 100x100 image - will take [45:55, 45:55] as estimate of PSF location)

PSF_init_col = (center_xy(1,1) - padding):(center_xy(1,1) + padding); % The 'x-direction' is an index for columns
PSF_init_row = (center_xy(1,2) - padding):(center_xy(1,2) + padding); % The 'y-direction' is an index for rows
PSF_init_location = zeros(nx, ny);
PSF_init_location(PSF_init_row,PSF_init_col) = 1; % Initial PSF estimate, represented by 1s (center of image + user-defined padding)

PSF_inverted = zeros(nx,ny,nf,nz); % New dataset with inversion considerations

for f = 1:nf
    for z = 1:nz

        temp_PSF_mag = abs(PSF(:,:,f,z)); % Magnitude of PSF at specific frequency and z-plane
        
        temp_PSF_mag_ROI = NaN(nx,ny); 
        temp_PSF_mag_ROI(PSF_init_location == 1) = temp_PSF_mag(PSF_init_location == 1); % Populate ROI with signal value, rest of matrix is NaN

        temp_PSF_exclusion = temp_PSF_mag; 
        temp_PSF_exclusion(PSF_init_location == 1) = NaN; % Populate exclusion zone with signal value in non-ROI region, ROI is NaN

        if mean2(temp_PSF_mag_ROI(~isnan(temp_PSF_mag_ROI))) < mean2(temp_PSF_exclusion(~isnan(temp_PSF_exclusion))) 
            PSF_inverted(:,:,f,z) = (abs(PSF(:,:,f,z))-1)*(-1); % Perform inversion on magnitude of PSF
        else
            PSF_inverted(:,:,f,z) = abs(PSF(:,:,f,z)); % Do not perform inversion - use original magnitude of PSF
        end
        
    end
end

%% (3A) CLUSTER: ROI 
% ---------------------------------------------------------------------- 
% The default ROI detection method. It detects the ROI by searching for a 
% center cluster of points with a value higher than the baseline 3dB loss 
% from 'maxima'. The center cluster must be within 25% of the center of the 
% image (ie. middle 50% of image) - the algorithm iterates in increments of 
% 5%, up to 50% (ie. 100% of image) until a center cluster is detected.
% ---------------------------------------------------------------------- 

sum_PSF_f = zeros(nx, ny, nz);
center_cluster_region_percent = 0.25:0.05:0.50;
percent_step_size = 0.05; 
circle_size = round(nx*percent_step_size); 

if ~var_input.maxima_ROI
    for z = 1:nz

        sum_PSF_f(:,:,z) = sum(PSF_inverted(:,:,:,z),3); % Summation across all frequencies (taking magnitude of complex dataset) for each z-plane

        baseline_threshold = (max(max(sum_PSF_f(:,:,z))) - min(min(sum_PSF_f(:,:,z))))*(10^-(cluster_threshold/10)); % -3dB = 10^-0.3 from maximum (default)
        threshold = max(max(sum_PSF_f(:,:,z))) - baseline_threshold;

        ROI_points = threshold < sum_PSF_f(:,:,z); % All points above 3dB threshold
        ROI_clusters = bwconncomp(ROI_points); % Identify all clusters above 3dB threshold
        ROI_clusters_centroids = regionprops(ROI_clusters,'Centroid'); % Calculate centroids of each cluster
        ROI_clusters_centroids = cell2mat(struct2cell(ROI_clusters_centroids)');

        clear ROI_clusters_centroids_region % Need to clear for each z-plane in order for conditional below to be valid at the start
        m = 0; 
        
        % Iterate until a cluster is found within the region 
        while ~exist('ROI_clusters_centroids_region','var') 

            m = m + 1;
            x_lower = center_xy(1,1) - floor(nx*center_cluster_region_percent(m)); % Updating search box for center cluster region
            y_lower = center_xy(1,2) - floor(ny*center_cluster_region_percent(m));

            x_upper = center_xy(1,1) + floor(nx*center_cluster_region_percent(m));
            y_upper = center_xy(1,2) + floor(ny*center_cluster_region_percent(m));

            c = 1;
            for k = 1:length(ROI_clusters_centroids(:,1))
                if (ROI_clusters_centroids(k,1) > x_lower && ROI_clusters_centroids(k,1) < x_upper) && (ROI_clusters_centroids(k,2) > y_lower && ROI_clusters_centroids(k,2) < y_upper)
                    ROI_clusters_centroids_region(c,:) = [ROI_clusters_centroids(k,1),ROI_clusters_centroids(k,2)]; % Holds all centroids that are within center cluster region search box
                    c = c + 1;
                end
            end
        end
        
        if m > 1 % Notifies user of successful search region (ie. if 50%, indicates that ROI is very off-centered)
            fprintf('Important: Center cluster for PSF was found within %d%% of center point of image! \n', center_cluster_region_percent(m)*100); 
        end
            
        center_cluster_i = knnsearch(ROI_clusters_centroids_region,center_xy); % Finding centroid closest to center of image
        center_cluster_i = find(ROI_clusters_centroids(:,1) == ROI_clusters_centroids_region(center_cluster_i,1)); % Finding index from original set of coordinates (looking at column 1 only)

        ROI_cluster_linear_index = ROI_clusters.PixelIdxList{center_cluster_i}; % Extract linear indices of center cluster (now ROI) from ROI_clusters struct
        [row_index_cluster, col_index_cluster] = ind2sub([nx, ny], ROI_cluster_linear_index); % Convert linear indices to their indices for row and column
        ROI_cluster = [row_index_cluster,col_index_cluster]; 

        % Creating a logical matrix with ROI = 1, non-ROI = 0 [x,y]
        determined_ROI = zeros(nx,ny); 
        num_pts_cluster = length(ROI_cluster(:,1));
        
        for k = 1:num_pts_cluster % # of rows
            determined_ROI(ROI_cluster(k,1),ROI_cluster(k,2)) = 1; 
        end
        
        % Defining a circular exclusion zone, using determined_ROI logical
        center_coords(z,:) = ROI_clusters_centroids(center_cluster_i,:); % Holds centroid coordinates (each z-plane)
        x_center = center_coords(z,1);
        y_center = center_coords(z,2);

        outer_boundary_cell = bwboundaries(determined_ROI);
        outer_boundary = outer_boundary_cell{1};
        boundary_distances = sqrt((outer_boundary(:,2) - x_center).^2 + (outer_boundary(:,1) - y_center).^2);

        [~, max_radius_i] = max(boundary_distances); % Determines maximum radius out of all boundary coordinates, and finds index of coordinates
        max_point_coord = outer_boundary(max_radius_i,:); % Boundaries for max radius
        
        radius_init(z,1) = ceil(sqrt((max_point_coord(1,2) - x_center)^2 + (max_point_coord(1,1) - y_center)^2)); % Use centroid and max_boundary coords to determine radius of circular mask
        
        % Check size of radius (to avoid using a very abnormal cluster)
        iteration_limit_init = nx - x_center; % left
        if x_center < iteration_limit_init 
            iteration_limit_init = x_center; % right
        end
        
        iteration_limit_init_y = ny - y_center; % top
        if y_center < iteration_limit_init_y 
            iteration_limit_init_y = y_center; % bottom
        end
        
        if iteration_limit_init_y < iteration_limit_init
            iteration_limit_init = iteration_limit_init_y;
        end

        % If cluster is an abnormal shape, could lead to radius that goes over image dimensions - a flag to switch to maxima method
        if radius_init(z,1) > iteration_limit_init 
            flag_maxima = 1;
            fprintf('Warning: could not find suitable ROI cluster and defaulted to maxima ROI method. Dataset may be very noisy.\n\n');
            break % Switching to maxima method
        end
        
        [xx,yy] = ndgrid((1:nx) - y_center,(1:ny) - x_center);
        circle_ez_logical = (xx.^2 + yy.^2) < radius_init(z,1)^2;

        for f = 1:nf % Mapping ROI and Exclusion Zone to each frequency
            ROI_mask{f,z} = determined_ROI; % Complete matrix of logical mask for ROI (does not vary w frequency, just recopied)
            ez_circle_initial{f,z} = circle_ez_logical;
        end
    end

end

%% (3B) MAXIMA: ROI
% ---------------------------------------------------------------------- 
% This is the alternative method of ROI detection, which can be selected by
% the user for use. It detects the maxima of the frequency-summed PSF and
% uses that as a center-point for a box, whose default size is that of the
% scattering probe [10x10mm]. 
% ---------------------------------------------------------------------- 

if var_input.maxima_ROI || flag_maxima 
    
    if flag_maxima % % Executing maxima method because cluster method had flag raised
        clearvars determined_ROI center_coords radius_init ROI_mask Ez_circle_initial % Clear variables that may have been used by cluster method
    end
    
    ROI_fixed_size = var_input.ROI_fixed_size; % As specified by user, otherwise default 10x10 mm [row,col > Y, X]
    ROI_fixed_size = floor(ROI_fixed_size/2); % 1px = 2mm
    ROI_fixed_step_size = floor(ROI_fixed_size/2); 
    
    for z = 1:nz
        sum_PSF_f(:,:,z) = sum(PSF_inverted(:,:,:,z),3); % Summation across all frequencies (taking magnitude of complex dataset) for each z-plane
        PSF_initial_estimate(:,:,z) = NaN([nx ny]); % Initialize PSF-sized matrix with NaN values
        PSF_initial_estimate(PSF_init_row,PSF_init_col,z) = sum_PSF_f(PSF_init_row,PSF_init_col,z); % Holds values within initial PSF estimate
        sum_PSF_f_max = find(PSF_initial_estimate(:,:,z) == max(max(PSF_initial_estimate(:,:,z))));
        [row_index_maxima, col_index_maxima] = ind2sub([nx, ny], sum_PSF_f_max); % X and Y coordinates of summation maxima
        
        ROI_X_maxima = (row_index_maxima - ROI_fixed_step_size(1,1)):(row_index_maxima + ROI_fixed_step_size(1,1)); % X values
        ROI_Y_maxima = (col_index_maxima - ROI_fixed_step_size(1,2)):(col_index_maxima + ROI_fixed_step_size(1,2)); % Y values
        determined_ROI = zeros(nx,ny); % Creating a logical matrix with ROI = 1, non-ROI = 0
        determined_ROI(ROI_X_maxima, ROI_Y_maxima) = 1;
        
        % Defining a circular exclusion zone, using determined_ROI logical
        center_coords(z,:) = [col_index_maxima, row_index_maxima]; % Holds centroid coordinates (each z-plane)
        x_center = center_coords(z,1);
        y_center = center_coords(z,2);

        outer_boundary_cell = bwboundaries(determined_ROI);
        outer_boundary = outer_boundary_cell{1};
        boundary_distances = sqrt((outer_boundary(:,2) - x_center).^2 + (outer_boundary(:,1) - y_center).^2);

        [~, max_radius_i] = max(boundary_distances); % Determines maximum radius out of all boundary coordinates, and finds index of coordinates
        max_point_coord = outer_boundary(max_radius_i,:); % Boundaries for max radius

        radius_init(z,1) = ceil(sqrt((max_point_coord(1,2) - x_center)^2 + (max_point_coord(1,1) - y_center)^2)); % Use centroid and max_boundary coords to determine radius of circular mask
        
        [xx,yy] = ndgrid((1:nx) - y_center,(1:ny) - x_center);
        circle_ez_logical = (xx.^2 + yy.^2) < radius_init(z,1)^2;
        
        for f = 1:nf
            ROI_mask{f,z} = determined_ROI; % Complete matrix of logical mask for ROI (does not vary w frequency, just recopied)
            ez_circle_initial{f,z} = circle_ez_logical; 
        end
    end

end


%% (4) PLOT ROI AGAINST FREQUENCY SUMMATION (OPTIONAL)
% ---------------------------------------------------------------------- 

if var_input.plot_ROI_summation % Plot frequency summation of PSF with automated ROI 
    
    figure; 
    for z = 1:nz
        ROI_xy_z = cell2mat(ROI_mask(1,z));
        
        if nz > 1 % Plot all z-planes in one figure
            subplot(round(nz/2),2,z);
        end

        imagesc(sum_PSF_f(:,:,z)); colorbar; 
        hold on;
        mask_ROI = cat(3, zeros(nx, ny), zeros(nx, ny), zeros(nx, ny));
        h = imagesc(mask_ROI); 
        viscircles(center_coords(z,:),radius_init(z,1),'Color','w', 'LineWidth', 0.1); 
        hold off;
        set(h, 'AlphaData', ~ROI_xy_z*0.25); % Automated ROI will appear as bright set of pixels
        title(['ROI with initial exclusion zone against frequency summation, z-plane: ' num2str(z)]);
        ax = gca; ax.FontName = 'Times New Roman';
        set(gca,'YDir','normal');
    end
    
end

%% (5) Exclusion Zone, Sizing 
% ---------------------------------------------------------------------- 
% Refining definition of initial exclusion zone by iterating through
% different exclusion zones (increments of 5% length-wise). If change in variance
% between exclusion zones is less than 4%, take the smaller exclusion zone.
% ---------------------------------------------------------------------- 

for f = 1:nf
    for z = 1:nz
        radius = radius_init(z,1); % Start circle with radius of rho_init, from initial sizing

        ROI_xy_f_z = cell2mat(ROI_mask(f,z)); % Extracting ROI for specific frequency and z-plane into logical matrix (row 1 = x, row 2 = y)
        PSF_mag_f_z = abs(PSF_inverted(:,:,f,z)); % Magnitude dataset, afrer inversion was performed
        
        ROI_mag = NaN(nx); % Initialize as NaN matrix
        ROI_mag(ROI_xy_f_z == 1) = PSF_mag_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
        
        ez_circle_f_z = cell2mat(ez_circle_initial(f,z));

        x_center = center_coords(z,1);
        y_center = center_coords(z,2);
        [xx,yy] = ndgrid((1:nx) - y_center,(1:ny) - x_center);

        if (nx - (x_center + radius)) > (x_center - radius)
            iteration_limit = round((x_center - radius)/circle_size) - 2; % Limiting side is the left
        else
            iteration_limit = round((nx - (x_center+radius))/circle_size) - 2; % Limiting side is the right
        end
        
        if (ny - (y_center + radius)) > (y_center - radius)
           iteration_limit_y = round((y_center - radius)/circle_size) - 2; % Limiting side is the bottom
        else
           iteration_limit_y = round((ny - (y_center + radius))/circle_size) - 2; % Limiting side is the top
        end

        if (iteration_limit_y < iteration_limit) % Select smaller range for x,y as iteration restraint
            iteration_limit = iteration_limit_y;
        end
        
        if iteration_limit < 1 % Temp fix for rytov bowtie
            iteration_limit = 1;
        end
        
        n_conditional = 1; % Count for incremental conditional
        for k_increment = 1:iteration_limit

            if (k_increment == 1)
                exclusion_zone_mag = PSF_mag_f_z; % Magnitude dataset for exclusion zone
                exclusion_zone_mag(ez_circle_f_z == 1) = NaN; % Exclusion zone is NaN in ROI + padding
                
                variance_iterations(f,k_increment,z) = (std2(abs(exclusion_zone_mag(~isnan(exclusion_zone_mag)))))^2;
            else
                radius = radius_2; 
            end
            
            radius_2 = radius + circle_size;
            Ez_circle_f_z_2 = (xx.^2 + yy.^2) < radius_2^2; % Next exclusion zone iteration

            exclusion_zone_mag_2 = PSF_mag_f_z;
            exclusion_zone_mag_2(Ez_circle_f_z_2 == 1) = NaN;
            
            variance_iterations(f,k_increment+1,z) = (std2(abs(exclusion_zone_mag_2(~isnan(exclusion_zone_mag_2)))))^2;
            delta_percent_var(f,k_increment,z) = abs(variance_iterations(f,k_increment+1,z) - variance_iterations(f,k_increment,z))/(variance_iterations(f,k_increment,z))*100;

            if delta_percent_var(f,k_increment,z) < 4.0 % Increment counter for satisfied conditional of change in variance < 4% between iterations
                n_conditional = n_conditional + 1;
            end

            if n_conditional == 2 % Conditional of change in variance < 4% between iterations only needs to be met once - then break out of loop
                break
            end
        end
        
        if (k_increment == iteration_limit) && (n_conditional == 1) % Conditional never met in all iterations, so selecting index with smallest change in variance to avoid extremely large exclusion zones
            percent_var = delta_percent_var(f,:,z);
            k_increment = find(delta_percent_var(f,:,z) == min(percent_var(percent_var > 0)));
            k_increment = k_increment - 1;
        end
       
        radius = radius_init(z,1) + circle_size*(k_increment-n_conditional+1); % Final radius determined for f,z
        determined_ez = (xx.^2 + yy.^2) < radius^2; % Next exclusion zone iteration
        
        radius_refined(f,z) = radius;
        ez_mask{f,z} = determined_ez; % Complete matrix of logical mask for final exclusion zones 
    end
end 

%% (6) CNR and SNR CALCULATIONS
% ---------------------------------------------------------------------- 
% CNR and SNR is calculated for each frequency and z-plane. The type of CNR 
% and SNR calculation dictates what is fed into the general equation. 

% Note: To calculate CNR, the standard equation used is | mean(ROI) -
% mean(background) | divided by standard dev(background). To
% calculate SNR, the standard equation used is mean(ROI)/standard 
% dev(background).
% ---------------------------------------------------------------------- 

% CNR for each frequency and z-plane
CNR_mag = zeros(nf, nz); 
CNR_phase = zeros(nf, nz);
CNR_real = zeros(nf, nz); CNR_imag = zeros(nf, nz);
CNR_real_imag = zeros(nf, nz); % Average of CNR_real and CNR_imag values
CNR_complex_mag = zeros(nf, nz);

% SNR for each frequency and z-plane
SNR_mag = zeros(nf, nz); 
SNR_phase = zeros(nf, nz);
SNR_real = zeros(nf, nz); SNR_imag = zeros(nf, nz);
SNR_real_imag = zeros(nf, nz); % Average of SNR_real and SNR_imag values
SNR_complex_mag = zeros(nf, nz);

for f = 1:nf
    for z = 1:nz
        PSF_phase(:,:,f,z) = phase_unwrap(angle(PSF(:,:,f,z))); % Phase of CO; need to unwrap (both x and y); using custom unwrap function 
    end
end

for f = 1:nf
    for z = 1:nz 

        ROI_xy_f_z = cell2mat(ROI_mask(f,z)); % Extracting ROI for specific frequency and z-plane into logical matrix (row 1 = x, row 2 = y)
        Ez_logical_f = cell2mat(ez_mask(f,z)); % Extracting exclusion zone for specific frequency and z-plane into logical matrix 
        
        % -------- Magnitude calculations
        
        if rho_mag_CNR ~= 0 || rho_mag_SNR ~= 0
            PSF_mag_f_z = abs(PSF_inverted(:,:,f,z)); % Magnitude dataset, inverted
            
            ROI_mag = NaN(nx); % Initialize as NaN matrix
            ROI_mag(ROI_xy_f_z == 1) = PSF_mag_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
            
            exclusion_zone_mag = PSF_mag_f_z; % Magnitude dataset for exclusion zone
            exclusion_zone_mag(Ez_logical_f == 1) = NaN; % Exclusion zone is NaN in ROI + padding
            
            if var_input.peak_signal    % Peak SNR calculation method
                if rho_mag_CNR ~=0
                    CNR_mag(f,z) = abs( max(max(ROI_mag(~isnan(ROI_mag)))) - max(max(exclusion_zone_mag(~isnan(exclusion_zone_mag)))) )/std2(abs(exclusion_zone_mag(~isnan(exclusion_zone_mag))));
                end
                if rho_mag_SNR ~= 0
                    SNR_mag(f,z) = max(max(ROI_mag(~isnan(ROI_mag))))/std2(abs(exclusion_zone_mag(~isnan(exclusion_zone_mag))));
                end
            else                        % Mean SNR calculation method
                if rho_mag_CNR ~= 0
                    CNR_mag(f,z) = abs( mean2(ROI_mag(~isnan(ROI_mag))) - mean2(exclusion_zone_mag(~isnan(exclusion_zone_mag))) )/std2(abs(exclusion_zone_mag(~isnan(exclusion_zone_mag))));
                end
                
                if rho_mag_SNR ~= 0
                    SNR_mag(f,z) = mean2(ROI_mag(~isnan(ROI_mag)))/std2(abs(exclusion_zone_mag(~isnan(exclusion_zone_mag))));
                end
            end
        end
        
        % -------- Phase calculations
        
        if rho_phase_CNR ~= 0 || rho_phase_SNR ~= 0
            PSF_phase_f_z = PSF_phase(:,:,f,z) - min(min(PSF_phase(:,:,f,z))); % Phase dataset, baseline to 0
            
            ROI_phase = NaN(nx); % Initialize as NaN matrix
            ROI_phase(ROI_xy_f_z == 1) = PSF_phase_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
            
            exclusion_zone_phase = PSF_phase_f_z; % Phase dataset for exclusion zone
            exclusion_zone_phase(Ez_logical_f == 1) = NaN; % Exclusion zone is NaN in ROI + padding
            
            if var_input.peak_signal    % Peak SNR calculation method
                if rho_phase_CNR ~= 0
                    CNR_phase(f,z) = abs( max(max(ROI_phase(~isnan(ROI_phase)))) - max(max(exclusion_zone_phase(~isnan(exclusion_zone_phase)))) )/std2(exclusion_zone_phase(~isnan(exclusion_zone_phase)));
                end
                if rho_phase_SNR ~=0 
                    SNR_phase(f,z) = max(max(ROI_phase(~isnan(ROI_phase))))/std2(exclusion_zone_phase(~isnan(exclusion_zone_phase)));
                end
            else                        % Mean SNR calculation method
                if rho_phase_CNR ~= 0
                    CNR_phase(f,z) = abs( mean2(ROI_phase(~isnan(ROI_phase))) - mean2(exclusion_zone_phase(~isnan(exclusion_zone_phase))) )/std2(exclusion_zone_phase(~isnan(exclusion_zone_phase)));
                end
                
                if rho_phase_SNR ~= 0
                    SNR_phase(f,z) = mean2(ROI_phase(~isnan(ROI_phase)))/std2(exclusion_zone_phase(~isnan(exclusion_zone_phase)));
                end
            end
        end
        
        % -------- Real and Imaginary calculations
        
        if rho_real_imag_CNR ~= 0 || rho_real_imag_SNR ~= 0
            PSF_complex_f_z = PSF(:,:,f,z); % Complex dataset
            PSF_real_f_z = real(PSF_complex_f_z) - min(min(real(PSF_complex_f_z))); % Real dataset, baseline to 0
            
            ROI_real = NaN(nx); % Initialize as NaN matrix
            ROI_real(ROI_xy_f_z == 1) = PSF_real_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
            
            exclusion_zone_real = PSF_real_f_z;
            exclusion_zone_real(Ez_logical_f == 1) = NaN;

            PSF_imag_f_z = imag(PSF_complex_f_z) - min(min(imag(PSF_complex_f_z))); % Imaginary dataset, baseline to 0
            
            ROI_imag = NaN(nx); % Initialize as NaN matrix
            ROI_imag(ROI_xy_f_z == 1) = PSF_imag_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
            
            exclusion_zone_imag = PSF_imag_f_z;
            exclusion_zone_imag(Ez_logical_f == 1) = NaN;
            
            if var_input.peak_signal    % Peak SNR calculation method
                if rho_real_imag_CNR ~= 0
                    CNR_real(f,z) = abs( max(max(ROI_real(~isnan(ROI_real)))) - max(max(exclusion_zone_real(~isnan(exclusion_zone_real)))) )/std2(exclusion_zone_real(~isnan(exclusion_zone_real)));
                    CNR_imag(f,z) = abs( max(max(ROI_imag(~isnan(ROI_imag)))) - max(max(exclusion_zone_imag(~isnan(exclusion_zone_imag)))) )/std2(exclusion_zone_imag(~isnan(exclusion_zone_imag)));
                end
                if rho_real_imag_SNR ~= 0
                    SNR_real(f,z) = max(max(ROI_real(~isnan(ROI_real))))/std2(exclusion_zone_real(~isnan(exclusion_zone_real)));
                    SNR_imag(f,z) = max(max(ROI_imag(~isnan(ROI_imag))))/std2(exclusion_zone_imag(~isnan(exclusion_zone_imag)));
                end
                
            else                        % Mean SNR calculation method
                if rho_real_imag_CNR ~= 0
                    CNR_real(f,z) = abs( mean2(ROI_real(~isnan(ROI_real))) - mean2(exclusion_zone_real(~isnan(exclusion_zone_real))) )/std2(exclusion_zone_real(~isnan(exclusion_zone_real)));
                    CNR_imag(f,z) = abs( mean2(ROI_imag(~isnan(ROI_imag))) - mean2(exclusion_zone_imag(~isnan(exclusion_zone_imag))) )/std2(exclusion_zone_imag(~isnan(exclusion_zone_imag)));
                end
                
                if rho_real_imag_SNR ~= 0
                    SNR_real(f,z) = mean2(ROI_real(~isnan(ROI_real)))/std2(exclusion_zone_real(~isnan(exclusion_zone_real)));
                    SNR_imag(f,z) = mean2(ROI_imag(~isnan(ROI_imag)))/std2(exclusion_zone_imag(~isnan(exclusion_zone_imag)));
                end
            end
            
            CNR_real_imag(f,z) = (CNR_imag(f,z) + CNR_real(f,z))/2;
            SNR_real_imag(f,z) = (SNR_imag(f,z) + SNR_real(f,z))/2;
        end
        
        % -------- Complex Magnitude calculations
        
        if rho_complex_mag_CNR ~= 0 || rho_complex_mag_SNR ~= 0
            PSF_complex_f_z = PSF(:,:,f,z); % Complex dataset
            
            ROI_complex = NaN(nx); % Initialize as NaN matrix
            ROI_complex(ROI_xy_f_z == 1) = PSF_complex_f_z(ROI_xy_f_z == 1); % Fill in ROI with signal, and the rest as NaN
            
            exclusion_zone_complex = PSF_complex_f_z; % Complex dataset for exclusion zone
            exclusion_zone_complex(Ez_logical_f == 1) = NaN;
            
            if var_input.peak_signal    % Peak SNR calculation method
                if rho_complex_mag_CNR ~= 0
                    CNR_complex_mag(f,z) = abs( max(max(ROI_complex(~isnan(ROI_complex)))) - max(max(exclusion_zone_complex(~isnan(exclusion_zone_complex)))) )/std2(exclusion_zone_complex(~isnan(exclusion_zone_complex)));
                end
                if rho_complex_mag_SNR ~= 0
                    SNR_complex_mag(f,z) = abs( max(max(ROI_complex(~isnan(ROI_complex)))) )/std2(exclusion_zone_complex(~isnan(exclusion_zone_complex)));
                end
                
            else                        % Mean SNR calculation method
                if rho_complex_mag_CNR ~= 0
                    CNR_complex_mag(f,z) = abs( abs(mean2(ROI_complex(~isnan(ROI_complex)))) - abs(mean2(exclusion_zone_complex(~isnan(exclusion_zone_complex)))) )/(std2(exclusion_zone_complex(~isnan(exclusion_zone_complex))));
                end
                
                if rho_complex_mag_SNR ~= 0
                    SNR_complex_mag(f,z) = abs( mean2(ROI_complex(~isnan(ROI_complex))))/std2(exclusion_zone_complex(~isnan(exclusion_zone_complex)));
                end
            end

        end
    end
end

%% (7) PSF QUALITY CLASSIFICATION USING THRESHOLDS
% ---------------------------------------------------------------------- 
% Output the frequency indices of 'good' PSFs by using thresholds. The CNR
% or SNR evaluated against the thresholds is the average across z-planes.
% ---------------------------------------------------------------------- 

for f = 1:nf
   CNR_mag_avg(f) = mean(CNR_mag(f,:));
   CNR_phase_avg(f) = mean(CNR_phase(f,:));
   CNR_real_imag_avg(f) = mean(CNR_real_imag(f,:));
   CNR_complex_mag_avg(f) = mean(CNR_complex_mag(f,:));
   SNR_mag_avg(f) = mean(SNR_mag(f,:));
   SNR_phase_avg(f) = mean(SNR_phase(f,:));
   SNR_real_imag_avg(f) = mean(SNR_real_imag(f,:));
   SNR_complex_mag_avg(f) = mean(SNR_complex_mag(f,:));
end

% Store 'good' frequencies in the freq_good structure for e/ classification
freq_good = struct('CNR_mag_i',{},'SNR_mag_i',{},'CNR_phase_i',{},'SNR_phase_i',{},'CNR_real_imag_i',{},'SNR_real_imag_i',{},'CNR_complex_mag_i',{},'SNR_complex_mag_i',{});

if rho_mag_CNR
    freq_good(1).CNR_mag_i = find(CNR_mag_avg > rho_mag_CNR)'; % Indices at which mag CNR is higher than user threshold
end

if rho_mag_SNR
    freq_good(1).SNR_mag_i = find(SNR_mag_avg > rho_mag_SNR)'; % Indices at which mag SNR is higher than user threshold
end

if rho_phase_CNR
    freq_good(1).CNR_phase_i = find(CNR_phase_avg > rho_phase_CNR)';
end

if rho_phase_SNR
    freq_good(1).SNR_phase_i = find(SNR_phase_avg > rho_phase_SNR)';
end

if rho_real_imag_CNR
    freq_good(1).CNR_real_imag_i = find(CNR_real_imag_avg > rho_real_imag_CNR)';
end

if rho_real_imag_SNR
   freq_good(1).SNR_real_imag_i = find(SNR_real_imag_avg > rho_real_imag_SNR)';
end

if rho_complex_mag_CNR
    freq_good(1).CNR_complex_mag_i = find(CNR_complex_mag_avg > rho_complex_mag_CNR)';
end

if rho_complex_mag_SNR
    freq_good(1).SNR_complex_mag_i = find(SNR_complex_mag_avg > rho_complex_mag_SNR)';
end

%% (8) PLOT SNR VS. FREQUENCY INDEX (OPTIONAL)
% ---------------------------------------------------------------------- 

if var_input.plot_CNR_SNR
    
    % -------- Magnitude plots
    
    if rho_mag_CNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z); 
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(CNR_mag(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(CNR_mag_avg),'r--');
                line([ftable(1),ftable(end)],[10*log10(rho_mag_CNR),10*log10(rho_mag_CNR)]); hold off; % Horizontal line indicating magnitude SNR threshold
                ylabel('CNR (dB)'); 
            else % linear scale
                plot(ftable,CNR_mag(:,z),'k-'); hold on; 
                plot(ftable,CNR_mag_avg,'r--');
                line([ftable(1),ftable(end)],[rho_mag_CNR,rho_mag_CNR]); hold off; % Horizontal line indicating magnitude SNR threshold
                ylabel('CNR'); 
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Magnitude CNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('CNR','CNR avg. (across z)','Threshold CNR');
        lgd.FontName = 'Times New Roman';
    end
        
    if rho_mag_SNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z); 
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(SNR_mag(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(SNR_mag_avg),'r--'); 
                line([ftable(1),ftable(end)],[10*log10(rho_mag_SNR),10*log10(rho_mag_SNR)]); hold off; % Horizontal line indicating magnitude SNR threshold
                ylabel('SNR (dB)'); 
            else % linear scale
                plot(ftable,SNR_mag(:,z),'k-'); hold on; 
                plot(ftable,SNR_mag_avg,'r--'); 
                line([ftable(1),ftable(end)],[rho_mag_SNR,rho_mag_SNR]); hold off; % Horizontal line indicating magnitude SNR threshold
                ylabel('SNR'); 
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Magnitude SNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('SNR','SNR avg. (across z)','Threshold SNR');
        lgd.FontName = 'Times New Roman';
    end
        
    % -------- Phase plots
    
    if rho_phase_CNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale 
                plot(ftable,10*log10(CNR_phase(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(CNR_phase_avg),'r--');
                line([ftable(1),ftable(end)],[10*log10(rho_phase_CNR),10*log10(rho_phase_CNR)]); hold off; % Horizontal line indicating phase SNR threshold
                ylabel('CNR (dB)'); 
            else % linear scale 
                plot(ftable,CNR_phase(:,z),'k-'); hold on; 
                plot(ftable,CNR_phase_avg,'r--');
                line([ftable(1),ftable(end)],[rho_phase_CNR,rho_phase_CNR]); hold off; % Horizontal line indicating phase SNR threshold
                ylabel('CNR'); 
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Phase CNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('CNR','CNR avg. (across z)','Threshold CNR');
        lgd.FontName = 'Times New Roman';
    end
       
    if rho_phase_SNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(SNR_phase(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(SNR_phase_avg),'r--');
                line([ftable(1),ftable(end)],[10*log10(rho_phase_SNR),10*log10(rho_phase_SNR)]); hold off; % Horizontal line indicating phase SNR threshold
                ylabel('SNR (dB)'); 
            else % linear scale
                plot(ftable,SNR_phase(:,z),'k-'); hold on; 
                plot(ftable,SNR_phase_avg,'r--');
                line([ftable(1),ftable(end)],[rho_phase_SNR,rho_phase_SNR]); hold off; % Horizontal line indicating phase SNR threshold
                ylabel('SNR'); 
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Phase SNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('SNR','SNR avg. (across z)','Threshold SNR');
        lgd.FontName = 'Times New Roman';
    end
    
    % -------- Real/Imaginary plots
    
    if rho_real_imag_CNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(CNR_real_imag(:,z)),'k-','linewidth',0.75); hold on;
                plot(ftable,10*log10(CNR_real_imag_avg),'r--','linewidth',0.75); 
                plot(ftable,10*log10(CNR_imag(:,z)),'b-.');
                plot(ftable,10*log10(CNR_real(:,z)),'g-.');
                line([ftable(1),ftable(end)],[10*log10(rho_real_imag_CNR),10*log10(rho_real_imag_CNR)]); hold off; % Horizontal line indicating real/imaginary SNR threshold
                ylabel('CNR (dB)'); 
            else % linear scale
                plot(ftable,CNR_real_imag(:,z),'k-','linewidth',0.75); hold on;
                plot(ftable,CNR_real_imag_avg,'r--','linewidth',0.75); 
                plot(ftable,CNR_imag(:,z),'b-.');
                plot(ftable,CNR_real(:,z),'g-.');
                line([ftable(1),ftable(end)],[rho_real_imag_CNR,rho_real_imag_CNR]); hold off; % Horizontal line indicating real/imaginary SNR threshold
                ylabel('CNR'); 
            end
            xlim([ftable(1),ftable(end)]);  grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Real and Imaginary CNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('CNR','CNR avg. (across z)','CNR Imag', 'CNR Real','Threshold CNR');
        lgd.FontName = 'Times New Roman';
    end
        
    if rho_real_imag_SNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(SNR_real_imag(:,z)),'k-','linewidth',0.75); hold on;
                plot(ftable,10*log10(SNR_real_imag_avg),'r--','linewidth',0.75); 
                plot(ftable,10*log10(SNR_imag(:,z)),'b-.');
                plot(ftable,10*log10(SNR_real(:,z)),'g-.');
                line([ftable(1),ftable(end)],[10*log10(rho_real_imag_SNR),10*log10(rho_real_imag_SNR)]); hold off; % Horizontal line indicating real/imaginary SNR threshold
                ylabel('SNR (dB)'); 
            else
                plot(ftable,SNR_real_imag(:,z),'k-','linewidth',0.75); hold on;
                plot(ftable,SNR_real_imag_avg,'r--','linewidth',0.75); 
                plot(ftable,SNR_imag(:,z),'b-.');
                plot(ftable,SNR_real(:,z),'g-.');
                line([ftable(1),ftable(end)],[rho_real_imag_SNR,rho_real_imag_SNR]); hold off; % Horizontal line indicating real/imaginary SNR threshold
                ylabel('SNR'); 
            end
            xlim([ftable(1),ftable(end)]);  grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Real and Imaginary SNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('SNR','SNR avg. (across z)','SNR Imag', 'SNR Real','Threshold SNR');
        lgd.FontName = 'Times New Roman';
    end
    
    % -------- Complex Magnitude plots
    
    if rho_complex_mag_CNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(CNR_complex_mag(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(CNR_complex_mag_avg),'r--');
                line([ftable(1),ftable(end)],[10*log10(rho_complex_mag_CNR),10*log10(rho_complex_mag_CNR)]); hold off; % Horizontal line indicating complex magnitude SNR threshold
                ylabel('CNR (dB)');
            else % linear scale
                plot(ftable,CNR_complex_mag(:,z),'k-'); hold on; 
                plot(ftable,CNR_complex_mag_avg,'r--');
                line([ftable(1),ftable(end)],[rho_complex_mag_CNR,rho_complex_mag_CNR]); hold off; % Horizontal line indicating complex magnitude SNR threshold
                ylabel('CNR');
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            ylim([-25 15]);
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Complex Magnitude CNR, z-plane: ' num2str(z)]);
            
        end
        lgd = legend('CNR','CNR avg. (across z)','Threshold CNR');
        lgd.FontName = 'Times New Roman';
    end
        
    if rho_complex_mag_SNR ~= 0
        figure; clf; 
        for z = 1:nz
            subplot(nz,1,z);
            if var_input.dB_scale == 1 % dB scale
                plot(ftable,10*log10(SNR_complex_mag(:,z)),'k-'); hold on; 
                plot(ftable,10*log10(SNR_complex_mag_avg),'r--');
                line([ftable(1),ftable(end)],[10*log10(rho_complex_mag_SNR),10*log10(rho_complex_mag_SNR)]); hold off; % Horizontal line indicating complex magnitude SNR threshold
                ylabel('SNR (dB)');
            else % linear scale
                plot(ftable,SNR_complex_mag(:,z),'k-'); hold on; 
                plot(ftable,SNR_complex_mag_avg,'r--');
                line([ftable(1),ftable(end)],[rho_complex_mag_SNR,rho_complex_mag_SNR]); hold off; % Horizontal line indicating complex magnitude SNR threshold
                ylabel('SNR');
            end
            xlim([ftable(1),ftable(end)]); grid minor;
            xlabel(f_plot_title); 
            ax = gca; ax.FontName = 'Times New Roman';
            title(['Complex Magnitude SNR, z-plane: ' num2str(z)]);
        end
        lgd = legend('SNR','SNR avg. (across z)','Threshold SNR');
        lgd.FontName = 'Times New Roman';
    end

end

%% (9) Setting variable output (for development, plotting purposes)
% ---------------------------------------------------------------------- 
% User can use these in the function call if needed for plotting the ROI
% and/or exclusion zone. Otherwise, can omit.
% ---------------------------------------------------------------------- 

exclusion_zone = struct('ez_mask', {}, 'center_coords', {}, 'radius', {});
exclusion_zone(1).ez_mask = ez_mask; 
exclusion_zone(1).center_coords = center_coords;
exclusion_zone(1).radius = radius_refined;

varargout{1,:} = ROI_mask; % Logical mask for ROI
varargout{2,:} = exclusion_zone; % Structure containing logical mask for exclusion zone + center coordinates and radius that define the circle
varargout{3,:} = PSF_inverted; % Initial magnitude CO dataset with inversion performed