function O=twoImageStitch2(I1,I2)

% I1 = (imread('/Users/JonBumstead/Desktop/gigaTest/test3/DSC_0128.JPG')); % image one
% I2 = (imread('/Users/JonBumstead/Desktop/gigaTest/test3/DSC_0127.JPG')); % image two
scaler = 0.4; % down sample to make the stitching faster

G1 = double(imresize(rgb2gray(I1),scaler)); % grayscale of image one
G1 = G1-min(G1(:));
G1 = G1./max(G1(:));

G2 = double(imresize(rgb2gray(I2),scaler)); % grayscale of image two
G2 = G2-min(G2(:));
G2 = G2./max(G2(:));

%%
offset = 65; % expected offshift of images
rangeShift = 30; % shift through this range of pixels 
offArrayX = round(linspace(offset-rangeShift/2,offset+rangeShift/2,rangeShift)); % array of corresponding pixel offsets
trace1 = G1(:,end); % last column in image one
clear r
% find shift in x-direction;
for p=1:rangeShift

    trace2 = G2(:,offArrayX(p));
    r(p) = mean(abs(trace1-trace2)); % quantify difference in columns of two images
    
end

[~,indX]=min(r); % location of column in image 2 that best matches the last column in image 1
pixShiftX = offArrayX(indX); % corresponding pixel 
trace2 = G2(:,pixShiftX); % column of image that best matches last column in image 1

r2=xcorr(trace1,trace2); % find shift in y-direction using cross correlation
Nr = size(r2,1); 
offArrayY = round(linspace(-Nr/2,Nr/2,Nr));
[~,indY]=max(r2);

pixShiftXo = round(offArrayX(indX)/scaler); % scale up pixel shift to size of original image
pixShiftYo = round(offArrayY(indY)/scaler); % scale up pixel shift to size of original image


%% Create master image out;
[nPy1,nPx1,~]=size(I1);
[~,nPx2,~]=size(I2);
bufferSize = abs(pixShiftYo); % buffer size around image to accommodate shifting image 2
pixShiftYo
O = zeros(nPy1,nPx1); % initiate master output image 
if pixShiftYo>=0
    O = padarray(O,bufferSize,0,'post'); % pad master output image
else
    O = padarray(O,bufferSize,0,'pre'); % pad master output image
end

O = uint8(cat(3,O,O,O)); % make output image RGB

if pixShiftYo>=0
    
    indArrayY = 1+pixShiftYo:nPy1;
    arrY=bufferSize+1:nPy1+bufferSize;

elseif pixShiftYo<0
    
    indArrayY = 1:nPy1+pixShiftYo;
    arrY=bufferSize+1:nPy1+bufferSize;
    
end


B1 = double(I1(indArrayY,end-pixShiftXo+1:end,:)); % Overlapping part in image 1
W1 = linspace(1,0,round(size(B1,2))); % linear weight for merging overlapping region
W1=repmat(W1,size(B1,1),1);

B2 = double(I2(indArrayY,1:pixShiftXo,:)); % overlapping part in image 2
W2 = fliplr(W1); % linear weight in opposite direction 

for p=1:3
    B1(:,:,p) = B1(:,:,p).*W1; % multiple overlapping region in image one by weight function
    B2(:,:,p) = B2(:,:,p).*W2; % multiple overlapping region in image two by weight function
end

B = uint8(sum(cat(4,B1,B2),4)); % merged overlap: a linear combination of image one and image two in overlapping region

I1(indArrayY,end-pixShiftXo+1:end,:)=B; % set overlapping region in image 1 to merged overlap

arrX=1:nPx1;
O(arrY,arrX,:)=I1; % put image 1 into output image

I2(indArrayY,1:pixShiftXo,:)=B; % set overlapping region in image 2 to merged overlap

arrX = (1:nPx2)-pixShiftXo+nPx1;
arrY = arrY+pixShiftYo;
O(arrY,arrX,:)=I2; % put image 2 into output image
