function [ R ] = radoNan( I, theta )
%RADONAN Summary of this function goes here
%   Detailed explanation goes here

if nargin<2, theta = 1:359; end

if any(isnan(I(:)))
%     tic;    
    % irregular interpolation
    [m,n,b] = size(I); 
    [Mf,Nf] = ndgrid(1:m,1:n); % full numbers
    Out = isnan(I); % exclude NaN's from the analysis
    if sum(Out(:)) > numel(I)-4 % to many Out
        R = zeros(length(1:.5:m),length(theta));
    else
        F = scatteredInterpolant(Mf(~Out),Nf(~Out),I(~Out), ...
            'linear', ...   % interpolation method
            'none');        % extrapolation method

        [M,N] = ndgrid(1:.5:m,1:.5:n);
        center = [size(M,1)./2 size(N,2)./2];

        % construct circle
        mC = center(1).* cosd(theta) + center(1);
        nC = center(2).* sind(theta) + center(2);
        In = poly2mask(mC,nC,size(M,1),size(N,2));
        theta = theta(1:180);

        Mc = M(In); Nc = N(In); % sample grid
        center = [M(end)./2 N(end)./2];    
        R = zeros(max((Nc.*2)-1),length(theta));
    
        for i = theta
            % new rotated coordinate frame
            Mnew = (sind(i).*(Mc-center(1)) - cosd(i).*(Nc-center(2))) +center(1);
            Nnew = (-cosd(i).*(Mc-center(1)) + sind(i).*(Nc-center(2))) +center(2);

            Inew = F(Mnew,Nnew);        % interpolation information
            if numel(Inew)==0
                R(:,i) = 0;
%                 keyboard;
            else

                try
                    Nnew = interp2(Nf,Mf,double(Out),Nnew,Mnew, ...
                        'nearest');             % interpolation NaN's
                catch
                    keyboard;
                end
                Inew(isnan(Nnew)|Nnew==1) = NaN;
                % radon insensative to NaN's
                try
                    R(:,i) = splitapply(@nanmean,Inew(:),(Nc(:).*2)-1);
                catch
                    keyboard;
                end
            end
        end
    end
%     toc;
else
    % build in matlab funtion is probably faster
%     tic;
    R = radon(I);
%     toc;
end

end

