Implementation Sobel operator in Matlab on YUV image and other extension image
Implementation
Sobel operator in Matlab on YUV image and other extension image
Today, we discuss Sobel operator and how to apply on image with step by step
discussion. If you are not familiar with the sobel operator or don’t know in
detail, don’t worry, we first discuss what is sobel operator followed by its
Matlab code.
Lets start our discussion.
First
question is what is Sobel filter/operator?
The algorithm we will
look at in this tutorial is an edge detection algorithm, specifically an edge
detection algorithm based on the Sobel operator. This algorithm works by
calculating the gradient of the intensity of the image at each point, finding
the direction of the change from light to dark and the magnitude of the change.
This magnitude corresponds to how sharp the edge is.
The Sobel Operator
To calculate the
gradient of each point in the image, the image is convolved with the Sobel
Kernel. Convolution is done by moving the kernel across the image, one pixel at
a time. At each pixel, the pixel and its neighbours are weighted by the
corresponding value in the kernel, and summed to produce a new value. This
operation is shown in the following diagram.
~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
k00 | k10 | k20 | | i00 | i10 | i20 | i30 | ..
|~~~~~+~~~~~+~~~~~+ |~~~~~+~~~~~+~~~~~+~~~~~+~~~
|
k01 | k11 | k21 | | i01 | i11
| i21 | i31 | ..
|~~~~~+~~~~~+~~~~~+ |~~~~~+~~~~~+~~~~~+~~~~~+~~~
|
k02 | k12 | k22 | | i02 | i12
| i22 | i32 | ..
|~~~~~+~~~~~+~~~~~+ |~~~~~+~~~~~+~~~~~+~~~~~+~~~
| i03 | i13 |
i23 | i33 | ..
Applying
the kernel to the
|~~~~~+~~~~~+~~~~~+~~~~~+~~~
to
the first pixel gives us | .. |
.. | .. | .. | ..
an
output value for that
pixel
|
o00 = k11 * i00 + k12 * i01 + k21 * i10 + k22 * i11
|
we treat it as though
the kernel has been overlaid onto the image, with the center pixel of the
kernel aligning we the first pixel in the image. Then we multiple each entry in
the kernel by the value beneath it, and sum them to produce a single single
output value from that pixel.
For the pixels on the
boundary, we just ignore any entries in the kernel that fall outside.
Edge detection using
the Sobel Operator applies two separate kernels to calculate the x and y
gradients in the image. The length of this gradient is then calculated and normalized
to produce a single intensity approximately equal to the sharpness of the edge
at that position.
The kernels used for
Sobel Edge Detection are shown below.
~~~~~~~~~~~~~
|-1 | 0 | 1 |
|~~~+~~~+~~~|
|-2 | 0 | 2 |
|~~~+~~~+~~~|
|-1 | 0 | 1 |
~~~~~~~~~~~~~
X gradient kernel
~~~~~~~~~~~~~
|-1 |-2 |-1 |
|~~~+~~~+~~~|
| 0 | 0 | 0 |
|~~~+~~~+~~~|
| 1 | 2 | 1 |
~~~~~~~~~~~~~
Y gradient kernel
|
The Algorithm
Now the algorithm can be broken down
into its constituent steps
1 Iterate over every pixel in the image
2 Apply the x gradient kernel
3 Apply the y gradient kernel
4 Find the length of the gradient using
pythagoras' theorem
5 Normalize the gradient length to the range
0-255
6 Set the pixels to the new values
This illustrates the main steps,
though we miss out some specifics such as what we do when we meet the
boundary of the image.
|
We need to do one tweak so that this will work
for our edge detection. When we are iterating over the image, we will be
changing value however when we move on to the next pixel, the just changed
value will still be in the area of influence of our kernel. For this reason we
need to make sure the image we read pixels from and the image we write pixels
to are different, and we can do this by just created an empty image, which we
will populate with the new values.
|
Now we will read pixels from the original image, but write the
edge values to
img_edge
Before we start fiddling with the pixel values, we need to
decide how we will handle the image boundaries. There are a few different ways
to do this, one way is to check each pixel to see if it is on the boundary of
the image, and ignore the non-existant pixels on that side. An alternative way
would be to wrap the pixels, so that a non-existant pixel on the left edge
actually maps to the pixel on the same row on the right edge. For simplicities
sake I am simply going to ignore the boundary pixels, instead choosing to start
one pixel in from the starting point, and 1 pixel away from the finishing edge.
To do this, change the iterators to start at 1, and finish at
img.width - 1 and img.height
- 1 respectively.
Now we can safely index adjacent pixels without falling off the edge of the
image.
Right, lets move on to
applying the gradients. For the x gradient, we need the column of pixels to the
left and right, for each of the 6 pixels we index we will add up the red, green
and blue intensities, multiply it by the kernel, and accumulated it into Gx. We
will also do the same for the y gradient. As Gx and Gy both use the same pixels
on occasions, it is best to apply the two kernels simultaneously.
Now that the kernels are
applied, Gx and Gy contain un-normalised values, which equal the relative
length of the gradient in their respective axis. We wish to calculate the
length of the gradient, and normalise it to a range suitable for displaying.
Lets calculate the
relative length of each gradient. We can do this with pythagoras' theorem:
|
|
This gives us an un-normalised length, to get the length
normalised, we need to know the range of
length
We know
that each pixel's intensity value has a range of 0 to 765, due to it being the
addition of 3 variables with range 0 to 255. Also, from looking at the kernels,
we can see the maximum amount you can accumulate is +/- 4 intensities.
Therefore the total range of Gx and Gy is between 0 and 3060. These values are
then squared, added and square rooted, giving us a final range of 0 to 4328 (
sqrt(2 * sqr(3060)) ). We can renormalise by dividing by the upper bounds and
multiplying by 255.
Right at the end of the code we need to change
displayImage to display our new image rather than
the original. Other than that, the script is finished, and we can test it to
see the results. The script should produce a new grayscale image, highlighting
edges from the original image.
This concludes this
tutorial. Below are the code listings for this tutorial. The C code will
actually continually capture images from the webcam, and you will find that on
the Raspberry Pi it is capable of displaying the edge detected image in
realtime.
Let’s
formulate the Formula;
So, we reach the required formula;
Here it is;
|
Gx = ( z7 + 2z8
+ z9) – ( z1 + 2Z2 + z3 )
&
Gy = (z3 + 2z6 +
z9) – ( z1 + 2Z4 + z7)
|
Now time to move ahead, now we are looking for Matlab
code, right?
Ok let’s start coding in Matlab;
For understand below code, you must be aware to the
following function’s.
Here is the complete Code; But this code have some Logical error. First you try to figure out
and in the next post we will discuss its logical error and correct them.
Incorrect
Solution
|
%Procedure for Sobel operator
Original_Image = imread('a.bmp');
Gray_Edge_Image = rgb2gray
(Original_Image);
Temp_Image = double (Gray_Edge_Image)
;
%Start iterate from first corner to
last corner
% size(Temp_Image,1) gives number of
rows; we subtract 2 bec. of leaving
boundry
for i = 1 : size(Temp_Image,1) - 2
for j = 1 : size (Temp_Image, 2) - 2
%calculating gradient in X direction
X_Gradient = ( ( Temp_Image ( i+2, j )
+ 2*Temp_Image( i+2, j+1 ) +
Temp_Image
( i+2 ,j+2 ) ) - ( Temp_Image( I , j ) +
2*Temp_Image(
I , j+1 ) + Temp_Image ( I , j+2 ) ) );
%calculating gradient in Y direction
Y_Gradient = ( ( Temp_Image( I , j+2 )
+ 2*Temp_Image( i+1, j+2 ) +
Temp_Image
( i+2, j+2 ) ) - ( Temp_Image( I , j )
+
2*Temp_Image( i+1, j ) + Temp_Image ( i+2, j ) ) );
%Subtracting 256 bec. of brightness ;
you may leave it
Gray_Edge_Image (i , j) = 256 - sqrt
( Y_Gradient * Y_Gradient
+ X_Gradient*X_Gradient );
end
end
figure, imshow(Original_Image), figure,
imshow(Gray_Edge_Image);
|
Are
you able to identify the logical error? No.
Ok let’s identify and correct them .
Identifying
logical error
We have already documents the Algorithm. In the
documented algorithm we missed one step and code it incorrectly, steps is
|
Normalize the gradient length
to the range 0-255
|
Just
subtracting the resultant value from 255 does not ensure that gradient length
to the range 0-255. Now, question is how to resolve it, Very easy, Here the
steps;
Steps to Resolve above logical
error
·
First find the maximum value of pixel and
store in Max_Value.
·
Afterword’s divide each pixel value by
Max_Value.
Now this step ensure that the pixel
values is in between 0 and 1, means in the range of [0 1].
·
Now, multiply the each pixel value to by 255
so that range become
255*[0 1] = [0 255]. This solve our
problem.
Above steps are easy to implement, here
is the complete code
Correct
Solution
|
%Modified code for Sobel filter
Original_Image = imread('a.bmp');
Gray_Edge_Image = rgb2gray (Original_Image);
Temp_Image = double (Gray_Edge_Image) ;
%variable for finding maximum value across Matrix
max_value = -1 ;
%Start iterate from first corner to last corner
% size(Temp_Image,1) gives number of rows; we
subtract 2 bec. of leaving boundry
for i = 2 : size(Temp_Image,1) - 1
for j =
2 : size (Temp_Image, 2) - 1
%calculating gradient in X direction
X_Gradient = ( ( Temp_Image( I + 1 ,
j – 1 ) + 2*Temp_Image(i+1,j)
+ Temp_Image
(i+1,j+1) ) - ( Temp_Image(i-1,j-1)
+
2*Temp_Image(i-1,j) + Temp_Image (i-1,j+1) ) );
%calculating gradient in Y direction
Y_Gradient = ( ( Temp_Image(i-1,j+1) + 2*Temp_Image(i,j+1)
+
Temp_Image (i+1,j+1) ) - ( Temp_Image(i-1,j-1)
+
2*Temp_Image(i,j-1) + Temp_Image (i+1,j-1)));
Gray_Edge_Image (i-1 , j-1) = sqrt ( Y_Gradient*Y_Gradient
+ X_Gradient*X_Gradient );
if
(max_value < Gray_Edge_Image ( i-1, j-1 ) )
max_value = Gray_Edge_Image ( i-1 ,j-1 ) ;
end
end
end
%converting whole image value in between 0 to
255; division makes it in range of [0,1]
for i = 2 : size(Temp_Image,1)-1
for j =
2 : size (Temp_Image, 2)-1
Gray_Edge_Image ( i-1
,j-1 ) =
((Gray_Edge_Image
( i-1 ,j-1 )/ max_value ) * 255 );
end
end
figure, imshow(Original_Image), figure,
imshow(Gray_Edge_Image);
|
 |
| Output |
Got Questions?
Feel free to ask me any
question because I'd be happy to walk you through step by step!
References and External Links
For Contact us….. Click on Contact us Tab
Posts
Seriously that little modification for taking a value for Gray_Edge_Image ( i-1 ,j-1 ) with respect to max_value for 255 is an most remarkable thing! that little modification changes the whole output deeply! thank u so so so much for such amazing work! keep it up! i have damn impressed with your way of thinking for extracting tiniest errors! your doing an amazing job sir!
ReplyDeleteHi i have a doubt. The max_value in this code never goes beyond 255, so i am not getting the idea of normalisation. Applying that changes the picture drastically.
ReplyDeleteWhy not jsut use uint8 conversion?
ReplyDeleteCanyou do an example but with prewitt edge detection???
ReplyDeleteI must say you had done A tremendous job, I appreciate all your efforts. Thanks a lot
ReplyDeleteWonderful article. I love your style of writing. It’s inspiring, and I am happy
ReplyDeleteWow, this post is pleasant, I will share this with my all friends and keep it up.
ReplyDeleteReally very interesting and very easy to understand, I really like this!
ReplyDelete