Order Statistics Filters

Here, we apply the median filter to the corrupted image by salt and peper. To get salt and paper I used matlab function I = imnoise(A,'salt & pepper',0.3) for 30% and J = imnoise(A,'salt & pepper',0.5) for 50% corruption. Again, the dollar image is used as the original image.

Implementation of the median filter is pretty much the same with spatial filtering. The changes are in zero padding, start and end points of the for loops because now we made our matrix size (row*4 X cols*4) instead of multiplying by 2 and in the helper “Filter” function. “Filter” is now named as “MyMedian” which actually takes the original image and the desired index that will be needed for finding the median of the 5X5 matrix. This function takes 5X5 neighboured indices, holds them in array and evaluate the median.

Let’s have look at the MedianFilter function. The actual work done is implemented as below after zero padding. Note that the values in the for loops are re arranged.

%zero padding
%first, new matrix is filled with 0
for k = 1: row + 4
    for l=1 : cols + 4
      C(k,l) = 0;
    end
end
%then values of the input is saved to output
%with the desired place to satisfy zero padding
for k = 1 : row
    for l = 1 : cols
      C(k+2,l+2) = A(k,l);
    end
end

%here, the actual work is done
%I wrote another function called 
%"MyMedian" that carries out median for
%only 5x5 part of the input matrix here 
%the 2 for loops make median filtering
%for all the possible non zero values of the C
for k = 3 : row + 2
    for l = 3 : cols + 2
      value = MyMedian(C,k,l);
      %a temp matrix is used not to override C and change its values
      %after a 5x5 convolution of it
      Temp(k,l) = value;
    end
end

Instead of equating value to Filter function now MyMedian function that evaluates the median is used. Theb the Temp is equated to C.

%temp is saved in C
for k = 3 : row + 2
    for l = 3 : cols + 2
      C(k,l) = Temp(k,l);
    end
end

Now, we can look at the MyMedian function.

for m = x-2: x + 2
    for n = y-2: y + 2
      k = k + 1;
      %the value of A saved in a matrix B
      B(k) = A(m,n);
    end
end

As said before it finds the median of any 5X5 matrix. To do this this matrix is saved in B then the median function of the matlab is used.

30 % Salt & Pepper Noised Image

Applying the Median Filter

Spatial Filtering

Spatial filtering is performed by convolving the image with a mask so that the new image g(x,y) is evaluatedas the following convolution formula.

In order to do that a matlab code called SpatialFilter is implemented. First, we need to do zero padding. I started to fill the output image with full of zeros as the following.

%C is the output image
%A is the original image
%B is the filter
function [C] =SpatialFilter(A,B)

for k = 1: row + 2
    for l=1 : cols +2
      C(k,l) = 0;
   end
end

Notice that 2 is added to end points and the new image’s size became row+2 X cols+2 where the row and cols is the input images’ row and column values. After that,the pixel values of the input image is saved on the output image C with the desired places which starts from the point (2,2) to (row,cols). Also notice that it isthe size of the input image.

for k = 1 : row
    for l = 1 : cols
      C(k+1,l+1) = A(k,l);
   end
end

So, our input image A is zero padded and saved in the output image C.

%here, the actual work is done

for k = 2 : row +1
    for l = 2 : cols +1
      value = Filter(C,B,k,l) / sumx;

      Temp(k,l) = value;
    end
end

Here, I used another function that I wrote called Filter. This function only does the convolution for the given index values that is a 3x3 matrix. After this partial convolution, the result is saved in a temporary matrix. The reason why I used a temporary matrix is that if I saved the new values in C then the old values of C would change and it would affect the result of the convolution. The explained producere is repeated for all the pixels in the zero padded input array A. See that I didn’t convolve the zero padded values in order to that for loops strats from 2. Finally, the Temp is saved in the output matrix C.

 
for k = 2 : row + 1
    for l = 2 : cols + 1
      C(k,l) = Temp(k,l);
    end
end   

Now, lets have look at Filter function. (The functionality of this function is explained above.)

 
% m an n is the neighboured indices of
%the given index so I traverse all
% the neighbours of the given
%index through the below loop
% k and l is used for indices of the filter

for m = x-1: x + 1
k = k +1;
    for n = y-1: y + 1
       l = l + 1;
      filter = filter +(A(m,n) * B(k,l));
    end
l = 0;
end

The for loops does the convolution of the wanted index with the filter by multiplying the given index and neighbour values with the filter and then summing them up. For filtering the dollor image is used. Results of using bluring and spatial filters are shown below.

Original 1 Dollar Image

Blurred 1 Dollar

After Laplacian Filter is Applied

Codes for spatial filter and filter are below.

%A is the image
%B is the filter
%x and y is the the indices that will convolve with the filter
%multiplies 3 x 3 neighbours of the given index of the image with 3 x 3 
%filter with corresponding indices and sum them
function [filter] = Filter(A,B,x,y)
filter = 0;
%used to take indices of the filter
k = 0;
l = 0;
% m  an n is the neighboured indices of the given index so I traverse all
% the neighbours of the given index through the below loop
% k and l is used for indices of the filter
for m = x-1: x + 1   
    k = k +1;
    for n = y-1: y + 1        
        l = l + 1;
        filter = filter +(A(m,n) * B(k,l));
    end
    l = 0;
