Sobel filter is a nice filter

I <3 Sobel filter

One of the first edge detection filters one encounters in image processing is Sobel filter. It is an adorable filter. It consists of the following two kernels
$$
S_x=\left[\begin{array}{ccc}-1 & -2 & -1 \\\,0 & \,0 &\,0 \\\,1 & \,2 & \,1\end{array}\right] \mbox{ and }S_y=\left[\begin{array}{ccc}-1 & 0 & 1 \\ -2 & 0 &2 \\-1 & 0 & 1\end{array}\right]
$$
(Remember our convention that the origin is at the top-left corner, $x$-axis is going down and $y$-axis is going to the right.)

You align the central zero of these operators and put them on top of each image pixel multiply the 9 pixels of the image with the 9 pixels of the operator and add the result. In other words convolve the image $u$ with $S_x$ and $S_y$. The results $u*S_x$ and $u*S_y$ give the location of the horizontal and vertical edge.

Another way to look at it is to observe the following:
\begin{align*}
&(u*S_x)(i, j)\\
&= u(i+1, j-1)-u(i-1, j-1)+2u(i+1, j)-2u(i-1, j)+u(i+1, j+1)-u(i-1, j+1)\\
&= 2\Big(\frac{u(i+1, j-1)-u(i-1, j-1)}{2}\Big)\\
& \, +4 \Big(\frac{u(i+1, j)-u(i-1, j)}{2}\Big)\\
& \, +2\Big(\frac{u(i+1, j+1)-u(i-1, j+1)}{2}\Big).\\
\end{align*}
It looks like there are three central difference schemes here.
Indeed,

$(u*S_x)(i, j)=2 D_xu(i, j-1)+4 D_xu(i, j) + 2 D_xu(i, j+1)$

and,

$(u*S_y)(i, j)=2 D_yu(i-1, j)+4 D_yu(i, j) + 2 D_yu(i+1, j),$

where $D_x$ and $D_y$ denote central difference schemes in the $x-$ and $y-$ directions respectively. I like to normalize theses operators by dividing by 8, which makes it a weighted average of three central difference schemes. i.e.
\begin{align*}
(u*S_x)(i, j)&\equiv \frac{D_xu(i, j-1)+2 D_xu(i, j) + D_xu(i, j+1)}{4}, \\
\\
(u*S_y)(i, j)&\equiv \frac{D_yu(i-1, j)+2 D_yu(i, j) + D_yu(i+1, j)}{4}.
\end{align*}
For coding Sobel filter one should be careful at the image boundary. I like to reflect the boundary. This is achieved by expanding the image u by one pixel and computing the convolution in the interior and then shrinking the image by one pixel.

So let's write the codes for expanding the image and shrinking the image. It would be useful to write these codes so that we could use them for expanding and shrinking by n number of pixels. Trust me, these codes will come handy later.

Shrinking an image

%==========================================================================

function u=shrink(U, n)

% Author: Prashant Athavale

% Date: 11/23/2011

% Please acknowledge my name, if you use this code.

% This program shrinks the image by n by removing its boundary n

% times.

%==========================================================================

if nargin==1

    n=1;

end

[M, N, P]=size(U);

if P==1

    u=U(n+1:M-n, n+1:N-n);

elseif P==3

    R=U(:,:,1); r=R(n+1:M-n, n+1:N-n);

    G=U(:,:,2); g=G(n+1:M-n, n+1:N-n);

    B=U(:,:,3); b=B(n+1:M-n, n+1:N-n);

    u=cat(3, r, g);

    u=cat(3, u, b);

end

end % end of the function shrink

%==========================================================================

Remarks

There are a few good things to note here. The line u=U(n+1:M-n, n+1:N-n); gets rid of the n-pixel boundaries of a grayscale image U.

Can you write a code for expanding the image?

If it is an RGB image then we have to do this separately for each channel, and then concatenate the channels. We could use the command cat for this. The number 3 in cat(3, r, g) gives the dimension in which the matrices are joined together. 


Sobel filter


Lastly, we code the Sobel filter itself.

%==========================================================================

function [A, ux, uy]=sobel_for(u, Abs)

% Author: Prashant Athavale

% Date: 11/23/2011

% Please acknowledge my name, if you use this code.

% This function performs sobel filtering.

%==========================================================================

if nargin==1

    Abs=0;

end

u=double(u);

% Let us expand the image u by 1 unit

u=expand(u);

[M, N, P]=size(u);

% preallocation of memory for faster implementation

ux=zeros(size(u));

uy=zeros(size(u));

