In this post I will explore a bit the question on how to calculate the discrete gradient and the discrete divergence of an image and a vector field, respectively.

Let be a discrete image, i.e. denotes the gray value at the -th pixel. The famous total variation denoising amounts to minimizing the functional

where the integral shall be understood as summation, stand for the gradient, i.e. the vector of the partial derivatives, and stands for the euclidean absolute value.

When using primal-dual methods, it is of crucial importance, that the used operators for the gradient and the divergence are *numerically* adjoint (up to the minus sign). That is, the numerical operations are adjoint in the following sense: If `grad`

is the operation for the gradient and `div`

is the operation for the divergence, then it should hold for any variables and of suitable size and with `gradu = grad(u)`

, `divv = div(v)`

that `sum(gradu(:).*v(:))`

and `-sum(u(:).*divv(:))`

are equal up to numerical precision. Due to the boundary treatment, the internal MATLAB operations `gradient`

and `divergence`

do not fulfill this requirement.

The most common discretization of the gradient uses discrete forward differences and a constant padding at the boundary (which means that Neumann boundary values are applied). In formula, this reads as

The respective adjoints are backward differences with zero boundary treatment (check it!). Apparently, there are many different ways to implement this routines (and the respective adjoints) in MATLAB. Here are four of them:

`For loops:`

Run through all pixels in two for-loops and assign the difference to the output. Of course, one should preallocate the output prior to the loop. But you may probably know the first two rules of MATLAB coding? If not, here they are: 1. Avoid for-loops. 2. Avoid for-loops, seriously. I put the routines into extra functions and created anonymous function to call the gradient and the divergence asgrad = @(u) cat(3,dxp(u),dyp(u)); div = @(V) dxm(V(:,:,1)) + dym(V(:,:,2));

`Shift and subtract:`

MATLAB is great in using vectors and matrices. And to avoid the for-loop one could also implement the forward difference in -direction by shifting the matrix and subtract the original one, i.e.`[u(:,2:end) u(:,end)] - u`

(and similarly for the other differences). Again, I wrote extra functions and used anonymous function as above.`Small sparse matrices from the left and from the right:`

MATLAB is also pretty good with sparse matrices. Since the derivatives in -direction only involve the subtraction of two elements in the same column, one can realize this by multiplying an image from the left with a sparse diagonal matrix with just two non-zero diagonals. Similarly, the derivative in -direction can be realized by multiplying from the right with a suitable matrix. More precisely, this approach is realized byDy = spdiags([-ones(M,1) ones(M,1)],[0 1],M,M); Dy(M,:) = 0; Dx = spdiags([-ones(N,1) ones(N,1)],[0 1],N,N); Dy(N,:) = 0; Dxu = Dx*u; Dyu = u*Dy';

(check it). Note that the adjoint of -derivative if simple the operation

`Dx'*u`

and the operation of the -derivative is`u*Dy`

