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