
%% Clear workspace
clc

%%     This script determines the rind thickness, major diameter, 
%      and minor diameter, of a stained hollow plant stalk cross-section.

% %     The script also creates Python files that contain the rind and pith 
%       spline coordinates. These Python files can be run in a third-party 
%       software (e.g, Abaqus/CAE) to create sketches.
      
%% %%%%% INPUTS %% %%%%%
inputImageFolderPath = 'C:\Users\odun3492\Documents\microscopeImages';      % insert folder path inside single quotes
spatialCalibrationFactor = 1/541;     % input spatial calibration factor (in mm/pixel) after the slash('/')

%% %%%%% OUTPUTS %% %%%%%
% %    A table that contains the determined phenotype measurements.
% %    Python files that contain the rind and hollow spline coordinates 
% %    A figure that shows the rind and hollow splines overlaid on the grayscale image of the analyzed cross-sectional input image.


%% Set working folder %%
identifyFiletypesToRead = fullfile(inputImageFolderPath,'*.tif');      % identify folder files with .tif extension
listOfStalkCrossSections = dir(identifyFiletypesToRead);      % create a list of folder files with .tif extension
mkdir Output\    % create new output folder within working folder 

%% Analyze each Stalk Cross-sectional Input Image %%

for i = 1:length(listOfStalkCrossSections)        % identify each captured stalk cross-section from working folder
    stalkCrossSectionName = listOfStalkCrossSections(i).name;     % specify stalk cross-section currently identified
    stalkCrossSectionFolder = fullfile(listOfStalkCrossSections(i).folder,stalkCrossSectionName);       % extract specified stalk sross-section from working folder
    
