
%% Clear workspace
clc

%%     This script determines the rind thickness, major diameter, minor diameter,
%       and number of vascular bundles of a stained pith-filled plant stalk
%       cross-section.

% %     The script also creates Python files that contain the rind, pith, 
%       and vascular bundles 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/164;     % 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, pith, and vascular bundles spline coordinates 
% %    A figure that shows the rind, pith, and vascular bundles splines overlaid on the grayscale image of the analyzed cross-sectional input image.

%% Set working folder %%
identifyFiletypesToRead = fullfile(inputImageFolderPath,'*.tif');      % identify files with .tif extension
listOfStalkCrossSections = dir(identifyFiletypesToRead);      % create a list of files with .tif extension
mkdir Output\    % create new 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
    thresholdValue = adaptthresh(convertedGrayscaleImage,0.8);     % find suitable threshold value to binarize image
    binarizedImage = imbinarize(convertedGrayscaleImage,thresholdValue);       % binarize image using threshold value
    
    segmentedRindArea = bwareafilt(binarizedImage,1);       % extract the rind(whole cross-sectional area) from binarized image
    segmentedRindArea = ~segmentedRindArea;        % invert stalk whole cross-sectional area
    segmentedRindArea = bwareafilt(segmentedRindArea,1);      % extract whole cross-sectional area again to make sure no pixel clusters (noise)
    rindArea = imfill(segmentedRindArea,'holes');        % remove any pixel clusters found in rind area
    rindBoundaryCoordinates = bwboundaries(rindArea);      % trace the boundaries of the rind and extract coordinates (in cell format)
    rindBoundaryCoordinates = rindBoundaryCoordinates{1};     % extract rind boundary coordinates from cell
    rindBoundaryCoordinates = unique(rindBoundaryCoordinates(:,1:2),'rows','stable');     % remove repeated rind boundary coordiantes
    rindBoundaryXcoordinates = rindBoundaryCoordinates(:,2);      % extract x-coordinates of rind boundary
    rindBoundaryYcoordinates = rindBoundaryCoordinates(:,1);      % extract y-coordinates of rind boundary
    rindSplineCoordinates = boundary(rindBoundaryXcoordinates,rindBoundaryYcoordinates);        % extract spline from rind boundary coordinates
    rindSplineXcoordinates = rindBoundaryXcoordinates(rindSplineCoordinates);     % extract x-coordinates of rind spline
    rindSplineYcoordinates = rindBoundaryYcoordinates(rindSplineCoordinates);     % extract y-coordinates of rind spline
    rindSplinePixelCoordinates = [rindSplineYcoordinates,rindSplineXcoordinates];       % combine rind x and y spline coordinates
    rindSpatialReferenceObject = imref2d(size(rindArea),spatialCalibrationFactor,spatialCalibrationFactor);       % create properties used to convert from pixel to mm
    for irind = 1:1:length(rindSplinePixelCoordinates)       % extract each row of rind spline coordinates
        [rindSplinemmCoordinates(irind,2),rindSplinemmCoordinates(irind,1)] = intrinsicToWorld(rindSpatialReferenceObject,rindSplinePixelCoordinates(irind,2),rindSplinePixelCoordinates(irind,1));     % convert pixels to mm
    end
    
    segmentedPithArea = imclearborder(binarizedImage);      % remove rind area from the binarized image
    segmentedPithArea = bwareafilt(segmentedPithArea,1);        % extract the pith from binarized image
    smoothingObject = strel('disk',10);        % create object to remove pixel clusters
    segmentedPithArea = imclose(segmentedPithArea,smoothingObject);       % remove pixel clusters from pixel area using smoothing object to dialate and erode image
    pithArea = imfill(segmentedPithArea,'holes');     % remove any extra pixel clusters from pith area
    pithBoundaryCoordinates = bwboundaries(pithArea);      % trace the boundaries of the pith and extract coordinates (in cell format)
    pithBoundaryCoordinates = pithBoundaryCoordinates{1};     % extract pith boundary coordinates from cell
    pithBoundaryCoordinates = unique(pithBoundaryCoordinates(:,1:2),'rows','stable');     % remove repeated pith boundary coordiantes
    pithBoundaryXcoordinates = pithBoundaryCoordinates(:,2);      % extract x-coordinates of pith boundary
    pithBoundaryYcoordinates = pithBoundaryCoordinates(:,1);      % extract y-coordinates of pith boundary
    pithSplineCoordinates = boundary(pithBoundaryXcoordinates,pithBoundaryYcoordinates);     % extract spline from pith boundary coordinates
    pithSplineXCoordinates = pithBoundaryXcoordinates(pithSplineCoordinates);       % extract x-coordinates of pith spline
    pithSplineYCoordinates = pithBoundaryYcoordinates(pithSplineCoordinates);       % extract y-coordinates of pith spline
    pithSplinePixelCoordinates = [pithSplineYCoordinates,pithSplineXCoordinates];     % combine pith x and y spline coordinates
    pithSpatialReferenceObject = imref2d(size(pithArea),spatialCalibrationFactor,spatialCalibrationFactor);       % create properties used to convert from pixel to mm
    for ipith = 1:1:length(pithSplinePixelCoordinates)        % extract each row of rind spline coordinates
        [pithSplinemmCoordinates(ipith,2),pithSplinemmCoordinates(ipith,1)] = intrinsicToWorld(pithSpatialReferenceObject,pithSplinePixelCoordinates(ipith,2),pithSplinePixelCoordinates(ipith,1));     % convert pixels to mm
    end
    
    rindSplinemmXCoordinates = rindSplinemmCoordinates(:,2);      % x-coordinates of rind spline in mm
    rindSplinemmYCoordinates = rindSplinemmCoordinates(:,1);      % y-coordinates of rind spline in mm
    pithSplinemmXCoordinates = pithSplinemmCoordinates(:,2);      % x-coordinates of pith spline in mm
    pithSplinemmYCoordinates = pithSplinemmCoordinates(:,1);      % y-coordinates of pith spline in mm
    
