function E2 = strain_calculations_BMSEED__20190605()
close all; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slightly edited for comments on Feb 1, 2017 by BMSEED
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Updated Dec 23, 2010 by BMIII
%
%     This function calculates lagrangian strain.  Specify the rest image
% as im1 and the maximum deformed image as im2.  You can paste the name and 
% path below.  It will display both images and allow you to select one 
% point at a time that is common to both images.  The number of points is
% set by the variable num.  FOUR is good in general.  Note that the
% deformed configuration is the LOWER CASE letters (abcd) and undeformed is
% UPPER CASE (ABCD)
%     It will then calculate Lagrangian strain three different ways:
% 1) it calculates the deformation gradient tensor by dividing the deformed 
%        by undeformed locations
% 2) by building a matrix of segment lengths and dividing
% 3) by dividing the deformed area by the undeformed area and turning into 
%        strain
% The first 2 methods give both normal strains and shear strain and
% represent the least squares fit of the overconstrained problem.
% The 3rd method assumes equal biaxial strain and zero shear strain
%     It will print all three on screen, but returns method 1 in variable 
% E2.  It also prints the location of the points on both images so that you
% can visually verify that you picked the same points in both.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% SPECIFYING IMAGES TO ANALYZE AND SETTING PARAMETERS

% Select undeformed ('initial') and deformed ('final') images
[filename1, pathname1] = uigetfile({'*.jpg';'*.jpeg';'*.png';'*.tif'}, ...
                'Select the image where the stretch is about to start');
im1name = fullfile(pathname1, filename1);
disp(['Undeformed file:   ' filename1]);
[filename2, pathname2] = uigetfile({'*.jpg';'*.jpeg';'*.png';'*.tif'}, ...
                'Select the image where strain appears to be at a max');
im2name = fullfile(pathname2, filename2);
disp(['Deformed file:     ' filename2]);
disp(' ');


% Select file directory
pathname3 = uigetdir;


