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