%% Error Check %%
    noOfRindSplinePixelCoordinates = length(rindSplinePixelCoordinates);           % count number of rind spline coordinates
    noOfRindSplinemmCoordinates = length(rindSplinemmCoordinates);        % count number of rind spline coordinates after conversion to mm
    rindConversionErrorCheck = noOfRindSplinemmCoordinates-noOfRindSplinePixelCoordinates;        % check that number of converted coordinates are similar to original
    rindConversionErrorCheck2 = noOfRindSplinemmCoordinates-(rindConversionErrorCheck-1);            % check if conversion increases coordinates
    if rindConversionErrorCheck~=0      % if number of converted coordinates are not equal to original
        rindSplinemmCoordinates(rindConversionErrorCheck2:end,:) = [];        % delete extra coordinates created from conversion
    end
    
    noOfPithSplinePixelCoordinates = length(pithSplinePixelCoordinates);      % count number of pith spline coordinates
    noOfPithSplinemmCoordinates = length(pithSplinemmCoordinates);        % count number of pith spline coordinates after conversion to mm
    pithConversionErrorCheck = noOfPithSplinemmCoordinates-noOfPithSplinePixelCoordinates;        % check that number of converted coordinates are similar to original
    pithConversionErrorCheck2 = noOfPithSplinemmCoordinates-(pithConversionErrorCheck-1);     % check if conversion increases coordinates
    if pithConversionErrorCheck~=0      % if number of converted coordinates are not equal to original
        pithSplinemmCoordinates(pithConversionErrorCheck2:end,:) = [];        % delete extra coordinates created from conversion
    end
    rindSplinemmXCoordinates = rindSplinemmCoordinates(:,2);      % x-coordinates of rind spline in mm
    rindSplinemmYCoordinates = rindSplinemmCoordinates(:,1);      % y-coordinates of rind spline in mm
    pithSplinemmXCoordinates = pithSplinemmCoordinates(:,2);      % x-coordinates of pith spline in mm
    pithSplinemmYCoordinates = pithSplinemmCoordinates(:,1);      % y-coordinates of pith spline in mm 

