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);

Cropping images is fun

Now that we know how to display images in Matlab using imshow and imagesc, let's see if we can write some basic codes involving manipulation of pixel values. 


The very first thing that I would like to try is to be able to crop an image... starting with the point (x_TopLeft, y_TopLeft) to (x_BottomRight, y_BottomRight).

Let's write the code below to crop the image: mri.jpg.

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

clear all; 

clc; 

close all;

A=imread('mri.jpg');

x_TopLeft=64;

y_TopLeft=64;

x_BottomRight=192; 

y_BottomRight=192;

M=x_BottomRight-x_TopLeft+1;

N=y_BottomRight-y_TopLeft+1; 

for i=1:M

    for j=1:N

        B(i, j, :)=A(x_TopLeft+i, y_TopLeft+j, :);

    end

end

figure;

imshow(A);

figure;

imshow(B);

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

If you write the above code in the editor file and save it as cropping.m and run it you will get two images, one is the original image and the other is the cropped image.

One thing new in the code is the mysterious appearance of a semicolon in B(i, j, :). It basically gives the element stored in B at the location (i, j) be it a grayscale or RGB. This comes in handy. 


Nevertheless, there are many problems with the above code. 

1. It assumes too much about the user's knowledge of the code. 

2. Every time if you want to crop an image you will essentially have to change the code. That is pain in the neck, and nobody likes it.

3. I used Matlab's imshow function, which I don't quite like in this case. I would like to not only crop the image, but also would like to see the cropped part closely.
4. What happens if the points for cropping give a trivial image?

... and so on and so on.

So, let's write a function instead which takes image or the name of the image. Below is the code for doing that.

Here, we are also writing another function new_imagesc, in the same file named cropping_fun.m

Note, it this case one must end both the functions with the keyword end.

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

function cropping_fun(image, x_TopLeft, y_TopLeft, x_BottomRight, y_BottomRight)

% Author: Prashant Athavale

% Date: 11/19/2011

% Please acknowledge my name if you use this code

% Function to crop images

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

clc;

close all;

% let us handle cases if the user inputs incorrect or insufficient inputs

if nargin>0

    if (ischar(image))&&(exist(image, 'file')==2)

        A=double(imread(image));

    elseif exist('image', 'var')

        A=double(image);

    end

    size_A=size(A);

end

if nargin==0

    A=double(imread('mri.jpg'));

    size_A=size(A);

    x_TopLeft=ceil(size_A(1)/4);

    y_TopLeft=ceil(size_A(2)/4);

    x_BottomRight=floor(3*size_A(1)/4);

    y_BottomRight=floor(3*size_A(2)/4);

elseif nargin==1

    x_TopLeft=ceil(size_A(1)/4);

    y_TopLeft=ceil(size_A(2)/4);

    x_BottomRight=floor(3*size_A(1)/4);

    y_BottomRight=floor(3*size_A(2)/4);

elseif nargin==2

    y_TopLeft=ceil(size_A(2)/4);

    x_BottomRight=floor(3*size_A(1)/4);

    y_BottomRight=floor(3*size_A(2)/4);

elseif nargin==3

    x_BottomRight=floor(3*size_A(1)/4);

    y_BottomRight=floor(3*size_A(2)/4);

elseif nargin==4

    y_BottomRight=floor(3*size_A(2)/4);

end

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

% size of the cropped image 

M=x_BottomRight-x_TopLeft+1;

N=y_BottomRight-y_TopLeft+1;

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

% Memory allocation for the output image for faster implementation of the

% code

image_dim=size(size_A);

if image_dim(2)==3

    P=3;

    B=zeros(M, N, P);

else

    B=zeros(M, N);

end

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

if ((M>0)&&(N>0))

    

    for i=1:M

        for j=1:N

            B(i, j, :)=A(x_TopLeft+i, y_TopLeft+j, :);

        end

    end

    new_imagesc(A);

    title('Original image');

    

    new_imagesc(B);

    title('Cropped image');

else

    error('Check the dimensions of the cropped image.');

end

end % End of the cropping_fun

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

function new_imagesc(A)

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

% Author: Prashant Athavale

% Date 11/19/2011

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

% This is a function to display images at the center of the screen using

