function [ theta ] = radonIce( I, alpha, levels)
%RADONICE Summary of this function goes here
%   Detailed explanation goes here

method = 'median'; % or "polynomial"

if nargin == 1
    alpha = 0.6;
    levels = 1;
    numtops = 3;
end
domain = 1:1:180; % angles
subdomain = .1:.1:180;

%% create bi-nomial filters
pascalTri = fliplr(pascal(2.*levels)); % triangle of Pascal

Rad = [];
for i = 1:2:2.*levels
    % binomial filter (1D)
    binom = diag(pascalTri,i);
    % laplacian filter (1D)
    lap = conv(binom,[1 -2 1]);
    lap = repmat(lap(:),1,length(lap));
    % laplacian filter (2D)
    lap2 = lap'+lap;
    % high pass filtering
    E = imfilter(I, lap2, 'symmetric'); % [-1 -2 -1; -2 12 -2; -1 -2 -1]
    %% randon transform
    RadonLevel = radoNan(E,domain);
%     RadonLevel = radon(E,domain);
    % stack results from multiple levels 
    Rad = cat(1, Rad, RadonLevel); 
end


%% look at statistics
stdRad = std(Rad,1); % standard deviation
medRad = median(cat(1, ... % running median filter
    circshift(stdRad,[0 -1]), stdRad, circshift(stdRad,[0 1])),1);

switch method
    case 'polynomial'
        %% estimate peak, through polynomial fit
        p = polyfit(domain,medRad,12);
        f = polyval(p,subdomain);
        [polyVal, polyIdx] = max(f);
        [pks,loc] = findpeaks(f,'sortstr', 'descend');
        
    otherwise
        %%  estimate peak through values
        [pks,loc] = findpeaks(medRad, 'minpeakdistance', 10, 'sortstr', 'descend');
        SNR = pks./std(stdRad);
end
% figure(); plot(domain,stdRad,'o',subdomain,f,'-')

%% culling
if ~isempty(pks)
    switch method
        case 'polynomial'
            if length(pks)~=1 % only one top
                if pks(2)./pks(1) < alpha % significant top
                    theta = NaN;
                else
                    theta = subdomain(polyIdx);
                end
            else
                theta = subdomain(polyIdx);
            end
        otherwise
            goodEs = max(find(SNR>alpha));
            if isempty(goodEs)
                theta = NaN.*ones(1,1,8);
            elseif goodEs>=4
                theta = permute(reshape([loc(1:4); pks(1:4)], 1, []), [1 3 2]);
            else
                theta = NaN.*ones(4,2);
                theta(1:goodEs,1:2) = [loc(1:goodEs); pks(1:goodEs)]';
                theta = permute(reshape(theta', 1, []), [1 3 2]);
            end
    end
else
    theta =NaN.*ones(1,1,8);
end

end