%% Pre-processing, Segmentation and Smoothing
    
    rawStalkCrossSectionImage = imread(stalkCrossSectionName);        % import raw cross-section image for processing
    convertedGrayscaleImage = rgb2gray(rawStalkCrossSectionImage);     % convert rgb image to grayscale image
    redLayerOfImage = rawStalkCrossSectionImage(:,:,1);     % segment red colour layer of image
    redLayerThreshold = graythresh(redLayerOfImage);        % find suitable threshold value to binarize image         
    binarizedRedLayer = imbinarize(redLayerOfImage,redLayerThreshold);      % binarize image using threshold value 
    greenLayerOfImage = rawStalkCrossSectionImage(:,:,2);       % segment green colour layer of image 
    greenLayerThreshold = graythresh(greenLayerOfImage);        % find suitable threshold value to binarize image
    binarizedGreenLayer = imbinarize(greenLayerOfImage,greenLayerThreshold);        % binarize image using threshold value
    blueLayerOfImage = rawStalkCrossSectionImage(:,:,3);        % segment blue colour layer of image 
    blueLayerThreshold = graythresh(blueLayerOfImage);      % find suitable threshold value to binarize image
    binarizedBlueLayer = imbinarize(blueLayerOfImage,blueLayerThreshold);       % binarize image using threshold value
    binarizedImage = binarizedRedLayer|binarizedGreenLayer|binarizedBlueLayer;      % combine all binaized image layers to get final binarized image
    
    segmentedHollowArea = imclearborder(binarizedImage);        % remove rind area from the binarized image 
    hollowArea = bwareafilt(segmentedHollowArea,1);     % extract the hollow area from binarized image
    hollowBoundaryCoordinates = bwboundaries(hollowArea);       % trace the boundaries of the hollow area and extract coordinates 
    [~,newHollowBoundaryCoordinates] = sort( cellfun( @length, hollowBoundaryCoordinates ), 'descend');     % sort obtained cell array of boundaries from largest to smallest using matrix dimension
    hollowBoundaryCoordinates =  hollowBoundaryCoordinates(newHollowBoundaryCoordinates);       % sort boundaries from largest to smallest using above order 
    hollowBoundaryCoord = hollowBoundaryCoordinates{1};     % extract coordinates of hollow boundary    
    hollowBoundaryXCoord = hollowBoundaryCoord(:,2);        % extract x-coordinates of hollow boundary
    hollowBoundaryYCoord = hollowBoundaryCoord(:,1);        % extract y-coordinates of hollow boundary 
    uniqueHollowBoundaryCoord = unique(hollowBoundaryCoord(:,1:2),'rows','stable');     % remove repeated hollow boundary coordiantes 
    hollowBoundaryXCoord = uniqueHollowBoundaryCoord(:,2);      % extract x-coordinates of hollow boundary
    hollowBoundaryYCoord = uniqueHollowBoundaryCoord(:,1);      % extract y-coordinates of hollow boundary 
    hollowSplineCoordinates = boundary(hollowBoundaryXCoord,hollowBoundaryYCoord);      % extract spline from hollow boundary coordinates 
    hollowSplineXCoord = hollowBoundaryXCoord(hollowSplineCoordinates);     % extract x-coordinates of hollow spline 
    hollowSplineYCoord = hollowBoundaryYCoord(hollowSplineCoordinates);     % extract y-coordinates of hollow spline 
    hollowSplinePixelCoord = [hollowSplineYCoord,hollowSplineXCoord];       % combine hollow x and y spline coordinates 
    hollowAreaSpatialReferenceObject = imref2d(size(hollowArea),spatialCalibrationFactor,spatialCalibrationFactor);     % create properties used to convert from pixel to mm
    for ihollow = 1:1:length(hollowSplinePixelCoord)        % extract each row of hollow spline coordinates
        [hollowSplinemmCoord(ihollow,2),hollowSplinemmCoord(ihollow,1)] = intrinsicToWorld(hollowAreaSpatialReferenceObject,hollowSplinePixelCoord(ihollow,2),hollowSplinePixelCoord(ihollow,1));       % convert pixels to mm
    end

    rindArea = binarizedImage;      % make rind area equal to binarized image whole area
    rindArea(hollowArea) = 0;       % extract hollow area from whole stalk area
    rindArea = bwareaopen(rindArea,5000);       % remove pixel clusters from whole stalk area image 
    rindArea = imcomplement(rindArea);      % invert stalk whole cross-sectional area
    rindArea = bwareafilt(rindArea,1);      % extract whole cross-sectional area again to make sure no pixel clusters (noise)
    rindBoundaryCoordinates = bwboundaries(rindArea);       % trace the boundaries of the rind and extract coordinates 
    rindBoundaryCoord = rindBoundaryCoordinates{1};     % extract rind boundary coordinates
    rindBoundaryXCoord = rindBoundaryCoord(:,2);        % extract x-coordinates of rind boundary
    rindBoundaryYCoord = rindBoundaryCoord(:,1);        % extract y-coordinates of rind boundary
    uniqueRindBoundaryCoord = unique(rindBoundaryCoord(:,1:2),'rows','stable');     % remove repeated rind boundary coordiantes 
    rindBoundaryXCoord = uniqueRindBoundaryCoord(:,2);      % extract x-coordinates of rind boundary
    rindBoundaryYCoord = uniqueRindBoundaryCoord(:,1);      % extract y-coordinates of rind boundary
    rindSplineCoordinates = boundary(rindBoundaryXCoord,rindBoundaryYCoord);        % extract spline from rind boundary coordinates 
    rindSplineXCoord = rindBoundaryXCoord(rindSplineCoordinates);       % extract x-coordinates of rind spline
    rindSplineYCoord = rindBoundaryYCoord(rindSplineCoordinates);       % extract y-coordinates of rind spline
    rindSplinePixelCoord = [rindSplineYCoord,rindSplineXCoord];     % combine rind x and y spline coordinates 
    rindAreaSpatialReferenceObject = imref2d(size(rindArea),spatialCalibrationFactor,spatialCalibrationFactor);     % create properties used to convert from pixel to mm
    for irind = 1:1:length(rindSplinePixelCoord)        % extract each row of rind spline coordinates
        [rindSplinemmCoord(irind,2),rindSplinemmCoord(irind,1)] = intrinsicToWorld(rindAreaSpatialReferenceObject,rindSplinePixelCoord(irind,2),rindSplinePixelCoord(irind,1));     % convert pixels to mm
    end
    
    rindSplinemmXCoord = rindSplinemmCoord(:,2);        % x-coordinates of rind spline in mm 
    rindSplinemmYCoord = rindSplinemmCoord(:,1);        % y-coordinates of rind spline in mm
    hollowSplinemmXCoord = hollowSplinemmCoord(:,2);        % x-coordinates of hollow spline in mm 
    hollowSplinemmYCoord = hollowSplinemmCoord(:,1);        % y-coordinates of hollow spline in mm 
    