end
end

%c is the output image
%A is the original image
%B is the filter
function [C] =SpatialFilter(A,B)
%input is converted to double
A = double(A);
%Initialize the row & columns
[row,cols] = size (A);
%sumx is the total sum of the filter
sumx = 0;

%zero padding
%first, new matrix is filled with 0
for k = 1: row + 2
    for l=1 : cols +2
    C(k,l) = 0;
    end
end
%then values of the input is saved to output 
%with the desired place to satisfy zero padding 
for k = 1 : row
    for l = 1 : cols    
        C(k+1,l+1) =  A(k,l);
    end
end

%sum of the filter index values is evaluated
for i = 1:3
    for j = 1:3
        sumx = sumx + B(i,j);
    end
end

%here, the actual work is done
%I wrote another function called "Filter" that carries out convolution for 
%only 3x3 part of the input matrix here the 2 for loops make this filtering
%for all the possible non zero values of the C 
for k = 2 : row +1
    for l = 2 : cols +1   
        value = Filter(C,B,k,l) / sumx;
        %!!!don't divide sumx when evaluating laplacian!!!
        %value = Filter(C,B,k,l);
        %a temp matrix is used not to override C and change its values
        %after a 3x3 convolution of it
        Temp(k,l) = value;        
    end
end
%temp is saved in C
for k = 2 : row + 1
    for l = 2 : cols + 1       
        C(k,l) = Temp(k,l);       
    end
end
end

Histogram & Equalization

In the following code we calculate the normalized histogram with the histogram equlization of an image. 

Before showing the whole code let's look at the content of the code.

Normalized Histogram

p(rk) = nk / n where n is the total numbers of pixel and nk is the number of pixels having gray level rk. In our gray level images number of gray levels is taken as 255 and probability of this gray levels frequencies is evaluated.

% for loops to take all indicies of input A

for i=1:rows

    for j=1:cols

      %value holds the matrix value at i,j indicies

      value=A(i,j);

      %occur holds the occurrrnce of each value in the input

      %since it is a matrix and the index of the matrix is equal to the

      %value, it calculates the occurances through the for loops

      %1 is added to index because 0 can't be an index

      occur(value+1)=occur(value+1)+1;

      %calculates the normalized histogram by dividing the occurrence of

      %each value by the pixel number

      %1 is added to index because 0 can't be an index

      histogram(value+1)=occur(value+1)/pixelnumbers;

  end

end

In matlab the histogram matrix is evaluated with the above two for loops which travels all the indexes of the input matrix A. Then it calculates the occurrences of each gray levels by simply doing occur(value+1)=occur(value+1)+1 where occur(10) = 100 means that the gray level value 10 has occurred 100 times. So, I get the occurrences of the pixels. Finally, histogram matrix is evaluated by dividing occurrences to number of pixels which gives the normalized histogram. You can see the normlized histogram of the below image.

Normalized histogram of above image

Histogram Equalization

In histogram equalization we need to find the cumulative distribution function of pixel occurrence. In the solution 9, I already found the occurrences. Now, I use occurrences and the histogram function to obtain equalized histogram function. Then, I multiply the new function with the number of bins that is 255. So, I get the equalized histogram.

% for loops to take all indicies of the histogram for given input
for i=1:size(histogram)
sum = sum+occur(i);
%repeatedly we add the number of occurrences to ongoing ones
cum(i)=sum;
%histogram equalization is done here
equalized_histogram(i)=round((cum(i)/pixelnumbers)*255);
end

After that I need to get new image that is histogram equalized. So, I replace the values in the input array from the values from equalized histogram and save it to the output. For example, say that I have the value x in the input array A (1,1) and the cdf of A (1,1) is y. Then, I make A(1,1)’s value y. The following for loops do that.

% for loops to take all indicies of input A

for i=1:rows

    for j=1:cols

      % replace the values of the input with the values from histogram

      % equalization and save it to the B 

      %1 is added to index because 0 can't be an index

     B(i,j)=equalized_histogram(A(i,j)+1);

   end

end

The whole code of the program is below. You can use the code using bar command to draw histograms as the following. Comments are included for better understanding.

 X = imread('Image name');

[b,e,h] = MyHistogram(X);

bar(e);

%Normalized Histogram && Histogram Equalization

%Takes the image as matrix
%Outputs the Equalized Image,Equalized Histogram and Normalized Histogram
%as matrices

function [B,equalized_histogram,histogram] = MyHistogram(A)

%input is converted to double
A = double(A);
%row & columns of A
[rows,cols] = size(A);
%pixel number is found
pixelnumbers = rows*cols;
%Initialize the histogram equalized output image with full of unit8 zeros
%otherwise we need to use imread(matrix,[]) instead of imread([])
B=uint8(zeros(rows,cols));
%occur matrix is used to hold occurrences of each value
occur=zeros(256,1);

