Image Sharpening using second order derivative –(Laplacian)

Prerequisite: Read EdgeDetection- fundamentals

The derivative operator Laplacian for an Image is defined as

For X-direction,

For Y-direction,

By substituting, Equations in Fig.B and Fig.C in Fig.A, we obtain the following equation

The equation  represented in terms of Mask:

0	1	0
1	-4	1
0	1	0

When the diagonals also considered then the equation becomes,

The Mask representation of the above equation,

1	1	1
1	-8	1
1	1	1

Now let’s discuss further how image sharpening is done using Laplacian.

Equation:

Where f(x,y)  is the input image

              g(x,y) is the sharpened image and

               c= -1 for the above mentioned filter masks.(fig.D and fig.E)

MATLAB CODE:

%Input Image

A=imread('coins.png');

figure,imshow(A);

%Preallocate the matrices with zeros

I1=A;

I=zeros(size(A));

I2=zeros(size(A));

%Filter Masks

F1=[0 1 0;1 -4 1; 0 1 0];

F2=[1 1 1;1 -8 1; 1 1 1];

%Padarray with zeros

A=padarray(A,[1,1]);

A=double(A);

%Implementation of the equation in Fig.D

for i=1:size(A,1)-2

    for j=1:size(A,2)-2

       

        I(i,j)=sum(sum(F1.*A(i:i+2,j:j+2)));

       

    end

end

I=uint8(I);

figure,imshow(I);title('Filtered Image');

The Laplacian derivative equation has produced grayish edge lines and other areas are made dark(background)

%Sharpenend Image

%Refer Equation in Fig.F

B=I1-I;

figure,imshow(B);title('Sharpened Image');

The Filter Image is combined with the Original input image thus the background is preserved and the sharpened image is obtained .

For a Filter Mask that includes Diagonal,

1	1	1
1	-8	1
1	1	1

Filtered Image:

Sharpened Image:

Edge Detection-Fundamentals

The derivatives of a digital function are defined in terms of differences.

The above statement made me to analyze about derivatives and how it is used for edge detection.  The first time when I came across the edge detection operation [Example: edge(Image,’sobel’)], I wondered how it worked. 

Consider a single dimensional array,

A =

5

4

3

2

2

2

2

8

8

8

6

6

5

4

0

MATLAB CODE:

x=1:15;

y=[5 4 3 2 2 2 2 8 8 8 6 6 5 4 0];

figure,

plot(x,y,'-o','LineWidth',3,'MarkerEdgeColor','k','Color','y');

title('Input Array');

First-order Derivative for one dimensional function f(x):

MATLAB  CODE:

x1=1:14;

y1=diff(y,1);

figure,

plot(x1,y1,'-o','LineWidth',3,'MarkerEdgeColor','k','Color','r');

-1

-1

-1

0

0

0

6

0

0

-2

0

-1

-1

-4

NOTE: The contiguous values are zero. Since the values are nonzero for non-contiguous values, the result will be thick edges.

The first-order derivative produces thicker edges.

Second-order Derivative for one dimensional function f(x):

MATLAB CODE:

x2=1:13;

y2=diff(y,2);

figure,

plot(x2,y2,'-o','LineWidth',3,'MarkerEdgeColor','k','Color','g');

0

0

1

0

0

6

-6

0

-2

2

-1

0

-3

The Second-order derivative gives finer result compared to first-order derivative. It gives fine detailed thin lines and isolated points.  Let’s see how the second-order derivative used for Image sharpening (Laplacian) in my upcoming post.

Convert HSI Image to RGB Image

CONVERT HSI IMAGE TO RGB IMAGE:

Refer How to convert the RGB image to HSI Image

