function filteredImg = EPFunc( img, type, w, sigma_d, sigma_r, bShift )
img = double(img);
% Description:
%In each subwindow, the mean and
% variance are computed. The output value (located at the center of the
% window) is set to the mean of the subwindow with the smallest
% variance.
% a ab b bc c
% ah abh abc bcd cd
% h agh X cde d
% gh fgh fge fde de
% g fg fg ef e
if strcmp(type, 'kuwahara filter')
PADDING = w;
%WINSZ = w*2+1;
NoiseImgpad = padarray(img,[PADDING PADDING],'replicate');
% NoiseImgpad = [1, 2, 3, 4, 5;
% 6, 7, 8, 9, 10;
% 11, 2, 12, 14, 15;
% 16 , 17, 18 , 19, 20;
% 21 , 22 , 23, 24,35;];
[padRows,padCols] = size(NoiseImgpad);
filteredImg = zeros(size(img));
h=waitbar(0,'Please wait..');
for i = (PADDING+1):(padRows-PADDING)
waitbar(i/padRows);
for j = (PADDING+1):(padCols-PADDING)
Wwn(1:w*3+1) = NoiseImgpad(i, j);
Wn(1:w*3+1) = NoiseImgpad(i, j);
Wne(1:w*3+1) = NoiseImgpad(i, j);
We(1:w*3+1) = NoiseImgpad(i, j);
Wse(1:w*3+1) = NoiseImgpad(i, j);
Ws(1:w*3+1) = NoiseImgpad(i, j);
Wsw(1:w*3+1) = NoiseImgpad(i, j);
Ww(1:w*3+1) = NoiseImgpad(i, j);
for k = 1:w
if(k>1)
kk = (k-1)*3+1;
else
kk = k;
end
for step = 1:w-1
Wwn(1, kk:kk+2) = [NoiseImgpad(i-k, j-k),NoiseImgpad(i-k+step, j-k),NoiseImgpad(i-k, j-k+step)];
Wn(1, kk:kk+2) = [NoiseImgpad(i-k, j), NoiseImgpad(i-k, j-step), NoiseImgpad(i-k, j+step)];
Wne(1, kk:kk+2) = [NoiseImgpad(i-k, j+k), NoiseImgpad(i-k+step, j+k), NoiseImgpad(i-k, j+k-step)];
We(1, kk:kk+2) = [NoiseImgpad(i, j+k), NoiseImgpad(i-step, j+k), NoiseImgpad(i+step, j+k)];
Wse(1, kk:kk+2) = [NoiseImgpad(i+k, j+k), NoiseImgpad(i+k-step, j+k), NoiseImgpad(i+k, j+k-step)];
Ws(1, kk:kk+2) = [NoiseImgpad(i+k, j), NoiseImgpad(i+k, j-step), NoiseImgpad(i+k, j+step)];
Wsw(1, kk:kk+2) = [NoiseImgpad(i+k, j-k), NoiseImgpad(i+k, j-k+step), NoiseImgpad(i+k-step, j-k)];
Ww(1, kk:kk+2) =[NoiseImgpad(i, j-k), NoiseImgpad(i-step, j-k), NoiseImgpad(i+step, j-k)];
end
end
s = var([Wwn(:), Wn(:), Wne(:), We(:), Wse(:), Ws(:), Wsw(:), Ww(:)]);
m = mean([Wwn(:), Wn(:), Wne(:), We(:), Wse(:), Ws(:), Wsw(:), Ww(:)]);
[y,k] = min(s);%get pos in matrix
% assign the mean of the subwindow with the least variance to
% the center pixel
filteredImg(i,j) = m(k);
end
end
close(h);
elseif strcmp(type, 'bilateral filter')
img = img/255;
[X,Y] = meshgrid(-w:w,-w:w);%set position of in mask
G = exp(-(X.^2+Y.^2)/(2*sigma_d^2));%e^[-(x^2-y^2)/2*sigma_d^2], set pos value in the full mask 11*11
% Apply bilateral filter.
dim = size(img);
filteredImg = zeros(dim);
h=waitbar(0,'Please wait..');
for i = 1:dim(1)
waitbar(i/dim(1));
for j = 1:dim(2)
sumCol = 0;
num = 0;
iters = 0;
oriI = i;
oriJ = j;
iters = 0;
% if i==50 && j==100
% a =1;
% end
while 1
% Extract local region.
iMin = max(i-w,1);
iMax = min(i+w,dim(1));
jMin = max(j-w,1);
jMax = min(j+w,dim(2));
I = img(iMin:iMax,jMin:jMax);
% Compute Gaussian intensity weights.
H = exp(-(I-img(i,j)).^2/(2*sigma_r^2));%e^[-(x^2-y^2)/2*sigma_r^2]
%get small part of the marks
C = G((iMin:iMax)-i+(w+1),(jMin:jMax)-j+(w+1));
F = H.*C;
sumCol = sumCol + sum(F(:).*I(:))/sum(F(:));
num = num + 1;
[nextJ, nextI] = meshgrid(jMin:jMax,iMin:iMax);%newI column the same, newJ row the same
nextI = sum(F(:).*nextI(:))/sum(F(:));
nextJ = sum(F(:).*nextJ(:))/sum(F(:));
nextI = floor(nextI + 0.5);
nextJ = floor(nextJ + 0.5);
if bShift == 0 || (i == nextI && j == nextJ) || iters > 100
i = oriI;
j = oriJ;
filteredImg(i,j) = sumCol / num;
break;
end
i = nextI;
j = nextJ;
iters = iters + 1;
end
end
end
close(h);
filteredImg = filteredImg*255;
elseif strcmp(type, 'Conditional mean filter')
filteredImg = img;
dim = size(filteredImg);
maxShift = 3;
iterTimes = 100;
colRange = 25;
h=waitbar(0,'Please wait..');
for y = 1:dim(1)
waitbar(y/dim(1));
for x = 1:dim(2)
xc = x;
yc = y;
PVc = filteredImg(y,x);
iters = 0;
while 1
xcOld = xc;
ycOld = yc;
YcOld = PVc;
sumPV = 0;
sumX = 0;
sumY = 0;
num = 0;
for ry = -w:w
y2 = yc + ry;
if y2 >= 1 && y2 <= dim(1)
for rx = -w:w
x2 = xc + rx;
if x2 >= 1 && x2 <= dim(2)
if ry*ry + rx*rx <= w^w
PV = filteredImg(y2,x2);
dPV = PVc - PV;
if sqrt(dPV*dPV) <= colRange
num = num + 1;
sumPV = sumPV + PV;
sumX = sumX + x2;
sumY = sumY + y2;
end
end
end
end
end
end
PVc = sumPV/num;
xc = floor(sumX/num + 0.5);
yc = floor(sumY/num + 0.5);
dx = xc-xcOld;
dy = yc-ycOld;
dY = PVc-YcOld;
shift = dx*dx+dy*dy+dY*dY;
iters = iters + 1;
if bShift == 0 || shift < maxShift || iters > iterTimes %for better smoothing
break;
end
end%wh