function [R1,R2,off] = viprectification_rectifywarp(I1,I2,H1,H2)
%Stereo rectification warp.
%  [R1,R2,OFF] = VIPRECTIFICATION_RECTIFYWARP(I1,I2,H1,H2) warps images I1
%  and I2 by the 3-by-3 transforms H1 and H2, respectively, padding both
%  resulting images such that their relative repositioning by the
%  transforms is preserved.

%   Copyright 2009 The MathWorks, Inc.

% Compute where corner points for both images appear under reprojection.
t1 = H1 * [1 1 1; size(I1,2) 1 1; 1 size(I1,1) 1; size(I1,2) size(I1,1) 1]';
t1 = bsxfun(@rdivide,t1,t1(3,:));
t2 = H2 * [1 1 1; size(I2,2) 1 1; 1 size(I2,1) 1; size(I2,2) size(I2,1) 1]';
t2 = bsxfun(@rdivide,t2,t2(3,:));
% Take bounds of above corner point location
box = [min(t1(1:2,:),[],2)' max(t1(1:2,:),[],2)';
       min(t2(1:2,:),[],2)' max(t2(1:2,:),[],2)'];
boxt = [floor(min(box(:,1:2),[],1)) ceil(max(box(:,3:4),[],1))];
% Compute output image size and origin.
sz = fliplr(1+(boxt(3:4)-boxt(1:2)));
off = 1-boxt(1:2); % Offset to return for point repositioning.
rectifiedOrigin = fliplr(boxt(1:2));

hGeoTx = vision.GeometricTransformer(...
    'OutputImagePositionSource','Property',...
    'OutputImagePosition',[rectifiedOrigin sz]);
R1 = step(hGeoTx,I1,single(H1([2 1 3],[2 1 3])));
R2 = step(hGeoTx,I2,single(H2([2 1 3],[2 1 3])));