. Together, the calculation of the gradient and the divergence was done by the anonymous functionsgrad = @(u) cat(3,u*Dx',Dy*u); div = @(V) V(:,:,1)*Dx + Dy'*V(:,:,2);

`Large sparse matrices:`

One could think about the following: Vectorize the image by`U = u(:)`

(which amounts to stacking the columns above each other). Then assemble a large sparse matrix which has just two non-zero diagonals to do the forward (and other) differences. More precisely this can be done by

(with`Dx`

and`Dy`

from above)DX = kron(Dx,speye(M)); DY = kron(speye(N),Dy); DxU = DX*U; DyU = DY*U;

Here, it is clear that the respective adjoint are just the multiplication with the transposed matrices. Here, the anonymous functions are

grad = @(u) [DX*u DY*u]; div = @(V) DX'*V(:,1) + DY'*V(:,2);

The different approaches have different pros and cons. Well, the for-loop only has cons: It is presumably slow and uses a lot indexing which easily leads to bugs. The shift-and-subtract method should go into an extra function to make it easy to use – but this is not necessarily a drawback. For the multiplication with the large matrix, one has to vectorize the image first and every time one wants to look at the result and need to do `reshape(U,N,M)`

. But let’s focus on speed: I implemented all methods, and let the run on square images of different sizes for 50 times (after a warmup) and measures the execution times with `tic`

and `toc`

. The assembly of the matrices did not enter the timing. Moreover, no parallel computation was used – just a single core. Finally, memory was not an issue since the larges matrices (of size ) only need roughly 100MB of memory.

Here is the table with the average times for one calculation of the gradient (in seconds):

For-loop | Shift-subtract | left-right | left | |

100 | 0.0004 | 0.0002 | 0.0001 | 0.0003 |

200 | 0.0018 | 0.0010 | 0.0005 | 0.0015 |

300 | 0.0057 | 0.0011 | 0.0014 | 0.0020 |

400 | 0.0096 | 0.0031 | 0.0022 | 0.0035 |

500 | 0.0178 | 0.0035 | 0.0030 | 0.0054 |

750 | 0.0449 | 0.0114 | 0.0097 | 0.0123 |

1000 | 0.0737 | 0.0189 | 0.0128 | 0.0212 |

1500 | 0.2055 | 0.0576 | 0.0379 | 0.0601 |

2000 | 0.3942 | 0.0915 | 0.0671 | 0.1136 |

2500 | 0.6719 | 0.1571 | 0.1068 | 0.1788 |

and here is the one for the divergences:

For-loop | Shift-subtract | left-right | left | |

100 | 0.0004 | 0.0003 | 0.0002 | 0.0002 |

200 | 0.0018 | 0.0015 | 0.0005 | 0.0008 |

300 | 0.0048 | 0.0016 | 0.0015 | 0.0012 |

400 | 0.0090 | 0.0036 | 0.0020 | 0.0022 |

500 | 0.0158 | 0.0057 | 0.0027 | 0.0035 |

750 | 0.0409 | 0.0132 | 0.0073 | 0.0069 |

1000 | 0.0708 | 0.0238 | 0.0130 | 0.0125 |

1500 | 0.2008 | 0.0654 | 0.0344 | 0.0370 |

2000 | 0.3886 | 0.1285 | 0.0622 | 0.0671 |

2500 | 0.6627 | 0.2512 | 0.1084 | 0.1361 |

As expected, the for-loop is clearly slower and also as expected, all methods basically scale quadratically (doubling amounts to multiplying the running time by four) since the work per pixel is constant. A little surprisingly, the multiplication from the left and from the right is fastest and also consistently a little faster then multiplication from the left with the larger sparse matrix. I don’t know, why the results are that different for the gradient and for the divergence. Maybe this is related to my use on anonymous functions or the allocation of memory?

August 23, 2013 at 4:31 pm

There are some extra tricks I noticed for the shift-and-subtract method, due to the fact that Matlab uses column-major format. Basically, the horizontal and vertical operations can have very different speeds but Matlab has very fast transpose operations, so you can circumvent most issues. Here’s an example.

Here’s an operator for horizontal differences. You can also use “diff”, but I wanted to later compose this (to put it into a Laplacian) so the output had to be the right size, and this also gives me control over boundary conditions.

diff_h = @(X) [zeros(n1,1),X(:,1:end-1)] – [X(:,1:end-1),zeros(n1,1) ]; % fast

Now, here’s the analogous operation for vertical differences. But the problem is that we want to augment a matrix by adding a row of zeros, not a column of zeros, and this is very slow due to the column-major format:

diff_v = @(X) [zeros(1,n2);X(1:end-1,:)] – [X(1:end-1,:);zeros(1,n2) ]; % slow (6x slower)

So instead, use this:

diff_v_t = @(Xt) ([zeros(n2,1),Xt(:,1:end-1)] – [Xt(:,1:end-1),zeros(n2,1) ])’; % fast

and then instead of calling diff_v(X), call diff_v(X’). It requires two extra transposes, but I found that it was still much faster.

For other boundary conditions, you can do similar tricks.

August 23, 2013 at 8:17 pm

Good to know that! I remember that I already heard that MATLAB uses column-major format some time ago but I didn’t know how that effected the concatenation in different directions. I try to add your proposal to the comparison soon.

September 2, 2013 at 9:04 pm

I just tried to reproduce your observations on my laptop but apparently there was not much of a difference between vertical and horizontal derivative by shift-and-subtract (vertical derivatives were about six to seven percent slower). Doing the vertical derivatives with two transposes and a horizontal derivative was about two times slower. Strange. I do not have an idea why that is…

September 23, 2013 at 10:05 am

Hi

I am a novice in matlab. I have an image and I want to compute the gradient and divergence of the image. I was till recently using the gradient command. However I found that it is painfully slow. So I read your blog. Can you tell me the gradient command in matlab does forward difference, backward difference or central difference ?

An can you give me an idea regarding what should I use instead of that so as to make it faster ?

Thanks in advance.

Devraj

September 23, 2013 at 7:37 pm

I would suggest to look at the code (type “edit gradient”). As far as I remember, it uses central differences except at the boundary (where either forward or backward differences are used). Concerning speed: The shift-and-subtract method as described in 2. is very simple to implement, but, as I wrote, it seems that, at least for me, “multiply by small sparse matrices” as described in 3. is even faster.

March 3, 2014 at 9:35 am

In the second Equation (D_xu)_{ij} -> (D_yu)_{ij}

?

March 3, 2014 at 2:27 pm

Of course – thanks!

March 3, 2014 at 9:51 am

Thank you Dirk. I think from this post and from comments, one can see here how to estimate tensor ranks of grad( u) as well as div(u) if the tensor rank of u is given.

July 15, 2015 at 3:00 pm

Hi Dirk,

I have a 3d Image volume(MRI) with a label image(1’s and 0’s). I need to calculte the gradient in and outside the object using Neumann boundary conditions. Would your approach work if use the voxels which are the most outer of my label? And to what value would i have to set them then? I am experiencing a lot of difficulties with this part of my implementation of a registration algorithm. If you could help me out a bit that would make my world a whole lot easier!

thanks in advance!!

Rob

July 15, 2015 at 3:12 pm

I guess it would be easiest to first label the boundary coefficients (detected, e.g., by morphological filtering, i.e by dilation of the labelled region). I think here you should take care in which direction take the derivative, i.e. if you take difference to the right and downwarts (in 2d), you only need boundary coefficients right and below. Then set the value at such a boundary pixel as the mean value of the values above and left the pixel (if both above and left pixel are inside) or just the value above or left (if that is the only neighbor inside). In 3D things get a bit more messy…

July 15, 2015 at 3:42 pm

Thanks i am going to try something like that. It’s quite a messy problem!

November 20, 2017 at 2:38 pm

Thanks for your explanation. However when I implement the algorithm, there is something wrong. I cannot figure out what is the problem. Can u help me out?

Following is my code:

function test

u = im2double(imread(‘imgs/Bishapur.jpg’));

if size(u, 3) == 3

u = rgb2gray(u);

end

grad = @(u) cat(3,dxp(u),dyp(u));

div = @(V) dxm(V(:,:,1)) + dym(V(:,:,2));

gradu = grad(u);

divv = div(gradu);

diff = sum(gradu(:).*gradu(:)) + sum(u(:).*divv(:)) ;

msg = sprintf(‘diff: %f’, diff);

disp(msg);

end

function fx = dxp(f)

fx = [f(:,2:end) f(:,end)] – f;

end

function fy = dyp(f)

fy = [f(2:end,:); f(end,:)] – f;

end

function fx = dxm(f)

fx = f – [f(:,1) f(:,1:end-1)];

end

function fy = dym(f)

fy = f – [f(1,:); f(1:end-1,:)];

end

November 20, 2017 at 2:54 pm

As far as I see, you got the boundary conditions wrong (any gradient has to do something on the boundary and the adjoint has to do the right thing, then).

November 20, 2017 at 3:00 pm

Could u please release your code as a clue? That will be very helpful for me.

November 20, 2017 at 3:10 pm

There is working code here for that task – if you want me to debug your code, then no, I don’t do that…

November 21, 2017 at 4:05 am

Sorry, I didn’t mean to let you debug my code. I just want to know how did you implement `dxp`, `dxm`etc, as I haven’t found source code for these functions in your blog.

After a few more search, I found an implementation which should be correct.

function fx = dxm(f)

fx = zeros(size(f));

fx(:, 2:end-1) = f(:,2:end-1) – f(:, 1:end-2);

fx(:, 1) = f(:, 1);

fx(:, end) = -f(:, end-1);

end

function fy = dym(f)

fy = zeros(size(f));

fy(2:end-1,:) = f(2:end-1,:) – f(1:end-2,:);

fy(1,:) = f(1,:);

fy(end,:) = -f(end-1,:);

end

Thanks again for your effort.