

function [ f ] = pocs(g,h,t,n)

f = g;

[height, width] = size(f);

m = (size(h,1)-1)/2;
h_norm = sum(sum(h.*h));

for i = 1:n
    for x = 2+m:width-m-1
        for y = 2+m:height-m-1
            
            HF = sum(sum(h.*f(y-m:y+m,x-m:x+m)));
                        
            if HF > (g(y,x)+t)
                f(y-m:y+m,x-m:x+m) = f(y-m:y+m,x-m:x+m) + h.*(g(y,x)-HF+t)/h_norm;
            elseif HF < (g(y,x)-t)
                f(y-m:y+m,x-m:x+m) = f(y-m:y+m,x-m:x+m) + h.*(g(y,x)-HF-t)/h_norm;
            end
        end
    end
    f(f>255) = 255;
    f(f<0) = 0;
    %f = double(uint8(f));
end

