% SedCT MATLAB GUI for sediment CT scans collected on a medical CT scanner
% Developed 2014, Brendan T. Reilly, Oregon State University
% Version 1.08, last modified July 21, 2022
% Please cite:
% Reilly, B.T., Stoner, J.S., Wiest, J., (2017) 
% SedCT: MATLAB™ tools for standardized and quantitative 
% processing of sediment core computed tomography (CT) data 
% collected using a medical CT scanner. Geochem. Geophys. Geosyst.,
% doi:10.1002/2017GC006884


function varargout = SedCT(varargin)
% SedCT MATLAB code for SedCT.fig
%      SedCT, by itself, creates a new SedCT or raises the existing
%      singleton*.
%

% Edit the above text to modify the response to help SedCT

% Last Modified by GUIDE v2.5 21-Jul-2022 13:45:32

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @SedCT_OpeningFcn, ...
                   'gui_OutputFcn',  @SedCT_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% Opening Function
function SedCT_OpeningFcn(hObject, eventdata, handles, varargin)
% Have SedCT open and display Benny the Beaver image.  Go Beavers!
RadioPlotNOPE(1, 2)
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);
% UIWAIT makes SedCT wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = SedCT_OutputFcn(hObject, eventdata, handles) 
varargout{1} = handles.output;

%-----------------------------
%GUI Functions
%-----------------------------
%Inputs Left Panel

%Input Folder------------------------------------
function Input_folder_Callback(hObject, eventdata, handles)
handles.metricdata.InputDICOM = 'No Path Selected';
InputDICOM = uigetdir('', 'Pick DICOM Directory');
if isnan(InputDICOM)
    set(hObject, 'String', NaN);
    errordlg('Select an input DICOM folder','Error');
end

handles.metricdata.InputDICOM = InputDICOM;
set(handles.path_display, 'String', handles.metricdata.InputDICOM);
guidata(hObject,handles)

%Coronal vs. Axial ----------------------
function coronal_radio_Callback(hObject, eventdata, handles)
set(handles.axial_radio, 'Value', 0);
set(handles.coronal_radio, 'Value', 1);

function axial_radio_Callback(hObject, eventdata, handles)
set(handles.axial_radio, 'Value', 1);
set(handles.coronal_radio, 'Value', 0);
set(handles.straighten, 'Value', 0);

%Straighten Check Box ------------------------------
function straighten_Callback(hObject, eventdata, handles)
if get(handles.axial_radio, 'Value') == 1;
    set(handles.straighten, 'Value', 0);
end

%Enter number of Pixels to sample -------------------
function pixel_sample_entry_Callback(hObject, eventdata, handles)
sa = str2double(get(hObject, 'String'));
if isnan(sa)
    set(hObject, 'String', 5000);
    errordlg('Input must be a number','Error');
end
%  File the new volume value
handles.metricdata.sa = sa;
guidata(hObject,handles)

function pixel_sample_entry_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


function ImageLocation_Callback(hObject, eventdata, handles)
imlocation = str2double(get(hObject, 'String'));
if isnan(imlocation)
    set(hObject, 'String', 50);
    errordlg('Input must be a number','Error');
end
%  File the new volume value
handles.metricdata.imlocation = imlocation;
guidata(hObject,handles)

function ImageLocation_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Inputs---------------------------------------
%Background value -- Added to deal with DICOM files generated on CT
%scanners that don't use HU values.  Default value is 0, which is water in HU units.
function background_value_Callback(hObject, eventdata, handles)
bv = str2double(get(hObject, 'String'));
if isnan(bv)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end
%  File the new volume value
handles.metricdata.bv = bv;
guidata(hObject,handles)

function background_value_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Max Slope value - Core specific, but helps isolate the values that are the
%most representative of sediment.  I find in nicely cores organic rich lake
%cores, this is better as a low value, but in heavily disturbed, high
%density sediments, sometime sthis values needs to be pretty high (5+).
function slope_entry_Callback(hObject, eventdata, handles)
MaxSlope = str2double(get(hObject, 'String'));
if isnan(MaxSlope)
    set(hObject, 'String', 1);
    errordlg('Input must be a number','Error');
end
%  File the new volume value
handles.metricdata.MaxSlope = MaxSlope;
guidata(hObject,handles)

function slope_entry_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
% function Input_slope_Callback(hObject, eventdata, handles)
% 
% MaxSlope = str2double(get(hObject, 'String'));
% if isnan(MaxSlope)
%     set(hObject, 'String', 1);
%     errordlg('Input must be a number','Error');
% end
% 
% % File the new volume value
% handles.metricdata.MaxSlope = MaxSlope;
% guidata(hObject,handles)
% function Input_slope_CreateFcn(hObject, eventdata, handles)
% 
% if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
%     set(hObject,'BackgroundColor','white');
% end

% Trim value parameter -- removes from the edges of the distribution.
% Again subjective based on the kind of sediment.
function Input_trim_Callback(hObject, eventdata, handles)

Trim = str2double(get(hObject, 'String'));
if isnan(Trim)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

% File the new volume value
handles.metricdata.Trim = Trim;
guidata(hObject,handles)
function Input_trim_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% Min Pixel Value -- Usually this is fine at zero, but sometimes not.  Gets
% rid of values if SedCT only isolates a few pixels it thinks is
% representative.
function MinPix_Callback(hObject, eventdata, handles)

Minpx = str2double(get(hObject, 'String'));
if isnan(Minpx)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

handles.metricdata.Minpx = Minpx;
guidata(hObject,handles)
function MinPix_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function TMask_Callback(hObject, eventdata, handles)


%The Mask parameters - most CT scans include volume above and below that
%are not related to the sediment.  This gets rid of those values.

TopMask = str2double(get(hObject, 'String'));
if isnan(TopMask)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

handles.metricdata.TopMask = TopMask;
guidata(hObject,handles)
function TMask_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
function BMask_Callback(hObject, eventdata, handles)

BotMask = str2double(get(hObject, 'String'));
if isnan(BotMask)
    set(hObject, 'String', 0);
    errordlg('Input must be a number','Error');
end

handles.metricdata.BotMask = BotMask;
guidata(hObject,handles)
function BMask_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%Load Button-----------------------------------------
%Load in the DICOM data and put it in a format that SedCT can use.
function Load_data_Callback(hObject, eventdata, handles)

RadioPlotNOPE(1, 2)

% Determine which section is to be loaded (up to four sections can be used
% at once).
if get(handles.Sec1, 'Value') == 1;
    ssc = 1;
