63 lines
1.8 KiB
Matlab
63 lines
1.8 KiB
Matlab
function [U,rmse] = integration_SCS(p,q,gt)
|
|
% Ref : Direct Analytical Methods for Solving Poisson Equation in
|
|
% Computer Vision problems - PAMI 1990
|
|
%
|
|
% example :
|
|
% p = ones(100,100);
|
|
% q = ones(100,100);
|
|
% u = simchony(p,q);
|
|
%
|
|
% This performs the least square solution to \nabla u = [p,q], i.e. :
|
|
% min \int_\Omega \| \nablua U - [p,q] \|^2
|
|
% where \Omega is square and the natural Neumann boundary condition
|
|
% \mu \cdotp (\nabla u -[p,q]) = 0 is used (2nd order derivatives
|
|
% are neglected), where \mu is the outer
|
|
% normal to \Omega.
|
|
%
|
|
% % Axis : O->y
|
|
% |
|
|
% x
|
|
%
|
|
% The problem comes to solve a linear system Ax=b, where A is a bloc
|
|
% Toplitz matrix. Fast solution is provided by Discrete Cosine Transform
|
|
|
|
% Compute div(p,q):
|
|
px = 0.5*(p([2:end end],:)-p([1 1:end-1],:));
|
|
qy = 0.5*(q(:,[2:end end])-q(:,[1 1:end-1]));
|
|
|
|
% Div(p,q) + Boundary Condition:
|
|
f = px+qy;
|
|
f(1,2:end-1) = 0.5*(p(1,2:end-1)+p(2,2:end-1));
|
|
f(end,2:end-1) = 0.5*(-p(end,2:end-1)-p(end-1,2:end-1));
|
|
f(2:end-1,1) = 0.5*(q(2:end-1,1)+q(2:end-1,2));
|
|
f(2:end-1,end) = 0.5*(-q(2:end-1,end)-q(2:end-1,end-1));
|
|
|
|
f(1,1) = 0.5*(p(1,1)+p(2,1)+q(1,1)+q(1,2));
|
|
f(end,1) = 0.5*(-p(end,1)-p(end-1,1)+q(end,1)+q(end,2));
|
|
f(1,end) = 0.5*(p(1,end)+p(2,end)-q(1,end)-q(1,end-1));
|
|
f(end,end) = 0.5*(-p(end,end)-p(end-1,end)-q(end,end)-q(end,end-1));
|
|
|
|
% Cosine Transform of f:
|
|
fsin = dct2(f);
|
|
|
|
% Denominator:
|
|
[x,y] = meshgrid(0:size(p,2)-1,0:size(p,1)-1);
|
|
denom = (2*cos(pi*x/(size(p,2)))-2) + (2*cos(pi*y/(size(p,1))) - 2);
|
|
Z = fsin./(denom);
|
|
Z(1,1) = 0.5*Z(1,2)+0.5*Z(2,1); % Or whatever...
|
|
|
|
% Inverse Cosine Transform:
|
|
U = idct2(Z);
|
|
|
|
% Ground truth available:
|
|
if (nargin>2)
|
|
moyenne_ecarts = mean(U(:)-gt(:));
|
|
U = U-moyenne_ecarts;
|
|
npix = size(p,1)*size(p,2);
|
|
rmse = sqrt((sum((U(:)-gt(:)).^2))/npix);
|
|
else
|
|
U = U-min(U(:));
|
|
end
|
|
|
|
end
|