MATLAB CODE:

 figure,imshow(HSI);title('HSI Image');  
 %Obtain the Hue, Saturation and Intensity components  
 H1=HSI(:,:,1);  
 S1=HSI(:,:,2);  
 I1=HSI(:,:,3);  
    
   

 %Multiply Hue by 360 to represent in the range [0 360]  
 H1=H1*360;                                               
    


 %Preallocate the R,G and B components  
 R1=zeros(size(H1));  
 G1=zeros(size(H1));  
 B1=zeros(size(H1));  
 RGB1=zeros([size(H1),3]);  
    


 %RG Sector(0<=H<120)  
 %When H is in the above sector, the RGB components equations are  


    
 B1(H1<120)=I1(H1<120).*(1-S1(H1<120));  
 R1(H1<120)=I1(H1<120).*(1+((S1(H1<120).*cosd(H1(H1<120)))./cosd(60-H1(H1<120))));  
 G1(H1<120)=3.*I1(H1<120)-(R1(H1<120)+B1(H1<120));  


    
 %GB Sector(120<=H<240)  
 %When H is in the above sector, the RGB components equations are  


    
 %Subtract 120 from Hue  
 H2=H1-120;  


    
 R1(H1>=120&H1<240)=I1(H1>=120&H1<240).*(1-S1(H1>=120&H1<240));  
 G1(H1>=120&H1<240)=I1(H1>=120&H1<240).*(1+((S1(H1>=120&H1<240).*cosd(H2(H1>=120&H1<240)))./cosd(60-H2(H1>=120&H1<240))));  
 B1(H1>=120&H1<240)=3.*I1(H1>=120&H1<240)-(R1(H1>=120&H1<240)+G1(H1>=120&H1<240));  


    
 %BR Sector(240<=H<=360)  
 %When H is in the above sector, the RGB components equations are  


    
 %Subtract 240 from Hue  
 H2=H1-240;  


    
 G1(H1>=240&H1<=360)=I1(H1>=240&H1<=360).*(1-S1(H1>=240&H1<=360));  
 B1(H1>=240&H1<=360)=I1(H1>=240&H1<=360).*(1+((S1(H1>=240&H1<=360).*cosd(H2(H1>=240&H1<=360)))./cosd(60-H2(H1>=240&H1<=360))));  
 R1(H1>=240&H1<=360)=3.*I1(H1>=240&H1<=360)-(G1(H1>=240&H1<=360)+B1(H1>=240&H1<=360));  


    
 %Form RGB Image  
 RGB1(:,:,1)=R1;  
 RGB1(:,:,2)=G1;  
 RGB1(:,:,3)=B1;  


    
 %Represent the image in the range [0 255]  
 RGB1=im2uint8(RGB1);  
 figure,imshow(RGB1);title('RGB Image');  
    
   

  

Explanation

Converting RGB Image to HSI

Converting RGB Image to HSI

H stands for Hue, S for Saturation and I for Intensity.

MATLAB CODE:

Read a RGB Image

A=imread('peppers.png');

figure,imshow(A);title('RGB Image');

%Represent the RGB image in [0 1] range

I=double(A)/255;

R=I(:,:,1);

G=I(:,:,2);

B=I(:,:,3);

%Hue

numi=1/2*((R-G)+(R-B));

denom=((R-G).^2+((R-B).*(G-B))).^0.5;

%To avoid divide by zero exception add a small number in the denominator

H=acosd(numi./(denom+0.000001));

%If B>G then H= 360-Theta

H(B>G)=360-H(B>G);

%Normalize to the range [0 1]

H=H/360;

%Saturation

S=1- (3./(sum(I,3)+0.000001)).*min(I,[],3);

%Intensity

I=sum(I,3)./3;

%HSI

HSI=zeros(size(A));

HSI(:,:,1)=H;

HSI(:,:,2)=S;

HSI(:,:,3)=I;

figure,imshow(HSI);title('HSI Image');

Explanation:

1.    Read a RGB image using ‘imread’ function.

2.    Each RGB component will be in the range of [0 255].  Represent the image in [0 1] range by dividing the image by 255.

3.    Find the theta value. If B<=G then H= theta. If B>G then H= 360-theta

4.    Use ‘acosd’ function to find inverse cosine and obtain the result in degrees.

5.    Divide the hue component by 360 to represent in the range [0 1]

6.    Similarly, find the saturation and the intensity components.

7.    Display the image.