else if get(handles.Sec2, 'Value') == 1;
        ssc = 2;
    else if get(handles.Sec3, 'Value') == 1;
            ssc = 3;
        else if get(handles.Sec4, 'Value') == 1;
                ssc = 4;
            else errordlg('Select a Section Number to Process','Error');
            end
        end
    end
end

% Determine if the CT scans are in axial or coronal format.  Log this value
% as orient.
if get(handles.coronal_radio, 'Value') == 1
    orient = 0;
else orient = 1;
end


%Parameters
%Defaults

% Use 5000 as number of pixels to sample, unless the user has entered
% another amount
sa = 5000;

try
    sa = handles.metricdata.sa; % number of pixels to sample
end
% try
    inputlocation = handles.metricdata.InputDICOM;
% catch
%     errordlg('Select an input folder','Error');
% end

% Sample central slice, unless user chooses another horizon
imlocationz = 50;
try
    imlocationz = handles.metricdata.imlocation; % percentage of way through CT scan to use slice image
end
imlocation = imlocationz/100;
if imlocation > 1;
    imlocation = 1;
end
if imlocation <= 0;
    imlocation = .01;
end


%Bring in DICOM file, make double precision, and extract information

% addpath(inputlocation);
l = dir(inputlocation);  % generates list of files/folders in the selected directory
d = l(3:end); % gets rid of the paths
dn = length(d); % number of files in the folder
dm = ceil(dn * imlocation); % use for whole rounds, this is the CT image slice
if isempty(dm)
    errordlg('Make sure to select a directory with DICOM files.  (Or make sure that the MATLAB image processing toolbox is installed).','Error')
    return
end

