function xreg = DepthReg(imgD, imgS, xs, xe, ys, ye, p, P, H)

[m, n] = size(imgD);
imgS = double(imgS);
  
fguess = 720; % max(xdim,ydim);
Aguess = [fguess 0 n/2; 0 fguess m/2; 0 0 1];

%----------------------------
% Normalize the measurements

for i = xs: xe
    for j = ys:ye
        x1 = inv(Aguess)*[j; i; 1];
        IH = inv(H);
        lambda = (IH(4, 4) - IH(1:3, 4)'*p)/((IH(1:3, 1:3)*x1)'*p-IH(4,1:3)*x1);
        X1 = [lambda*x1; 1]; 
%        X1 = inv(H(:,:, 1))*X1;
%        X1 = [X1(1)/X1(4); X1(2)/X1(4);X1(3)/X1(4);1];
%         X2 = [R T; 0 0 0 1]*X1;
        X2 = X1;
%        x2=P(22:24,:)*(inv(H(:,:,8))*X2);
        x2 = P*X2;
        x2 = x2/x2(3);
        x2 = Aguess *x2;
        x = x2(2);
        y = x2(1);
        fx = floor(x);
        cx = ceil(x);
        fy = floor(y);
        cy = ceil(y);
        if (fx < 1 || cx > m || fy < 1 || cy > n)
            imgD(i, j) = 0;
            continue
        else
            imgD(i,j) = imgS(fx, fy);
%             imgD(i, j) = (x - fx)*(y - fy)*imgS(fx, fy) ...
%                 + (x - fx)*(cy - y)*imgS(fx, cy) ...
%                 + (cx - x)*(y - fy)*imgS(cx, fy) ...
%                 + (cx - x)*(cy - y)*imgS(cx, cy);
        end
    end
end

xreg = imgD;