Monday, November 21, 2011

Life is short, and we are all going to die

Life is short, and we are all going to die...

So we need to utilize our time as effectively as possible.
Fortunately, in Matlab it is easy to keep track of time using the command t=cputime
for example:

t=cputime;
% your code
e=cputime-t

This e gives the elapsed time while Matlab was implementing your code.
Last time we wrote codes for  finding partial derivatives in Matlab.
Let's see how much time we spent doing it. Following is the code that gives me the elapsed time.


>> u=imread('lenna3.png');
>> ts=cputime;ux=Dx_plus(u);es=cputime-ts


es =


    0.3900

So, it takes 0.39 seconds to take the x-derivative.

Recall we used nested for loops in our codes before. Unfortunately, it is not the fasted way of doing things in Matlab. As a rule of thumb, avoid for loops as much as possible in Matlab. They slow down your computations.


Let us see if we can avoid the for loops in the previous code. You, being a smart reader must have observed that all we were doing effectively in the Dx_plus code was
(i) shifting all the rows of the input matrix u up by 1 unit and 
(ii) subtracting u from this shifted matrix.

If we can do this we can avoid for loops. In matlab there is a cute little command called circshift. This command shifts the matrix entries. e.g. u_shifted=circshift(u, [-1 0]);shifts  the first dimension of the matrix u  up by 1, while keeping the second dimension unchanged. Using this we can write the following (hopefully) fast function to approximate the partial derivative of u in the x-direction.


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

function dx_plus=Dfx_plus(u, h)

% Author: Prashant Athavale
% Date 11/20/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,

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

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, ~, ~]=size(u);
    if strcmp(class(u), 'double')==0
        u=double(u);
    end
    
    u_iplus1=circshift(u, [-1 0]);
    % this shifts the first dimension back by 1
    dx_plus=(u_iplus1-u)/h;
    dx_plus(M, :, :)=0;
end


end % end of the function Dfx_plus

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


Remarks
Note the mysterious appearance of "~" (tilde) in [M, ~, ~]=size(u); it essentially means that we only need the first dimension of the size of u. 
dx_plus(M, :, :)=0; indicates that the last row of the result is forced to be zero. 

Let's run the code: 

>> tf=cputime;ufx=Dfx_plus(u);ef=cputime-tf

ef =

    0.0100

This is incredible. It runs the code in 0.01 seconds, as opposed to 0.39 seconds. i.e. it is 39 times faster that the previous code. 

You can see the result using the following code in the command window:
>> figure; imshow(ufx/255 +0.5);

I got the following result:


Partial x derivative of Lenna


Q: Can you write a fast code to find the partial derivative in the y-direction?

No comments:

Post a Comment