Constrained Weighted Least Squares (Image Interpolation)

4 次查看（过去 30 天）
Royi Avital
Royi Avital 2014-4-16
回答： Matt J 2014-4-17
Hello,
I'd like your assistance with solving a constrained weighted least squares problem.
I can't use `lsqlin` since the size of the images creates matrices beyond the abilities of my computer (Or any PC for that matter).
We're given a matrix I of M by N. On this matrix set of L pixels S is set to be anchors. All the other pixels should be recalculated based on that pixels using the following constrains:
Which is equivalent of:
Where the matrix A just "Selects" the anchor pixels from the interpolated vector.
Now, the trick for most pixels, outside the relevant pixel neighborhood, the weights are zeros. Namely, the equation above could be written as:
Now, assume the neighborhood is something like K x K neighborhood. I think this should allow moving into a sparse problem.
Assuming we have all the W_rs needed, how can I built this problem in a manner MATLAB's can solve it?
Thank You.
Giving an example for w_rs:
Where I is the original values of the pixels at pixel r and s respectively.
Of course it is normalized per pixel.
  4 个评论
显示 2更早的评论
Royi Avital
Royi Avital 2014-4-16
编辑：Royi Avital 2014-4-16
Hi Matt,
Yes, W_rs = W_sr.
Yet it is spatially variant.
I added a definition (One of many) at the question itself.
Thank You.
Matt J
Matt J 2014-4-16
OK. Just to be clear, I didn't ask if W_rs=W_sr. I asked if it has the separable form
W_rs= f_s *f_r
Your recently posted example for W_rs does not have this form.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2014-4-16
编辑：Matt J 2014-4-16
Here's what I'd propose. Note that your constraints are highly trivial and shouldn't require an iterative solver. Just treat the E(p)=I(p) as knowns (which is what they are) and apply the solver exclusively to the unknowns E(~p).
%%Fake initial data
M=200;
N=M;
K=3;
I=rand(M,N);
sig=1;
p=false(1,M*N);
p(1:10:end)=1;
tic; %%Set up equation coefficient matrix
xx=1:M;
yy=1:N;
delta=floor((-K:K)/2); %assume K odd
[x,xdelta,y,ydelta]=ndgrid(xx,delta,yy,delta);
xd=x+xdelta;
yd=y+ydelta;
idx=(xd>=1 & xd<=M & yd>=1 & yd<=N); %legitimate indices
x=x(idx); xd=xd(idx);
y=y(idx); yd=yd(idx);
r=sub2ind([M,N],x,y); r=r(:);
s=sub2ind([M,N],xd,yd); s=s(:);
W=spdiags(exp(-(I(r)-I(s))/2/sig^2),0, length(r),length(r));
A=W*(sparse(1:length(r),r,1)-sparse(1:length(s),s,1));
toc;
tic; %%Solve equations
E=zeros(M*N,1);
E(p)=I(p);
E(~p) = -A(:,~p)\(A*E); %solve
toc;
  9 个评论
显示 7更早的评论
Royi Avital
Royi Avital 2014-4-17
Since I don't really understand how did you achieve this solution, How can I incorporate the normalization?
Thank You for your assistance.
Moreover, How come this solution solves the other problem I linked to?
Matt J
Matt J 2014-4-17
编辑：Matt J 2014-4-17
I think the best thing would be to compute W_rs as a sparse matrix
W=sparse(r,s,exp(-(I(r)-I(s).^2)/2/sig^2, N*M,N*M);
Now you can just normalize the rows to sum to 1
sW=sum(W,2);
nzRows=sW>0;
W(nzRows,:)=bsxfun(@rdivide,W(nzRows,:),sW(nzrows));
and now
[r,s,w]=find(W);
Although, I think that normally all(nzRows)=true.
Note that once you have the data triples r,s,w and/or the matrix W, then either of your objective functions can be expressed almost as written and fed to fmincon,e.g.,
function f=objfun1(E,r,s,w)
f=sum(w.*(E(r)-E(s)).^2);
end
function f=objfun2(E,W)
f=norm(E-W*E).^2;
end

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2014-4-17
Maybe a more transparent way to see the problem is to notice that it can be expressed in the equivalent matrix/vector form,
min E.'*(I-W)*E (*)
ignoring your constraints for a moment. To see the above, take the given form in your post (divided by two) and do the following massaging, where W is the sparse matrix with components W(r,s).
f/2=sum_rs(w_rs.*(E(r)-E(s)).^2)/2;
=sum_rs(w_rs.*(E(r)^2)/2 +sum_rs(w_rs*E(s).^2)/2-sum(w_rs*E(r)*E(s))
=norm(E).^2 - E.'*W*E
=E.'(I-W)*E
which is exactly (*). Thus, once you've generated W using the sparse command, everything is expressible as a simple quadratic program.
Similarly the other version of the problem can be expressed
min. f=norm(E-W*E)^2
or equivalently,
min E*((I-W)^2)*E
So, the only difference between the two problems is that the Hessian 2*(I-W) in the first gets squared in the second 2*(I-W)^2.
To deal with "constraints" in either problem, write the problem generically as
min f(E) = E.'*Q*E/2 (**)
and rewrite E in the partitioned form E=[Eu;Ep] where Eu are the unknown E(r) and Ep are the known E(p). Similarly, parition Q as
Q=[Quu, Qup;Qup.' Qpp]
Substituting into () will lead to an unconstrained quadratic program in Eu
min f(Eu) = Eu.'*Quu*Eu/2 -(Qup*Ep).'*Eu
whcih can be minimized analytically
Eu=Quu\(Qup*Ep)
How computationally hard this all is, depends on the sparsity of Q and hence the sparsity of W. The number of bytes needed to hold the sparse matrix W will be approximately 24*N*M*K^2.
For N-M=4000 and K=7, this comes out to about 17GB. So, you'll need to cut back somewhere, I think.

类别

在 Help Center 和 File Exchange 中查找有关 Interpolation 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by