% implementation of the for loop

t=cputime;

for i=2:M-1

    for j=2:N-1

        ux(i, j)=(u(i+1, j-1)-u(i-1, j-1)...

            +2*u(i+1, j)-2*u(i-1, j)...

            +u(i+1, j+1)-u(i-1, j+1))/8;

        uy(i, j)=(u(i-1, j+1)-u(i-1, j-1)...

            +2*u(i, j+1)-2*u(i, j-1)...

            +u(i+1, j+1)-u(i+1, j-1))/8;

    end

end

% Let us shrink the images back to the original size

ux=shrink(ux);

uy=shrink(uy);

if Abs~=0

    if P==1

        A=(abs(ux)+abs(uy))/2;

    else

        A=zeros(size(ux));

        for i=1:P

            A(:,:,i)=(abs(ux(:,:,i))+abs(uy(:,:,i)))/2;

        end

    end

else

    A=sqrt(ux.*ux+uy.*uy);

end

end % end of the function sobel_for

%==========================================================================

Ok, so I named this code sobel_for as I am using for loop. Can you code this without using them? Does that make the code faster?


Results

For this filtering I used the following image:

Aishwarya Rai 

For the following code ...

u=imread('eyes.jpg');
A=sobel_for(u, 1);
figure; imshow(A/92);


 I got the following result ... 

Result of Sobel filtering

Our eyes are most sensitive to gradients

Our eyes, and our sense in general, are most sensitive to gradients. May be that is the reason when we are asked to draw a picture from memory, we tend to draw the sketch involving primarily the edges. It is also the reason why it is difficult for a person from California to get adjusted to the winter in Canada.

Anyway we are ready to talk about gradient (u_x, u_y) of an image u. We already have the code Dfx_plus for computing the partial derivative in the x-direction.

Image gradient 

Now let's compute the gradient and store it in matrices Gx and Gy. Here is the trivial code to do that.

%==========================================================================

function [Gx, Gy]=gradient_plus(u)

% Author: Prashant Athavale

% Date: 1122/2011

% Please acknowledge my name if you use this code, thank you.

% This function computes the gradient f the image u.

% It uses forward approximation of the derivatives and thus the name

% gradient_plus

%==========================================================================

if strcmp(class(u),'uint8')

    u=double(u);

end

Gx=Dfx_plus(u);

Gy=Dfy_plus(u);

end

%==========================================================================

Magnitude of the image gradient

Let us now write another code for computing the magnitude of the gradient.

%==========================================================================

function Gm=grad_plus_magnitude(u)

% Author: Prashant Athavale

% Date: 11/22/2011

% Please acknowledge my name if you use this code, thank you.

% This function takes the image stored in u and gives the magnitude of the

% gradient

%==========================================================================

if strcmp(class(u),'uint8')

    u=double(u);

end

[Gx, Gy]=gradient_plus(u);

Gm=sqrt(Gx.*Gx+Gy.*Gy);

end

%==========================================================================

Remark 1: Pointwise operation

Note in the above code you see the line Gx.*Gx. Here the .* gives point wise multiplication of the entries in the matrix Gx. This is useful as one more time we could avoid a for loop, making things faster.

Let's implement the following in the command window.

>> u=imread('barbara.png');

>> Gm=grad_plus_magnitude(u);

>> figure; imshow(Gm/255+0.5);

>> edges=(Gm>15);

>> figure; imshow(edges);

Remark 2: Logical assignments

Note the line edges=(Gm>15). Here is another instance we avoided a for loop. What it does is, it creates a matrix edges of class logical. The value edges(i, j)=1 if Gm(i, j) > 15. Write a code using a nested for loop to do this and see which one is faster.

Results

Following are the images I got for Gm and edges as a result of the above code. This is how a gradient looks like. I don't know about you, but I was very happy when I saw a gradient of an image for the first time.

Magnitude of the gradient in the image Barbara

Edges, defined as the points where gradient is more than 15

Partial derivatives with backward difference scheme


The partial derivatives with backward difference scheme are given by

u_x:=(u(x, y)-u(x-h, y))/h  and u_y:=(u(x, y)-u(x, y-h))/h .

I recommend that you code these as a practice.

