function [uN,U] = laplace_deriv(im);

%laplace_deriv   -- compute derivatives: U_x, U_y, V_x, V_y
%                   also U, V and V coefficients(DST2D).
%                   (The Model: Image = U + V)
%
% Inputs:
%    im   -- M x N matrix of given image
%
% Outputs:
%    ux,uy -- (M-2) x (N-2) matrix of derivatives in x, y direction
%             for part U of the above model.
%    vx,vy -- (M-2) x (N-2) matrix of derivatives in x, y direction
%             for part V of the above model.
%    u, v  -- (M-2) x (N-2) matrix. By the model: image = u + v
%    vsc   -- DST 2D coefficients of v.

if(ndims(im)~=2)
  error('IM has to be a 2D array');
end

[m,n] = size(im);

dx = 1/(n-1);
dy = 1/(m-1);

[x,y] = meshgrid([1:n-2]*dx,[1:m-2]*dy);
u = zeros(m-2,n-2);
v = zeros(m-2,n-2);



b1 = im(2:end-1,1);
b2 = im(end,2:end-1);
b3 = im(2:end-1,end);
b4 = im(1,2:end-1);

d = im(1,1);
b = im(1,n) - d;
c = im(m,1) - d;
a = im(m,n) - b - c - d;

u = a*x.*y + b.*x + c.*y + d;
ux= a.*y + b;
uy= a.*x + c;
u1 = a.*[dy:dy:1-dy] + b;
u3 = u1;
u2 = a.*[dx:dx:1-dx] + b;
u4 = u2;

b1 = b1 - c*[1:m-2]'*dy - d;
b4 = b4 - b*[1:n-2]*dx - d;
b3 = b3 - b - (a+c)*[1:m-2]'*dy - d;
b2 = b2 - c - (a+b)*[1:n-2]*dx - d;

% Now for the actual "solution" of the poisson equation
db1 = idst(b1);
db2 = idst(b2);
db3 = idst(b3);
db4 = idst(b4);
lambdan = pi*[1:n-2];
lambdam = pi*[1:m-2]';

%if x>=20, then (sinh(x)-exp(x)/2)/exp(x) is below the double precision.
%we will use this threshold to switch from the sinh value to the
%exponential approximation

p = repmat(lambdam,1,n-2);
q = repmat(lambdan,m-2,1);