% Matlab's imagesc function

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

figure;

screen_size=get(0,'screensize');

screen_width=screen_size(3);

screen_height=screen_size(4);

figure_width=screen_size(4)/2; % width of the figure

figure_height=screen_size(4)/2; % height of the figure

x_lowerbottom=(screen_width-figure_width)/2; 

y_lowerbottom=(screen_height-figure_height)/2;

set(gcf,'position',[x_lowerbottom y_lowerbottom figure_width figure_height]);

image_dim=size(size(A));

% What to do if the image is color/RGB or grayscale

if image_dim(2)==3

    M=max(max(max(A)));

elseif image_dim(2)==2

    M=max(max(A));

end

if M>1

    M=1;

else

    M=255;

end

% imagesc likes the images from 0 to 255

% if they are not then we make them so.

A=double(A);

imagesc(uint8(A*M));

if image_dim(2)==3

    colormap('colormap');

else

    colormap('gray');

end

end % end of function new_imagesc

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

Here is the screenshot of the output of the function cropping_fun:

Creating awesome functions in Matlab is easy

When you open Matlab you will see Editor window, Command window, Workspace, Command history etc. See the screenshot attached here.

Functions


A function is a set of commands. A function is written in an editor window. For our first function let's write a protocol to display an image.

function practical(image)

A=imread(image);

figure;

imshow(A);

end

Here the name of the function is "practical" the input argument is "image".
The name of the file must match the name of the function. Thus you save this file as practical.m
The function is terminated with the keyword "end". This is optional if you have only one function in the file.


Running the function


To run this function write the following in the command window.

>> practical('mri.jpg')

Now, this assumes that the image mri.jpg is in your working directory.


Making your code childproof 


You see, the code that you write is not for you. Someone else is going to use it. The end user is most likely not a programmer, he is a normal person... Usually the functions that you would create will be pretty long. The user doesn't have the time to read what the function does. How it works, how to run it, etc. As a good programmer your goal should be to write a code such that anyone can run it. One can not make assumptions about the knowledge of the user.


Let's see what more the above code assumes.


1. First it assumes that image is the name of the image that you want to display. Well, what if it is already a variable that some other code has stored the image in?


2. It assumes that the user knows how to run the code in Matlab, i.e. by writing practical('mri.jpg') in the command window.  The user might just try to run the code by hitting the "run" button in the editor window.


3. It assumes that the image is uint8 or else the imshow will not work well.  As we noted before to work with images effectively I like to convert them to class double. But if you convert the image with class double, the imshow will give a silly output as it wants the function to range from 0 to 1. 


4. It assumes that the file 'mri.jpg' exists in the Matlab's path.


Childproof function


Below is a completely childproof code. It will work no matter what!

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

function practical(image)

% Author: Prashant Athavale

% This functions displays an image

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

clc; % this clears the Command Window

close all; % this closes previous figures if open

if nargin==0 % nargin gives the number of input argument

    image='mri.jpg'; % if the user gives no input arguments then

end

% let's see if the input is name of a file or is it a variable

if (ischar(image)&&(exist(image, 'file')==2))||(exist('image', 'var')==1)

    

    if (ischar(image)&&(exist(image, 'file')==2))

        A=double(imread(image));

    elseif (exist('image', 'var')==1)&&(ischar(image)==0)

        A=double(image);

    end

    

    if exist('A', 'var')

        size_of_the_image=size(A);

        size_of_size=size(size_of_the_image);

    end

end

if exist('A', 'var')&&(size_of_size(2)==3) % i.e. if the image is color, or RGB image

    if max(max(max(A)))>1

        M=255;

    else

        M=1;

    end

    figure;

    imshow(A/M);

elseif exist('A', 'var')&&(size_of_size(2)==2) % i.e. if the image is  grayscale image

    if max(max(A))>1

        M=255;

    else

        M=1;

    end

    figure;

    imshow(A, [0, M]);

end

if exist('A', 'var')==0

    if (ischar(image))

        text=['The input file: ' image ' does not exist in current directory.'];

        errordlg(text);

    else

        text='The input does not exist in current directory.';

        errordlg(text);

    end

end

end % end of the function practical

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