function [ux,uy] = image_correlation_example(ref,cur,distperpix,gridx,gridy,size_ref)
%    Simple code illustrating how to use MATLAB to measure displacements
%    on the surface of a specimen by correlating images of its surface
%    before and after deformation
%
%    ref - reference image
%    cur - deformed image
%    distperpix - physical distance per pixel
%    gridx - horizontal position of grid points for displacement output
%    (physical coords, with origin at bottom left of image)
%    gridy - vertical position of grid points for displacement output
%    (physical coords)
%    The grid is generated in ROW MAJOR order - i.e. coordinates 
%    of points are (x1,y1), (x1,y2), (x1,y3)... (x1,y_nrows); (x2,y1) etc
%    This is the transpose of what is used, eg, in contours and gradient
%    size_ref Size (in pixels) of small square in ref. image that will be used to track
%    displacements


[n,m] = size(ref);

size_cur = 5*size_ref;  % size of cropped region in current image, pixels

%
refDisplayGrid = ref;    % Reference image with superimposed grid
curDisplayGrid = cur;    % Current image with superimposed grid

% Convert grid to pixel coords
  hpoints = floor(gridx/distperpix);   
  vpoints = floor(gridy/distperpix); 
  vpoints = n-vpoints;                 % Flip columns to put origin @ top left
  [hgrid,vgrid] = ndgrid(hpoints,vpoints);  % Grid pts; row major order
   
   count = 0;
%%
for j=1:size(hgrid,2)
    for i = 1:size(hgrid,1)
        vpoint = vgrid(i,j);
        hpoint = hgrid(i,j);
        
        %   Crop a small piece of the current image
        vlo = max(vpoint-size_cur,1);
        vhi = min(vpoint+size_cur,size(cur,1));
        hlo = max(hpoint-size_cur,1);
        hhi = min(hpoint+size_cur,size(cur,2));
        cursubimage = cur(vlo:vhi,hlo:hhi);
        %   Crop a small piece of the reference image
        vlo = max(vpoint-size_ref,1);
        vhi = min(vpoint+size_ref,size(ref,1));
        hlo = max(hpoint-size_ref,1);
        hhi = min(hpoint+size_ref,size(ref,2));
        refsubimage = ref(vlo:vhi,hlo:hhi);
        %   Correlate
        c = normxcorr2(refsubimage,cursubimage);
        [vpeak,hpeak] = find(c==max(c(:)));
        %   The coords of the top left corner of the current image in the ref
        voffset = vpeak-size(refsubimage,1);
        hoffset = hpeak-size(refsubimage,2);
        vdisp(i,j) = (voffset-(size(cursubimage,2)/2-size(refsubimage,2)/2));
        hdisp(i,j) = (hoffset-(size(cursubimage,1)/2-size(refsubimage,1)/2));
        %   Superimpose rectangles on the two images to illustrate the correlation
        ref_rectangle = [hpoint-size(refsubimage,2)/2,vpoint-size(refsubimage,1)/2,size(refsubimage,2),size(refsubimage,1)];
        refDisplayGrid = insertShape(refDisplayGrid, 'Rectangle', ref_rectangle, 'LineWidth', 2, 'Color',[1,0,0]);
        cur_rectangle = [hoffset+hpoint-size(cursubimage,2)/2,voffset+vpoint-size(cursubimage,1)/2,size(refsubimage,2),size(refsubimage,1)];
        curDisplayGrid = insertShape(curDisplayGrid, 'Rectangle', cur_rectangle, 'LineWidth', 2, 'Color',[1,0,0]);
                    curDisplayGrid = insertShape(curDisplayGrid, 'Rectangle', ref_rectangle, 'LineWidth', 1, 'Color',[0,0,1]);
%       Store pixel coords of each point before and after deformation
        count = count + 1;
        refpoints(count,2) = double(vpoint);
        refpoints(count,1) = double(hpoint);
        defpoints(count,2) = double(vpoint + vdisp(i,j));
        defpoints(count,1) = double(hpoint + hdisp(i,j));
    end
end

%% Plot the grids and displacements
% figure
% imshow(refDisplayGrid)
figure
imshow(curDisplayGrid)


%% Improve the displacement resolution using cpcorr
adjustedpoints = cpcorr(defpoints,refpoints,cur,ref);

displacements = adjustedpoints-refpoints;
n = length(hgrid);
m = length(vgrid);
% Create matrices of displacements (physical coords; origin at bottom left
% of image
ux = distperpix*reshape(displacements(:,1),n,m);
uy = -distperpix*reshape(displacements(:,2),n,m); 
     
     

end