%% Error Check %% 
    noOfRindSplinePixelCoord = length(rindSplinePixelCoord);        % count number of rind spline coordinates
    noOfRindSplinemmCoord = length(rindSplinemmCoord);      % count number of rind spline coordinates after conversion to mm 
    rindConversionErrorCheck = noOfRindSplinemmCoord - noOfRindSplinePixelCoord;         % check that number of converted coordinates are similar to original
    rindConversionErrorCheck2 = noOfRindSplinemmCoord - (rindConversionErrorCheck-1);       % check if conversion increases coordinates 
    if rindConversionErrorCheck~=0      % if number of converted coordinates are not equal to original
        rindSplinemmCoord(rindConversionErrorCheck2:end,:) = [];        % delete extra coordinates created from conversion 
    end
    
    noOfHollowSplinePixelCoord = length(hollowSplinePixelCoord);        % count number of hollow spline coordinates
    noOfHollowSplinemmCoord = length(hollowSplinemmCoord);      % count number of hollow spline coordinates after conversion to mm
    hollowAreaConversionErrorCheck = noOfHollowSplinemmCoord - noOfHollowSplinePixelCoord;      % check that number of converted coordinates are similar to original
    hollowAreaConversionErrorCheck2 = noOfHollowSplinemmCoord - (hollowAreaConversionErrorCheck-1);     % check if conversion increases coordinates
    if hollowAreaConversionErrorCheck~=0        % if number of converted coordinates are not equal to original
        hollowSplinemmCoord(hollowAreaConversionErrorCheck2:end,:) = [];        % delete extra coordinates created from conversion
    end
    rindSplinemmXCoord = rindSplinemmCoord(:,2);        % x-coordinates of rind spline in mm
    rindSplinemmYCoord = rindSplinemmCoord(:,1);        % y-coordinates of rind spline in mm 
    hollowSplinemmXCoord = hollowSplinemmCoord(:,2);        % x-coordinates of hollow spline in mm 
    hollowSplinemmYCoord = hollowSplinemmCoord(:,1);        % y-coordinates of hollow spline in mm 
    
%% Plot Rind and Pith Splines %% 

    imshow(convertedGrayscaleImage)     % display grayscale to overlay splines on
    hold on
    plot(rindSplineXCoord,rindSplineYCoord,'g','Linewidth',2)       % plot rind splines 
    plot(hollowSplineXCoord,hollowSplineYCoord,'r','Linewidth',2)       % plot hollow splines 
    set(gca,'YDir','reverse');      % reverse y-axis to get correct orientation of splines
    
%% Determine the Rind Thickness %%

    [pairwiseDistance,S]=pdist2(hollowSplinemmCoord,rindSplinemmCoord,'euclidean','Smallest',1);        % find the shortest pairwise distance between cooridinates of the hollow and rind boundary and their indices
    rindThickness=mean(pairwiseDistance);       % calculate the mean of the shortest pairwise distances to determine rind thickness 

