Help With Adaptive Median Filter

Could anyone help with this; I am trying to implement adaptive median filter ,and my code isn't generating the right result .

I= imread('Obama.jpg');
I=double(I);
I=I(:,:,1);
p3=0.05; %default
p4=0.95;
b = I;
x2 = rand(size(b));
d = find(x2 < p3/2);
b(d) = 0; % Minimum value
d = find(x2 >= p4);
b(d) = 1; % Maximum (saturated) value

   A=b+I;
  smax= 9;

K= zeros(size(A));
A=double(A);
[nrows ncols] = size(A);

for rows= smax:nrows-smax
for cols= smax:ncols-smax
for s=3:2:smax

inc=1;

ul=round(s/2);
for r= -ul:ul
for c= -ul:ul

region(inc)= A(rows+r,cols+c);

inc=inc+1;
end
end

kount= sort(region);
rmin= kount(1);
rmax= kount(inc-1);
rmed=kount(inc/2);
A1= rmed-rmin;
A2=rmed-rmax;

if A1>0 && A2<0  %%%go to stage B
B1= A(rows,cols)- rmin;
B2= A(rows,cols)-rmax;%
if B1>0 && B2<0
J= A(rows,cols);
else
J= rmed;
end
else
J=rmed;
end
end
%end
%end

K(rows,cols)=J;

end
end

figure (2),imshow (uint8(K))
figure (1),imshow (uint8(A))

[Image Analyst]

You forgot to attach the image. But why aren't you using medfilt2()?

[Masoni]

The idea is not to use Matlab function ,but to implement this Algorithm

    Zmin = minimum gray level value in Sxy
  Zmax = maximum gray level value in Sxy
  Zmed = median of gray levels in Sxy
  Zxy = gray level at coordinates (x, y)
  Smax = maximum allowed size of Sxy

● Algorithm Level A: A1 = Zmed - Zmin A2 = Zmed - Zmax if A1 > 0 AND A2 < 0, go to level B else increase the window size if window size < Smax, repeat level A else output Zxy Level B: B1 = Zxy - Zmin B2 = Zxy - Zmax if B1 > 0 AND B2 < 0, output Zxy else output Zmed

[Image Analyst]

So have you stepped through the code in the debugger to figure out where the wrong result is happening?

[Masoni]

Well that the point , I can't figure out why my image isn't filtered. When I reduce the code to normal median filter the code work great .

Answer 1

[Image Analyst]

You have 3 too many "for" loops in there, and too many "if" statements. Also the alpahbet soup naming of the variables didn't help - it was so confusing. I gave descriptive variable names and got rid of excess for loops and bogus if statements. It's a lot simpler and clearer now:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

I= imread('peppers.png');
I=double(I);
I=I(:,:,1);
p3=0.05; %default
p4=0.95;
b = I;
x2 = rand(size(b));
d = find(x2 < p3/2);
b(d) = 0; % Minimum value
d = find(x2 >= p4);
b(d) = 1; % Maximum (saturated) value

originalImage=b+I;
windowWidth= 9;
halfWidth = floor(windowWidth/2);
middleIndex = ceil(windowWidth^2/2);

filteredImage= zeros(size(originalImage));
originalImage=double(originalImage);
[nrows, ncols] = size(originalImage);
threshold = 5; % graylevels - whatever.

for rows= windowWidth:nrows-windowWidth
  for cols= windowWidth:ncols-windowWidth
    thisWindow= originalImage(rows-halfWidth : rows + halfWidth, cols-halfWidth : cols + halfWidth);
    sortedValuesInWindow = sort(thisWindow(:), 'Ascend');
    minValue = sortedValuesInWindow(1);
    maxValue = sortedValuesInWindow(end);
    medianValue = sortedValuesInWindow(middleIndex);
    diffFromMin = abs(medianValue - minValue);
    diffFromMax = abs(medianValue - maxValue);
    if diffFromMin> threshold || diffFromMax > threshold
      % It's noise.  Replace by the median.
      newValue = medianValue;
    else
      % It's not noise.  Use the original value.
      newValue = originalImage(rows,cols);
    end
    filteredImage(rows,cols)=newValue;
  end
end

subplot(1,2,1);
imshow (uint8(originalImage))
title('Original Image', 'FontSize', fontSize);

subplot(1, 2, 2);
imshow (uint8(filteredImage))
title('Filtered Image', 'FontSize', fontSize);

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Masoni]

Thanks ,I appreciate the help. I know I need to work on my coding skill

[sindhu k]

Can we apply this on dicom images what would be the syntax

[Image Analyst]