%% Plot Rind and Pith Splines %%
    imshow(convertedGrayscaleImage)     % display grayscale to overlay splines on
    hold on
    plot(rindSplineXcoordinates,rindSplineYcoordinates,'g','Linewidth',2)       % plot rind splines
    plot(pithSplineXCoordinates,pithSplineYCoordinates, 'r','Linewidth',2)       % plot pith splines
    set(gca,'YDir','reverse');      % reverse y-axis to get correct orientation of splines
    
%% Determine the Rind Thickness %%

    [pairwiseDistance,S] = pdist2(pithSplinemmCoordinates,rindSplinemmCoordinates,'euclidean','Smallest',1);      % find the shortest pairwise distance between cooridinates of the pith and rind boundary and their indices
    rindThickness = mean(pairwiseDistance);      % calculate the mean of the shortest pairwise distances to determine rind thickness
    
%% Determine Major Diameter (to inside and outside of ear groove) and Minor 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);       % x-coordinate of rind area centroid
    [centroidXmmCoordinate,centroidYmmCoordinate] = intrinsicToWorld(rindSpatialReferenceObject,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(rindSpatialReferenceObject,centroidTopVertCoord,centroidBottomVertCoord);        % convert centroid coordinates to mm
    [leftCoordHorzLine,rightCoordHorzLine] = intrinsicToWorld(rindSpatialReferenceObject,centroidleftHorzCoord,centroidrightHorzCoord);       % convert centroid coordinates to mm
    
    % % Determine Major Diameter to End of Ear Groove % %
    rindVertMaxandMin = any(rindArea,2);      % find the first non-zero array in rind area matrix
    rindVertMax = find(rindVertMaxandMin,1,'first');       % find the maximum top coordinate of rind area    
    rindVertMin = find(rindVertMaxandMin,1,'last');        % find the minimum bottom coordinate of rind area
    [rindVertMaxmm,rindVertMinmm] = intrinsicToWorld(rindSpatialReferenceObject,rindVertMax,rindVertMin);       % convert 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
    
    MajorDiaToInsideOfGroove = bottomCoordVertLine-topCoordVertLine;      % determine major diameter to middle of ear groove
    MinorDia = rightCoordHorzLine - leftCoordHorzLine;      % determine minor diameter
    MajorDiaToOutsideOfGroove = rindVertMinmm - rindVertMaxmm;       % determine major diameter to end of ear groove

%% Determine the Number of Vascular Bundles %%

    smoothingObject = strel('disk',13);     % create object to remove pixel clusters
    pithVascularBundles = imclose(segmentedPithArea,smoothingObject);     % remove pixel clusters from pith area using smoothing object to dialate and erode image
    invertedvascularBundles = ~pithVascularBundles;      % invert vascular bundle areas
    vascularBundles = bwareafilt(invertedvascularBundles,1);        % extract pith area from binarized image
    vascularBundles = invertedvascularBundles&~vascularBundles;     % extract vascular bundles from pith 
    vascularBundles = bwareafilt(vascularBundles,[2000 10000]);     % extract vascular bundles without pixel clusters from binarized image
    vascularBundleSplines = bwboundaries(vascularBundles);       % trace the boundaries of the vascular bundles and extract coordinates (in cell format)
    noOfVascularBundles = length(vascularBundleSplines);     % count number of vascular bundles extracted
    hold on
    for ivb = 1:length(vascularBundleSplines)      % for each vascular bundle boundary
        bundle = vascularBundleSplines{ivb};       % extract coordinates of each vascular bundle boundary to plot
        plot((bundle(:,2)),(bundle(:,1)),'r','Linewidth',2)       % plot vascular bundle boundaries
    end
    
%% Save Outputs in New Output 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,6];      % specify table size
    tableVarTypes = {'cell','double','double','double','double','double'};     % specify table column variables 
    tableVarNames = {'fileName','rindThickness','majorDiameterToInsideOfGroove','minorDiameter','majorDiameterToOutsideOfGroove','numberOfVascularBundles'};     % 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).majorDiameterToInsideOfGroove = MajorDiaToInsideOfGroove;        % fill third column of new row with major diameter to groove middle
    structnewRow(i,1).minorDiameter = MinorDia;        % fill fourth coulmn of new row with minor diameter
    structnewRow(i,1).majorDiameterToOutsideOfGroove = MajorDiaToOutsideOfGroove;      % fill fifth column of new row with major diameter to groove end
    structnewRow(i,1).numberOfVascularBundles = noOfVascularBundles;        % fill sixth column of new row with total vascular bundle count
    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 pith boundary lines?', 'Accept or Reject rind and pith 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
            
            %% Select, Delete and/or Add Vascular Bundles %%

            vascularBundlesPrompt = questdlg('Would you like to manually delete and/or Add Vascular Bundles', 'Manually Select Vascular Bundles','Yes','No','No');        % create question dialogue box for user to either delete and/or add or save current vascular bundles. 
            % response after selecting an option from question dialogue
            switch vascularBundlesPrompt
                case 'Yes'      % response if user wants to delete and/or add vascular bundles
                    
                    [identifiedBundles,bundleCount] = bwlabel(vascularBundles);       % identify and count the number of vascular bundles extracted from pith area
                    bundleBoundaries = bwboundaries(identifiedBundles);        % trace the boundaries of each vascular bundle and extract coordinates
                    imshow(convertedGrayscaleImage)     % display grayscale image of stalk scross-section to overlay splines on
                    hold on
                    for ibundle = 1:length(bundleBoundaries) % run through each row of vascular bundle boundary coordinates
                        bundleBoundary = bundleBoundaries{ibundle};     % extract boundary coordinates of each vascular bundle
                        plot(bundleBoundary(:,2),bundleBoundary(:,1),'r')       % overlay each vascular bundle boundary on grayscale image
                    end
                    uiwait(msgbox('STEP 1: Select Vascular Bundles you would like to delete. Once complete press Enter. (You must click OK before selecting vascular bundles. If there are no vascular bundles you would like to delete just press enter.)','Delete Vascular Bundles'));        % display message for user to select vascular bundles to be deleted (user needs to select ok for the script to continue running)
                    [deletePointX,deletePointY] = getpts;     % use mouse to select vascular bundle boundaries to delete on displayed grayscale image
                    close all 
                    imshow(identifiedBundles)        % display binarized vascular bundles image to enable user re-select vascular bundles
                    hold on
                    plot(deletePointX,deletePointY,'+');    % plot crosses on selected vascular bundles so user can easily identify them 
                    selectedBundles = bwselect(identifiedBundles,deletePointX,deletePointY);        % create a new binarized image of selected vascular bundles 
                    identifiedBundles(selectedBundles) = 0;       % delete selected vascular bundle splines
                    [newIdentifiedBundles,oldBundleCount] = bwlabel(identifiedBundles);        % identify and count vascular bundles after deletion
                    newBundleBoundaries = bwboundaries(newIdentifiedBundles);      % trace boundaries of each vascular bundle and extract coordinates
                    imshow(convertedGrayscaleImage)     % display grayscale image of stalk scross-section to select vascular bundles to be counted  
                    hold on
                    plot(rindSplineXcoordinates,rindSplineYcoordinates,'g','Linewidth',2)       % overlay rind and pith splines on grayscale image
                    plot(pithSplineXCoordinates,pithSplineYCoordinates, 'r','Linewidth',2)       % overlay rind and pith splines on grayscale image
                    set(gca,'YDir','reverse');      % reverse y-axis to get correct orientation of splines
                    hold on
                    for inewBundle = 1:length(newBundleBoundaries)       % run through each row of vascular bundle boundary coordinates
                        newBundleBoundary = newBundleBoundaries{inewBundle};       % extract boundary coordinates of each vascular bundle
                        plot(newBundleBoundary(:,2),newBundleBoundary(:,1),'r','Linewidth',2)       % overlay each vascular bundle boundary on grayscale image
                    end
                    uiwait(msgbox('STEP 2: Select Vascular Bundles you would like to add and press Enter when complete. (You must click OK before selecting vascular bundles. If there are no Vascular Bundles you would like to add just press enter.)','Add Vascular Bundles'));      % display message dialogue box for user to add vascular bundles to count (user needs to select ok for the script to continue running)
                    [addPointX,addPointY] = getpts;       % use mouse to select vascular bundles to be added
                    circleCentre = [addPointX addPointY];        % use seleccted points to create a center point for plotting circles on selected vascular bunudles
                    newBundleCount = length(addPointX);      % count number of newly selected vascular bundles
                    radii = 30*ones(newBundleCount,1);     % specify radius of circle to be plotted around selected vascular bundles
                    plotCircles = viscircles(circleCentre,radii,'Color','g');       % plot circles around selected vascular bundles
                    totalBundleCount = oldBundleCount + newBundleCount;      % calculate total vascular bundle count

                    %% Replace Number of Vascular Bundles in Table %% 
                    newTable{i,6} = totalBundleCount;       % specify cell to replace total vascular bundle count
                    
                    %% Save Output Image %% 
                    saveas(gcf,fullfile(outputFolderPath,outputImageFileName),'tif')        % save output image in output folder
                   
                    %% Write Python File to Sketch Vascular Bundle Splines %%

                    vbBoundaries = cell(size(newBundleBoundaries));     % set boundary coordinates
                    splinename = 'VascularBundleSplines';     % append suffix to python filename
                    filename = sprintf('%s%s.py',inputImageFileName,splinename);       % set filename for vascular bundle splines
                    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
                    for ivbboundary = 1:length(newBundleBoundaries)        % run through each row of vascular bundle boundary coordinates
                        eachVbBoundary = newBundleBoundaries{ivbboundary};     % extract boundary coordinates of each vascular bundle
                        eachVbBoundary = unique(eachVbBoundary(:,1:2),'rows','stable');        % remove any repeated coordinates in boundary coordinates
                        eachVbBoundaryX = eachVbBoundary(:,2);        % x-coordinates of vascular bundle boundary
                        eachVbBoundaryY = eachVbBoundary(:,1);        % y-coordinates of vascular bundle boundary
                        vbSpline = boundary(eachVbBoundaryX,eachVbBoundaryY);       % extract splines from vascular bundle boundaries
                        vbSplineX = eachVbBoundaryX(vbSpline);        % x-coordinates of vascular bundle spline
                        vbSplineY = eachVbBoundaryY(vbSpline);        % y-coordinates of vascular bundle spline 
                        vbSplines = [vbSplineY*spatialCalibrationFactor,vbSplineX*spatialCalibrationFactor];       % convert pixel to mm
                        for ivbboundary2 = 1:1:length(vbBoundaries)     % run through each row of vascular bundle boundary coordinates 
                            vbBoundaries{ivbboundary} = vbSplines;     % replace each cell containing boundary coordinates with spline coordinates
                        end
                    end
                    vbBoundaries = vbBoundaries(~cellfun('isempty',vbBoundaries));      % delete empty cells from cell array
                    for ivbsplines = 1:length(vbBoundaries)     % run through each row of vascular bundle spline coordinates
                        eachVbSpline = vbBoundaries{ivbsplines};        % extract coordinates of each vascular bundle spline 
                        fprintf(fileID,'mdb.models[''Model-1''].sketches[''__profile__''].Spline(points=(\n\t'); % print Abaqus statement to sketch splines
                            for iprintspline = 1:1:size(eachVbSpline)      % run through each vascular bundle spline coordinates
                                formatSpec = '(%4.4f, %4.4f), ';        % set the amount of significant figures for spline coordinates
                                fprintf(fileID,formatSpec,eachVbSpline(iprintspline,1),eachVbSpline(iprintspline,2));       % print spline coordinates to file
                            end
                        fprintf(fileID,'\t))\n');
                    end
                    fprintf(fileID,'mdb.models[''Model-1''].sketches.changeKey(fromName=''__profile__'', toName=''vbspline'')\n');      % print Abaqus statement to save sketch


                case 'No'
                    
                    %% Save Output Image %% 

                    saveas(gcf,fullfile(outputFolderPath,outputImageFileName),'tif')        % save output image in output folder

                    %% Write Python File to Sketch Vascular Bundle Splines %%

                    vbBoundaries2 = cell(size(vascularBundleSplines));     % set boundary coordinates
                    splinename = 'VascularBundleSplines';     % suffix for filename
                    filename = sprintf('%s%s.py',inputImageFileName,splinename);       % set filename for vascular bundle splines
                    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
                    for ivbboundary = 1:length(vascularBundleSplines)        % run through each row of vascular bundle boundary coordinates
                        eachVbBoundary = vascularBundleSplines{ivbboundary};     % extract boundary coordinates of each vascular bundle
                        eachVbBoundary = unique(eachVbBoundary(:,1:2),'rows','stable');        % remove any repeated coordinates in boundary coordinates
                        eachVbBoundaryX = eachVbBoundary(:,2);        % x-coordinates of vascular bundle boundary
                        eachVbBoundaryY = eachVbBoundary(:,1);        % y-coordinates of vascular bundle boundary
                        vbSpline = boundary(eachVbBoundaryX,eachVbBoundaryY);       % extract splines from vascular bundle boundaries
                        vbSplineX = eachVbBoundaryX(vbSpline);        % x-coordinates of vascular bundle spline
                        vbSplineY = eachVbBoundaryY(vbSpline);        % y-coordinates of vascular bundle spline 
                        vbSplines = [vbSplineY*spatialCalibrationFactor,vbSplineX*spatialCalibrationFactor];       % convert pixel to mm
                        for ivbboundary2 = 1:1:length(vbBoundaries2)     % run through each row of vascular bundle boundary coordinates 
                            vbBoundaries2{ivbboundary} = vbSplines;     % replace each cell containing boundary coordinates with spline coordinates
                        end
                    end
                    vbBoundaries2 = vbBoundaries2(~cellfun('isempty',vbBoundaries2));      % delete empty cells from cell array
                    for ivbsplines = 1:length(vbBoundaries2)     % run through each row of vascular bundle spline coordinates
                        eachVbSpline = vbBoundaries2{ivbsplines};        % extract coordinates of each vascular bundle spline 
                        fprintf(fileID,'mdb.models[''Model-1''].sketches[''__profile__''].Spline(points=(\n\t'); % print Abaqus statement to sketch splines
                            for iprintspline = 1:1:size(eachVbSpline)      % run through each vascular bundle spline coordinates
                                formatSpec = '(%4.4f, %4.4f), ';        % set the amount of significant figures for spline coordinates
                                fprintf(fileID,formatSpec,eachVbSpline(iprintspline,1),eachVbSpline(iprintspline,2));       % print spline coordinates to file
                            end
                        fprintf(fileID,'\t))\n');
                    end
                    fprintf(fileID,'mdb.models[''Model-1''].sketches.changeKey(fromName=''__profile__'', toName=''vbspline'')\n');      % print Abaqus statement to save sketch
 
            end
            
            
