function [reconDist1, reconDist2, reconValidRange, z] = exp_synDistorted(img1, img2, blockx, blocky, b, deltapsReals)
% Created by Yuandong Tian, Oct. 5, 2009
% Email: yuandong@andrew.cmu.edu
% surface reconstruction using coefficients.

[m, n] = size(img1);
[warp.Xs, warp.Ys] = getWarping(m, n, blockx, blocky, b, deltapsReals);

% try a projection step
z = zeros(size(warp.Xs));
for i = 1:size(warp.Xs, 3)
    z(:, :, i) = frankotchellappa(warp.Xs(:, :, i), warp.Ys(:, :, i));
    % and then take the derivatives
    [xs, ys] = gradient(z(:, :, i));
    
    % that is a scaling factor, for some reasons...so fit the energy them
    energySrc = sum(sum(warp.Xs(:, :, i).^2)) + sum(sum(warp.Ys(:, :, i).^2));
    energyDst = sum(xs(:).^2) + sum(ys(:).^2);
    ratio = sqrt(energySrc / energyDst);
    
    xs = xs * ratio; 
    ys = ys * ratio;
    
    warp.Xs(:, :, i) = xs;
    warp.Ys(:, :, i) = ys;
    
    z(:, :, i) = z(:, :, i) * ratio;
end;

reconDist1 = simulate(img1, warp);
reconDist2 = simulate(img2, warp);
reconValidRange = [blockx(1), blocky(1), blockx(end)-1, blocky(end)-1];
