Prefer Your Language

Search This Blog

Matlab code for Sobel operator, Implementation Sobel operator in Matlab on YUV video File



Implementation Sobel operator in Matlab on YUV video File

Today, we discuss Sobel operator and how to apply on YUV video file 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 frame 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.

If you are not aware how to apply sobel operator on image, please refer this link for that.
The Sobel Operator
To calculate the gradient of each point in the frame, the frame is convolved with the Sobel Kernel. Convolution is done by moving the kernel across the frame, 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 frame, with the center pixel of the kernel aligning we the first pixel in the frame. 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 frame. 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.      1.  Load YUV File.
2.    Load Frame of YUV video file.
3.    Load Y , U and V Plane separately
4.    Iterate over every pixel in the Plane for each Plane for each Frame
5.    Apply the x gradient kernel
6.    Apply the y gradient kernel
7.    Find the length of the gradient using pythagoras' theorem
8.    Normalize the gradient length to the range 0-255
9.    Set the pixels to the new values
10. Combined all three Plane and make a frame
11. Write this frame to video file yuv format
12. Save the YUV file.

This illustrates the main steps, though we miss out some specifics such as what we do when we meet the boundary of the frame.

Formula for calculating gradient in X and Y direction ;
Let say (i,j) is the middle pixel of 3*3 grid then we have the following formula


               
Gx =  - ( - Z1 + Z3 - 2*Z4 + 2*Z6 – Z7 + Z9 )  / 8

&

Gy = -(-Z1 – 2*Z2 – Z3 + Z7 + 2*Z8 + Z9 ) /8


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;
Code have Three procedure.
Save below procedure in new M file with name Sobel_Operator.m





%Procedure; Sobel operator procedure; it takes one frame, apply sobel operator and return it
function[ Gray_Edge_Image] = Sobel_operator( orig_frame , a)


        gh_m = zeros(size(orig_frame,1), size(orig_frame,2)) ;
        gv_m = zeros(size(orig_frame,1), size(orig_frame,2)) ;


           
            Gray_Edge_Image = zeros (size(orig_frame,1) , size(orig_frame,2)) ;
      
           
            %make a temporary gray scale image
            Temp_Image = double (orig_frame) ; 
           
           
           
            %Apply sobel operator
           
            %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
                   
                    Gradient_v = -(-Temp_Image (i-1,j-1) + Temp_Image (i-1,j+1) - 2*Temp_Image (i,j-1) + 2*Temp_Image (i,j+1) - Temp_Image (i+1,j-1) + Temp_Image (i+1,j+1) ) /8 ;
                        
                    Gradient_h = -(-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) ) /8 ;
                   
                    %obtaining Gradient for a pixel
                    Gray_Edge_Image (i-1 , j-1) = sqrt (  Gradient_v * Gradient_v  + Gradient_h*Gradient_h );
                   
                    Gray_Edge_Image ( i-1 , j-1 ) = Thrash_Hold (Gray_Edge_Image ( i-1 , j-1 ), a ) ;
                    
                end
            end
end

And here the main code; save this code in new M file with name 
Thrash_Hold.m

%Procedure for Thrashold value
%here threshold i defined threshold empirecaly adaptive thresholding


function[ Value ] = Thrash_Hold ( Z , a )


            if( Z >= a)
                    Value = 255;
            else
                    Value = 0;
            end
end

And save this file as Sobel.m
%Reading and writing YUV video file frame by frame


% close all;
% clear all;


%all global variable(s)
height = 480    ;  %height of the video
width = 640;   %Width of the video
num_frames = 1;   %number of framse in video
THRESH_HOLD_Y = 20 ;
THRESH_HOLD_U = 2 ;
THRESH_HOLD_V = 2;


filename = Input_file_name.yuv';  %file name of video
filename_copy1 = OutputFilename_matlab.yuv' ; % output obtained by our code
filename_copy2 = 'Output_stand.yuv' ; %Output obtained by Standard edge method




%open video file of .yuv format
fid_read = fopen (filename, 'rb');


fid_write1 = fopen (filename_copy1,'wb') ;


fid_write2 = fopen (filename_copy2,'wb') ;


