Create The Image Laplacian Matrix Effectively

Hello,

I want to build the Spatial Laplacian of a given operation on an image.

The Matrix is given by:

The matrix Dx / Dy is the forward difference operator -> Hence its transpose is the backward difference operator.

The matrix Ax / Ay is diagonal matrix with weights which are function of the gradient of the image.

It is defined by:

Where Ix(i) is the horizontal gradient of the input image at the i-th pixel.

As said above Ax(i, j) = 0, i ~= j.

It is the same for Ay with the direction modification.

Assuming input Image G -> g = vec(G) = G(:).

I want to find and image U -> u = vec(U) = U(:) s.t.:

How can I solve it most efficiently in MATLAB?

How should I build the sparse Matrices?

Thank You.

[Matt J]

It looks like deconvreg in the Image Processing Toolbox does the above (or something similar), but without linear algebraic methods, probably.

[Royi Avital]

I don't think deconvreg is related.

Answer 1

[Matt J]

 [M,N]=size(inputImage);

 g=inputImage(:);

 Dx=diff(speye(N),1,1);
  Dx=kron(Dx,speye(M));

 Dy=diff(speye(M),1,1);
  Dy=kron(speye(N),Dy);

    sp=@(V) spdiags(V(:),0,numel(V),numel(V));

 Ax=sp(Dx*g);
 Ay=sp(Dy*g);

 Lg=Dx.'*Ax*Dx + Dy.'*Ay*Dy;

 u=(speye(size(Lg))+lambda*Lg)\g;

[Royi Avital]

Hi Mat,

Thank you for your answer.

Few remarks:

1. Wouldn't be faster to use 'diff' and then put it into the sparse matrix instead of doing it by matrix multiplication?
2. Let's say the weights are given by exponent weight (See my update to the question). How would you do it?

[Matt J]

Hi Royi,

1. Yes, probably.
2. Ax=sp(exp(-(Dx*g)/2/alpha^2)). Or implement Dx*g using diff() as you mentioned.