function [u,v,reL2] = homoDecom(img,r,p,choice)
[m,n] = size(img);
u = zeros(m,n);
v = u;
lm = floor((m-1)/2^p);
ln = floor((n-1)/2^p);
uNx = zeros(m,n);
uNy =uNx;
if choice >= 0,
[uNx,uNy] =NormDev(img,.5);
end
for j = 1:2^p,
    for k = 1:2^p,
        [uN,u0] = laplace_deriv(img(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1)));
        Dx = uNx(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1));
        Dy = uNy(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1));
        [mDx,nDx] = size(Dx);
        D = zeros(mDx,nDx); D([1 mDx],:) = Dx([1 mDx],:); D(:,[1 nDx]) = Dy(:,[1 nDx]);
        if choice >= 0,
            u(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1)) = u0 + biharmNormal(D-uN);
        else
            u(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1)) = u0;
        end
        v1 = img(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1)) - u(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1));
        fv = fftshift(fft2(v1));
        cm = floor((lm-r*lm)/2)+1;
        cn = floor((ln-r*ln)/2)+1;
        fv(1:cm,:) = 0; fv(end-cm+1:end,:) = 0; fv(:,1:cn) = 0; fv(:,end-cn+1:end) = 0;
        v(((j-1)*lm+1):(j*lm+1),((k-1)*ln+1):(k*ln+1)) = real(ifft2(ifftshift(fv)));
    end
end

diff = img - u - v;
reL2 = norm(diff(:))/norm(img(:));
    
