Image processing: Minimizing function (regularized least square problem)

Hello,

I'm trying to minimize this function (by A):

argmin A (|L(A)|^2 + a*||A-B||^2*)

where:

• A is a MxN image
• L is the Laplacian Operator
• .|| is the usual norm operator
• a is a weight parameter
• B is a matrix of size (M+2*k)xN where k is an integer parameter.
• * indicates that we just consider the pixels in the boundary (we want to preserve in A the pixels in the boundary of B).

Maybe the problem has a trivial solution, but I'm absolutely blocked.

If you need more details, it's (4) equation in this paper.

I will be very grateful for any help provided.

Answer 1

[Matt J]

Below is a direct linear algebraic solution. It works by rewriting the problem as

 min ||  H*A(:)-y ||^2

for an appropriate matrix H and vector y. Whether it will be practical for you depends on the magnitudes of M,N and the power of your computer. Also, you will have to modify sqrt(a)*Imn and B(:) to include only the boundary pixels (just delete the rows corresponding to those pixels).

[M,N]=size(A);

Im=speye(M);
In=speye(N);
Imn=speye(M*N);

Dx=diff(Im,2,1);
Dy=diff(In,2,1);
LA=kron(Dx,In)+kron(Im,Dy); %Laplacian operator

H=[LA; sqrt(a)*Imn];
y=[zeros(size(LA,1),1); B(:) ];

A_solution=reshape( H\y, M,N);

[gui_tech]

Thank you for your answer. Actually, I don't understand how do you construct the H matrix. Can you explain the process a little bit?

Thank you again

[Matt J]

Not sure how detailed you want me to be. In your original least squares function, the first term has a Laplacian operator applied to A(:). The second term has sqrt(a)*eye(M*N) applied to A(:). The H matrix is just the vertical stacking of these 2 operators.

BTW, I've been assuming that the "usual norm operator" that you're using is the Frobenius norm.

[gui_tech]

I don't see how do you rewrite the problem. Maybe I have a theorical lack.

By the way, I'm working on your solution, I've seen there is a mistake in my original post: B has size (M+2*k)xN where k is an integer parameter

Yes, I'm using the Frobenius norm.

[Matt J]

B has size (M+2*k)xN where k is an integer parameter

If A and B are of different sizes, how is A-B defined?

[gui_tech]

In the original paper they define the 'boundary condition' as a regularization term. In the original post I wrote this in * remark. A-B are just defined in the boundary pixels (only in top and bottom), we want: A(1:k,:)=B(1:k,:) and A((M-k):M,:)=B(M:M+k,:)

[Matt J]

I don't see how do you rewrite the problem. Maybe I have a theorical lack.

Forget about A being a matrix for a moment. Suppose you just had a vector unknown, x, and a two-term lsq cost function like this:

argmin_x  norm(P*x-p)^2 +norm(Q*x-q)^2

Can you see that this is equivalent to

argmin_x norm([P;Q]*x - [p;q])^2

and can be solved as [P;Q]\[p;q] ? That's all we're doing here, except in your case the vector is A(:) instead of x.

A-B are just defined in the boundary pixels (only in top and bottom), we want: A(1:k,:)=B(1:k,:) and A((M-k):M,:)=B(M:M+k,:)

Then you would need to modify the code as follows

Amask=false(size(A)); 
Amask(1:k,M-k:M],:)=true;

Bmask=false(size(B)); 
Bmask(1:k,M:M+k],:)=true;

H=[LA; sqrt(a)*Imn(Amask,:)]; %definition of H

y=[zeros(size(LA,1),1); B(Bmask) ]; %definition of y

[gui_tech]

Thank you for your answer. With your last reply I understand better what we're doing. I've been testing the algorithm and the results are not as I expected. For a 8bit image (the same that they use in the article, Cameraman) I get bizarre intensity values of the order of -1*10^3. Maybe the error (if there is an error) is in the laplacian operator? I expected the size of the laplacian operator to be M^2xN^2 and it's not. I'm testing different implementations of the operator but I've not obtained anything good yet.

[Matt J]

You need to convert the 8-bit image to double so that MATLAB will use double precision math (not integer math) to do the calculations.

I don't know why you expect the Laplacian operator to be M^2xN^2. The vector A(:) has M*N elements and so the matrix that pre-multiplies it must have M*N columns, not N^2 columns. Similarly, the number of rows should be approximately M*N, since there are M*N pixel locations where the Laplacian gets evaluated. It won't be exactly M*N, because the DIFF command reduces the length of the vector, e.g.

>> diff([1,3,4])

ans =

       2     1

[Matt J]

One possibility, I edited a line of the code from

y=[zeros(size(LA,1),1); B(Bmask,:) ]; %definition of y

to this

y=[zeros(size(LA,1),1); B(Bmask) ]; %remove colon

Make sure you have the 2nd version.

[gui_tech]

I had the good code and I already converted the int8 to double. I am not sure if the results are good (I'll see soon), but your help has been crucial to advance in my program and I'm sure it's the line to follow. I'm very grateful for your patience and your help!

[ZHUN]

Matt J, kron(Dx,In) and kron(Im,Dy) have different size, how did you add them together please?

[OJ27]

if have the following min problem.

I am using young's inequality to say |k|_L1 times(|p|_L2)-||tc||+theta||grad(p)||_L1

I need some help coding this one

[OJ27]

in this case, when I create the H, I would replace sqrt(a)*Imn by k, however I don't know what form I have to put k since they are different dimensions. Should I make k sparse? also, since they are different norms, will that change the implementation

[Matt J]

OLGA, there is no algebraic solution to the the problem you've written. You will have to use an iterative solver, probably ga(). For ga(), you don't need matrix implementations of the different operators. Just implement the objective function in function form and pass an appropriate function handle to the solver.

Answer 2

[Laurentiu]

gui_tech: I am working on the implementation of the paper as well. Were you able to fix the problem with the low intensities in the solution ? I am also getting these -1*10^3 intensities in A.

Any help would be appreciated.

[Matt J]

Possibly, my implementation of the Laplace operator is not what the paper expects. It might be worth trying to use

http://www.mathworks.com/matlabcentral/fileexchange/44669-func2mat-convert-linear-function-to-matrix

to convert whatever routine you normally use to compute L(A) to matrix form.