function [output_vals,FT_vals,FFT_bins_out] = fft_circle(x0,y0,r,img,FFT_bins)
% developed by Nils Grosse Hokamp (University Hospital Cologne, Germany) and 
% Brendan Eck (Case Western Reserve University, Cleveland, OH, USA). 
% Please cite corresponding publications when using
% x0 = column index of center point (in pixels)
% y0 = row index of center point (in pixels)
% r = radius of circle (in pixels)
% img = CT image
% FFT_bins = 2d array with the first column indicating start index of a bin
%           and second column indicating end index of a bin. Number of rows
%           is the number of bins to be output. If second index =0, then
%           algorithm will use the maximum frequency value index. If both
%           indexes are 100, 200, or 400, it will partition frequency bins
%           into one bin, 2 bins, or 4 bins, respectively.
% FFT_bins_out = 2d array of bins that were used to calculate output_vals

thetas = [0:.01:360]*pi/180; % very small theta increments just to be sure no values are missed
x = r*cos(thetas);
y = r*sin(thetas);

x = round(x0 + x); % just keep pixel values; no interpolation
y = round(y0 + y); % just keep pixel values; no interpolation

linearIDX = unique(sub2ind(size(img),y,x),'stable'); % keep only each unique voxel
vals = img(linearIDX);
FT_vals = fft(vals);
FT_vals = abs(FT_vals); % just look at magnitude, not phase
count = 1;
for i=1:size(FFT_bins,1)
    start_idx = FFT_bins(i,1);
    end_idx = FFT_bins(i,2);
    if start_idx>=100&&end_idx>=100
        n_partitions = start_idx/100;
        max_idx = floor(length(FT_vals)/2);
        for kk=1:n_partitions
            start_idx = (kk-1)*floor(max_idx/n_partitions)+1;
            end_idx = (kk)*floor(max_idx/n_partitions);
            output_vals(count) = sum(FT_vals(start_idx:end_idx));
            FFT_bins_out(count,1) = start_idx;
            FFT_bins_out(count,2) = end_idx;
            count = count+1;
        end
        continue
    end
    if end_idx==0
        end_idx = floor(length(FT_vals)/2);
    end
    output_vals(count) = sum(FT_vals(start_idx:end_idx));
    FFT_bins_out(count,1) = start_idx;
    FFT_bins_out(count,2) = end_idx;
    count = count+1;
end



end