%% Determine Minor Diameter and Major Diameter %%
    
    rindAreaCentroid=regionprops(rindArea,'Centroid');      % find the centroid of the rind area 
    rindAreaCentroid=cat(1,rindAreaCentroid.Centroid);      % convert the centroid into a matrix 
    centroidXCoordinate=rindAreaCentroid(:,1);      % x-coordinate of rind area centroid 
    centroidYCoordinate=rindAreaCentroid(:,2);      % y-coordinate of rind area centroid 
    [centroidXmmCoordinate,centroidYmmCoordinate]=intrinsicToWorld(rindAreaSpatialReferenceObject,centroidXCoordinate,centroidYCoordinate);     % convert centroid coordinates to mm
    centroidXCoordinate=int16(centroidXCoordinate);     % convert x coordinate of cenroid to whole number
    centroidYCoordinate=int16(centroidYCoordinate);     % convert y coordinate of cenroid to whole number 
    centroidVertLine=rindArea(:,centroidXCoordinate);       % extract vertical centerline from rind area matrix using centroid x-coordinate  
    centroidHorzLine=rindArea(centroidYCoordinate,:);       % extract horizontal centerline from pith area matrix using centroid y-coordinate 
    centroidTopVertCoord=find(centroidVertLine,1,'first');      % find top coordinate of vertical centre line      
    centroidBottomVertCoord=find(centroidVertLine,1,'last');        % find bottom coordinate of vertical centre line 
    centroidleftHorzCoord=find(centroidHorzLine,1,'first');     % find first(left) coordinate of horizontal centre line 
    centroidrightHorzCoord=find(centroidHorzLine,1,'last');     % find last (right) coordinate of horizontal centre line 
    [topCoordVertLine,bottomCoordVertLine]=intrinsicToWorld(rindAreaSpatialReferenceObject,centroidTopVertCoord,centroidBottomVertCoord);       % convert centroid coordinates to mm
    [leftCoordHorzLine,rightCoordHorzLine]=intrinsicToWorld(rindAreaSpatialReferenceObject,centroidleftHorzCoord,centroidrightHorzCoord);       % convert centroid coordinates to mm
    hold on
    plot([centroidXCoordinate,centroidXCoordinate],[centroidTopVertCoord,centroidBottomVertCoord],'Linewidth',1.5);     % plot vertical centre line
    plot([centroidleftHorzCoord,centroidrightHorzCoord],[centroidYCoordinate,centroidYCoordinate],'Linewidth',1.5);     % plot horizontal centre line
    
    MajorDia = bottomCoordVertLine - topCoordVertLine;      % determine major diameter
    MinorDia = rightCoordHorzLine - leftCoordHorzLine;      % determine minor diameter
  
%% Save Outputs in New Output Folder %%
%     mkdir Output\    % create new folder within working folder 
    newFolderName = '\Output\';      % string used to create output folder path 
    outputFolderPath = strcat(inputImageFolderPath,newFolderName);     % create path for output folder
    [~,inputImageFileName,ext] = fileparts(stalkCrossSectionName);       % extract parts of the file name
    outputImageFileName = strcat(inputImageFileName,'OutputImage');     % append OutputImage to output image file name

%% Create Table Containing Determined Phenotype Values %%
    tableSize = [0,4];      % specify table size
    tableVarTypes = {'cell','double','double','double'};     % specify table column variables 
    tableVarNames = {'fileName','rindThickness','majorDiameter','minorDiameter'};     % specify table column headers
    phenotypesMeasured = table('Size',tableSize,'VariableTypes',tableVarTypes,'VariableNames',tableVarNames);      % create table using specifications
    structnewRow(i,1).fileName = inputImageFileName;     % fill first column of new row with plant stalk name
    structnewRow(i,1).rindThickness = rindThickness;      % fill second column of new row with rind thickness
    structnewRow(i,1).majorDiameter = MajorDia;        % fill third column of new row with major diameter 
    structnewRow(i,1).minorDiameter = MinorDia;        % fill fourth coulmn of new row with minor diameter
    newTable = [phenotypesMeasured;struct2table(structnewRow)];   % add new row to existing table