% Load in Dicom Files ---------------------
%Load in Dicom info
h = waitbar(0,'Loading DICOM File'); % use a waitbar, because this can take a while on slow computers
ord = nan(dn, 1); % preallocate the ord variable -- this variable is used to make sure DICOMS are imported in the correct order, regardless of file name.
check_w = ord;
check_h = ord;
count = 0;
for n = 1:dn % do this for every file in the directory
    % try function is used in case the user has non-dicom files in the
    % folder.  Prevents errors when using the dicominfo function
    try 
        info = dicominfo([inputlocation '/' d(n).name]); % requires image toolbox
        count = count + 1;
        if isfield(info, 'ImageComments')
            if strcmp(info.ImageComments, 'Coronal.Ref') % Check for these pesky Coronal Reference Files
                count = count-1;
            end
        end
        if ~isfield(info, 'SliceThickness') % Add a catch if the metadata is missing Slice Thickness info
            info.SliceThickness = 0.5;
        end
    if dn == 1 | count == 2
        ID = []; % Do this so that scanners with slightly different metadata don't cause errors
        try
            Name1 = info.PatientID; % At OSU, core info is usually stored in these three variables.
            ID = [ID Name1]; % This will be the ID used from here onwards
        end
        try
            Name = info.PatientName.FamilyName;
            ID = [ID ' ' Name]; % This will be the ID used from here onwards
        end
        try
            Section = info.ImageComments;
            ID = [ID ' ' Section]; % This will be the ID used from here onwards
        end
        if isempty(ID)
            ID = 'SedCT could not read metadata';
        end
        
        try
            rsslope = info.RescaleSlope;
            rsintercept = info.RescaleIntercept;
        catch
            rsslope = 1;
            rsintercept = 0;
        end
        
        try
            pxl2 = info.PixelSpacing(1)/10; % And this defines the size of the pixels in two dimensions (This will be different for Axial vs. Coronal, hence the if statement
                if orient == 1
                    pxl1 = info.SliceThickness/10;
                    pxl2 = pxl2;
                else pxl1 = pxl2;
                end
        catch
             errordlg('DICOM depth metadata cannot be read.  Please contact btreilly@ucsd.edu','Error')
             close (h)
             return
        end
        
        
        
    end
    
    if isfield(info, 'ImageComments')
        if ~strcmp(info.ImageComments, 'Coronal.Ref') %Only count if not a Coronal Refernce File
            ord(n) = info.InstanceNumber; % Instance number is used to make sure that DICOM files are loaded in the correct order.
            check_w(n) = info.Width;
            check_h(n) = info.Height;
        end
    else
        try
            ord(n) = info.InstanceNumber; % Instance number is used to make sure that DICOM files are loaded in the correct order.
        catch
            ord(n) = n;
        end
        check_w(n) = info.Width;
        check_h(n) = info.Height;
    end
    end
    waitbar(n/dn,h); 
end
if isempty(find(~isnan(ord) == 1))
    errordlg('Make sure to select a directory with DICOM files.  (Or make sure that the MATLAB cessing toolbox is installed).','Error')
    close(h)
    return
end

% Make sure files are in correct order
i = check_w ~= mode(check_w) | check_h ~= mode(check_h);
ord(i) = NaN;
i = ~isnan(ord); % if ord is a NaN, it is because there was a non-DICOM file in the directory
d = d(i); % resorting the directory
[~, i] = sort(ord(i)); % resorting the directory
d = d(i); % resorting the directory
dn = length(d); % Number of DICOM files in the directory
dm = ceil(dn  * imlocation); % Central Slice
% dm = ceil(dn * 1/4); %use for buttons! (you can make your image whatever slice you want)
if isempty(dm) | dm == 0;
    errordlg('Make sure to select a directory with DICOM files.  (Or make sure that the MATLAB image processing toolbox is installed).','Error')
        close(h)
    return
end

%Load in Dicom Data
if orient == 0 %for coronal/sagittal data
    ctimage = double(dicomread([inputlocation '/' d(dm).name]))*rsslope+rsintercept; %reference image from the dm or central slice & rescale based on metadata
    [l1, l2] = size(ctimage); % getting slice size, typically 1624 x 1624 for OSU, but different for other scanners
    i = ctimage < -500; % if the HU value is less than -500, lets get rid of it.  This is way less than the density of water.  Adios!
    keep = sum(i) ~= l2;  % defines the rows we want to keep (i.e. if all the HU values are less than -500, get rid of it to reduce data size).
    if isempty(find(keep==1))
        keep = ones(size(keep)); % But don't throw away all the data if nothing is there...
        errordlg('Looks like your image does not contain any mud.  Try adjusting the location of your image slice.','Error')
    end
    ctimage = ctimage(:, keep); %throw away no data regions
    % If straightening is selected
    if get(handles.straighten, 'Value') == 1
        dx = diff(ctimage')'; % all this does is try to fins where the HU values increase the most -- the hope is it is the transition between air and core liner
        temp = dx > 100;
        [~, mi] = max(temp');
        use = [ceil(.25*l1):1:l1-ceil(.25*l1)];
        l = fitlm(use, mi(use)); % Figure out the slope of the first high rates of change relative to the DICOM file
        errordlg('You need the statistics toolbox to use the "Vertically Straighten" function.  If you are missing this tool box, just unclick the checkmark box and you should be able to proceed without a problem.','Error');
        s = l.Coefficients.Estimate(2);
        rot = rad2deg(atan(s));
        xr = imrotate(ctimage, rot/2, 'nearest', 'crop'); %rotate the image to straighten
        ctimage = xr;
%         i = ctimage < -500;
%         keep = sum(i) ~= l1;
%         ctimage = ctimage(:, keep) %throw away no data regions
    end
    dic = zeros(l1, length(find(keep == 1)), dn); %preallocate matrix to load in files
    for n = 1:dn % for each file
        temp = double(dicomread([inputlocation '/' d(n).name]))*rsslope+rsintercept; % And make it double and rescale
        temp = temp(:, keep); % again only keeping what we care about
        if get(handles.straighten, 'Value') == 1
            temp = imrotate(temp, -rot, 'nearest', 'crop'); % roatate as above
        end
        dic(:, :, n) = temp; %store data in matrix
        waitbar(n/dn,h); 
    end
    a = size(dic, 2) * size(dic, 3); % how many pixels are there?
    if a < sa %make sure you have enough pixels to meet the use defined number
        sa = a; % if not? Just use what you have.
    end
    y = zeros(sa, dn); %preallocate the matrix that will hold the sampled pixels
    for n = 1:l1 % do it for the number of pixel rows in the DICOM files
        i = randperm(a, sa); %index for random selection
        y(:, n) = dic(n, i); %data to keep
    end
    clear dic % this can be a lot of data, so lets get rid of it
    
    %Make Cross Section for Choosing Slice Image
    imlocSlice = zeros(dn, l2);
    keep = ceil(l1/2); %keep center line only for cross section;
    for n = 1:dn
        temp = double(dicomread([inputlocation '/' d(n).name]))*rsslope+rsintercept; %read data
        imlocSlice(n, :) = temp(keep, :);
    end
    temp = zeros(size(imlocSlice));
    i = imlocSlice > 0;
    temp(i) = 1;
    keep = sum(temp)>0;
    imlocSlice = imlocSlice(:, keep);
    
    %Axial Slices
else [l1, l2] = size(dicomread([inputlocation '/' d(dm).name])); %for axial data - which is a little easier to work with
    a = l1*l2;
    keep = ceil(l2 * imlocation); %keep center line only
    ctimage = zeros(dn, l1); % preallocate
    if a < sa %make sure you have enough pixels
        sa = a;
    end
    y = zeros(sa, dn); %preallocate
    for n = 1:dn
        temp = double(dicomread([inputlocation '/' d(n).name]))*rsslope+rsintercept; %read data
        ctimage(n, :) = temp(keep, :); %make coronal image
        i = randperm(a, sa); %index for random selection
        y(:, n) = temp(i); %keep only part of the data
        if n == ceil(dn/2)
            imlocSlice = temp;
        end
        waitbar(n/dn,h); 
    end
end

close(h)
subplot(3, 3, 5)
hold on
image(imlocSlice)
if orient == 0
    daspect([info.SliceThickness/10 pxl1 1])
else
    daspect([pxl1 pxl1 1])
end
colormap(gray)
l2p = length(imlocSlice(:, 1));
ti = [0:l2p/10:l2p];
plot([1 length(imlocSlice(1, :))], imlocation * [l2p l2p], 'r', 'linewidth', 3)
for n = 1:length(ti);
    plot([1 length(imlocSlice(1, :))], ti(n)*[1 1], 'r--')
end
title('Horizon for Slice Image')
ylabel('Image Location (%)')
set(gca, 'xtick', [], 'ytick', ti, 'yticklabels', [0:10:100], 'ydir', 'reverse')
ylim([0 l2p])
xlim([0 length(imlocSlice(1, :))])



set(handles.data_load_display, 'String', ID);
%save data for each section to retreive later
handles.metricdata.section(ssc).name = ID;
handles.metricdata.section(ssc).pixelsize = [pxl1 pxl2];
handles.metricdata.section(ssc).imagein = ctimage;
handles.metricdata.section(ssc).rawdata = y;
guidata(hObject, handles);


%Process CT Data -----------------------------------
function process_data_Callback(hObject, eventdata, handles)

% Check to see what section we are working with
if get(handles.Sec1, 'Value') == 1;
    ssc = 1;
else if get(handles.Sec2, 'Value') == 1;
        ssc = 2;
    else if get(handles.Sec3, 'Value') == 1;
            ssc = 3;
        else if get(handles.Sec4, 'Value') == 1;
                ssc = 4;
            else errordlg('Select a Section Number to Process','Error');
            end
        end
    end
end
if get(handles.coronal_radio, 'Value') == 1
    orient = 0;
else orient = 1;
end


%Parameters
%Defaults (these are described above, in the user guise, and in the
%publication)
MaxSlope = 1;
TrimRange = 500;
manpx = 0;
BotMask = 0;
TopMask = 0;
backv = 0;

% See if the user has enetered values other than the defaults
try
    MaxSlope = handles.metricdata.MaxSlope; %maximum rate of change allowed
end
try
    TrimRange = handles.metricdata.Trim; % px trimmed from either end
end
try
    manpx = handles.metricdata.Minpx; %Number of pixels allowed
end
try
    TopMask = handles.metricdata.TopMask;
end
try
    BotMask = handles.metricdata.BotMask;
end
try
    backv = handles.metricdata.bv;
end

% Retreive the data logged while importing the DICOM files
try
    ID = handles.metricdata.section(ssc).name;
    pxl1 = handles.metricdata.section(ssc).pixelsize(1);
    pxl2 = handles.metricdata.section(ssc).pixelsize(2);
    ctimage = handles.metricdata.section(ssc).imagein;
    y = handles.metricdata.section(ssc).rawdata;
catch
    errordlg('Load a DICOM Directory','Error');
end


pxcut = 10; %Standard Deviations from the mean for min px count

h = waitbar(0,'Processing CT Data');
z = sort(y); % sort the randomly sampled pixels
i = z <= backv; % get rid of anything that is above the background
% if isempty(find(i==1))
%     errordlg('No Background Detected.  Increase Background Value.','Error');
% end
z(i) = NaN;
i = sum(z, 2, 'omitnan') > 0;
z = z(i, :);
% Now we isolate the longest continuous selection of sorted pixels, below a
% user defined slope
dz = diff(z);
z = z(1:end-1, :);
i = dz > MaxSlope;
z(i) = NaN;
[l1, l2] = size(z);
z = [nan(1, l2); z];
% Replace the unrepresentative values with NaNs
for n = 1:l2
    F = find(isnan(z(:, n)));
    D = diff(F);
    [M, L] = max(D);
%     i1n(n) = sum(D(1:L-1))+2;
%     i2n(n) = sum(D(1:L-1))+2+M;
    i1 = sum(D(1:L-1))+2 + TrimRange;
    i2 = sum(D(1:L-1))+2+M - TrimRange;
    if i2 < 1;
        i2 = 1;
    end
    z(1:i1, n) = NaN;
    z(i2:end, n) = NaN;
    waitbar(n/l2,h);
end
% i = find(sum(z, 'omitnan') > 0);
% startcore = i(1);
startcore = 1; % startcore for now is 1, to keep things consistent
z = sort(z);
i = sum(z, 2, 'omitnan') ~= 0;
z = z(i, :);  % shed some of the NaN values
[ww, w] = size(z);

%Calcualte mean and standard deviation of the isolated values
gs = mean(z, 'omitnan');
sd = std(z, 'omitnan');
%count pixels
px = sum(~isnan(z));
meanpx = mean(px, 'omitnan');
maxpx = max(px);
stdpx = std(px, 'omitnan');
i = px < meanpx - pxcut * stdpx;
gs(i) = NaN;
sd(i) = NaN;
i = find(px < manpx);
gs(i) = NaN;
sd(i) = NaN;

midgs = mean(gs, 'omitnan');
stdgs = std(gs, 'omitnan');
i = gs <= 0;
gs(i) = NaN;
sd(i) = NaN;

%Dimensions on depth scale
depth = ((pxl1 * w)/w:(pxl1 * w)/w :pxl1 *w) - startcore*(pxl1);

%Mask
upmask = depth < TopMask;
gs(upmask) = NaN;

downmask = find(depth - BotMask + startcore*(.025) <= 0);
gs(end-downmask+1) = NaN;


%interpolate to common depthscale
t = timeseries([gs', sd', px'], depth);
t1 = resample(t, min(depth):.025:max(depth));
depthi = t1.time;
gsi = t1.data(:, 1);
sdi = t1.data(:, 2);
pxi = t1.data(:, 3);
imr = imresize(ctimage, [length(depthi) pxl2*40*size(ctimage, 2)]);
if isempty(z)
    errordlg('No solution.  Try changing paramters, such as decreasing the trim value.  Check that data was loaded correctly as axial or sagittal','Error');
    close(h)
    return
end
z1 = imresize(z, [size(z, 1), length(depthi)]);

close(h)

%Visualize plot
subplot(1, 11, 1:2)
top = (midgs + 3*stdgs);
imshow(imr)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
set(gca, 'xtick', [])

subplot(1, 11, 3:5)
plot(gsi, depthi, 'k')
hold on
set(gca,'YDir','reverse');
ylim([min(depthi), max(depthi)])
xlabel('CT #')
xlim([min(gsi) max(gsi)])
ylabel('Depth (cm)')
hold off

subplot(1, 11, 6:7)
plot(sdi, depthi, 'k')
set(gca,'YDir','reverse');
ylim([min(depthi), max(depthi)])
xlabel('Standard Deviation')
title(ID, 'Interpreter', 'none')
xlim([min(sdi) max(sdi)])
set(gca, 'YTicklabel', [])


subplot(1, 11, (8:10))
mesh([1:1:ww], depth, z')
view(0, 90)
set(gca,'YDir','reverse');
ylim([min(depth), max(depth)])
xlim([0, max(px)])
colorbar
xlabel('Values Used')
if (midgs - 3*stdgs) > 0
    caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
set(gca, 'YTicklabel', [])


%Create output
output = [depth' gs' sd' px'];
outputi = [depthi gsi sdi pxi];

handles.metricdata.section(ssc).DICOM = ctimage;
handles.metricdata.section(ssc).pixel = pxl1;
handles.metricdata.section(ssc).resultsnointerp = output;
handles.metricdata.section(ssc).results = outputi;
handles.metricdata.section(ssc).meshplot = z';
handles.metricdata.section(ssc).CTresize = z1';
handles.metricdata.section(ssc).name = ID;
handles.metricdata.section(ssc).startcore = startcore;
handles.metricdata.section(ssc).CT = imr;
guidata(hObject, handles);


clear l w midgs stdgs output outputlocation pixellength px pxcut s savelocation sd startcore stdpx w y
clear ID MaxSlope TrimRange depth di dm filename gs z

%----------------------------
%Right Panel
%Radiobutton Functions
function RadioSelect(hObject, eventdata, handles)
set(handles.Sec1, 'Value', get(handles.Sec1, 'Min'));
set(handles.Sec2, 'Value', get(handles.Sec2, 'Min'));
set(handles.Sec3, 'Value', get(handles.Sec3, 'Min'));
set(handles.Sec4, 'Value', get(handles.Sec4, 'Min'));
set(handles.st12, 'Value', get(handles.st12, 'Min'));
set(handles.st23, 'Value', get(handles.st23, 'Min'));
set(handles.st34, 'Value', get(handles.st34, 'Min'));
set(handles.comALL, 'Value', get(handles.comALL, 'Min'));

function RadioPlot(CTscan, ResultsInput, MeshPlot, Name)
midgs = mean(ResultsInput(:, 2), 'omitnan');
stdgs = std(ResultsInput(:, 2), 'omitnan');
ID = Name;
depthi = ResultsInput(:, 1);

subplot(1, 11, 1:2)
imshow(CTscan)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
set(gca, 'xtick', [])

subplot(1, 11, 3:5)
plot(ResultsInput(:,2), ResultsInput(:,1), 'k')
hold on
set(gca,'YDir','reverse');
ylim([min(depthi), max(depthi)])
xlabel('CT #')
xlim([min(ResultsInput(:,2)) max(ResultsInput(:,2))])
ylabel('Depth (cm)')
hold off

subplot(1, 11, 6:7)
plot(ResultsInput(:,3), ResultsInput(:,1), 'k')
set(gca,'YDir','reverse');
ylim([min(depthi), max(depthi)])
xlabel('Standard Deviation')
title(ID, 'Interpreter', 'none')
xlim([min(ResultsInput(:,3)) max(ResultsInput(:,3))])
set(gca, 'YTicklabel', [])


subplot(1, 11, (8:10))
mesh([1:1:size(MeshPlot, 2)], depthi, MeshPlot)
view(0, 90)
set(gca,'YDir','reverse');
ylim([min(depthi), max(depthi)])
xlim([0, max(ResultsInput(:, 4))])
colorbar
xlabel('Values Used')
if (midgs - 3*stdgs) > 0
    caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
set(gca, 'YTicklabel', [])


function RadioStitchPlot(CTscan1, Data1, CTscan2, Data2, Stitch, StitchCT, sover)
midgs = mean(Stitch(:,2), 'omitnan');
stdgs = std(Stitch(:,2), 'omitnan');

subplot(1, 11, 1)
imshow(CTscan1)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
title('Top')

subplot(1, 11, 2)
imshow(CTscan2)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
title('Bottom')

subplot(1, 11, 8)
imshow(StitchCT)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
title('Stitched')

subplot(1, 11, 3:5)
imshow(StitchCT)
colormap('gray')
if (midgs - 3*stdgs) > 0
   caxis([(midgs - 3*stdgs) (midgs + 3*stdgs)])
else caxis([0 (midgs + 3*stdgs)])   
end
title('Stitch Section')
ylim([(sover-5)*40 (max(Data1(:,1))+5)*40])

subplot(1, 11, 6:7)
plot(Data1(:,2), Data1(:,1), 'r')
hold on
plot(Data2(:,2), Data2(:,1), 'b')
ylim([(sover-5) (max(Data1(:,1))+5)])
set(gca,'YDir','reverse');
title('Stitch Section')
xlabel('CT #')
hold off

subplot(1, 11, 9:10)
plot(Data1(:,2), Data1(:,1), 'r', 'LineWidth', 2)
hold on
plot(Data2(:,2), Data2(:,1), 'b', 'LineWidth', 2)
plot(Stitch(:,2), Stitch(:,1), 'k')
ylim([min(Stitch(:,1)), max(Stitch(:,1))])
xlim([min(Stitch(:, 2)), max(Stitch(:, 2))])
set(gca,'YDir','reverse');
title('Complete Stitch')
xlabel('CT #')
hold off

function RadioPlotNOPE(input1, input2)
subplot(1, 11, 1:10)
plot(input1, input2);
subplot(1, 11, 3:9);
x = imread('blank_screen1.tif');
imshow(x)

%Radiobuttons
function Sec1_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.Sec1, 'Value', 1)
try RadioPlot(handles.metricdata.section(1).CT,... 
    handles.metricdata.section(1).results,...
    handles.metricdata.section(1).CTresize,...
    handles.metricdata.section(1).name)
catch
    RadioPlotNOPE(1,2)
end
function Sec2_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.Sec2, 'Value', 1)
try RadioPlot(handles.metricdata.section(2).CT,... 
    handles.metricdata.section(2).results,...
    handles.metricdata.section(2).CTresize,...
    handles.metricdata.section(2).name)
catch
    RadioPlotNOPE(1,2)
end
function Sec3_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.Sec3, 'Value', 1)
try RadioPlot(handles.metricdata.section(3).CT,... 
    handles.metricdata.section(3).results,...
    handles.metricdata.section(3).CTresize,...
    handles.metricdata.section(3).name)
catch
    RadioPlotNOPE(1,2)
end
function Sec4_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.Sec4, 'Value', 1)
try RadioPlot(handles.metricdata.section(4).CT,... 
    handles.metricdata.section(4).results,...
    handles.metricdata.section(4).CTresize,...
    handles.metricdata.section(4).name)
catch
    RadioPlotNOPE(1,2)
end

function st12_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.st12, 'Value', 1)
try RadioStitchPlot(handles.metricdata.section(1).CT,...
    handles.metricdata.stitch(1).top,...
    handles.metricdata.section(2).CT,...
    handles.metricdata.stitch(1).bot,...
    handles.metricdata.stitch(1).stitch,...
    handles.metricdata.stitch(1).stitchCT,...
    handles.metricdata.stitch(1).startoverlap)
catch
    RadioPlotNOPE(1,2)
end

function st23_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.st23, 'Value', 1)
try RadioStitchPlot(handles.metricdata.section(2).CT,...
    handles.metricdata.stitch(2).top,...
    handles.metricdata.section(3).CT,...
    handles.metricdata.stitch(2).bot,...
    handles.metricdata.stitch(2).stitch,...
    handles.metricdata.stitch(2).stitchCT,...
    handles.metricdata.stitch(2).startoverlap)