% Select Excel file to export numbers to
prompt = {'Enter spreadsheet filename (exclude ".xlsx"):'};
dlg_title = 'Spreadsheet filename';
defaultans = {'temp'};
options.Resize = 'on';
answers = inputdlg(prompt, dlg_title, [1 65], defaultans, options);
file = [pathname3 '\' answers{1} '.xlsx'];  % converting to string variable
        % % [filename3, pathname3] = uigetfile({'*.xlsx';'*.xls'}, ...
        % %                     'Select Excel file to export results to');
        % % file = fullfile(pathname3, filename3);
        % % disp(['Results saved to:  ' filename3]);
        % % disp(' ');


% Specifying number of slices, runs per slice, and number ID for each image
prompt1 = 'How many slices / areas of interest?:  ';
sli = input(prompt1);
prompt2 = 'How many iterations (of 4 points) to run for each slice?:  ';
runs = input(prompt2);
disp(' ');


% Calculating the time elapsed between images
prompt3  = 'Frame rate of video (fps):  ';  % video frame rate (fps)
video_fps = input(prompt3);
prompt4  = 'Frame # of 1st image:  ';       % initial, undeformed frame
ini       = input(prompt4);
prompt5  = 'Frame # of 2nd image:  ';       % final, max deformed frame
fin       = input(prompt5);

fr_elapsed = fin - ini;                     % frames elapsed
t = fr_elapsed/video_fps;                   % time elapsed
disp(['Frames elapsed: ' num2str(fr_elapsed)]);
disp(['Time elapsed:   ' num2str(t)   ' sec']);
disp(' ');



%% PREPARATION PRIOR TO DATA ANALYSES

% Specifying parameters for data analysis and recording
num = 4;            % number of points to take from each image
FALSE = 0;
TRUE = 1;
repeat = FALSE;     % if TRUE, takes one extra measurement 
                    %     ("repeats" the function)

% Matrix with RGB values for different iterations                    
color_rgb = [   0   0   1 ;   1   0   0 ;    0   1   0 ; ... 
                1 0.5   0 ; 0.5   0   1 ;    0   0   0 ; ...
                0   1   1 ;   1   0   1 ;  0.4 0.2   0 ; ...
              0.5 0.5 0.5  ];
              % blue,  red,      green      % orange,  purple,  black,
              % cyan,  magenta,  brown      % gray

% Initializing variables where results will be stored 
out_E2 = [];    out_E3 = [];    out_E4 = [];
out_pts = [];   id_col = [];



%% CALCULATING STRAIN USING THREE DIFFERENT METHODS

% Read in both image files and display them in separate figures
figure(1);  im1=imread(im1name);
imshow(im1);    title('Undeformed Image');
figure(2);  im2=imread(im2name);
imshow(im2);    title('Deformed Image');

        
% Nested for loops that will repeat for each slice and for each iteration
for i = 1:sli           % Loop for the # of slices/areas of interest
    
    j = 0;
    while j < runs      % Loop that will run until all iter. are complete
                
        % Adding an empty line at the beginning of the first run
        if (i == 1   &&   j == 0   &&   repeat == FALSE)
            disp(' '); 
        end

        % Getting pre- and post- coordinates (user input) for each image
        out = zeros(4,4);   % pre-allocating for speed
        for k = 1:num
            figure(1);          % pull points from image 1 - undeformed
            [x,y]=ginputc(1, 'Color', 'r', 'LineWidth', 1, ...'
                                           'LineStyle', ':');
            figure(2);          % pull points from image 2 - deformed
            [xpost,ypost]=ginputc(1, 'Color', 'c', 'LineWidth', 1, ...
                                                   'LineStyle', ':');
            out(k,:)=[x,y,xpost,ypost]; % compiling coordinates
        end

        % Splitting x-, y-coordinates of undeformed and deformed states
        A = out(1,1:2)';    B = out(2,1:2)';    % UNDEFORMED config.
        C = out(3,1:2)';    D = out(4,1:2)';
        a = out(1,3:4)';    b = out(2,3:4)';    % DEFORMED configuration
        c = out(3,3:4)';    d = out(4,3:4)';


% -- METHOD 1 --

        % This sets up the problem to give both deformation gradient and
        % rigid body displacement
        undef2 = [A B C D ; 1 1 1 1]; 
        def2   = [a b c d];
        Fp2 = def2/undef2;      % least squares pseudo inverse 
                                % for overconstrained
        F2 = Fp2(:,1:2);                % deformation gradient tensor
        % % % p2 = Fp2(:,3);                  % rigid body displacement
        E2 = 0.5*(F2'*F2-eye(2));       % calculate Lagrangian strain

        % Calculations of strain rates
        srx = E2(1,1)/t;    % strain rate in x (strain_x/time)... [1/s]
        sry = E2(2,2)/t;    % strain rate in y (strain_y/time)... [1/s]

        % Depending on the calculated strains, points might be re-selected
        if (E2(1,1) > 0) && (E2(2,2) > 0)   % x and y are both positive
            j = j + 1;
            repeat = FALSE;

        else    % either x or y is negative
            if repeat == FALSE  % function repeats one more time for iter.
                % j = j;
                repeat = TRUE;
                disp(['    Warning: Slice #' num2str(i) ', Iteration #' ...
                        num2str(j+1) ' must be repeated ... ' ...
                        'negative strain calculated']);
                continue
            else
                j = j + 1;
                repeat = FALSE;
                disp(['    Warning: A negative strain was calculated '...
                           'in the re-do of this iteration']);
            end
        end
 

% -- METHOD 2 --           
        % Solving for Lagrangian strain tensor (from Holmes' HW2 solution)

        % Build a matrix of the length of line segments in both 
        % configurations
        S1=B-A;     S2=C-A;     S3=D-A;     S4=C-B;     S5=D-B;     S6=D-C;
        s1=b-a;     s2=c-a;     s3=d-a;     s4=c-b;     s5=d-b;     s6=d-c;

        % Deformed configuration vector
        def3(1,1) = s1'*s1 - S1'*S1;
        def3(2,1) = s2'*s2 - S2'*S2;
        def3(3,1) = s3'*s3 - S3'*S3;
        def3(4,1) = s4'*s4 - S4'*S4;
        def3(5,1) = s5'*s5 - S5'*S5;
        def3(6,1) = s6'*s6 - S6'*S6;
        def3      = 0.5*def3;
    
        % Undeformed configuration matrix
        undef3(1,:) = [ S1(1)^2   S1(2)^2   2*S1(1)*S1(2) ];
        undef3(2,:) = [ S2(1)^2   S2(2)^2   2*S2(1)*S2(2) ];
        undef3(3,:) = [ S3(1)^2   S3(2)^2   2*S3(1)*S3(2) ];
        undef3(4,:) = [ S4(1)^2   S4(2)^2   2*S4(1)*S4(2) ];
        undef3(5,:) = [ S5(1)^2   S5(2)^2   2*S5(1)*S5(2) ];
        undef3(6,:) = [ S6(1)^2   S6(2)^2   2*S6(1)*S6(2) ];

        temp = undef3\def3;
        E3 = [temp(1), temp(3); temp(3), temp(2)];
 

% -- METHOD 3 --           
        % Strain is determined using the undeformed and deformed 
        % configuration areas.
        % Assumptions:  1) Equibiaxial
        %               2) No shear

        % Norm of a vector gives its length
        areapost  = 0.5*(norm(a-c) .* norm(b-d));   % post-strain area
        areapre   = 0.5*(norm(A-C) .* norm(B-D));   % pre-strain area
        arearatio = areapost/areapre;
        E4 = 0.5*(arearatio - 1);

    
% -- Plots of four selected points in images --

        % Plots of points on undeformed and deformed images
        post = [a, b, c, d, a];
        pre  = [A, B, C, D, A];
        figure(1); hold on;
        plot( pre(1,:),  pre(2,:), 'bo', 'MarkerSize', 15);
        hold off
        figure(2); hold on;
        plot(post(1,:), post(2,:), 'ro', 'MarkerSize', 15);
        hold off

        % Plotting the points of different iterations into one figure
        if      mod(j, 10) ~= 0
                        plot_ind = mod(j, 10);  % to specify line color
                        set_num = floor(j/10);  % to specify figure number
        else
                        plot_ind = 10;
                        set_num = floor(j/10)-1;
        end
        figure(10*i + set_num); 
        hold on; box on;
        title(['Slice #' num2str(i) ...
                    ' - Set of points #' num2str(set_num+1) ]);
        plot(  pre(1,:),  pre(2,:), 'Color', color_rgb(plot_ind,:), ...
                    'Marker', 'o', 'MarkerSize', 5, 'LineStyle', '-', ...
                    'DisplayName', 'pre-'  ); % '.', 15
        plot( post(1,:), post(2,:), 'Color', color_rgb(plot_ind,:), ...
                    'Marker', 'o', 'MarkerSize',  5, 'LineStyle', ':', ...
                    'DisplayName', 'post-' );
        legend( 'show');
        set(gca, 'Ydir', 'reverse');    % to match the image files

    
% -- Repeating an iteration of measurements --

        % Asking user about re-doing an iteration through a pop-up window
        prompt_redo = ['Re-do Slice #' num2str(i) ...
                            ', Iteration #' num2str(j) '?'];
        choice = questdlg(prompt_redo, 'Verification of Points', ...
                                    'Yes', 'No', 'No');
        
        % Removing four clicked points, regardless of user input
        figure(1); h1 = findobj(gca,'type','line');
        delete(h1(1));
        figure(2); h2 = findobj(gca,'type','line');
        delete(h2(1));
                
        % Switch to take into account user input
        switch choice
            case 'Yes'      % iteration will be repeated
                disp(['    Slice #' num2str(i) ', Iteration #' ...
                            num2str(j) ' will be repeated.']);
                j = j-1;
                figure(10*i + set_num);     % plotted line will be deleted
                hh = findobj(gca,'type','line');
                delete(hh(1:2));
            	continue
            case 'No'       % code resumes, user does not exit while loop
                            % (points and results will be stored)
        end
        disp(['Slice #' num2str(i) ', Iteration #' ...
                            num2str(j) ' has been stored.']);


% -- Storing all results an selected points into variables --

        id_col = [id_col; {['S' num2str(i) ', I' num2str(j)]}];

        % Putting the results for each iteration together
        out_E2 = [ out_E2;    E2(1,1),  E2(2,2),  srx,  sry  ];
        out_E3 = [ out_E3;    E3(1,1),  E3(2,2),    0,    0  ];
        out_E4 = [ out_E4;         E4,        0,    0,    0  ];

        % Compiling the points for each iteration
        pts = [ -1, A(1), B(1), C(1), D(1), a(1), b(1), c(1), d(1); ...
                -2, A(2), B(2), C(2), D(2), a(2), b(2), c(2), d(2); ...
                 0,    0,    0,    0,    0,    0,   0 ,   0 ,    0  ];
        out_pts = [ out_pts; pts ];

    end
    
end



%% ADDING AVERAGES TO THE RESULTS TABLE

% Adding two extra rows to the ID and 'results' variables
id_col = [id_col; {'Mean'; 'St_dev'}];
out_E2 = [ out_E2; zeros(2,4) ];
out_E3 = [ out_E3; zeros(2,4) ];
out_E4 = [ out_E4; zeros(2,4) ];

% Calculating the mean of the strain and strain rate measurements
for k = 1:4
    % Mean of iteration results
    out_E2(end-1,k) = mean( out_E2(1:end-2,k) );
    out_E3(end-1,k) = mean( out_E3(1:end-2,k) );
    out_E4(end-1,k) = mean( out_E4(1:end-2,k) );
    
    % Standard deviation of iteration results
    out_E2(end,k)   =  std( out_E2(1:end-2,k) );
    out_E3(end,k)   =  std( out_E3(1:end-2,k) );
    out_E4(end,k)   =  std( out_E4(1:end-2,k) );    
end

% Creating ID columns for each of the three methods
for k = 1:length( out_E2(:,1) )
    id_col_E2{k,1} = strcat(id_col{k,1}, ' - E2');
    id_col_E3{k,1} = strcat(id_col{k,1}, ' - E3');
    id_col_E4{k,1} = strcat(id_col{k,1}, ' - E4');
end

% Putting all the results (strains, points) together
out_Es  = [out_E2;   zeros(1,4);   out_E3;   zeros(1,4);   out_E4];
out_pts = [out_pts; zeros( length(out_Es(:,1))-length(out_pts(:,1)), 9)];
out_f   = [out_Es,  zeros( length(out_Es(:,1)), 1),  out_pts];

% Creating the first column for the output table
first_col = [id_col_E2;  {''};  id_col_E3;  {''};  id_col_E4];



%% OVERLAY PLOTS OF SELECTED POINTS OVER BOTH IMAGES

for i = 1:sli   % loop that runs for every slice / area of interest
    
    xtra_figs_per_sli = floor(runs/10) + 1;
    for j = 1:xtra_figs_per_sli

        % Selecting the axes of the undeformed and deformed photos
        figure(10*i + 100 + j-1);             % copy of undeformed figure 
        s1 = subplot(1,2,1);    imshow(im1);    
        hh1_undef = gca;
        figure(10*i + 100 + j-1);             % copy of deformed figure  
        s2 = subplot(1,2,2);    imshow(im2);  
        hh2_def   = gca;

        % Extracting the plotted lines for both configurations
        figure(10*i + j-1);
        points_set = findobj(gca,'type','line');
        hh1_points_set = points_set(2:2:end);
        hh2_points_set = points_set(1:2:end);

        % Copying the plotted lines onto the photos
        copyobj( hh1_points_set, hh1_undef );
        copyobj( hh2_points_set, hh2_def   );

        % Zooming into area of interest
        x_min = min([hh1_points_set.XData hh2_points_set.XData]) - 10;
        x_max = max([hh1_points_set.XData hh2_points_set.XData]) + 10;
        y_min = min([hh1_points_set.YData hh2_points_set.YData]) - 10;
        y_max = max([hh1_points_set.YData hh2_points_set.YData]) + 10;
        figure(10*i + 100 + j-1);
        subplot(1,2,1);     xlim([x_min x_max]);    ylim([y_min y_max]);
        title(['Slice #' num2str(i) ' - Set of points #' num2str(j) ]);
        subplot(1,2,2);     xlim([x_min x_max]);    ylim([y_min y_max]);

        % Reducing the space between subplots
        d = 0.05;       % distance between images
        g1 = get(s1,'position');
        set(s1,'position',[0   0   0.5   1]);
        set(s2,'position',[g1(1)+g1(3)+d   0   0.5   1]);
        
    end
end



%% COMPILING RESULTS AND CREATING A 'CLEAN' TABLE

% Header for table of results and creating the variable with the table info
first_row = { '',   'strain_x'      ,  'strain_y'       ,  ...
                    'strain_rate_x' ,  'strain_rate_y'  ,  '', ...
              '',   'Pt1_i'         ,  'Pt2_i'          ,  ...
                    'Pt3_i'         ,  'Pt4_i'          ,  ...
                    'Pt1_f'         ,  'Pt2_f'          ,  ...
                    'Pt3_f'         ,  'Pt4_f'      };
out_table = [first_row; first_col num2cell(out_f)];


% Cleaning up the table by converting cells with 0's to blank cells 
for i = 1:length(out_table(:,1))
    for j = 1:length(out_table(1,:))
        if out_table{i,j} ==  0,    out_table{i,j} =  '';    end
        if out_table{i,j} == -1,    out_table{i,j} = 'x';    end
        if out_table{i,j} == -2,    out_table{i,j} = 'y';    end
    end
end
writecell(out_table, file);
disp(' ');
disp(out_table(:,1:5));


% Saving results (i.e. data and figures)
prompt_save_data = 'Save data and figures from this run?';
choice_save_data = questdlg(prompt_save_data, 'Saving results', ...
                            'Yes', 'No', 'No');
if strcmp(choice_save_data, 'Yes')       
    % saving the results in a mat file
    save([pathname3 '\' answers{1} '.mat'], 'out_f', 'out_table');
    
    % saving the figures as png's and fig's
    figHandles = get(groot, 'Children');
    num_figs = length(figHandles);
    for k = 1:(num_figs-2)
        figlabel = num2str(k);
        savefig(figHandles(num_figs-k-1), [pathname3 '\' figlabel '.fig']);
        saveas( figHandles(num_figs-k-1), [pathname3 '\' figlabel '.png']);
    end
end



%% OUTPUT
% out


end