%% Write Python File to Sketch Rind Spline %%

            rindSpline = rindSplinemmCoordinates;       % 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
            formatSpec = '(%4.4f, %4.4f)\n';        
            fprintf(fileID,formatSpec,rindSpline(1,1),rindSpline(1,2));     % re-print first spline coordinate to close the spline
            fprintf(fileID,'\t))\n');       
            fprintf(fileID,'mdb.models[''Model-1''].sketches.changeKey(fromName=''__profile__'', toName=''rindspline'')\n');        % print Abaqus statement to save sketch

%% Write Python File to Sketch Pith Spline %%

            pithSpline = pithSplinemmCoordinates;       % set spline coordinates
            splineName='PithSpline';      % suffix for filename
            filename = sprintf('%s%s.py',inputImageFileName,splineName);       % set filename for pith 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 ipithspline = 1:1:size(pithSpline,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,pithSpline(ipithspline,1),pithSpline(ipithspline,2));       % print spline coordinates to file
            end
            formatSpec = '(%4.4f, %4.4f)\n';
            fprintf(fileID,formatSpec,pithSpline(1,1),pithSpline(1,2));       % re-print first spline coordinate to close the spline
            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).majorDiameterToInsideOfGroove = 0;        % fill third column of new row with major diameter to groove middle
            structnewRow(i,1).minorDiameter = 0;        % fill fourth coulmn of new row with minor diameter
            structnewRow(i,1).majorDiameterToOutsideOfGroove = 0;      % fill fifth column of new row with major diameter to groove end
            structnewRow(i,1).numberOfVascularBundles = 0;        % fill sixth column of new row with total vascular bundle count
            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