function uOut = apaf(u,T, alfa,Wmax, M)
 
 if ~exist('T')
     T=15;
 end
 
 if ~exist('alfa')
     alfa=0.9;
 end
 
 if ~exist('Wmax')
     Wmax=17;
 end
 
 if ~exist('M')
     M=13;
 end
 
zero = uint8(0);

[height width] = size(u);
%H = fspecial('average',M);
uAvg = u;%imfilter(u,H);

uOut = uAvg;
Wxy(Wmax, Wmax) = zero;
sizeWxy = numel(Wxy);

centreW = (Wmax/2) + 0.5;
offsetW = centreW - 1;

sizeC = 3:2:Wmax;
offsetC = uint8((sizeC./2) - 0.5);
countC = numel(offsetC);
bordC = sizeC * 4 - 4;

lbC = (centreW - offsetC);
ubC = (centreW + offsetC);

col1 = zeros(1, countC);
colMax = zeros(1, countC);
row1 = zeros(1, countC);
rowMax = zeros(1, countC);

pZerWek = zeros(1, countC);

for x=centreW:(height-centreW)
    for y=centreW:(width-centreW)        
        Wxy = uAvg((x - offsetW:x + offsetW), (y - offsetW:y + offsetW));
        
        I = Wxy(centreW,centreW);
        if (nnz(Wxy == I) == sizeWxy)
           continue
        end
        
        WxyM = ~(((Wxy<(I-T)) | (Wxy>(I+T))));
        
        for i=1:countC
            C = WxyM (lbC(i):ubC(i),lbC(i):ubC(i));
            
            col1(i) = sum(C(:,1));
            colMax(i) =sum( C(:,end));
            row1(i) = sum(C(1,(2:end-1)));
            rowMax(i) = sum(C(end,(2:end-1)));
           
        end
        
        pZerWek = (bordC - col1 - colMax - row1 - rowMax)./bordC;
        pZerWek(pZerWek>alfa) = NaN;
        
        if (sum(sum(isnan(pZerWek))) == countC)
            continue;
        else
            [maxPZer index] = max(pZerWek);

            subWxy = Wxy(lbC(index):ubC(index),lbC(index):ubC(index));
            maskC = WxyM(lbC(index):ubC(index),lbC(index):ubC(index));

            uOut(x,y) = mean(subWxy(maskC));
        end
    end
end