Taking partial derivatives is easy in Matlab ( also why I don't like the class uint8 )

Can we take partial derivatives of the images? It is really easy. Grayscale digital images can be considered as 2D sampled points of a graph of a function u(x, y) where the domain of the function is the area of the image.

Here, we try to find the partial derivative of u with respect to x. There are many ways to approximate this, one of the ways is to look at the forward difference. i.e. $$\frac{\partial u}{\partial x}(x_0,\, y_0)\approx \frac{u(x_0+h,\, y_0)-u(x_0,\,y_0)}{h}=\frac{u(i+1,\, j)-u(i,\, j)}{h}$$
Note that we are taking the origin at the top left corner, and the sampling distance, i.e. the distance between the two pixels is h. I like to take h=1. Some of my friends like to take h=1/M if the image is MxM image.

Ok, so let's write the code to find the partial derivative of u in x-direction.

The first attempt 

clear all;

u=imread('lenna1.png'); 

clc;

close all;

[M, N, P]=size(u);

h=1.0;

dx_plus=zeros([M, N, P]); % the partial derivative in x direction

for i=1:M-1

    for j=1:N

        dx_plus(i, j, :)=(u(i+h, j, :)-u(i, j, :))/h;

    end

end

figure;

imshow(dx_plus/255+0.5);

This seems like a reasonable code. My claim is that it is wrong. Can you find out the mistake?
.
.
.
.
.
.
.
.
.
The problem here is the line u=imread('lenna1.png');
This gives stores the image lenna.png in the variable u of class uint8.
The class uint8 stores only integers. And the operations between two uint8 is a uint8. So if dx_plus has a negative value, it is stored as zero, i.e. we miss out all negative derivatives.

The corrected code

Try the following code yourself and see the difference.

%==========================================================================

clear all;

A=imread('lenna1.png');

clc;

u=A;

U=double(u); % U is now of class double

[M, N, P]=size(u);

h=1.0;

dx_plus=zeros([M, N, P]);

dx_plus2=zeros([M, N, P]);

for i=1:M-1

    for j=1:N

        dx_plus(i, j, :)=(u(i+h, j, :)-u(i, j, :))/h;

        dx_plus2(i, j, :)=(U(i+h, j, :)-U(i, j, :))/h;

    end

end

difference=abs(dx_plus2-dx_plus);

figure;

imshow(difference/255+0.3);

%==========================================================================

Function to find partial derivative in the x-direction

Finding partial derivatives is fun and it is something that you will be doing many times. So why not write a function for that?

%==========================================================================

function dx_plus=Dx_plus(u, h)

% Author: Prashant Athavale

% Date 11/19/2011

% Please acknowledge my name if you use this code, thank you.

% This function takes a grayscale image variable u and gives its derivative

% The derivative is the forward x-derivative,

% i.e. we try to approximate the derivative with (u(x+h)-u(x))/h

% but note that the x-axis is taken going down, top left is the origin

% and the distance between the two pixels, h, is taken as 1-unit

% Some people like to take the distance h=1/M for square images where M is the dimension of the image.

%==========================================================================

if nargin==0;

    error('At least one input is needed');

elseif nargin>0

    if nargin==1

        h=1; % the distance between two pixels is defined as 1

    end

    [M, N, P]=size(u);

    if strcmp(class(u), 'double')==0

        u=double(u);

    end

    dx_plus=zeros([M, N, P]);

    for i=1:M-1

        for j=1:N

            dx_plus(i, j, :)=(u(i+1, j, :)-u(i, j, :))/h;

        end

    end

end

end % end of the function Dx_plus

%==========================================================================

How about the partial derivative in the y-direction ?

Here is a function for that...

%==========================================================================

function dy_plus=Dy_plus(u, h)

% Author: Prashant Athavale

% Date 11/19/2011

% This function takes a grayscale image variable u and gives its derivative

% The derivative is the forward x-derivative,

% i.e. we try to approximate the derivative with (u(y+h)-u(y))/h

% but note that the y-axis is taken going to the right 

% the origin is at the top left corner.

% and the distance between the two pixels, h, is taken as 1-unit

% Some people like to take the distance h=1/M for square images where M is the dimension of the image.

%==========================================================================

if nargin==0;

    error('At least one input is needed');

elseif nargin>0

    if nargin==1

        h=1; % the distance between two pixels is defined as 1

    end

    [M, N, P]=size(u);

    if strcmp(class(u), 'double')==0

        u=double(u);

    end

    dy_plus=zeros([M, N, P]);

    for i=1:M

        for j=1:N-1

            dy_plus(i, j, :)=(u(i, j+1, :)-u(i, j, :))/h;

        end

    end

end

end % end of the function Dy_plus

%==========================================================================

My result


Well here is the x-partial derivative I got  with the following commands:

>> u=imread('lenna3.png');
>> ux=Dx_plus(u);
>> figure; imshow(ux/255 +0.5);