Adaptive median filter

clc;
clear all;
close all;
a=imread('cameraman.tif');
figure(1),imshow(a),title('origanal image');
b=imnoise(a,'salt & pepper',.02)
figure(2),imshow(b),title('noisy image');
Smax=9;

for i=1:254
    for j=1:254
        n=b(i:i+2,j:j+2)
        Zmin=min(n(:));
        Zmax=max(n(:));
        Zmed=median(n(:));
        sx=3;
        sy=3;

                A1=Zmed-Zmin;
                A2=Zmed-Zmax;
              if (A1>0) && (A2<0)
                  B1 = Zxy-Zmin;

                  B2 = Zxy-Zmax;
                  if (B1>0) && (B2<0)
                        b(i:i+2,j:j+2) = n(i, j);
                      break;
                  else
                       b(i:i+2,j:j+2) = Zmed;
                      break;
                  end 
              else 
                  sx=sx+2;
                  sy=sy+2;
                  if (sx > Smax) && (sy > Smax)
                      b(i:i+2,j:j+2) = n(i, j);
                  end 
                  end 
      end 
      end 

figure(3),imshow(b),title('denoised image');

i am trying to implement the following algorithm but where the above pro. is wrong i am unable to find ..

please help me.......

Adaptive median filter changes size of Sxy (the size of the neighborhood) during operation.

● Notation
  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

[Ian Wildemann]

I figuired it out. After the last if block, add an else block that has b(i,j) = zmed. Worked for me.

[ciprian goia]

can you explain how?

[DGM]

For a clearer picture of the algo/pseudocode posted above, see this PDF.

http://site.iugaza.edu.ps/ayash/files/Lab9_Adaptive-Filters-Color-Image-Processing.pdf

For other implementations of varying efficiency, see:

https://www.mathworks.com/matlabcentral/fileexchange/30068-adaptive-median-filter-matlab-code (no border handling)

https://github.com/Ali478/Matlab-Adaptive-Median-Filter (includes border handling)

https://github.com/mortezamg63/Adaptive-Median-Filter (slow, less effective, performs poorly with small smax)

[pranali ukey]

It looks like in the above code Zxy is undefined. I tried to run this code and the filtered image is still having more noise content in it.

[DGM]

The above code is incomplete and nonfunctioning.

I posted information about a complete implementation based on the same pseudocode in the answer below.

[pranali ukey]

Hi, I'm trying above code on the cameraman image but the quality with this filter is not so good when compared with other mean filter. It looks like the image is not getting filtered but it is showing the noised image at the output. If anyone has tried that with any image please answer.

[DGM]

I already did answer.

The above code is incomplete and nonfunctioning. It doesn't do anything.

a = imread('cameraman.tif'); % this is uint8
b=imnoise(a,'salt & pepper',.02);

c = b; % if you overwrite the input, then you can't compare anything
Smax=9;
for i=1:254 % don't embed image sizes in the code
    for j=1:254 % just don't
        n=b(i:i+2,j:j+2); % sample the input image, not the output image!
        Zmin=min(n(:)); % this inherits its class from b
        Zmax=max(n(:)); % this inherits its class from b
        Zmed=median(n(:)); % this inherits its class from b
        sx=3;
        sy=3;
        % these also inherit their class from b
        % if b is uint8, A2 will always be non-negative
        A1=Zmed-Zmin; % should always be >= 0
        A2=Zmed-Zmax; % should always be <= 0, but will actually always be exactly 0
        if (A1>0) && (A2<0) % if the input is uint8, this will never be true
            % noise discrimination
            % this will never be entered
            % if it does, it will immediately throw an error
            % if the error is fixed, it does nothing.
            B1 = Zxy-Zmin;
            B2 = Zxy-Zmax; % this will be truncated for uint8 inputs for the same reason
            if (B1>0) && (B2<0) % if the input is uint8, this will never be true
                % the LHS of this assignment should be a scalar
                % even if this could be reached, it would cause an error
                % since i,j will be out-of-range subscripts for the size of n.
                % even if that were changed to c(i,j) = n(1,1)
                % the assignment would be completely redundant
                % since the output is already a copy of the input
                c(i:i+2,j:j+2) = n(i, j);
                break; % this would cause entire chunks of rows to get skipped
            else
                % the LHS of this assignment should be a scalar
                c(i:i+2,j:j+2) = Zmed;
                break; % this would cause entire chunks of rows to get skipped
            end
        else
            % this is ostensibly the window resizing routine
            % but the window is never actually resized
            % and the region stats are never recalculated
            sx=sx+2; % this will always be 5
            sy=sy+2; % this will always be 5
            if (sx > Smax) && (sy > Smax) % this will never be true
                % even if this could be reached, it would cause an error
                % for the same reason as before
                % this assignment should still be prefaced by a noise 
                % discrimination test anyway
                c(i:i+2,j:j+2) = n(i, j);
            end
        end
    end
end

% the array c was never changed
% the output is identical to the input
% simply because that's how c was created in the first place
isequal(b,c)