Instead of imread() use dicomread()

I= dicomread('ct_slice.dcm');

[DGM]

The noise generation portion from the original script is still broken.

% this is OP's attempt to add noise, but it's messed up
% the result is erroneously weak and asymmetric noise
% while easy to filter, the filtered image gets grossly truncated
p3=0.05; %default
p4=0.95;

b = I; % both I and b are uint8-scale double
x2 = rand(size(b));
d = find(x2 < p3/2); % it makes no sense to apply half if parameter is one-sided
b(d) = 0;
d = find(x2 >= p4);
b(d) = 1; % but the noise values are unit-scale
originalImage = b+I; % and the result doubles the uint8-scaling
% there's no reason to add b+I, since b is largely a copy of I

% ...
imshow(uint8(originalImage)) % and everything gets truncated
% ...
imshow(uint8(filteredImage)) % and everything gets truncated

Consider the test image:

This would yield the corresponding noisy image (originalImage) and filtered image. Here, both are scaled by half to illustrate that the "salt" noise isn't merely concealed by display truncation.

While this is obviously an incorrect noise generation routine, the filtered result has lost detail in non-noise regions. So there are multiple problems, but let's solve one at a time.

Since we're avoiding IPT, we can fix the noise generation in lieu of using imnoise().

% a image
I = imread('circuit.tif');
I = im2double(I);

% noise parameters
noisedensity = [0.05 0.05]; % noise density [pepper salt]
noisevalues = [0 1]; % noise intensity [pepper salt]

% generate noisy image
originalImage = I;
x2 = rand(size(originalImage));
mask = x2 <= noisedensity(1);
originalImage(mask) = noisevalues(1);
mask = x2 >= (1-noisedensity(2));
originalImage(mask) = noisevalues(2);

% this is basically unchanged
windowWidth = 9;
halfWidth = floor(windowWidth/2);
middleIndex = ceil(windowWidth^2/2);

filteredImage= zeros(size(originalImage));
originalImage=double(originalImage);

[nrows, ncols] = size(originalImage);
threshold = 5/255; % rescaled to unit-scale

for rows= windowWidth:nrows-windowWidth
    for cols= windowWidth:ncols-windowWidth
        thisWindow= originalImage(rows-halfWidth : rows + halfWidth, cols-halfWidth : cols + halfWidth);
        sortedValuesInWindow = sort(thisWindow(:), 'Ascend');
        minValue = sortedValuesInWindow(1);
        maxValue = sortedValuesInWindow(end);
        medianValue = sortedValuesInWindow(middleIndex);
        diffFromMin = abs(medianValue - minValue);
        diffFromMax = abs(medianValue - maxValue);
        if diffFromMin> threshold || diffFromMax > threshold
            % It's noise.  Replace by the median.
            newValue = medianValue;
        else
            % It's not noise.  Use the original value.
            newValue = originalImage(rows,cols);
        end
        filteredImage(rows,cols)=newValue;
    end
end

subplot(1,2,1);
imshow(originalImage)
subplot(1,2,2);
imshow(filteredImage)

So we've generated appropriate impulse noise, and the filter removed it. While we can ignore the lack of edge handling for a simple homework assignment, this isn't really any better than a plain median filter without any noise discrimination.

imshow(medfilt2(originalImage,[1 1]*windowWidth))

As to why the filter doesn't work, I think the purpose of this test is being confused. There are two tests in the described median filter, one which tests to see if the local median can be distinguished from the local extrema. The second test checks the current pixel against the local extrema. The first test is used to determine when to change the window size. The second test is what's used to determine whether the current pixel is noise.

diffFromMin = abs(medianValue - minValue);
diffFromMax = abs(medianValue - maxValue);

if diffFromMin> threshold || diffFromMax > threshold
    % It's noise.  Replace by the median.
    newValue = medianValue;
else
    % It's not noise.  Use the original value.
    newValue = originalImage(rows,cols);
end

The end result is that the wrong conditions are being used in the test. We're using the window-control test as the pixel-replacement test. Either way, it uses a fixed window, so it's not an adaptive filter in the manner that's suggested by these assignments.

The attached demos in @Image Analyst's other answer are working examples of fixed-window median filters, and would suffice for this example. This is the prior noisy image as processed by those demos.

This is the prior noisy image as processed using MIMT amedfilt(), which is an adaptive filter similar to the text that the question references.

op = amedfilt(inpict,5);

At this noise level, the difference between the adaptive and fixed filter is minor, and the fixed filter will be much faster.

Answer 2

[Image Analyst]

Attached are my adaptive median filter demos for removing salt and pepper noise.