%% Before starting this routine

% You must use matlab version 2016b or later or the file won't run. You
% will also need to have access to the statistics package in matlab
    
% This file analyzes the time, displacement, and load data in a .csv
% file from a puncture test. Sample raw data .csv files are provided with 
% this matlab file which demonstrate file naming conventions and formating.
% This MATLAB file should be placed in a root directory together with sub
% folders named "Data" (Where the raw data files should be placed) and "Graphs"
% (Where plots of the load-extension curve will be output). After this
% file has finished running, the four output variables (diameter, rind
% thickness, peak force magnitude, and Integrative Puncture Score) will
% be available in the workspace.

%% OPEN RAW DATA FILES AND COMPUTE KEY MEASURMENTS FROM EACH GRAPH
    
%Clear workspace
clear all
close all

%Save all names of raw data files
global starting_folder;
starting_folder = cd();  %save the file path to the current folder as "starting_folder"
dinfo = dir('Data');   %look in the 'Data' folder to get all of the file names

cd('Data') %browse to 'data' folder
   
%analyze each test file
for i = 3:length(dinfo(:,1)) % Start at three because first two entries are an artifact of using the dir command 
    fileName{i} = dinfo(i).name;   % file naming structure = "set_stalk_internode.csv"
    name1=split(fileName(i),"_"); 
    set(i) = str2num(name1{1});             
    name2=split(name1(3),".");
    stalk(i) = str2num(name1{2});
    internode(i) = str2num(name2{1});
 
%run the fucntion 'process data' to analyze each raw data file and return key
%measurments from each graph to the Table T. The function
%definition is provided at the end of this file. Note if using a
%matlabl version prior to 2016b the function definition must be in
%a seperate file.
    if exist(fileName{i}, 'file')==2 %only open the file if it exist
        forceDispData = xlsread(fileName{i},'A3:C2000'); %read data from cell A3 to C2000 and store in 'forceDispData'
        forceDispData(:,3)=abs(forceDispData(:,3));
        %call function to analyse test_data
        [               
            diameter(i),rnd_thickness(i),peaks_magnitude1(i),Integrative_Puncture_Score(i),...
        ]...
        = process_data(forceDispData,i,fileName(i));
    end
cd(starting_folder);
cd('Data');
clear('forceDispData')
end
                    
cd(starting_folder);                   
clear('num', 'num2', 'num3', 'num4','i','j','dinfo','folder_names','starting_folder','test_data','test_names')
    
%% ==========================================================================================================================================
% Fuction Definition
%============================================================================================================================================ 
%pass a three column matrix into the function.  first column = time, second
%column = extension(mm), third column = load (kN).  The function calls this
%matrix the variable 'num'.  The word "index" used below refers to the index location within
%this matrix of load and extension data.

function [diameter, rnd_thickness, peaks_magnitude1, Integrative_Puncture_Score...
         ]=process_data(dataArray,ii,name)
global starting_folder;

%initialize function outputs to NaN incase function terminates before values are assinged
[diameter, rnd_thickness, peaks_magnitude1, Integrative_Puncture_Score...
] = deal(NaN);
 
%set variables used to determine key points in force-disp graphs 
maxExtension = 4;           %maximum extension before probe probe returns to starting point
contactDerivative = 0.5;    %derivative of load data used to find first contact with stalk
contactForce = 1.8;         %force threshold in Newtons used to find first contact with stalk
exitDerivative = - 0.1;     %derivative of load data used to find probe exit (should be negative)
midDerivative = 0.075;      %derivative used to find exit of first rind and entry of second rind
midOffset = 65;
minPeakProminence = 5;      %minimum peak prominence used to find rind tissue
minPeakHeight = 10;         %minimum force (N) to be considered rind tissue
ticker=0;                   %Used to find rind thickness line 174
check=0;                    %Used to find rind thickness line 173
rind_threshold=10;          %change this variable to alter the point selected for calculating rind thickness


%Trim data (eliminate data collected while probe returning to starting position)
%find where data is greater than an extension of 'maxExtension' and delete the rest    
maxExtensionIndex = find(dataArray(:,2) > maxExtension, 1);  %index location max_extension
dataArray(maxExtensionIndex:length(dataArray),:)=[]; %eliminate data from the probe returning
            
%find index of zero extension (this is the location of the backside of the
%stalk, i.e., the side of stalk resting against the flat horizontal plate during puncture tests)
zero_extension_index = find(dataArray(:,2) > 0, 1); 
                
%Find stalk Diameter 
%find first contact using derivatives   
diameterIndex_1 = find(diff(dataArray(:,3))>contactDerivative,1,'first')-3; %find first data point with derivative greater than 'contactDerivative'
diameter_1 = -(dataArray(diameterIndex_1,2)); %'diameter' of stalk = distance from first contanct to zero extension
%find first contact using force threshold
diameterIndex_2 = find(dataArray(:,3)>contactForce,1,'first')-3; %find first data point with force greater than 'contactForce'
diameter_2 = -(dataArray(diameterIndex_2,2)); %'diameter' of stalk = distance from first contanct to zero extension
%Compare diameter measurments
diameterIndex = min([diameterIndex_1,diameterIndex_2]);
diameter = max([diameter_1,diameter_2]);
%error Checking
if diameter < 5 || diameterIndex < 5
    diameter = NaN;
return                
end
            
%Find the two peaks with the highest prominence
[peaks_magnitude,peaks_index,~,peaks_prom] = findpeaks(dataArray(:,3),'MinPeakHeight', minPeakHeight,'MinPeakProminence',minPeakProminence);
 A=[peaks_magnitude,peaks_index,peaks_prom];
            
%error Checking
size(A);
if ans(1)<2
    [peaks_magnitude,peaks_index,peaks_prom]=deal(NaN);
    return