%% Accept or Reject Analyzed Plant Stalk Cross-section
  
    userPrompt = questdlg('Would you like to accept these rind and hollow boundary lines?', 'Accept or Reject rind and hollow boundary lines','Accept','Reject','Accept');     % create question dialogue box for user to accept or reject splines
    %  response after selecting an option from question dialogue
    switch userPrompt
        case 'Accept'   % response if splines are accepted
            
            %% Save Output Image %% 
            
            saveas(gcf,fullfile(outputFolderPath,outputImageFileName),'tif')        % save output image in output folder
            
            %% Write Python File to Sketch Rind Spline %%

            rindSpline = rindSplinemmCoord;     % set spline coordinates 
            splineName='RindSpline';      % suffix for filename
            filename = sprintf('%s%s.py',inputImageFileName,splineName);      % set filename for rind spline     
            filename = strcat(outputFolderPath,filename);       % concatenate folder path to save file
            fileID = fopen(filename,'w');       % start a new text file
            fprintf(fileID,'from sketch import *\n\n');     % print Abaqus statement to import sketch module 
            fprintf(fileID,'mdb.models[''Model-1''].ConstrainedSketch(name=''__profile__'', sheetSize=200.0)\n');       % print Abaqus statement to create sketch and set sheet size
            fprintf(fileID,'mdb.models[''Model-1''].sketches[''__profile__''].Spline(points=(\n\t');        % print Abaqus statement to sketch splines
            for irindspline = 1:1:size(rindSpline,1)        % run through each row of spline coordinates
                formatSpec = '(%4.4f, %4.4f), ';        % set the amount of significant figures for spline coordinates
                fprintf(fileID,formatSpec,rindSpline(irindspline,1),rindSpline(irindspline,2));     % print spline coordinates to file
            end
            fprintf(fileID,'\t))\n');
            fprintf(fileID,'mdb.models[''Model-1''].sketches.changeKey(fromName=''__profile__'', toName=''rindspline'')\n');
            
            %% Write Python File to Sketch Hollow Spline %%

            hollowSpline = hollowSplinemmCoord;     % set spline coordinates 
            splineName='HollowSpline';        % suffix for filename
            filename = sprintf('%s%s.py',inputImageFileName,splineName);      % set filename for hollow spline
            filename = strcat(outputFolderPath,filename);       % concatenate folder path to save file
            fileID = fopen(filename,'w');       % start a new text file
            fprintf(fileID,'from sketch import *\n\n');     % print Abaqus statement to import sketch module 
            fprintf(fileID,'mdb.models[''Model-1''].ConstrainedSketch(name=''__profile__'', sheetSize=200.0)\n');       % print Abaqus statement to create sketch and set sheet size
            fprintf(fileID,'mdb.models[''Model-1''].sketches[''__profile__''].Spline(points=(\n\t');        % print Abaqus statement to sketch splines
            for ihollowspline = 1:1:size(hollowSpline,1)        % run through each row of spline coordinates
                formatSpec = '(%4.4f, %4.4f), ';        % set the amount of significant figures for spline coordinates
                fprintf(fileID,formatSpec,hollowSpline(ihollowspline,1),hollowSpline(ihollowspline,2));     % print spline coordinates to file
            end
            fprintf(fileID,'\t))\n');
            fprintf(fileID,'mdb.models[''Model-1''].sketches.changeKey(fromName=''__profile__'', toName=''pithspline'')\n');        % print Abaqus statement to save sketch
            
        case 'Reject'       % response if splines are rejected
            
            %% Save Output Image and Update Table %% 
            
            saveas(gcf,fullfile(outputFolderPath,outputImageFileName),'tif')        % save output image in output folder

            structnewRow(i,1).fileName = inputImageFileName;     % fill first column of new row with plant stalk name
            structnewRow(i,1).rindThickness = 0;      % fill second column of new row with rind thickness
            structnewRow(i,1).majorDiameter = 0;        % fill third column of new row with major diameter 
            structnewRow(i,1).minorDiameter = 0;        % fill fourth coulmn of new row with minor diameter
            newTable = [phenotypesMeasured;struct2table(structnewRow)];   % add new row to existing table


    end
    
%%
close all
writetable(newTable,'Output\PhenotypeMeasurements.xls')        % save table in an excel file     

end

clc 
clear all