Table Tennis Ball Detection-MATLAB CODE:

Using the concept of roundness to detect a circular object, the table tennis ball in mid air or in palm is detected.

FLOW CHART:

1.     Read input image

·        Read a RGB image with ball in mid air or ball in palm.

                   MATLAB CODE:

                       %Read the input image

         Img=imread(Filename);

         axes('Position',[0 .1 .74 .8],'xtick',[],'ytick',[]);

         imshow(Img);title('Original Image');

2.     Image Pre-processing

·        Convert the RGB image to grayscale image.

·        Apply median filter

·        Adjust the brightness and contrast of the image using ‘imadjust’ function.

                      MATLAB CODE:

           I=rgb2gray(Img);     % Converting RGB Image to

                                % Gray Scale Image

           I=im2double(I);      % Converting Gray scale Image

                                % to Double type

           J = medfilt2(I,[3 3]); % Median Filter , 

                                  % 3x3 Convolution

                                  % on Image

           I2 = imadjust(J);     % Improve to quality of Image

                                 % and adjusting

                                 % contrast and brightness values

3.     Threshold the image

·        The pixel value greater than the threshold value is converted to one else zero.

                     MATLAB CODE:

           Ib = I2> 0.9627; 

4.     Image Labeling

·        Label the connected components using ‘bwlabel’ function

·        Remove components that are smaller in size.

                     MATLAB CODE:

           %Labelling

           [Label,total] = bwlabel(Ib,4); % Indexing segments by

                                          % binary label function

            %Remove components that is small and tiny

            for i=1:total

                if(sum(sum(Label==i)) < 500 )

                    Label(Label==i)=0;

  

                end

            end

5.     Find the image properties: Area, Perimeter and Centroid

·        Using ‘regionprops’ function, find the Area, Perimeter, Bounding Box and Centroid.

                     MATLAB CODE:

            %Find the properties of the image

             Sdata = regionprops(Label,'all'); 

6.     Calculate the Roundness

·        Roundness = 4*PI*A/P^2

           MATLAB CODE:

                %Find the components number

          Un=unique(Label);

          my_max=0.0;

            %Check the Roundness metrics

            %Roundness=4*PI*Area/Perimeter.^2

            for i=2:numel(Un)

               Roundness=(4*pi*Sdata(Un(i)).Area)/Sdata(Un(i)).Perimeter.^2;

               my_max=max(my_max,Roundness);

                if(Roundness==my_max)

                   ele=Un(i);

                end

            end

7.     Find the component with the maximum roundness value

·        Find the max of the Roundness value for all the labeled components

8.     Show the detected table tennis ball

·        Use the ‘BoundingBox’ values to plot rectangle around the ball

·        Mark the centroid of the ball

                     MATLAB CODE:

          %Draw the box around the ball

           box=Sdata(ele).BoundingBox;

           box(1,1:2)=box(1,1:2)-15;

           box(1,3:4)=box(1,3)+25;

          %Crop the image

           C=imcrop(Img,box);

          %Find the centroid

          cen=Sdata(ele).Centroid;

          %Display the image

           axes('Position',[0 .1 .74 .8],'xtick',[],'ytick',[])

           imshow(Img);

           hold on

           plot(cen(1,1),cen(1,2),'rx');%Mark the centroid

9.     Generate report

·        Find the radius using the Equidiameter obtained using ‘regionprops’ function.

·        Display the radius,Area,Perimeter and Centroid of the ball.

·        Show the Binary and Original image of the cropped ball.

MATLAB CODE:

        Rad=(sdata(ele).EquivDiameter)/2;

        Rad=strcat('Radius of the Ball :',num2str(Rad));

     

        Area=sdata(ele).Area;

        Area=strcat('Area of the ball:',num2str(Area));

       

        Pmt=sdata(ele).Perimeter;

        Pmt=strcat('Perimeter of the ball:',num2str(Pmt));

       

        Cen=sdata(ele).Centroid;

        Cent=strcat('Centroid:',num2str(Cen(1,1)),',',num2str(Cen(1,2)));

BALL IN MID AIR: