function [spec, r_grid_bin] = radSpec(inputImage, n_bins)

%RADSPEC radial spectrum
%
%   [radSpec] = radSpec(im, n_bins)
%
%   Returns amplitude spectrum of image in frequency domain using circular
%   bounds for calculation
%
%   Input(s):
%   'im'        input image N1xN2xN3
%   'n_bins'    number of bins for histogram calculation
%
%   Output(s):
%   'radSpec'	spectrum response
%
%   See also DEBLURRING_KERNEL_ESTIMATION
%
%
%   Copyright (c) 2019 Mahdi S. Hosseini
%
%   University of Toronto
%   The Edward S. Rogers Sr. Department of,
%   Electrical and Computer Engineering (ECE)
%   Toronto, ON, M5S3G4, Canada
%   Tel: +1 (416) 978 6845
%   email: mahdi.hosseini@mail.utoronto.ca

[N_1, N_2, ~] = size(inputImage);
if mod(N_1, 2) == 0
    inputImage = padarray(inputImage, [1, 0], 'pre');
end
if mod(N_2, 2) == 0
    inputImage = padarray(inputImage, [0, 1], 'pre');
end
[N_1, N_2, N_3] = size(inputImage);
r_max = floor(max([N_1, N_2])/2);

for channel = 1: N_3
    I = inputImage(:,:,channel);
    FT = fftshift(fft2(I, 2*r_max+1, 2*r_max+1));
    frequency_spectrum = abs(FT);
    if false
        frequency_spectrum(:, floor(N_2/2)+1) = 0;
        frequency_spectrum(floor(N_1/2)+1, :) = 0;
    end
    h_bin = 1/n_bins;
    r_grid_bin = [h_bin:h_bin:1]'*pi;
    r_grid = [-r_max: r_max]/r_max*pi;
    [X, Y] = meshgrid(r_grid);
    Y = -Y;
    A = sqrt(X.^2 + Y.^2);    
    circle_mask = A <= pi;
    pixel_area = pi^3/sum(circle_mask(:));
    for iteration = 1: n_bins
        mask_out =  A <= r_grid_bin(iteration);
        if iteration ~= 1
            mask_in = A <= r_grid_bin(iteration-1);
            mask = xor(mask_out, mask_in);
        else
            mask = mask_out;
        end
        S(iteration) = pixel_area * sum(mask(:));
        spec(iteration, channel) = ...
            sum(frequency_spectrum(mask))/S(iteration);
    end
    spec(:, channel) = spec(:, channel)/spec(1, channel);
end