%h1 = sinh(p.*(1-x))./sinh(p);
h1 = exp(-p.*x) .* ( 1 - exp(-2*p.*(1-x)))./(1-exp(-2*p));
%h1x = -p.*cosh(p.*(1-x))./sinh(p);
h1x = -p.*exp(-p.*x).*(1+exp(-2*p.*(1-x)))./(1-exp(-2*p));
u11 = -lambdam'.*(1+exp(-2*lambdam'))./(1-exp(-2*lambdam'));
u13 = -lambdam'.*exp(-lambdam').*2./(1-exp(-2*lambdam'));
tmp1 = repmat(pi*[1:m-2],m-2,1).*repmat(y(:,1),1,m-2);
%tmpu4 = lambdam'.*zeros(1,m-2);
tmpu4 = zeros(1,m-2);
tmpu2 = lambdam'.*ones(1,m-2);
tmp1 = repmat(pi*[1:m-2],m-2,1).*cos(tmp1);
tmpu4 = lambdam';
tmpu2 = lambdam'.*cos(tmpu2);
h1 = repmat(db1,1,n-2).*h1;
h1y= tmp1 * h1;
u14 = tmpu4 * h1;
u12 = tmpu2 * h1;
%h3 = sinh(p.*x)./sinh(p);
h3 = exp(p.*(x-1)) .* ( 1 - exp(-2*p.*x))./(1-exp(-2*p));
%h3x = p.*cosh(p.*x)./sinh(p);
h3x = p.*exp(p.*(x-1)).*(1+exp(-2*p.*x))./(1-exp(-2*p));
u31 = lambdam'.*exp(-lambdam').*2./(1-exp(-2*lambdam'));
u33 = lambdam'.*(1+exp(-2.*lambdam'))./(1-exp(-2*lambdam'));

h3 = repmat(db3,1,n-2).*h3;
h3y= tmp1 * h3;

u34 = tmpu4 * h3;
u32 = tmpu2 * h3;
%h2 = sinh(q.*y)./sinh(q);
h2 = exp(q.*(y-1)) .* ( 1 - exp(-2*q.*y))./(1-exp(-2*q));
%h2y = q.*cosh(q.*y)./sinh(q);
h2y = q.*exp(q.*(y-1)).*(1+exp(-2*q.*y))./(1-exp(-2*q));
u22 = lambdan.*(1+exp(-2*lambdan))./(1-exp(-2*lambdan));
u24 = lambdan.*exp(-lambdan).*2./(1-exp(-2*lambdan));
tmp2 = repmat(pi*[1:n-2],n-2,1).*repmat(x(1,:)',1,n-2);
tmpu1 = zeros(1,n-2);
tmpu3 = lambdan.*ones(1,n-2);
tmp2 = repmat(pi*[1:n-2],n-2,1).*cos(tmp2);
tmpu1 = lambdan;
tmpu3 = lambdan.*cos(tmpu3);
h2 = repmat(db2,m-2,1).*h2;
h2x = h2 * tmp2';
u21 = h2 * tmpu1';
u23 = h2 * tmpu3';
%h4 = sinh(q.*(1-y))./sinh(q);
h4 = exp(-q.*y) .* ( 1 - exp(-2*q.*(1-y)))./(1-exp(-2*q));
%h4y = -q.*cosh(q.*(1-y))./sinh(q);
h4y = -q.*exp(-q.*y).*(1+exp(-2*q.*(1-y)))./(1-exp(-2*q));
u44 = -lambdan.*(1+exp(-2*lambdan))./(1-exp(-2*lambdan));
u42 = -lambdan.*exp(-lambdan).*2./(1-exp(-2*lambdan));
h4 = repmat(db4,m-2,1).*h4;
h4x = h4 * tmp2';
u41 = h4 * tmpu1';
u43 = h4 * tmpu3';

h1x = repmat(db1,1,n-2).*h1x;
u11 = db1.*u11';
u13 = db1.*u13';
h3x = repmat(db3,1,n-2).*h3x;
u31 = db3.*u31';
u33 = db3.*u33';
h2y = repmat(db2,m-2,1).*h2y;
u22 = db2.*u22;
u24 = db2.*u24;
h4y = repmat(db4,m-2,1).*h4y;
u42 = db4.*u42;
u44 = db4.*u44;

ux = ux + dst(h1x+h3x) + h2x + h4x;
u1 = u1' + dst(u11+u31) + u21+u41;
u3 = u3' + dst(u13+u33) + u23+u43;

uy = uy + dst(h2y'+h4y')' + h1y + h3y;
u2 = u2 + dst(u22 + u42) + u12 +u32;
u4 = u4 + dst(u24 + u44) + u14 +u34;

uN = zeros(m,n);
uN(2:end-1,1) = u1; uN(2:end-1,end) =u3;
uN(1,2:end-1) = u4; uN(end,2:end-1) =u2;

u = u + dst(h1+h3) + dst(h2'+h4')';
U = im;
U(2:end-1,2:end-1)=u;
% Compute v, vsc
tmp=im;
tmp(:,1)=[];
tmp(:,end)=[];
tmp(1,:)=[];
tmp(end,:)=[];
v=tmp-u;
clear im tmp;

vsc=idst(idst(v)')';

% Compute vx and vy
tmpv1 = repmat(pi*[1:n-2]',1,n-2).*cos(repmat(pi*[1:n-2]',1,n-2).*repmat(x(1,:),n-2,1));
vx = dst(vsc * tmpv1);
tmpv2 = repmat(pi*[1:m-2],m-2,1).*cos(repmat(pi*[1:m-2],m-2,1).*repmat(y(:,1),1,m-2));
vy = dst((tmpv2 * vsc)')';