Time = cputime ;
%loop over every frames of video
for index= 1 : num_frames
 
    %reading first frame ; deg will be the vector of frames and num will
    %contain number of data read successfully ; In YUV 4:2:0; each frame is of size = 1.5*width*height

    [deg,num]=  fread( fid_read, 1.5*width*height, 'uchar'); %read frame
 
    %if all frame have been read than print msg and exit
    if (num == 0)
                fprintf('read %d frames and exiting',index);
   
                break;
    end


    %reading frame for the component Y,U and V; seprately, starts from here
   
    %since Y frame; from 1 to width*Height (of size M*N)
    deg_y = deg(1:width*height);
   
    %U frame; from 1+width*height (since, we read Y frame previously from 1 ) to width*Height + width*Height / 4 (of size M*N/4)
    deg_u = deg(1+(width*height):(width*height)+(width*height/4));
   
    %V frame; from 1 + 1.25*width*height (since, we read U frame previously upto 1.25*width*height ) to end (of size M*N/4)
    deg_v = deg(1+(1.25*width*height) : width*height*1.5 );




   
    %reshaping the Y, U and V from to its size M*N , (M/2)(N/2) and (M/2)(N/2) respectively;  reshape will convert 1D array to 2D array
    %Reshaping Y frame 
    orig_y_imgtop1 = reshape(deg_y,width,height)'; % here ' means transpose;
   
    %Reshaping U frame
    orig_u_imgtop1 = reshape(deg_u,width/2,height/2)';
   
    %Reshaping V frame
    orig_v_imgtop1 = reshape(deg_v,width/2,height/2)';
   
    %applying sobel operator on each frame and store result back to the
    %respective arrays
   
    %output obtained by standard sobel operator
    [orig_y_imgtop_st, Y_Matlab] = edge (orig_y_imgtop1, 'sobel',THRESH_HOLD_Y , 'nothinning'  );
    [orig_u_imgtop_st, U_Matlab] = edge (orig_u_imgtop1, 'sobel',THRESH_HOLD_U , 'nothinning' );
    [orig_v_imgtop_st, V_Matlab] = edge (orig_v_imgtop1, 'sobel',THRESH_HOLD_V, 'nothinning' );
      
    orig_y_imgtop_st = 255 * orig_y_imgtop_st ;
    orig_u_imgtop_st = 255 * orig_u_imgtop_st ;
    orig_v_imgtop_st = 255 * orig_v_imgtop_st ;
       
    %output obtained by my sobel operator
    orig_y_imgtop_m  = Sobel_operator (orig_y_imgtop1,THRESH_HOLD_Y );
    orig_u_imgtop_m =  Sobel_operator (orig_u_imgtop1, THRESH_HOLD_U  );
    orig_v_imgtop_m =  Sobel_operator (orig_v_imgtop1, THRESH_HOLD_V  );
  
  
    %wrting back to new video file
   
    %take the transpose
    orig_y_imgtop1 = orig_y_imgtop_m';
   
    orig_y_imgtop2 = orig_y_imgtop_st';
   
    %convert 2D array to 1D array using (:)
    deg_y_new1 = orig_y_imgtop1(:);
    deg_y_new2 = orig_y_imgtop2(:);
   
   
    %take the transpose
    orig_u_imgtop1 = orig_u_imgtop_m';
    
    orig_u_imgtop2 = orig_u_imgtop_st';
   
    %convert 2D array to 1D array using (:)
    deg_u_new1 = orig_u_imgtop1(:);
   
    deg_u_new2 = orig_u_imgtop2(:);
   
    %take the transpose
    orig_v_imgtop1 = orig_v_imgtop_m';
    orig_v_imgtop2 = orig_v_imgtop_st';
   
    %convert 2D array to 1D array using (:)
    deg_v_new1 = orig_v_imgtop1(:);
    deg_v_new2 = orig_v_imgtop2(:);
   
    %lets concatenate each frame into single vector
    deg_new1 = [deg_y_new1; deg_u_new1; deg_v_new1];
    deg_new2 = [deg_y_new2; deg_u_new2; deg_v_new2];
   
    %now write into another yuv file
    fwrite( fid_write1, deg_new1, 'uchar');
    fwrite( fid_write2, deg_new2, 'uchar');
   
   
end


Elapsed_time = cputime - Time 


%close everyting
fclose (fid_read) ;
fclose (fid_write1) ;
fclose (fid_write2) ;
 



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

1 comment:

Help us to improve our quality and become contributor to our blog

Become a contributor to this blog. Click on contact us tab
Blogger Template by Clairvo