%histogram matrix is used for normalized histogram of the input image
histogram =zeros(256,1);
%equalized_histogram matrix is used for equalized histogram of the input
equalized_histogram=zeros(256,1);
%cum is used to hold cumalitives for equalized histogram
cum=zeros(256,1);
%sum is set to 0
sum=0;

% for loops to take all indicies of input A
for i=1:rows
    for j=1:cols
        %value holds the matrix value at i,j indicies
        value=A(i,j);
        %occur holds the occurance of each value in the input
        %since it is a matrix and the index of the matrix is equal to the
        %value, it calculates the occurances through the for loops
        %1 is added to index because 0 can't be an index 
        occur(value+1)=occur(value+1)+1;
        %calculates the normalized histogram by dividing the occurance of
        %each value by the pixel number
        %1 is added to index because 0 can't be an index 
        histogram(value+1)=occur(value+1)/pixelnumbers;

    end
end

%bar(histogram);

% for loops to take all indicies of the histogram for given input
for i=1:size(histogram)
   sum = sum+occur(i);
   %repeatedly we add the number of occurrences to ongoing ones
   cum(i)=sum; 
   %histogram equalization is done here
   equalized_histogram(i)=round((cum(i)/pixelnumbers)*255);
end

% for loops to take all indicies of input A
for i=1:rows
    for j=1:cols 
        % replace the values of the input with the values from histogram
        % equalization and save it to the B
        %1 is added to index because 0 can't be an index 
        B(i,j)=equalized_histogram(A(i,j)+1);        
    end
end

%bar(equalized_histogram);

end

Bit Plane Slicing

Intensities is formed by 8 bits so we have 256 levels of gray. We will use this bits for bit slicing so that the image will be sliced according to selected bit. For example, say that I want to slice the 8th bit then, I will take 1000000 (128 in decimal) and make the bitwise operation with the corresponding image. If 1st bit as 00000001 (1 in decimal) is taken than the resulting image will be full of dots because I take the least significant bit in the bitwise operation and the most significant bit when I choose 10000000 in bit slicing. So 8th bit will give me an image close to original one. See the following code.

 %Takes the image and the bit in decimal as arguments
%Decimal arguments for the bit numbers:1,2,4,8,16,35,64,128
%Insert 192 for 11000000, 8th and 7th bits. 224 for 11100000 and 240 for 11110000.
 for i=1:row
   for j=1:cols
   %Bitwise operation is handled
   B(i,j) = double(bitand(A(i,j),bit));
   end
end

Here, matlab bitwise operation, which takes its arguments in decimal, is used to take the desired bit, and that bit is assigned to output image. See the BitSlice function. Resulting image on 1 dollar image is as the following.

1st bit is sliced

4th bit is sliced

6th bit is sliced

8th bit is sliced

So, as you can see whenever I slice the more significant bit, the image is about to disappear.

Contrast Stretching

In contrast strertching we expand the range of intesity levels so that it spans the full intensity range by applying the following.

function [B] = ContrastStretching(A,min,max)
%input is converted to double
A = double(A);
%Initialize the row & columns
[row,cols] = size (A);
for k = 1 : row
    for l = 1 : cols
      %Do the following according to comparission
      %Apply the equations
      if A(k,l) <= min
      B(k,l) = 0;
        else if  A(k,l) > min && A(k,l) < max
        B(k,l) = (A(k,l)-min) / (max -min);
          else if A(k,l) >= max
          B(k,l) =1;
          end
        end
      end            
    end
end
end

Washed out pollen image

Corrected pollen image

Gamma Transformation and Correction

For this image enhancement we use the following game transformation s = c r^γ where c and γ are positive constants. As the formula implies, we hava a straight line when γ is 1. When γ is less than 1 we have high intensity level for the corresponding input intensity level which means we have a lighter image and vice versa when γ is greater than 1.

As we did in previous example, we change the inside of nested for loops so that it becomes B(k,l) = c*(A(k,l).^y); .

when y = 0,6

   

when y = 0,3

Suppose I have this washed out image here. To make it more visible, I can use gamma transformation so that I get darker images with y = 4 and y =5.

y = 4

Log Transformation

We, now try to apply s = c log(1+r) transformation to an image. It's smillar to negative of an image. The only difference is in the for loops. 

We should replace line B(k,l) = (L-A(k,l)); with B(k,l) = c*log(1 + (A(k,l))); which applies the log transform that converts narrow range of low intensity values to a wider range intensity values. Notice that in the log function I have high values of output for the low values of the input and the oppsite for inverse log function.

function [B] =LogTransform(A)
A = double(A);
[row,cols] = size (A);
%Function coefficient
c = 1;
  for k = 1 : row
    for l = 1 : cols
        %Outputs from the log function is
        %assigned to output array
        B(k,l) = c*log(1 + (A(k,l)));
    end
  end
end

Before Log Transform

After Log Transform