function HW3_2019

% Solution to ENGN1750 HW3 Problem 2 AY 2019-20
% 
 close all
 
% Problem 2.1
% Create the 701x701 reference image
ref_im = randi(100,701,701)/100;
%imshow(ref_im)

% Problem 2.2
% Translate the image 10 pixels vertically and horizontally
 disp(1:701,1:701,1) = -10;
 disp(1:701,1:701,2) = 10;
 def_im = imwarp(ref_im,disp);
 imshow(def_im)
 title('Translated Image for prob 2.1')
 
 
 % Problem 2.3
 % Use image correlation code to find the displacements
  distperpix = 1;
  gridx = distperpix*(100:50:600);
  gridy = distperpix*(100:50:600);
  template_size_pixels = 20;
  [ux,uy] = image_correlation_example(ref_im,def_im,distperpix,gridx,gridy,template_size_pixels);
% 
  [refx,refy] = ndgrid(gridx,gridy);   % Create grid of undeformed coords for contour plot
 figure('position',[200 200 400 800]);
 subplot(2,1,1), contour(refx,refy,ux,'fill','on');
 ylabel({'y (pixels)'},'FontSize',16);
 xlabel({'x (pixels)'},'FontSize',16);
 title({'x displacement'},'FontSize',18);
 colorbar('peer',gca);
 subplot(2,1,2), contour(refx,refy,uy,'fill','on');
 ylabel({'y (pixels)'},'FontSize',16);
 xlabel({'x (pixels)'},'FontSize',16);
 title({'y displacement'},'FontSize',18);
 colorbar('peer',gca);
 

% Problem 2.2 use the image correlation code
distperpix = 1;
gridx = distperpix*(100:50:600);
gridy = distperpix*(100:50:600);
template_size_pixels = 20;
% 
% % Problem 2.3
  [n,m] = size(ref_im);
 eyy = 0.07;
 for i=1:n
    for j=1:m
        y = (n-i);
        disp(i,j,1) = 0;
        disp(i,j,2) = eyy*y;
    end
 end
 def_im = imwarp(ref_im,disp);
 [ux,uy] = image_correlation_example(ref_im,def_im,distperpix,gridx,gridy,template_size_pixels);
 [n,m] = size(ux);
 hspacing = gridx(2)-gridx(1);
 vspacing = gridy(2)-gridy(1);
 [uyx,uyy] = gradient(uy',vspacing,vspacing);

 figure
 contour(refx,refy,uyy','fill','on')
 ylabel({'y (pixels)'},'FontSize',16);
 xlabel({'x (pixels)'},'FontSize',16);
 title({'eyy'},'FontSize',18);
 colorbar('peer',gca);
 
 [n,m] = size(ref_im);
 exy = 0.07;
 for i=1:n
    for j=1:m
        y = (n-i);
        x = j;
        disp(i,j,1) = -exy*y;
        disp(i,j,2) = exy*x;
    end
 end
 def_im = imwarp(ref_im,disp);
 [ux,uy] = image_correlation_example(ref_im,def_im,distperpix,gridx,gridy,template_size_pixels);
 [n,m] = size(ux);
 hspacing = gridx(2)-gridx(1);
 vspacing = gridy(2)-gridy(1);
 [uxx,uxy] = gradient(ux',hspacing,vspacing);
 [uyx,uyy] = gradient(uy',hspacing,vspacing);

 figure
 contour(refx,refy,(uxy'+uyx')/2.,'fill','on')
 ylabel({'y (pixels)'},'FontSize',16);
 xlabel({'x (pixels)'},'FontSize',16);
 title({'exy'},'FontSize',18);
 colorbar('peer',gca); 
 
% Problem 2.3
nu = 0.3;
[n,m] = size(ref_im);
for i=1:n
    for j=1:m
        x = -0.5+(j-1)/(m-1);
        y = (n-i)/(n-1);
        r = sqrt(x^2+y^2);
        if (r==0) 
            theta = 0;
        else
           theta = acos(x/r);
        end
        u1 = -30*sqrt(r)*(1-2*nu+sin(theta/2)^2)*cos(theta/2);
        u2 = 30*sqrt(r)*(2-2*nu-cos(theta/2)^2)*sin(theta/2);
        disp(i,j,1) = u1;
        disp(i,j,2) = u2;
    end
end

def_im = imwarp(ref_im,disp);
[ux,uy] = image_correlation_example(ref_im,def_im,distperpix,gridx,gridy,template_size_pixels);

 [n,m] = size(ux);
 hspacing = gridx(2)-gridx(1);
 vspacing = gridy(2)-gridy(1);
 [uxx,uxy] = gradient(ux',hspacing,vspacing);
 [uyx,uyy] = gradient(uy',hspacing,vspacing);
 for i=1:m
     for j=1:n
         ematrix = [uxx(i,j),(uxy(i,j)+uyx(i,j))/2.;...
                   (uxy(i,j)+uyx(i,j))/2.,uyy(i,j)];
         pvals = eig(ematrix);
         pstrain(i,j) = max(pvals);
     end
 end

 figure
 contour(refx,refy,pstrain','fill','on')
 ylabel({'y (pixels)'},'FontSize',16);
 xlabel({'x (pixels)'},'FontSize',16);
 title({'Max Principal Strain'},'FontSize',18);
 colorbar('peer',gca); 

end

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
%    hpoints - horizontal position of grid points for displacement output
%    (physical coords, with origin at bottom left of image)
%    vpoints - 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


[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