catch
    RadioPlotNOPE(1,2)
end

function st34_Callback(hObject, eventdata, handles)
RadioSelect(hObject, eventdata, handles)
set(handles.st34, 'Value', 1)
try RadioStitchPlot(handles.metricdata.section(3).CT,...
    handles.metricdata.stitch(3).top,...
    handles.metricdata.section(4).CT,...
    handles.metricdata.stitch(3).bot,...
    handles.metricdata.stitch(3).stitch,...
    handles.metricdata.stitch(3).stitchCT,...
    handles.metricdata.stitch(3).startoverlap)
catch
    RadioPlotNOPE(1,2)
end

function comALL_Callback(hObject, eventdata, handles)
RadioPlotNOPE(1, 2)

RadioSelect(hObject, eventdata, handles)
set(handles.comALL, 'Value', 1)

if get(handles.secch1, 'Value') == 0;
    errordlg('Use check boxes to select sections, starting with section 1.','Error');
end

if get(handles.secch1, 'Value') == 1;
    num = 1;
end
if num == 1 & get(handles.secch2, 'Value') == 1;
    num = 2;
end
if num == 2 & get(handles.secch3, 'Value') == 1;
    num = 3;
end
if num == 3 & get(handles.secch4, 'Value') == 1;
    num = 4;
end

if num == 1
    CTout = handles.metricdata.section(1).CT;
    RESout = handles.metricdata.section(1).results;
    MESHout = handles.metricdata.section(1).CTresize;
    IDout = handles.metricdata.section(1).name;
end

if num >= 2
    CTout = handles.metricdata.stitch(1).stitchCT;
    RESout = handles.metricdata.stitch(1).stitch;
    IDout = handles.metricdata.section(1).name;
    
    MESH1 = handles.metricdata.section(1).CTresize;
    MESH2 = handles.metricdata.section(2).CTresize;
    overlap = handles.metricdata.stitch(1).overlap;
    
    i = sum(CTout, 'omitnan') == 0;
    CTout = CTout(:, ~i);
    
    [mtl, mtw] = size(MESH1);
    [mbl, mbw] = size(MESH2);
    
    MESH1a = zeros(mtl, max([mtw, mbw]));
    i = MESH1a == 0;
    MESH1a(i) = NaN;
    MESH1a(1:mtl, 1:mtw) = MESH1;
    MESH2a = zeros(mbl, max([mtw, mbw]));
    i = MESH2a == 0;
    MESH2a(i) = NaN;
    MESH2a(1:mbl, 1:mbw) = MESH2;
    
    mstitch = [MESH1a(1:end-floor(overlap/2), :); MESH2a(ceil(overlap/2):end, :)];
    
    MESHout = mstitch;
end

if num == 3 | num == 4
    CToutt = handles.metricdata.stitch(1).stitchCT;
    RESoutt = handles.metricdata.stitch(1).stitch;
    MESHoutt = MESHout;
    IDoutt = handles.metricdata.section(1).name;
    
    pixellength2 = handles.metricdata.section(3).pixel;
    CToutb = handles.metricdata.section(3).CT;
    RESoutb = handles.metricdata.section(3).results;
    MESHoutb = handles.metricdata.section(3).CTresize;
    IDoutb = handles.metricdata.section(3).name;

    top = RESoutt;
    bot = RESoutb;

    startcore = handles.metricdata.section(1).startcore + handles.metricdata.section(2).startcore;
    % forcalc = tl - Val + m;
    cordepth = handles.metricdata.stitch(1).startoverlap + handles.metricdata.stitch(2).startoverlap;
    overlap = handles.metricdata.stitch(2).overlap;

    bot(:, 1) = bot(:, 1) + cordepth;

    stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];
    
    [mtl, mtw] = size(MESHoutt);
    [mbl, mbw] = size(MESHoutb);
    
    MESH1a = zeros(mtl, max([mtw, mbw]));
    i = MESH1a == 0;
    MESH1a(i) = NaN;
    MESH1a(1:mtl, 1:mtw) = MESHoutt;
    MESH2a = zeros(mbl, max([mtw, mbw]));
    i = MESH2a == 0;
    MESH2a(i) = NaN;
    MESH2a(1:mbl, 1:mbw) = MESHoutb;
    
    mstitch = [MESH1a(1:end-floor(overlap/2), :); MESH2a(ceil(overlap/2):end, :)];
    
    % Now stitch the image
    shift = handles.metricdata.stitch(1).shift + handles.metricdata.stitch(2).shift;
    topCT = CToutt;
    botCT = CToutb;
    
    
    [tl, tw] = size(topCT);
    [bl, bw] = size(botCT);
    [Vaw, WW] = min([tw, bw]);
    
    
    stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
    i = stitchCT == 0;
    stitchCT(i) = NaN;


    if shift >= 0
    %     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
        stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
    %     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
        catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
        end
    end
    if shift < 0
        stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    %     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
        try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        end
    %     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end

    i = stitchCT == 0;
    stitchCT(i) = NaN;

    CTout = stitchCT;
    RESout = stitch;
    MESHout = mstitch;
    IDout = handles.metricdata.section(1).name;
    
    i = sum(CTout, 'omitnan') == 0;
    CTout = CTout(:, ~i);
   
end

if num == 4
    CToutt = CTout;
    RESoutt = RESout;
    MESHoutt = MESHout;
    IDoutt = IDout;
    clear CTout RESout MESHout IDout
    
    pixellength2 = handles.metricdata.section(4).pixel;
    CToutb = handles.metricdata.section(4).CT;
    RESoutb = handles.metricdata.section(4).results;
    MESHoutb = handles.metricdata.section(4).CTresize;
    IDoutb = handles.metricdata.section(4).name;

    top = RESoutt;
    bot = RESoutb;

    startcore = handles.metricdata.section(1).startcore + handles.metricdata.section(2).startcore...
        + handles.metricdata.section(3).startcore;
    % forcalc = tl - Val + m;
    cordepth = handles.metricdata.stitch(1).startoverlap + handles.metricdata.stitch(2).startoverlap...
        + handles.metricdata.stitch(3).startoverlap;
    overlap = handles.metricdata.stitch(3).overlap;

    bot(:, 1) = bot(:, 1) + cordepth;

    stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];
    
    [mtl, mtw] = size(MESHoutt);
    [mbl, mbw] = size(MESHoutb);
    
    MESH1a = zeros(mtl, max([mtw, mbw]));
    i = MESH1a == 0;
    MESH1a(i) = NaN;
    MESH1a(1:mtl, 1:mtw) = MESHoutt;
    MESH2a = zeros(mbl, max([mtw, mbw]));
    i = MESH2a == 0;
    MESH2a(i) = NaN;
    MESH2a(1:mbl, 1:mbw) = MESHoutb;
    
    mstitch = [MESH1a(1:end-floor(overlap/2), :); MESH2a(ceil(overlap/2):end, :)];
    
    
    % Now stitch the image
    shift = handles.metricdata.stitch(1).shift + handles.metricdata.stitch(2).shift...
        + handles.metricdata.stitch(3).shift;
    topCT = CToutt;
    botCT = CToutb;
    
    
    [tl, tw] = size(topCT);
    [bl, bw] = size(botCT);
    [Vaw, WW] = min([tw, bw]);
    
    
    stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw tw]));
    i = stitchCT == 0;
    stitchCT(i) = NaN;


    if shift >= 0
    %     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
        stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
    %     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
        catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
        end
    end
    if shift < 0
        stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    %     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
        try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
        end
    %     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end

    i = stitchCT == 0;
    stitchCT(i) = NaN;

    CTout = stitchCT;
    RESout = stitch;
    MESHout = mstitch;
    IDout = handles.metricdata.section(1).name;
    
    i = sum(CTout, 'omitnan') == 0;
    CTout = CTout(:, ~i);
end


RadioPlot(CTout,...
    RESout,...
    MESHout,...
    IDout)

handles.metricdata.out.CTout = CTout;
handles.metricdata.out.RESout = RESout;
guidata(hObject, handles);

 
%Checkboxes
function secch1_Callback(hObject, eventdata, handles)
function secch2_Callback(hObject, eventdata, handles)
function secch3_Callback(hObject, eventdata, handles)
function secch4_Callback(hObject, eventdata, handles)


%Stitch Commands
function ST_Stitch_Callback(hObject, eventdata, handles)
RadioPlotNOPE(1, 2)

if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end

try top = handles.metricdata.section(stt).results;
catch errordlg('Process a DICOM for top section','Error');
end
try bot = handles.metricdata.section(stt+1).results;
catch errordlg('Process a DICOM for bottom section','Error');
end

[tl, tw] = size(top);
[bl, bw] = size(bot);
[Val, LL] = min([tl, bl]);

for n = 1:Val
    in = tl - Val + n;
    i2 = (top(in:end, 2) - bot(1:length(top(in:end, 2)), 2)).^2;
    i4(n) = sqrt(mean(abs(i2), 'omitnan'));
    clear i2
end

[lc, m] = min(i4);


startcore = handles.metricdata.section(stt).startcore;
forcalc = tl - Val + m;
cordepth = (tl - Val + m - startcore)/40;
overlap = abs(Val - m);
if overlap <= 0;
    overlap = 2;
end

bot(:, 1) = bot(:, 1) + cordepth;

stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];

% Now stitch the image
% pixellength = handles.metricdata.section(stt).pixel;
% pixellength2 = handles.metricdata.section(stt+1).pixel;
topCT = handles.metricdata.section(stt).CT;
botCT = handles.metricdata.section(stt+1).CT;


[tl, tw] = size(topCT);
[bl, bw] = size(botCT);
[Vaw, WW] = min([tw, bw]);

% calcst = zeros(20, (tw+bw));
calcst = zeros(1, (tw+bw));

i5 = zeros(length(Vaw), 1);

for n = 1:Vaw
    i = calcst > -9999;
    calcst(i) = 0;
%     calcst(1:20, ceil(Vaw/4):ceil(Vaw/4)+tw-1) = topCT(forcalc:forcalc+19, 1:tw);
    calcst(1, ceil(Vaw/4):ceil(Vaw/4)+tw-1) = topCT(forcalc, 1:tw);

%     calcst(:, n:bw+n-1) = (calcst(:, n:bw+n-1) - botCT(1:20, :)).^2;
    calcst(:, n:bw+n-1) = (calcst(:, n:bw+n-1) - botCT(1, :)).^2;
%     i5(n) = sqrt(mean(mean(calcst, 'omitnan'), 'omitnan'));
    i5(n) = sqrt((mean(calcst, 'omitnan')));
end
    
[lc, m] = min(i5);


shift = m - ceil(Vaw/4);


stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
i = stitchCT == 0;
stitchCT(i) = NaN;
id1t = length(botCT(floor(overlap/2)+2:end, 1));
id1c = length(botCT(floor(overlap/2)+1:end, 1));

overlap
bw
tw
size(stitchCT)
size(botCT)
size(topCT)


if shift >= 0
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    try   stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end
end
if shift < 0
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    end
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
end

i = stitchCT <= 0;
stitchCT(i) = NaN;


handles.metricdata.stitch(stt).stitch = stitch;
handles.metricdata.stitch(stt).stitchCT = stitchCT;
handles.metricdata.stitch(stt).topCT = topCT;
handles.metricdata.stitch(stt).botCT = botCT;
handles.metricdata.stitch(stt).top = top;
handles.metricdata.stitch(stt).bot = bot;
handles.metricdata.stitch(stt).startoverlap = cordepth;
handles.metricdata.stitch(stt).overlap = overlap;
handles.metricdata.stitch(stt).shift = shift;
guidata(hObject, handles);

RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)
function ST_Back_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end

stitchCT = handles.metricdata.stitch(stt).stitchCT;
topCT = handles.metricdata.stitch(stt).topCT;
overlap = handles.metricdata.stitch(stt).overlap;
cordepth = length(topCT(1:end-floor(overlap/2), 1))+1;

shiftn = handles.metricdata.stitch(stt).shift + 1;

[l, w] = size(stitchCT);
newp = zeros([l, w+1]);
newp(1:cordepth, 1:end-1) = stitchCT(1:cordepth, :);
newp(cordepth:end, 2:end) = stitchCT(cordepth:end, :);

i = newp == 0;
newp(i) = NaN;

handles.metricdata.stitch(stt).stitchCT = newp;
handles.metricdata.stitch(stt).shift = shiftn;
guidata(hObject, handles);

RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)

function ST_Forward_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end

stitchCT = handles.metricdata.stitch(stt).stitchCT;
topCT = handles.metricdata.stitch(stt).topCT;
overlap = handles.metricdata.stitch(stt).overlap;
cordepth = length(topCT(1:end-floor(overlap/2), 1))+1;

shiftn = handles.metricdata.stitch(stt).shift - 1;

[l, w] = size(stitchCT);
newp = zeros([l, w+1]);
newp(1:cordepth, 2:end) = stitchCT(1:cordepth, :);
newp(cordepth:end, 1:end-1) = stitchCT(cordepth:end, :);

i = newp == 0;
newp(i) = NaN;

handles.metricdata.stitch(stt).stitchCT = newp;
handles.metricdata.stitch(stt).shift = shiftn;
guidata(hObject, handles);

RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)
function ST_up_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end


top = handles.metricdata.stitch(stt).top;
bot = handles.metricdata.stitch(stt).bot;
cordepth = handles.metricdata.stitch(stt).startoverlap;
overlap = handles.metricdata.stitch(stt).overlap;

bot(:, 1) = bot(:, 1) +  1/40;
overlap = overlap - 1;

stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];

% Now stitch the image
topCT = handles.metricdata.stitch(stt).topCT;
botCT = handles.metricdata.stitch(stt).botCT;
shift = handles.metricdata.stitch(stt).shift;

[tl, tw] = size(topCT);
[bl, bw] = size(botCT);

stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
i = stitchCT == 0;
stitchCT(i) = NaN;


if shift >= 0
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end
end
if shift < 0
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    end
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
end

i = stitchCT == 0;
stitchCT(i) = NaN;

handles.metricdata.stitch(stt).top = top;
handles.metricdata.stitch(stt).bot = bot;
handles.metricdata.stitch(stt).stitch = stitch;
handles.metricdata.stitch(stt).stitchCT = stitchCT;
handles.metricdata.stitch(stt).overlap = overlap;
guidata(hObject, handles);


RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)
function ST_down_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end


top = handles.metricdata.stitch(stt).top;
bot = handles.metricdata.stitch(stt).bot;
cordepth = handles.metricdata.stitch(stt).startoverlap;
overlap = handles.metricdata.stitch(stt).overlap;

bot(:, 1) = bot(:, 1) -  1/40;
overlap = overlap + 1;

stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];

% Now stitch the image
topCT = handles.metricdata.stitch(stt).topCT;
botCT = handles.metricdata.stitch(stt).botCT;
shift = handles.metricdata.stitch(stt).shift;

[tl, tw] = size(topCT);
[bl, bw] = size(botCT);

stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
i = stitchCT == 0;
stitchCT(i) = NaN;


if shift >= 0
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end
end
if shift < 0
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    end
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
end

i = stitchCT == 0;
stitchCT(i) = NaN;

handles.metricdata.stitch(stt).top = top;
handles.metricdata.stitch(stt).bot = bot;
handles.metricdata.stitch(stt).stitch = stitch;
handles.metricdata.stitch(stt).stitchCT = stitchCT;
handles.metricdata.stitch(stt).overlap = overlap;
guidata(hObject, handles);


RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)


% --- Executes on button press in ST_upC.
function ST_upC_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end


top = handles.metricdata.stitch(stt).top;
bot = handles.metricdata.stitch(stt).bot;
cordepth = handles.metricdata.stitch(stt).startoverlap;
overlap = handles.metricdata.stitch(stt).overlap;

bot(:, 1) = bot(:, 1) +  1/2;
overlap = overlap - 20;

stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];

% Now stitch the image
topCT = handles.metricdata.stitch(stt).topCT;
botCT = handles.metricdata.stitch(stt).botCT;
shift = handles.metricdata.stitch(stt).shift;

[tl, tw] = size(topCT);
[bl, bw] = size(botCT);

stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
i = stitchCT == 0;
stitchCT(i) = NaN;


if shift >= 0
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end
end
if shift < 0
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    end
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
end

i = stitchCT == 0;
stitchCT(i) = NaN;

handles.metricdata.stitch(stt).top = top;
handles.metricdata.stitch(stt).bot = bot;
handles.metricdata.stitch(stt).stitch = stitch;
handles.metricdata.stitch(stt).stitchCT = stitchCT;
handles.metricdata.stitch(stt).overlap = overlap;
guidata(hObject, handles);


RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)


% --- Executes on button press in ST_downC.
function ST_downC_Callback(hObject, eventdata, handles)
if get(handles.st12, 'Value') == 1;
    stt = 1;
else if get(handles.st23, 'Value') == 1;
        stt = 2;
    else if get(handles.st34, 'Value') == 1;
            stt = 3;
            else errordlg('Select Compile Sections Option to Stitch','Error');
        end
    end
end


top = handles.metricdata.stitch(stt).top;
bot = handles.metricdata.stitch(stt).bot;
cordepth = handles.metricdata.stitch(stt).startoverlap;
overlap = handles.metricdata.stitch(stt).overlap;

bot(:, 1) = bot(:, 1) -  1/2;
overlap = overlap + 20;

stitch = [top(1:end-floor(overlap/2), :); bot(ceil(overlap/2):end, :)];

% Now stitch the image
topCT = handles.metricdata.stitch(stt).topCT;
botCT = handles.metricdata.stitch(stt).botCT;
shift = handles.metricdata.stitch(stt).shift;

[tl, tw] = size(topCT);
[bl, bw] = size(botCT);

stitchCT = zeros(length(topCT(1:end-floor(overlap/2), 1)) + length(botCT(ceil(overlap/2):end, 1)), abs(shift)+ max([bw, tw]));
i = stitchCT == 0;
stitchCT(i) = NaN;


if shift >= 0
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1+shift:tw+shift) = topCT(1:end-floor(overlap/2), :);    
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
    end
end
if shift < 0
    stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), end-tw+1:end) = topCT(1:end-floor(overlap/2), :);
%     stitchCT(1:length(topCT(1:end-floor(overlap/2), 1)), 1:tw) = topCT(1:end-floor(overlap/2), :);    
    try stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    catch stitchCT(end-length(botCT(floor(overlap/2)+1:end, 1)):end, 1:bw) = botCT(ceil(overlap/2):end, :);
    end
%     stitchCT(end-length(botCT(floor(overlap/2)+2:end, 1)):end, end-bw+1:end) = botCT(ceil(overlap/2):end, :);
end

i = stitchCT == 0;
stitchCT(i) = NaN;

handles.metricdata.stitch(stt).top = top;
handles.metricdata.stitch(stt).bot = bot;
handles.metricdata.stitch(stt).stitch = stitch;
handles.metricdata.stitch(stt).stitchCT = stitchCT;
handles.metricdata.stitch(stt).overlap = overlap;
guidata(hObject, handles);


RadioStitchPlot(handles.metricdata.section(stt).CT,...
    handles.metricdata.stitch(stt).top,...
    handles.metricdata.section(stt+1).CT,...
    handles.metricdata.stitch(stt).bot,...
    handles.metricdata.stitch(stt).stitch,...
    handles.metricdata.stitch(stt).stitchCT,...
    handles.metricdata.stitch(stt).startoverlap)


%Create Outputs
function Create_Out_Callback(hObject, eventdata, handles)
[FileName, PathName] = uiputfile('*.dpro');

data = handles.metricdata.out.RESout;

fid = fopen([PathName '/' FileName], 'w');
fprintf(fid, '%s,', ['Depth_cm']);
fprintf(fid, '%s,', ['HU_Number']);
fprintf(fid, '%s,', ['StDev']);
fprintf(fid, '%s \r\n', ['Pixel_Count']);
fprintf(fid, '%f, %f, %f, %f \r\n', data');
fclose(fid);

outim = uint16(handles.metricdata.out.CTout);
imwrite(outim, [PathName '/' FileName '.tiff'], 'tiff');

%Clear All
function Clear_All_Callback(hObject, eventdata, handles)

RadioPlotNOPE(1,2)

% set(handles.path_display, 'String', 'No Path Selected');
% set(handles.data_load_display, 'String', 'No Data Loaded')
try
    handles = rmfield(handles, 'metricdata');
end
guidata(hObject, handles);
close all
SedCT



%Menu Items
function File_Callback(hObject, eventdata, handles)
function Open_Callback(hObject, eventdata, handles)
Input_folder_Callback(hObject, eventdata, handles)
function output_Callback(hObject, eventdata, handles)
Create_Out_Callback(hObject, eventdata, handles)
function clear_Callback(hObject, eventdata, handles)
Clear_All_Callback(hObject, eventdata, handles)
function Exit_Callback(hObject, eventdata, handles)
close SedCT

function Processing_Callback(hObject, eventdata, handles)
function Load_Callback(hObject, eventdata, handles)
Load_data_Callback(hObject, eventdata, handles)
function Process_Callback(hObject, eventdata, handles)
process_data_Callback(hObject, eventdata, handles)

function Analysis_Callback(hObject, eventdata, handles)
function SedCTimage_Callback(hObject, eventdata, handles)
SedCTimage

%Mystery Buttons
function Done_Callback(hObject, eventdata, handles)





%display
function data_load_display_Callback(hObject, eventdata, handles)

function data_load_display_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