ans = logical
   1

I'm not going to rewrite old broken half-implemented code when I've already provided a working implementation of the same algorithm.

% the same image
a = imread('cameraman.tif');
b = imnoise(a,'salt & pepper',0.1);

c = amedfilt(b); % included in MIMT

imshow([b c],'border','tight')

[DGM]

If you want something divorced from MIMT conveniences and without the extra functionality I added, then:

% the same image
a = imread('cameraman.tif');
b = imnoise(a,'salt & pepper',0.1);

c = amedfiltsimple(b); % read the synopsis

imshow([b c],'border','tight')

Answer 1

[Image Analyst]

You need to take the contents of the "if" block and make it a function that takes the center location and the window size as input arguments. Then in the "else" block you need to increase the window size and call the function. Something like

% Determine window size.
if (A1>0) && (A2<0)
  windowSize = standardSize; % 3 or whatever
else
  windowSize = % your formula for bigger window
end
% Make the assignment.
newImage(x,y) = GetMedian(x,y,windowSize)

[Image Analyst]

You call checkMedian like this:

result = checkMedian;

but the definition of checkMedian is this

function  checkMedian()
% more code

Exactly how is checkMedian() supposed to return a value that you can stuff into result if your function definition says it does not return anything? It can't. Define it to return an output argument if you want something to be returned.

[Ian Wildemann]

How does making this function help? Isn't the original code already doing these things without the function?

Answer 2

[DGM]

Synchronicity is strange. After seeing this one question a month ago, I kept seeing similar threads in the following days. There are a lot of incomplete, poorly-working or poorly-documented implementations to be found. One of the best I found was this one. I figured I might as well put one together that works and is readable.

MIMT amedfilt() is similar to the concept described in the above pdf, though I've changed it to have broader utility. For the sake of comparison, I also implemented a more typical fixed-window median noise removal filter, fmedfilt().

Let's start with an extremely noisy image.

% an image
inpict0 = imread('circuit.tif');

% with a bunch of noise
inpict = imnoise(inpict0,'salt & pepper',0.6);

% use fixed-window filter and defaults
% here, the second parameter is the filter width
op1 = fmedfilt(inpict,11);

% use adaptive filter
% here, the second parameter is the maximum filter radius
op2 = amedfilt(inpict,5);

While this is an extreme case, we can see that both filters work. The adaptive filter tends to be better at retaining contrast. The fixed filter results are good, but a bit noisier, especially around edges.

The fixed filter's defaults are similar to other good examples shared on the forum. The noise discrimination is a simple test to see whether the current pixel is saturated (e.g. at either 0 or 255 for uint8). If for some reason, (e.g. some intermediate processing or compression), the noise pixels are not completely saturated, they will simply not be replaced. While this can be solved by adding some sort of tolerance parameter to the fixed filter, the adaptive filter is unaffected.

% make the noise pixels less extreme (by 1%)
mk = inpict~=inpict0;
inpict(mk) = imadjust(inpict(mk),[0 1],[0.01 0.99]);

% fails completely with tight tolerance
op1 = fmedfilt(inpict,11);

% wider tolerance catches noise
op2 = fmedfilt(inpict,11,0.015); 

% amedfilt() is unaffected
op3 = amedfilt(inpict,5);

As to the manner in which amedfilt() deviates from the prescribed algorithm or expectations, there are two things to mention. First is that since I've included these in MIMT, both of these tools are written to be independent of Image Processing Toolbox, but they may rely on MIMT tools instead. I leave it as an exercise to the student to figure that out.

The second alteration was a lazy attempt to make the filter more robust against lossy compression. It's not uncommon for test images to be transported in wholly inappropriate formats like JPG. If an image with this sort of impulse noise is compressed, it's quite a bit more difficult to filter than the prior example might suggest. This change introduces an extra parameter which controls which order statistic is used in the window adjustment and noise discrimination tests. By default, this is the local min/max, but can be adjusted. Since the window size is variable, these are specified in a modified relative manner.

% an image
inpict0 = imread('circuit.tif');

% add noise
inpict = imnoise(inpict0,'salt & pepper',0.1);

% subject the image to lossy compression
% this is a MIMT convenience tool that just saves the image and reads it back
inpict = jpegger(inpict,70);

% requires wide tolerance; noisy
op1 = fmedfilt(inpict,11,0.15);

% works, but performs poorly
op2 = amedfilt(inpict,5);

% improved
op3 = amedfilt(inpict,5,1);

% compare results against original
err1 = immse(op1,inpict0)
err2 = immse(op2,inpict0)
err3 = immse(op3,inpict0)

err1 =

140.1061

err2 =

329.9459

err3 =

43.4217

I'm sure there are other canonical approaches to handle such an image, but my interest in this additional feature was utilitarian and largely unrelated to the original motive. Once I had a working noise removal filter, I found myself wondering how long it would be before I needed to apply it to a screenshot that someone posted to the forum. For that, I'm sure my modification will suffice.