%Implements Enis' circular buffon needle idea. Input arguments are I
%(M-by-N grayscale input image), windowSize (scalar denoting the size of
%the quadratic search window), distThresh (scalar denoting the value of maximal
%difference between pixel values, when, if not exceeded, raises the
%counter), sampling size (2-element vector or scalar denoting the sampling interbal of the image. Optional, if not defined it is set to one in x- and y- direction.)
%Outputs are scalar count, denoting all directions that on average excceed
%the distance threshold and an output image, where pixelvalues deniote the
%number of exceeded distances on average over all admissable directions.
function [resI]=circularNeedleAnglesSameSize7(I,windowSize,distThresh,samplingSize)
if nargin==3
    samplingSize=[1 1];
end
orImSize=size(I);
if length(orImSize)~=2
    error('Input image is not of size M-by-N')
end
if length(samplingSize)==1
    samplingSize=[1 1]*samplingSize;
end
%error threshold for finding number of admissable angles
errorThresh=2;
%check if window size is odd
if rem(windowSize,2)~=1
    windowSize=windowSize+1;
end
radius=floor(windowSize/2);
%I=uint8(appendMirrorImage(double(I),windowSize));
I=appendMirrorImage(double(I),windowSize);
imSize=size(I);
%compute acceptable angles for the image sizes
%set initial value of smallest angle

alpha=pi/4;
flag=0;
newAlpha=alpha;
while flag==0
    oldAlpha=newAlpha;
    newAlpha=oldAlpha/2;
    %compute corrdinates of new angle
    [newX,newY]=pol2cart(newAlpha,radius);
    newX=round(newX);
    newY=round(newY);
    %compute corrdinates of old angle
    [oldX,oldY]=pol2cart(oldAlpha,radius);
    oldX=round(oldX);
    oldY=round(oldY);
    %if the difference between the new and old coordinates is smaller than errorThresh-many
    %pixels, abort while
    if abs(newX-oldX)<=errorThresh && abs(newY-oldY)<=errorThresh
        %set final iteration value to previous value
        alpha=oldAlpha;
        flag=1;
    end
end

angles=(0:alpha:2*pi-alpha);
diff=zeros(1,radius);

aa=1:length(angles);
rr=1:radius;
rrl = length(rr);
[colTable,rowTable] =  pol2cart(repmat(angles',[1,length(rr)]),repmat(rr,[length(angles),1]));
indx = mod(angles(aa)/(pi/4),2)==1;
colTable(indx,:)=sign(colTable(indx,:)).*repmat(rr,[sum(double(indx)),1]);
rowTable(indx,:)=sign(rowTable(indx,:)).*repmat(rr,[sum(double(indx)),1]);

%visualize results in image
%resI=zeros(length(1:samplingSize(1):orImSize(1)),length(1:samplingSize(2):orImSize(2)),length(angles));
resI=zeros(length(1:orImSize(1)),length(1:orImSize(2)),length(angles));

m2=0;
%for all rows of I

for mm=radius+1:samplingSize(1):imSize(1)-radius
    m2=m2+1;
    n2=0;
    %for all columns of I
    for nn=radius+1:samplingSize(2):imSize(2)-radius
        n2=n2+1;
        
        indx = (mm-round(rowTable))';
        indy = (nn+round(colTable))';
        indx = indx(:);
        indy = indy(:);
        
        diff1 = abs(I(mm,nn) - I(indx+imSize(1)*(indy-1)));
        diff1 = reshape(diff1,[radius,length(angles)])';
        if (m2*samplingSize(1)>orImSize(1)) || (n2*samplingSize(2)>orImSize(2))
            resI(((m2-1)*samplingSize(1))+1:orImSize(1),(n2-1)*samplingSize(2)+1:orImSize(2),aa)=sum((diff1>distThresh),2);
        else
            resI(((m2-1)*samplingSize(1))+1:(m2*samplingSize(1)),(n2-1)*samplingSize(2)+1:n2*samplingSize(2),aa)=sum((diff1>distThresh),2);
        end
    end
end
end

function outI=appendMirrorImage(inI,filtSize)
appendNo=floor(filtSize/2);
imSize=size(inI);
outI=zeros(imSize(1)+2*appendNo,imSize(2)+2*appendNo);
%filter image, horizontal
for mm=1:size(inI,1)
    outI(mm+appendNo,:)= cat(2,inI(mm,appendNo:-1:1),inI(mm,:),inI(mm,end:-1:end-appendNo+1)) ;
end
%filter image, vertical
for nn=1:size(inI,2)
    outI(:,nn+appendNo)=cat(1,inI(appendNo:-1:1,nn),inI(:,nn),inI(end:-1:end-appendNo+1,nn));
end
%for the outer squares (that are black after the previous code),
%compute their values as the average of two values
for ii=1:appendNo
    outI(ii,1:appendNo)=0.5*(inI(appendNo+1-ii,1)+inI(1,1:appendNo));
    outI(ii,end-appendNo+1:end)=0.5*(inI(appendNo+1-ii,end)+inI(1,end:-1:end-appendNo+1));
    outI(end-appendNo+ii,1:appendNo)=0.5*(inI(end+1-ii,1)+inI(end,1:appendNo));
    outI(end-appendNo+ii,end-appendNo+1:end)=0.5*(inI(end+1-ii,end)+inI(end,end:-1:end-appendNo+1));
end
end