end 
            
[~,I] = sort(A(:,3),'descend'); % sort by prominence and return indices 
B = A(I,:);                     % use column indicies to sort all columns of A
B = B(1:2,:);                   % Keep two most prominent peaks
[~,I] = sort(B(:,2));           % Sort B by 'peaks_index' first peak to appear in data appears first in B
C = B(I,:);                     % use column indices to sort other data in B
%save index location, magnitude and prominence of peaks
peaks_magnitude1 = C(1,1);
peaks_magnitude2 = C(2,1);  
peaks_index = C(:,2);
prom2 = C(2,3);
            
%Find point where probe exits stalk using derivatives marching forwards from second peak
exit_index = peaks_index(2) + midOffset + find( diff( dataArray(peaks_index(2)+midOffset:end,3) ) > exitDerivative, 1, 'first');
        
%error checking  
if isempty(exit_index) 
    return
end

%Find midpoint of stalk
midPointIndex = round(mean([diameterIndex,zero_extension_index]));

%Find trough of first peak marching forwards from first peak
%error checking
if isempty(find( diff( dataArray( peaks_index(1)+midOffset:peaks_index(2),3) ) > (-midDerivative), 1, 'first'))
    return
end
            
trough_index(1) = peaks_index(1) + midOffset + find( diff( dataArray( peaks_index(1)+midOffset:peaks_index(2),3) ) > (-midDerivative), 1, 'first');

%Find trough of second peak marching backwards from second peak
ForcederDifferential=diff(transpose(polyder(dataArray(:,3))));
ticker=round(length(ForcederDifferential)/1.8);
check=ForcederDifferential(ticker+1);

while  check < rind_threshold
    if ticker > zero_extension_index || ticker > peaks_index(2)
        ticker=round(length(ForcederDifferential)/2);
    end
    
    if ticker  > round(length(ForcederDifferential)-175)
        rind_threshold=rind_threshold-10;
        ticker=round(length(ForcederDifferential)/1.5);
    else       
        ticker=ticker+1;                                                       
        check=ForcederDifferential(ticker);
    end
  
end
    trough_index(2)=ticker;
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
%% %Plot graphs and key  points on graph and save png's 
figure
name=split(name,".");
name=name(1);
           
rnd_thickness = -(dataArray(trough_index(2),2) - dataArray(zero_extension_index,2));
         
if 1 > rnd_thickness || rnd_thickness > 2.5    
    plot(dataArray(:,2),dataArray(:,3),'color', [0.5,0.5,0.5], 'Linewidth', 3)    
else 
    plot(dataArray(:,2),dataArray(:,3), 'Linewidth', 3)
end

ylabel('Force (N)')
xlabel('Displacement (mm)')
title([name],'Interpreter','none');
hold on 
set(gcf, 'Position', [500, 500, 600, 300])
plot(dataArray(diameterIndex,2),dataArray(diameterIndex,3),'b*')
plot(dataArray(peaks_index(1),2),dataArray(peaks_index(1),3),'r*')
plot(dataArray(peaks_index(2),2),dataArray(peaks_index(2),3),'r*')
plot(dataArray(exit_index,2),dataArray(exit_index,3),'b*')
plot(dataArray(midPointIndex,2),dataArray(midPointIndex,3),'k*')
plot(dataArray(zero_extension_index,2),dataArray(zero_extension_index,3),'k*')
plot(dataArray(trough_index(1),2),dataArray(trough_index(1),3),'b+')
plot(dataArray(trough_index(2),2),dataArray(trough_index(2),3),'b+')
pause(0.5)
hold off
cd(starting_folder)
cd('Graphs')
print(strcat(name{1},'.png'),'-dpng')
cd(starting_folder);
            
   
%%
%Stretch rind thickness data. This code will normalize the location of the
%peak force to always occur at at zero point.
if peaks_index(2) > zero_extension_index 
    stretchedData2 = dataArray(trough_index(2):peaks_index(2),:);
else %(peak force occurs before zero plane)
    stretchedData2 = dataArray(trough_index(2):zero_extension_index,:);
end
            
l=length(stretchedData2);
stretchedX=linspace(stretchedData2(1,2),0,l);
stretchedData2(:,2)=stretchedX;
stretchedDatam2p=vertcat(dataArray(midPointIndex:trough_index(2),3),stretchedData2(:,3));
positionm2p=flipud(vertcat(dataArray(midPointIndex:trough_index(2),2),stretchedData2(:,2)));
            
           
            
%Calculate elliptical moment of inertia and section modulus
%moi of ellipse = pia^3b/8. We only measured one diameter
%on one axis and we are most interested in correlations so we
%can drop the pi/8.  Additinoal we now the major and minor
%axes of the stalk are highly correlated and are related to
%one another by a constant so droping this constant we get the
%elliptical moi can be estimated as raidus^4. Weight the moi by the 
%force the probe excerted on the stalk. Using the parallel axis theorem, 
%split the force displacement curve into mulitple slices, find 
%the moi of each slice, weight that slice by the force and them sum 
%all of slices. That is moi = r(i)^4*force(i)-r(i-1)^4*force(i-1)  
position = diameter / 2 + stretchedData2(:,2);
moi_m2p4=0;
Integrative_Puncture_Score=0;
if length(position)<10  %error checking
    return
else
    for j = 2:length(stretchedDatam2p) %just trying out weighting diameter to different powers
        moi_m2p4 = moi_m2p4 + abs( positionm2p(j)^4 * stretchedDatam2p(j) - positionm2p(j-1)^4 * stretchedDatam2p(j-1)); 
    end
end

Integrative_Puncture_Score=moi_m2p4/(diameter/2);
          
close all           
end 


