%% CSVtoThicknessContour.m

% This code accepts a CSV file from LJ-Navigator software as input and
% assumes that the laser scan values are first reported in column 17 of the
% imported CSV file (standard starting point).

% This code requires several user inputs:
   % 1. Line 20 - File name
   % 2. Line 25 - Indicate platform used for scanning
   % 3. Line 30 - Scanning frequency
   
% This code relies on the user to identify the region of interest.

% The user will first be asked to click and drag to draw a rectangle around
% the sample. As long as the full sample is inside the rectangle, no 
% information will be lost. A tighter rectangle is preferred.

% The user will then be asked to identify several points that should have a
% thickness of zero. Click the points with your mouse and then, once done
% clicking ALL OF YOUR POINTS, press Enter.

%% Housekeeping

clear; clear all; clc; close all force;

%% Load Data

% Import information about the scan

delimiterIn = ','; %specify that input is a CSV file
headerlinesIn = 1; %specify that headers are on the first row and first column
m = importdata('ExampleCSVFile.csv',delimiterIn,headerlinesIn); %change name of csv file here
m=struct2cell(m); %convert the imported data from struct to cell
Numbers = m{1,1}; %extract the numbers matrix
Headings = m{2,1}; %extract the headings matrix (note that width is also included in the heading matrix because it is in the first row of the output)
WasGlassUsed = 1; % if the sample was placed on glass, this should be 1; otherwise, 0

% Enter the sampling frequency

fs = 20; % [Hz] change this depending on the sampling frequency that the data was collected at
delt_t = 1/fs; % [s] calculate the time step

% Determine the size of the numbers and headings matrices

Begin = 17; % this should stay the same, but may change based on other outputs you've requested
Thickness = Numbers(:, Begin:1:end);

%% Plot Raw Data

% This section produces a plot of the raw, unfiltered, unedited data

figure; % call a new figure
surf(Thickness); % use surf to plot the thickness contour
view(0, 90); % reposition the view
shading interp; % remove lines from surface plot
title('Click and drag to draw a rectangle around your sample.'); % title
Box = drawrectangle; % requesting user input to define the region of interest

%% Interpreting Rectangle Input

% This section produces a smaller region of interest

XRange = round([Box.Vertices(1, 1), Box.Vertices(3, 1)]); % determining the columns contained by the rectangle
YRange = round([Box.Vertices(1, 2), Box.Vertices(3, 2)]); % determing the rows contained by the rectangle
XRange = sort(XRange); % sorting the columns in ascending order
YRange = sort(YRange); % sorting the rows in ascending order

Thickness = Thickness(YRange(1):1:YRange(2), XRange(1):1:XRange(2)); % extracting the region of interest

%% Adjusting for Laser Tilt (If Any)

% This section adjusts the thickness contour according to the slope from
% right to left (if there is any) of the first and last rows of the
% Thickness contour. These two rows should have zero thickness from side
% to side.

FirstRow = Thickness(1, :); % pulling the first row, which should be flat ideally
LastRow = Thickness(end, :); % pulling the last row, which should also be flat ideally
Slope(1) = FirstRow(end) - FirstRow(1); % determining the slope of the first row
Slope(2) = LastRow(end) - LastRow(1); % determining the slope of the second row
Slope = mean(Slope); % calculating the average slope
Adjustment = linspace(0, Slope, size(Thickness, 2)); % setting up an adjustment matrix based on the slope of the laser head

for i = 1:1:size(Thickness, 1) % for every row in the Thickness contour...
   Thickness(i, :) = Thickness(i, :) - Adjustment; % add the adjustment matrix to it
end

%% Finding What Zero Should Be

% This section uses user input to determine what "zero thickness" currently
% is and adjusts the entire matrix based on this offset

figure; % calling a new figure
surf(Thickness); % plotting the Thickness contour using the surf function
view(0, 90); % repositioning the view
title({'Please select three locations that should be zero.', 'Then press Enter.'}); % title
shading interp; % removing black lines from the surface plot

[ZeroX, ZeroY] = ginput; % requesting user input about locations with zero thickness
ZeroX = round(ZeroX); % rounding points to be indices
ZeroY = round(ZeroY); % rounding points to be indices
ZeroThickness(1:1:size(ZeroX, 1), 1:1:1) = zeros; % initialization

for i = 1:1:size(ZeroX, 1) % for every zero thickness point entered by the user...
   ZeroThickness(i, :) = Thickness(ZeroY(i), ZeroX(i)); % pull its associated raw thickness
end

AverageZero = mean(ZeroThickness); % find the average thickness of the user-selected zero points

Thickness = Thickness + abs(AverageZero); % adjust the Thickness contour based on the user-selected zero points

%% Removing Negative Spikes

% This section removes large negative spikes from the data

for j = 1:1:size(Thickness, 2) % for every column in the Thickness contour...
   for i = 1:1:size(Thickness, 1) % and for every row in the Thickness contour...
      if Thickness(i, j) < 0 % if the current entry is less than zero
         Thickness(i, j) = 0; % make it zero
      end
   end
end

%% Accounting for Increased Thickness From Scanning Stage (If Necessary)

% This section removes "extra thickness" from the sample, which occurs if
% the sample was scanned while resting on a glass insert. Our glass insert
% is 2.45 mm, so this value will need to be edited if this is not the case
% for your current scan.

if WasGlassUsed == 1 % if the sample was scanned using a glass platform...
   for j = 1:1:size(Thickness, 2) % for every column in the Thickness contour...
      for i = 1:1:size(Thickness, 1) % and for every row in the Thickness contour...
         if Thickness(i, j) >= 1 % if the thickness value belongs to the sample...
            Thickness(i, j) = Thickness(i, j) - 2.45; % subtract the thickness of the glass platform
         end
      end
   end
end 

%% Removing Positive Spikes

% This section removes large positive spikes from the Thickness contour and
% assumes that the sample is no thicker than 1.3 mm at any location. For
% this TissueMend sample, this is true; however, this will not always be
% the case for every soft tissue sample and this value will need to be
% adjsuted accordingly.

for j = 1:1:size(Thickness, 2) % for every column in the Thickness contour...
   for i = 1:1:size(Thickness, 1) % and for every row in the Thickness contour...
      if Thickness(i, j) >= 1.3 % if the thickness value is greater than the expected thickness of the sample....
         Thickness(i, j) = 0; % make it zero
      end
   end 
end

%% Zeroing Leftover Spikes

% This section removes any remaining negative spikes

for j = 1:1:size(Thickness, 2) % for every column in the Thickness contour...
   for i = 1:1:size(Thickness, 1) % and for every row in the Thickness contour...
      if Thickness(i ,j) < 0 % if there is still a negative thickness value...
         Thickness(i, j) = 0; % make it zero!
      end
   end 
end

%% Filtering and Plotting

% This section first filters the Thickness contour using a 2D Gaussian
% filter, then plots the filtered Thickness contour for visualization.

Filter = imgaussfilt(Thickness); % apply a 2D Gaussian filter to the sample

figure; % calling a new figure
surf(Filter); % plotting the Thickness contour using the surf function
c = colorbar ('southoutside'); %calling a colorbar and changing colobar location
c.Ruler.TickLabelFormat='%.1f mm'
%title('Thickness (mm)'); % adding a title
view(0, 90); % repositioning the view
shading interp; % removing black lines from the surface plot
set(gca, 'FontName', 'Times', 'FontWeight', 'bold', 'FontSize', 14); % formatting
axis off; % formatting
caxis([0 1.5]); % formatting



