function [array, lims] = normalize_array(array, output_range, input_range)
%NORMALIZE_ARRAY stretches an n-D array to the desired range
% array - n-D array in single or double
% optional:
% output_range - empty (output will then be [0.0, 1.0]), 
%                scalar (output range will be [0, output_range] in that case),
%                or 2-element vector specifying output number range
% input_range - empty or not provided (the range is calculated from input data)
%               scalar (the percentage of array elements outside specified range (outliers) - between 0 and 1)
%               2-element vector (the minimum and maximum of the grayscale range of the input)

% TODO - update description

% validate inputs
validateattributes(array, {'single', 'double'}, {'nonempty'})
if nargin < 3 || isempty(input_range) % range not provided or empty
    input_range = 0.0; % default is 0% outliers
end
if numel(input_range) == 1 % range provided as single percentage of outliers
    validateattributes(input_range, {'numeric'}, {'scalar', '>=', 0.0, '<=', 1.0})
    lo = quantile(array(:), input_range/2);
    hi = quantile(array(:), 1-input_range/2);
else % range provided as a 2-element vector of minimum and maximum value
    validateattributes(input_range, {'numeric'}, {'numel', 2, 'nonnan', 'finite'})
    lo = double(input_range(1));
    hi = double(input_range(2));
end
if nargin < 2 || isempty(output_range)
    output_range = [0.0, 1.0];
else
    validateattributes(output_range, {'single', 'double'}, {'vector'})
    if numel(output_range) == 1
        output_range = [0, output_range];
    else
        assert(numel(output_range)==2, ...
            'range should be empty, a scalar, or 2-element vector')
    end
end

% normalize and scale
array = output_range(1) + (output_range(2) - output_range(1)) .* ...
        (array - lo) ./ (hi - lo);
lims = [lo, hi];
end


