Contents

% To run this progam without GUI, uncomment next line:
 %S.GUI_flag = 0;

Clear Variables. Turn off Obnoxious warnings...

if S.GUI_flag == 1       % program was called from the GUI
else
    clear all;
    close all;
    S.GUI_flag = 0;
end

%Image Filename and path
%warning('off','MATLAB:MKDIR:DirectoryExists');                                                                                 % Elliminates warning that the directory being made already exists
warning('off','MATLAB:nearlySingularMatrix');                                                                                   % Your fit sucked.  Turns off the warning
warning('off','MATLAB:polyfit:PolyNotUnique');                                                                                  % Fit still sucked.  Turns off the warning
warning('off','Images:initSize:adjustingMag');                                                                                  % Image was too big to display on the screen. Matlab resized it.  Turns the warning off.
warning('off','curvefit:fit:equationBadlyConditioned');                                                                         % Your fit sucked.  Turns off the warning
warning('off','curvefit:prepareSurfaceData:removingNaNAndInf');                                                                 % NaN's in the thresholded data.  Turns off the warning
warning('off','curvefit:prepareSurfaceData:swapXAndYForTable');                                                                 % Image isn't square.  Turns off the warning
warning('off','MATLAB:MKDIR:DirectoryExists');
Undefined function or variable 'S'.

Error in ==> Nano6_Flatten_surf_poly_GUI at 6
if S.GUI_flag == 1       % program was called from the GUI

Set relevant input paramters. Scalebar sizes, fitting limits and so forth

ScaleBar                    = 1000;                                                                                             % ScaleBar Lenght in nm

if S.GUI_flag == 1
         FontSize           = 11;                                                                                              % Font Size for images
         FontName           = 'Arial';                                                                                          % Font for images
else
    FontSize                = 16;                                                                                              % Font Size for images
    FontName                = 'Times';                                                                                          % Font for images
end

print_dpi                   = 150;                                                                                              % DPI for images.  72 for presentations.  300 for publications

amplitude_fract             = 0.95;                                                                                             % Amplitude variability for fitting histogramed peaks to determine thresholds
crop_fraction               = 0.0;

channels_nb = 0;

if S.GUI_flag                       % program was called from the GUI
                                    % so read values from the GUI

    existingChannels = zeros(14, 1);
    for k = 1:7
            channel_name = get(S.channelProcessText(k), 'string');
            channel_name = channel_name{1};
        switch channel_name
                case 'Height'
                     existingChannels(1) = 1;
                case 'Deflection'
                     existingChannels(2) = 1;
                case 'Phase'
                     existingChannels(3) = 1;
                case 'Amplitude'
                     existingChannels(4) = 1;
                case 'Amplitude Error'
                     existingChannels(5) = 1;
                case 'Potential'
                     existingChannels(6) = 1;
                case 'Peak Force Error'
                     existingChannels(7) = 1;
                case 'DMTModulus'
                     existingChannels(8) = 1;
                case 'LogDMTModulus'
                     existingChannels(9) = 1;
                case 'Adhesion'
                     existingChannels(10) = 1;
                case 'Deformation'
                     existingChannels(11) = 1;
                case 'Dissipation'
                     existingChannels(12) = 1;
                case 'Input1'
                     existingChannels(13) = 1;
                case 'Input2'
                     existingChannels(14) = 1;
        end
    end

    channels_indices = find(existingChannels); % returns the index of the non-zeros elements of existingChannels

    % Get all the parameters for each channel:
    for k = 1:14
        Med_corr(k)                 = get(S.Med_corr(k), 'value');
        Med_corr_2(k)               = get(S.Med_corr_2(k), 'value');
        vert_med_corr(k)            = get(S.vert_med_corr(k), 'value');
        poly_order_horiz(k)         = get(S.poly_order_horiz(k), 'value')-1;
        poly_order_vert(k)          = get(S.poly_order_vert(k), 'value')-1;
        z_height_reference(k)       = str2num(get(S.z_height_reference(k), 'string'));
        min_spacing_stdev_user(k)   = str2num(get(S.min_spacing_stdev_user(k), 'string'));
        max_spacing_stdev_user(k)   = str2num(get(S.max_spacing_stdev_user(k), 'string'));
        min_z(k)                    = str2num(get(S.min_z(k), 'string'));
        max_z(k)                    = str2num(get(S.max_z(k), 'string'));
        peaks_fit(k)                = get(S.peaks_fit(k), 'value')-1;
        peak_reference_height(k)    = get(S.peak_reference_height(k), 'value')-1;
        peak_to_flatten(k)          = get(S.peak_to_flatten(k), 'value')-1;
        Calc_flatten_lim(k)         = get(S.Calc_flatten_lim(k), 'value');
        overwrite_xlim(k)           = get(S.overwrite_xlim(k), 'value');
        sort_peaks(k)               = get(S.sort_peaks(k), 'value');
    end

    % Remove the parameters for the channels that do not exist in the file
    Med_corr                 = Med_corr(existingChannels~=0);
    Med_corr_2               = Med_corr_2(existingChannels~=0);
    vert_med_corr            = vert_med_corr(existingChannels~=0);
    poly_order_horiz         = poly_order_horiz(existingChannels~=0);
    poly_order_vert          = poly_order_vert(existingChannels~=0);
    peaks_fit                = peaks_fit(existingChannels~=0);
    peak_reference_height    = peak_reference_height(existingChannels~=0);
    peak_to_flatten          = peak_to_flatten(existingChannels~=0);
    Calc_flatten_lim         = Calc_flatten_lim(existingChannels~=0);
    overwrite_xlim           = overwrite_xlim(existingChannels~=0);
    sort_peaks               = sort_peaks(existingChannels~=0);
    z_height_reference       = z_height_reference(existingChannels~=0);
    min_spacing_stdev_user   = min_spacing_stdev_user(existingChannels~=0);
    max_spacing_stdev_user   = max_spacing_stdev_user(existingChannels~=0);
    min_z                    = min_z(existingChannels~=0);
    max_z                    = max_z(existingChannels~=0);

    for k = 1:7
        channels_flatten(k)         = get(S.channelProcessSelect(k), 'value');
        if channels_flatten(k) == 1
            channels_nb = channels_nb + 1;
        end
    end

    % Remove the parameters for the channels that are not going to be processed:
    Med_corr                 = Med_corr(channels_flatten~=0);
    Med_corr_2               = Med_corr_2(channels_flatten~=0);
    vert_med_corr            = vert_med_corr(channels_flatten~=0);
    poly_order_horiz         = poly_order_horiz(channels_flatten~=0);
    poly_order_vert          = poly_order_vert(channels_flatten~=0);
    peaks_fit                = peaks_fit(channels_flatten~=0);
    peak_reference_height    = peak_reference_height(channels_flatten~=0);
    peak_to_flatten          = peak_to_flatten(channels_flatten~=0);
    Calc_flatten_lim         = Calc_flatten_lim(channels_flatten~=0);
    overwrite_xlim           = overwrite_xlim(channels_flatten~=0);
    sort_peaks               = sort_peaks(channels_flatten~=0);
    z_height_reference       = z_height_reference(channels_flatten~=0);
    min_spacing_stdev_user   = min_spacing_stdev_user(channels_flatten~=0);
    max_spacing_stdev_user   = max_spacing_stdev_user(channels_flatten~=0);
    min_z                    = min_z(channels_flatten~=0);
    max_z                    = max_z(channels_flatten~=0);
    channels_indices         = channels_indices(channels_flatten~=0);

%     channels_flatten = 1:channels_nb;


    Offset_data = zeros(channels_nb, 1);
    fit_hist = zeros(channels_nb, 1);
    for k = 1:channels_nb
         if z_height_reference(k) == 0
            Offset_data(k) = 0;
         else
            Offset_data(k) = 1;
         end
         if peaks_fit(k) == 0
             fit_hist(k) = 0;
         else
             fit_hist(k) = 1;
         end
    end

    histogram_bins              = str2num(get(S.histogram_bins, 'string'));
    smoothing_level             = str2num(get(S.smoothing_level, 'string'));
%     Calc_flatten_lim            = get(S.Calc_flatten_lim, 'value');
    thresh_range                = str2num(get(S.thresh_range, 'string'));
%     overwrite_xlim              = get(S.overwrite_xlim, 'value');
%     fit_hist                    = get(S.fit_hist, 'value');
    min_peak_height             = str2num(get(S.min_peak_height, 'string'));
    min_peak_distance           = str2num(get(S.min_peak_distance, 'string'));
    peak_threshold              = str2num(get(S.peak_threshold, 'string'));
%     sort_peaks                  = get(S.sort_peaks, 'value');


else    % program wasn't called from the GUI
    % 0 = 'N' (No), 1 = 'Y' (Yes)
    Med_corr                    = [1,0,0,1,1,1,1];                                                                    % For each channel, apply or not the median correction
    Med_corr_2                  = [1,0,0,1,1,1,1];
    vert_med_corr               = [0,0,0,0,0,0,0];

    initial_poly_order_horiz    = [2,         0,        0,        0,        0,        0,        0];
    initial_poly_order_vert     = [2,         0,        0,        0,        0,        0,        0];
    poly_order_horiz            = [3,         0,        0,        0,        0,        1,        1];                                     % Polynomial order for horizontal fitting
    poly_order_vert             = [3,         0,        0,        0,        0,        1,        1];                                     % Order for vertical fit

    Offset_data                 = [1,         0,        0,        0,        0,        0,        0];  %'Y';
    histogram_bins              = 256;
    smoothing_level             = 6;
%   z_off_thresh                = [ 2.0e-2,   8.0e-3,   1.3e-2,   1.0e-2,   5.0e-3,   2.0e-2,   0.0e-0];                                % Min peak threshold to set offset
    z_height_reference          = [ 5.0e-0,   0.0e-0,   0.0e-0,   0.0e-0,   0.0e-0,   0.0e-0,   0.0e-0];                                % Reference Z height in different units for each channel

    Calc_flatten_lim            = [1, 1, 1, 1, 1, 1, 1];%'Yes';
    thresh_range                = 3;  %2.5
    min_spacing_stdev_user      = [ 5.0e-2,   2.0e+0,   2.0e-1,   5.0e-2,   1.0e+1,   5.0e-3,   1.0e-3];                                % Min height stdev in different units for each channel
    max_spacing_stdev_user      = [ 5.0e-1,   2.0e+1,   2.0e+0,   1.0e-0,   4.0e+1,   2.0e-2,   5.0e-2];                                % Max height stdev in different units for each channel

    overwrite_xlim              = [1, 1, 1, 1, 1, 1, 1]; %'Yes';
    min_z                       = [0.0e+0,  -3.0e+1,  -5.0e+0,  -1.0e+1,  -1.0e+2,  -1.0e-2,  -1.0e-1];                                % Min channel display values in different units for each channel
    max_z                       = [ 10.0e+0,   3.0e+1,   5.0e+0,   1.0e+1,   1.0e+2,   3.0e-2,   1.0e-1];                                % Max channel display values in different units for each channel


    channels_flatten            = [1 0 0 0 0 0 0]; %[1:3];
    channels_indices            = find(channels_flatten); % indices of the non-zero elements of channels_flatten %[1:3];
    channels_nb                 = length(find(channels_flatten)); % number of non-zero elements of channels_flatten

    fit_hist                    = [1,         1,        1,        1,        1,        1,        1];  %'Yes';
    peaks_fit                   = [2,1,1,2,1,2,1];   %2

    peak_reference_height       = [2, 1, 1, 1, 1, 1, 1];    %2                                                                            % Peak for the offset shift (peaks sorted from smallest to highest position)

    peak_to_flatten             = [2, 1, 1, 1, 1, 1, 1]; %2                                                                               % Peak to be flattened for each channel; if sort_peaks = N, then peak_to_flatten = 1 for flattening the peak with highest amplitude

    min_peak_height             = 1e2;                                                                                                    % Min peak height for the peak finding
    min_peak_distance           = 5;                                                                                                 % Min distance between peaks for the peak finding
    peak_threshold              = 5;                                                                                                    % To find peaks that are at least greater than their neighbors by the threshold peak_threshold
    sort_peaks                  = [1, 1, 1, 1, 1, 1, 1]; %'Y';                                                                                              % if N, peaks are sorted from highest amplitude to smallest one,
                                                                                                                                        % and if Y, peaks are sorted from smallest index to highest index

end

disp_percent                = 1; %'Yes';
percent_search              = [1 10 90 99];

sequential_flattening       = 1; %'Y';                                                                                               % If Y, the flattening is done sequentially; if N, the final mask is directly applied to the inital image after first median correction

plot_intermediate_steps     = 0; %'N';                                                                                                  % To decide if the intermediate images and histograms at each step are plotted or not
plot_intermediate_image     = 0; %'N';                                                                                                  % To decide if the image and histogram after first median corr and global tilt corr are plotted or not

peak_pos_tolerance          = 5.0e-1;                                                                                               % Tolerance on the displacement of each peak from one step to the other for a good identification of each peak

Select Files to read

if S.GUI_flag                  % program was called from the GUI

    Filename = evalin('base', 'Filename'); % retrieves the variable Filename that had been saved from the GUI
    ImageDir = get(S.selectFolderText, 'string');
    ImageDir = [ImageDir,filesep];

    if S.processAll  == 1        % Process all images
            %Filename = get(S.selectFileText, 'string');
            Filename_length = length(Filename);
    else
              %Filename = {Filename, Filename};
              Filename = Filename{1};
              Filename = {Filename, Filename};
              Filename_length = 1;
    end
else        % program was not called from the GUI
        [Filename, ImageDir]        = uigetfile(    {'*.*','Nanoscope Files (*.*)'}, ...                                                % Allows for user selection of the images to process
                                                'Pick a file',...
                                                'C:\Users\LBNI\Desktop\Test ata\',...                            % Default starting directory.  Useful to change depending on where files are saved
                                                'MultiSelect','on');                                                                % Allows for multiple images to be selected.  At least 2 must be selected
        Filename_length = length(Filename);
end

Prep for file output

mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep]);                                                                                       % Makes a new directory for the final jpeg output
mkdir([ImageDir,'Matlab',filesep,'Figure Fig',filesep]);                                                                                        % Makes a new directory for the final fig output

Loop the import and processing routiene for the number of selected files

load('Sky','mycmap_sky')

    searchstring(1).label='@2:Z scale:';
    searchstring(2).label='Samps';
    searchstring(3).label='Lines:';
    searchstring(4).label='Data offset';
    searchstring(5).label='Image Data:';
    searchstring(6).label='Scan Size:';
    searchstring(7).label='@Sens. Zsens:';
    searchstring(8).label='@Z magnify:';

    searchstring(9).label='Scan Rate:';
    searchstring(10).label='@Sens. DeflSens: V';
    searchstring(11).label='@2:CantFrequency: V';
    searchstring(12).label='@2:CantDrive: V';
    searchstring(13).label='@2:TMSetAmplitude: V';
    Init_poly_coef      = cell(length(Filename),length(channels_flatten));

for index = 1:Filename_length                                                                                               % Loops through the number of files
   file_name = [ImageDir,Filename{index}];                                                                                     % Generates the filname to open for this particular iteration

Read in DI Image and image propreties This code is mostly from Jaco De Groot at University College London

    % openNano version 6.0
    % Last updated: 01-06-2006 by Macerano Blanco for correctly reading
    % Nanoscope 6 files
    %
    % [images Channel_Info ]=openNano(filename)
    %
    % Reads images from a Nanoscope 6 image files. The images are stored in
    % "images" and can be displayed with imshow(images(:,:,channel_number). The
    % relevant channel number can be found in the print output fromt this subroutine
    % Basic channel information can be found in Channel_Info that for each
    % channel contains, the name, trace/retrace, unit and x and y scan size
    %
    % Example:
    % [images channel_info]=openNano('C:\Briefcase\Parchment\AFM\Dimension\CR41\CR41.000');
    % The output of this command is:
    %  Channel                                     Unit
    %   1. Height                   Trace          um
    %   2. Deflection               Trace          V
    %
    % The following command images the trace height. The height is specified in
    % microns:
    % imshow(images(:,:,1),[]);
    % This image can be stored as a gray scale image with
    % imwrite(mat2gray(images(:,:,1)),'g:/traceheight.tif','tif');

    % This function/script is authorized for use in government and academic
    % research laboratories and non-profit institutions only. Though this
    % function has been tested prior to its posting, it may contain mistakes or
    % require improvements. In exchange for use of this free product, we
    % request that its use and any issues that may arise be reported to us.
    % Comments and suggestions are therefore welcome and should be sent to
    %
    % Jaco de Groot
    % Bone & Mineral Centre, Dept of Medicine
    % University College London
    % 5 University Street
    % The Rayne Building
    % London    WC1E 6JJ
    % Tel: 020-7679 6143
    % Mob: 07745-948245
    % Fax: 020-7679 6219
    % Email: jacodegroot@yahoo.se
    % http://www.ucl.ac.uk/medicine/bmc/


    % Define End of file identifier
        % Opend the file given in argument and reference as
        % fid.  Also if there was an error output error
        % number and error message to screen


    fid = fopen(file_name,'r');
        [message,errnum] = ferror(fid);
        if(errnum)
            fprintf(1,'I/O Error %d \t %s',[errnum,message]);
            %   break
        end
        header_end=0; eof = 0; counter = 1; byte_location = 0;
        nstrings=size(searchstring,2);
        parcounter = ones(nstrings,1);
        paramters = struct('trace',cell(1,nstrings),'values',cell(1,nstrings),'channel',cell(1,nstrings));

        for ij=1:nstrings
            parcounter(ij)=1;
            parameters(ij).trace=0;
        end;

        while( and( ~eof, ~header_end ) )

        byte_location = ftell(fid);
        line = fgets(fid);

        for ij=1:nstrings
            if strfind(line,searchstring(ij).label)
                if (extract_num(line))
                    b=strfind(line, 'LSB');
                     if (b>0)
                        parameters(ij).values(parcounter(ij))=extract_num(line(b(1):end));
                     else parameters(ij).values(parcounter(ij))=extract_num(line);
                    end;

                else
                    b= strfind(line,'"');
                    if (b>0)
                        parameters(ij).channel(parcounter(ij)).name=line(b(1)+1:b(2)-1);
                    else
                        if (strfind(line,'Trace')>0)
                            parameters(ij).trace=1;
                        end;
                    end;
                end;
                parcounter(ij)=parcounter(ij)+1;
            end
        end;

        if( (-1)==line )
            eof  = 1;
        end

        if ~isempty( strfind( line, '\*File list end' ) )
            header_end = 1;
        end

        counter=counter+1;
        end
    fclose(fid);

    param           = parameters;
    scaling         = param(1).values;                                                                                                                                                 % Scaling parameters
    spl             = param(2).values;                                                                                                                                                 % Samples per line
    linno           = param(3).values;                                                                                                                                                 % No of lines
    image_pos       = param(4).values;                                                                                                                                                 % Data position
    ScanSize        = param(6).values(1);
    Z_Sensitivity   = param(7).values;                                                                                                                                                 % Data position
    Z_Magnification = param(8).values;                                                                                                                                                 % Zmagnification

    Scan_Rate       = param(9).values(1);
    Defl_Sens       = param(10).values;
    Drive_Freq      = param(11).values;
    Drive_Amp       = param(12).values;
    Amp_Set_Point   = param(13).values;

    [tok remain]    = strtok(Filename{index},'.');
    file_number     = str2double(remain(2:4));

    Image_Data_File = [ImageDir,'Matlab',filesep,'Imaging_Parameters.txt'];
    if  index == 1
        fid = fopen(Image_Data_File, 'a');
        fprintf(fid, [  'File Number','\t',...
                        'Samples per Line','\t',...
                        'Lines per Image','\t',...
                        'Scan Rate (Hz)','\t',...
                        'Deflection Sensitivity (nm/V)','\t',...
                        'Drive Frequency (kHz)','\t',...
                        'Drive Amplitude (mV)','\t',...
                        'Amplitude Setpoint (mV)','\n'],...
                        'char');

        fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                        num2str(spl(1),'%0.0f'),'\t',...
                        num2str(linno,'%0.0f'),'\t',...
                        num2str(Scan_Rate,'%0.3f'),'\t',...
                        num2str(Defl_Sens,'%0.3f'),'\t',...
                        num2str(Drive_Freq,'%0.2f'),'\t',...
                        num2str(Drive_Amp,'%0.2f'),'\t',...
                        num2str(Amp_Set_Point,'%0.2f'),'\n'],...
                        'char');
        fclose(fid);
    else
        fid = fopen(Image_Data_File, 'a');
        fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                        num2str(spl(1),'%0.0f'),'\t',...
                        num2str(linno,'%0.0f'),'\t',...
                        num2str(Scan_Rate,'%0.3f'),'\t',...
                        num2str(Defl_Sens,'%0.3f'),'\t',...
                        num2str(Drive_Freq,'%0.2f'),'\t',...
                        num2str(Drive_Amp,'%0.2f'),'\t',...
                        num2str(Amp_Set_Point,'%0.2f'),'\n'],...
                        'char');
        fclose(fid);
    end


    L = length(image_pos);                                                                                                                                                             % Sets the number of image channels
    channel_info = struct('Trace',cell(L,1),'Name',cell(L,1),'Finalscaling',cell(L,1),'Unit',cell(L,1));
    finalscaling = zeros(L,1);
    for im=1:L                                                                                                                                                                         % Loops for the number of channels
        channel_info(im).Trace=char((param(6).trace)*'Trace  '+(1-param(6).trace)*'Retrace');                                                                                          % Determines if Channel is Trace or Retrace
        channel_info(im).Name=param(5).channel(im).name;                                                                                                                               % Name of channel 1
        channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);                                                                                                     % Computes some of the scaling paramters
        mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep]);
        mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep,'hist',filesep]);
        mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep,'grey',filesep]);

        switch channel_info(im).Name                                                                                                                                                   % Sets Units according to recorded channel
            case 'Height'                                                                                                                                                              % If height...
               channel_info(im).Unit='nm' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);                                                                                                                    % Sets bit scaling

            case 'Deflection'                                                                                                                                                          % If Deflection...
               channel_info(im).Unit='V';                                                                                                                                              % Units are V
               finalscaling(im)=(scaling(im))/(2^16);                                                                                                                                  % Sets bit scaling

            case 'Phase'                                                                                                                                                               % If Phase...
                channel_info(im).Unit='Degree';                                                                                                                                        % Units are degrees
                finalscaling(im)=(scaling(im))/(2^16);                                                                                                                                 % Sets bit scaling

            case 'Amplitude'                                                                                                                                                           % If Potential...
                channel_info(im).Unit='mV';                                                                                                                                             % Units are V
                channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);                                                                                             % Sets final scaling value
                finalscaling(im)=(scaling(im))/(2^16);

            case 'Amplitude Error'                                                                                                                                                           % If Potential...
                channel_info(im).Unit='mV';                                                                                                                                             % Units are V
                channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);                                                                                             % Sets final scaling value
                finalscaling(im)=(scaling(im))/(2^16);

            case 'Potential'                                                                                                                                                           % If Potential...
                channel_info(im).Unit='V';                                                                                                                                             % Units are V
                channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);                                                                                             % Sets final scaling value
                finalscaling(im)=(scaling(im))/(2^16);                                                                                                                                 % Sets bit scaling

            case 'Peak Force Error'                                                                                                                                                              % If height...
               channel_info(im).Unit='nN' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'DMTModulus'                                                                                                                                                              % If height...
               channel_info(im).Unit='MPa' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'LogDMTModulus'                                                                                                                                                              % If height...
               channel_info(im).Unit='log(MPa)' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'Adhesion'                                                                                                                                                              % If height...
               channel_info(im).Unit='nN  ' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'Deformation'                                                                                                                                                              % If height...
               channel_info(im).Unit='nm' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'Dissipation'                                                                                                                                                              % If height...
               channel_info(im).Unit='eV' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);

            case 'Input1'                                                                                                                                                              % If height...
               channel_info(im).Unit='nm' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity*Defl_Sens*1.50)/(2^16);


            case 'Input2'                                                                                                                                                              % If height...
               channel_info(im).Unit='nm' ;                                                                                                                                            % Units are nm
               channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);                                                                              % Sets final scaling value
               finalscaling(im)=(scaling(im)*Z_Sensitivity*22*Defl_Sens)/(2^16);

            otherwise                                                                                                                                                                  % If unknown/other
                channel_info(im).Unit='V';                                                                                                                                             % Units are V
                channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);                                                                                             % Sets final scaling value
                finalscaling(im)=(scaling(im))/(2^16);                                                                                                                                 % Sets bit scaling
        end;
    end;

    if S.GUI_flag                  % program was called from the GUI
        if S.processAll  == 0        % Process one image
            disp({['Image ',num2str(index),' of 1'];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});    % Displays progress and scaling
        else
           disp({['Image ',num2str(index),' of ',num2str(length(Filename))];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});    % Displays progress and scaling
        end
    else
        disp({['Image ',num2str(index),' of ',num2str(length(Filename))];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});    % Displays progress and scaling
    end

    fid = fopen(file_name,'r');

    images = zeros(linno, spl(1), L);
    for i = 1:L
       fseek(fid,image_pos(i),-1);
       A = fread(fid, [spl(i) linno],'int16');
       images(:,:,i) = rot90(finalscaling(i)*A);
    end;

Crop Images. Helps deal with some of the turn around ripple

images   = images(ceil(crop_fraction/2*size(images,1))+1:floor((1-crop_fraction/2)*size(images,1))-1,ceil(crop_fraction/2*size(images,2))+1:floor((1-crop_fraction/2)*size(images,2))-1,:);
images   = images(:,:, channels_flatten == 1);  % keep only the channels to be processed
channel_info   = channel_info(channels_flatten == 1);  % keep only the info about channels to be processed

This section actuall does all of the image processing. Different blocks will be explained below Code Loops through each of the channels and applies all of the corrections

    for j                                           = 1:channels_nb                                                                                                                  % Loops for the number of different channels in the image

        disp({'Channel:'; channel_info(j).Name})

     try
        if S.GUI_flag                       % program was called from the GUI
            if S.processAll  ==  0  && j == 1         % process only one image AND first channel
                raw_images = zeros(size(images,1), size(images,2), 14);    % Initialization
                med_corr_1_images = zeros(size(images,1), size(images,2), 14);    % Initialization
                poly_flat_1_images = zeros(size(images,1), size(images,2), 14);    % Initialization
                med_corr_2_images = zeros(size(images,1), size(images,2), 14);    % Initialization
                poly_flat_2_images = zeros(size(images,1), size(images,2), 14);    % Initialization
                final_images = zeros(size(images,1), size(images,2), 14);    % Initialization
            end
        end

        if S.GUI_flag                       % program was called from the GUI
            if S.processAll  ==  0          % process only one image
                raw_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
            end
        end

        % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' raw.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);


        max_peak_ampl               = 0;                                                                                                    % Variable to keep track of the maximum amplitude of the histogram for nice plotting

        % Plots the inital image and histogram
        ScanZ                                       = max_z(j)-min_z(j);

        if plot_intermediate_steps == 1
            if index == 1 && j == channels_flatten(1)
                summary_fig                             = figure('Name','Summary', 'OuterPosition', [1 1 1800 800]);
            else
                close(summary_fig)
                summary_fig                             = figure('Name','Summary', 'OuterPosition', [1 1 1800 800]);
                %figure(summary_fig);
            end
            % summary_fig                                 = figure('OuterPosition', [1 1 1800 800]);
            subplot(2, 7, 1);
            plot_nice_image;                                                                                                                       % Plots the current image
            title('Initial image')
            subplot(2, 7, 8);
            hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
            title ('Initial histogram')
        end

        [height_ini position_ini]                   = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);                                                                 % Calculates the inital histogram
        hist_std_ini                                = std(position_ini);                                                                                                                % Computes the inital standard deviation for improvement check

Does an initial Median correction to fix basic line skips. Don't use on grids or other similar structures

        saved_image_no_med_corr                                           = images(:,:,j);    % Saves the image before the median correction in order to use it later (to compare global tilt with or without median corr)

        if Med_corr(j) == 1                                                                                                                                                       % Allows user to turn median correction on or off

            disp('--- First median correction ----')
            [images(:,:,j), hist_std, improv_param]                   = medianCorrection(images(:,:,j), hist_std_ini, 4, histogram_bins);

            if S.GUI_flag                       % program was called from the GUI
                if S.processAll  ==  0          % process only one image
                    med_corr_1_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
                end
            end


            if plot_intermediate_steps == 1
                figure(summary_fig)
                subplot(2, 7, 2);
                plot_nice_image;
                title('After 1st median correction')
                subplot(2, 7, 9);
                hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
                title('Histogram')
                if (improv_param == 1)
                       xlabel('This step improved the image');
                else
                       xlabel('This step did NOT improve the image');
                end

                % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' 1st med.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);

            end
        else
            [height_ini position_ini]                   = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);                                                                 % Calculates the inital histogram
            hist_std                                    = std(position_ini);
        end

        if vert_med_corr(j) == 1                                                                                                                                                   % Allows user to turn median correction on or off
             [images(:,:,j), hist_std, improv_param]                   = medianCorrection(images(:,:,j)', hist_std, 4, histogram_bins);
             images(:,:,j)                                             = images(:,:,j)';
        end

        final_flattening_image                                        = images(:,:,j);                                                                                                % Saves the image to do a direct final flattening with the final mask

Does an initial tilt correction to fix some of the initial sample tilt. Needed before doing any histogrammed flattening % Determines the number of coeffecients for the user selected polynomial order

        saved_image                                 = images(:,:,j);

        % FIRST order for polynomial fitting (1st order for both horizontal
        % and vertical fit, WITH median correction
        initial_poly_order_horizontal               = 1;
        initial_poly_order_vertical                 = 1;
        coef_length                                 = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);                                                      % Determines the number of coeffecients for the user selected polynomial order
        [xInput, yInput, zOutput]                   = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));                                              % Preps vectors for surface fitting
        poly_type                                   = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );                                     % Sets the fit type to be an X-Y polynomial of user choice
        poly_opts                                   = fitoptions( poly_type );                                                                                                          % Sets the fit  construction
        if index == 1
            poly_opts.Lower                         = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                         =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                          = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
            poly_opt.Upper                          = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
        end
        [poly_result]                               = fit( [xInput, yInput], zOutput, poly_type, poly_opts );                                                                           % Fits the polynomial to the actual data
        images(:,:,j)                               = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Init_poly_coef{index,j}                     = poly_result;

        saved_image_1st_order_fit_med_cor     = images(:,:,j);

        images(:,:,j)                               = saved_image_no_med_corr;

        % FIRST order for polynomial fitting (1st order for both horizontal
        % and vertical fit, WITHOUT median correction
        initial_poly_order_horizontal               = 1;
        initial_poly_order_vertical                 = 1;
        coef_length                                 = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);                                                      % Determines the number of coeffecients for the user selected polynomial order
        [xInput, yInput, zOutput]                   = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));                                              % Preps vectors for surface fitting
        poly_type                                   = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );                                     % Sets the fit type to be an X-Y polynomial of user choice
        poly_opts                                   = fitoptions( poly_type );                                                                                                          % Sets the fit  construction
        if index == 1
            poly_opts.Lower                         = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                         =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                          = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
            poly_opt.Upper                          = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
        end
        [poly_result]                               = fit( [xInput, yInput], zOutput, poly_type, poly_opts );                                                                           % Fits the polynomial to the actual data
        images(:,:,j)                               = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Init_poly_coef{index,j}                     = poly_result;

        saved_image_1st_order_fit                   = images(:,:,j);

%         comparison_1st_2nd_order_fit=figure('Name', 'Comparison between 1st and 2nd order polynomial fit');
%         subplot(2, 2, 1)
%         plot_nice_image;
%         title('After inital tilt correction - First order')
%         subplot(2, 2, 3);
%         hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
%         title('Histogram')

        images(:,:,j)                               = saved_image;

        % SECOND order for polynomial fitting (2nd order for both horizontal
        % and vertical fit, WITH median correction
        initial_poly_order_horizontal               = 2;
        initial_poly_order_vertical                 = 2;
        coef_length                                 = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);                                                      % Determines the number of coeffecients for the user selected polynomial order
        [xInput, yInput, zOutput]                   = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));                                              % Preps vectors for surface fitting
        poly_type                                   = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );                                     % Sets the fit type to be an X-Y polynomial of user choice
        poly_opts                                   = fitoptions( poly_type );                                                                                                          % Sets the fit  construction
        if index == 1
            poly_opts.Lower                         = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                         =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                          = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
            poly_opt.Upper                          = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
        end
        [poly_result]                               = fit( [xInput, yInput], zOutput, poly_type, poly_opts );                                                                           % Fits the polynomial to the actual data
        images(:,:,j)                               = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Init_poly_coef{index,j}                     = poly_result;

        saved_image_2nd_order_fit_med_cor           = images(:,:,j);

        images(:,:,j)                               = saved_image_no_med_corr;

        % SECOND order for polynomial fitting (2nd order for both horizontal
        % and vertical fit, WITHOUT median correction
        initial_poly_order_horizontal               = 2;
        initial_poly_order_vertical                 = 2;
        coef_length                                 = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);                                                      % Determines the number of coeffecients for the user selected polynomial order
        [xInput, yInput, zOutput]                   = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));                                              % Preps vectors for surface fitting
        poly_type                                   = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );                                     % Sets the fit type to be an X-Y polynomial of user choice
        poly_opts                                   = fitoptions( poly_type );                                                                                                          % Sets the fit  construction
        if index == 1
            poly_opts.Lower                         = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                         =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                          = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
            poly_opt.Upper                          = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
        end
        [poly_result]                               = fit( [xInput, yInput], zOutput, poly_type, poly_opts );                                                                           % Fits the polynomial to the actual data
        images(:,:,j)                               = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Init_poly_coef{index,j}                     = poly_result;

        saved_image_2nd_order_fit                   = images(:,:,j);

         % In order to really find the best tilt correction, not take the
         % std of the whole image but the part of the image between cdf =
         % 5% and 95%:
         tilt_corr_names = ['saved_image_1st_order_fit        '; ...
                            'saved_image_2nd_order_fit        '; ...
                            'saved_image_1st_order_fit_med_cor'; ...
                            'saved_image_2nd_order_fit_med_cor'];
%           std_plot = figure;
%           hold on;
%           colors_plot = ['b', 'g', 'r', 'm'];
         for ii = 1:4   % for the four different tilt corrections
             eval(['images(:,:,j) = ', num2str(tilt_corr_names(ii,:)), ';']);
             [hcdf xcdf] = ecdf(reshape(images(:,:,j),size(images,1)*size(images,2),1));  % computes cumulative distribution function
             k = 1;
             while hcdf(k) <= 0.01
                 k = k+1;
             end
             x_1percent = xcdf(k);  % finds the position for which cdf = 1%
             index_1percent = floor((x_1percent-min(min(images(:,:,j))))/(max(max(images(:,:,j)))-min(min(images(:,:,j)))) * histogram_bins);
             index_1percent = max([1,index_1percent]);

             k = length(hcdf);
              while hcdf(k) >= 0.99
                 k = k-1;
              end
             x_99percent = xcdf(k); % finds the position for which cdf = 99%
             index_99percent = floor((x_99percent-min(min(images(:,:,j))))/(max(max(images(:,:,j)))-min(min(images(:,:,j)))) * histogram_bins);
             index_99percent = min([index_99percent,length(xcdf)]);

             [height position]       = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);                                                                 % Calculates the histogram of the corrected height data
             height                  = smooth(height,smoothing_level,'loess');
             [peak_height peak_index]= findpeaks(height,'MinpeakDistance', min_peak_distance, 'Sortstr','descend');                                            % Finds up to four peaks in the corrected image to identify different terraces
             peak_index              = peak_index(1:min(length(peak_index),peaks_fit(j)));
             peak_height             = height(peak_index);
             tilt_hist_std(ii)       = std(position(index_1percent:index_99percent));
%              figure(std_plot)
%              plot(position(index_1percent:index_99percent), height(index_1percent:index_99percent), num2str(colors_plot(ii)))
         end
%          hold off;

        if plot_intermediate_steps == 1
         figure(summary_fig)
        end


          disp('--- Tilt correction ----')

         if (min(tilt_hist_std) > hist_std)             % if none of the different tilt corrections improved the image
             % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' tilt fail.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
             images(:,:,j) = saved_image;            % Cancels the tilt correction
            disp('Neither the 1st order nor the 2nd order inital tilt correction did improve the image')
            if plot_intermediate_steps == 1
                subplot(2, 7, 10);
                xlabel('This step did NOT improve the image')
            end
         else
           [minimum_std best_tilt_corr] = min(tilt_hist_std);
           images_size = size(saved_image_1st_order_fit, 1);
           all_tilt_corr_im = [saved_image_1st_order_fit; saved_image_2nd_order_fit; saved_image_1st_order_fit_med_cor; saved_image_2nd_order_fit_med_cor];
           images(:,:,j) = all_tilt_corr_im(images_size*(best_tilt_corr-1)+1:images_size*best_tilt_corr, :);      % Puts in images(:,:,j) the best image of the median/tilt corrections step
           titles = ['tilt 1st order without med cor'; 'tilt 2nd order without med cor'; 'tilt 1st order  with  med  cor'; 'tilt 2nd order  with  med  cor'];
           disp(['best tilt correction: ', titles(best_tilt_corr, :)])
           if plot_intermediate_steps == 1
               subplot(2, 7, 10);
               xlabel(['best corr: ', titles(best_tilt_corr, :)])
           end
           % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index}, ' ', titles(best_tilt_corr, :), '.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
         end

        if plot_intermediate_steps == 1
            figure(summary_fig)
            subplot(2, 7, 3);
            plot_nice_image;
            title('After inital tilt correction')
            subplot(2, 7, 10);
            hold on;
        end

        % This subtracts the median of the image from the image and
        % calculates the range of the image values so that all histograms have the same range
        images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));
        min_hist = min(min(images(:,:,j))) + 0.05*min(min(images(:,:,j)));
        max_hist = max(max(images(:,:,j))) - 0.05*max(max(images(:,:,j)));

       if S.GUI_flag                       % program was called from the GUI
            if S.processAll  ==  0          % process only one image
                ini_tilt_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
            end
      end

%        close(std_plot)

Generates initial fitting bound estimates

        image_inside_range                      = images(:,:,j);
        image_inside_range((images(:,:,j) < min_hist))   = NaN;                                                                                                                              % Thresholds out all of the data below
        image_inside_range((images(:,:,j) > max_hist))   = NaN;

        [height_est position_est]               = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);                                                                 % Calculates the histogram of the corrected height data (First image only)
        height_est                              = smooth(height_est,smoothing_level,'loess');

        [peak_height_est peak_index_est]        = findpeaks(height_est,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                    % Finds up to four peaks in the corrected image to identify different terraces

        if sort_peaks(j) == 1
            peak_index_est                      = sort(peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)])),'ascend');
        else
            peak_index_est                      = peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)]));
        end
        peak_height_est                         = height_est(peak_index_est);

        clear std_est
        if      length(peak_index_est) == 1
            std_est(1)                          = std(reshape(images(:,:,j),size(images,1)*size(images,2),1));
        elseif  length(peak_index_est) == 2
            std_est(1)                          = std(  images(images(:,:,j) <= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
            std_est(2)                          = std(  images(images(:,:,j) >= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
        elseif  length(peak_index_est) == 3
            std_est(1)                          = std(  images(images(:,:,j) <= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
            std_est(2)                          = std(  images(position_est(round((peak_index_est(1)+peak_index_est(2))/2)) <= images(:,:,j) <= position_est(round((peak_index_est(2)+peak_index_est(3))/2))));
            std_est(3)                          = std(  images(images(:,:,j) >= position_est(round((peak_index_est(2)+peak_index_est(3))/2))));
        else
        end

       if Calc_flatten_lim(j) == 1
            max_spacing_stdev(j)                    = 1.2*max(std_est);
            min_spacing_stdev(j)                    = 0.01*min(std_est);
        else
            min_spacing_stdev                       = min_spacing_stdev_user;
            max_spacing_stdev                       = max_spacing_stdev_user;
       end


         if plot_intermediate_steps == 1
            figure(summary_fig)
            subplot(2, 7, 10);
            title('Histogram')
            image_inside_range                      = images(:,:,j);
            image_inside_range((images(:,:,j) < min_hist))   = NaN;
            image_inside_range((images(:,:,j) > max_hist))   = NaN;
            hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
         end

        if plot_intermediate_image == 1
            save_intermediate_image;
        end

%         min_peak_distance                           = floor(max(std_est) * histogram_bins / (max(position_est)-min(position_est)));%floor((max(std_est)+min(std_est))/2 * histogram_bins / (max(position_est)-min(position_est)))              % Estimates the minimum distance between peaks
%         if min_peak_distance < 1
%             min_peak_distance = 1;
%         end
%
        [peak_height_est peak_index_est]        = findpeaks(height_est,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                    % Finds up to four peaks in the corrected image to identify different terraces
        peak_index_est                      = peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)]));

        ref_peak_pos = position_est(peak_index_est);        % Peaks reference positions to find corresponding peaks
%         min_peak_distance                           = floor((max(std_est)+min(std_est))/2 * histogram_bins / (max(position_est)-min(position_est))); % = 2;            % Sets the minimum distance between peaks to a minimum now that the reference peaks positions have been found
%          if min_peak_distance < 1
%             min_peak_distance = 1;
%          end

         if plot_intermediate_steps == 1
            figure(summary_fig)
            subplot(2, 7, 10)
            plot(position_est(peak_index_est), height_est(peak_index_est), 'k^', 'MarkerEdgeColor','k', 'MarkerFaceColor', 'g')
            hold off;
         end

This section histograms the image, identifies up to 3 distinct levels in the historam, fits them to Gaussians, and uses this to determin an approriate threshold range to use for a masked polynomial approach. This leads to truly flat planes that are not biased by steps and terraces in the image data

        saved_image                                 = images(:,:,j);                                                                                                                % Saves the image to restore it if this step doesn't improve it

        image_inside_range                      = images(:,:,j);
        image_inside_range((images(:,:,j) < min_hist))   = NaN;
        image_inside_range((images(:,:,j) > max_hist))   = NaN;

        [height_1 position_1]                       = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);                         % Calculates the histogram of the corrected height data (First image only)
        height_1                                    = smooth(height_1,smoothing_level,'loess');
        [peak_height_1 peak_index_1]                = findpeaks(height_1,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                  % Finds up to four peaks in the corrected image to identify different terraces

        peak_index_1                                = find_closest(ref_peak_pos, peak_index_1(1:min(7, length(peak_index_1))), position_1); % peak_index_1(1:...): to keep only the 5 higher peaks

        if sort_peaks(j) == 1
            peak_index_1 = sort(peak_index_1(1:min(peaks_fit(j), length(peak_index_1))),'ascend');
        else
            peak_index_1 = peak_index_1(1:min(peaks_fit(j), length(peak_index_1)));
        end
        peak_height_1 = height_1(peak_index_1);

        func_form                                   = [];                                                                                                                               % Initializes the functional form to fit
        t_hist_lower_a                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the lower bounds for fitting
        t_hist_lower_s                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the lower bounds for fitting
        t_hist_lower_x                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the lower bounds for fitting

        t_hist_upper_a                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the upper bounds for fitting
        t_hist_upper_s                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the upper bounds for fitting
        t_hist_upper_x                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the upper bounds for fitting

        t_hist_start_a                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the start point for fitting
        t_hist_start_s                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the start point for fitting
        t_hist_start_x                              = zeros(1,min(3,length(peak_height_1)));                                                                                            % Initializes the start point for fitting

        for i                                       = 1:min(3,length(peak_height_1))                                                                                                    % Loops for the minimum of the number of peaks or 3
            t_hist_lower_a(i)                       = (1-amplitude_fract)*peak_height_1(i);                                                                                             % Sets the lower amplitude bound for the fit
            t_hist_lower_s(i)                       = min_spacing_stdev(j);                                                                                                             % Sets the lower standard deviation
            t_hist_lower_x(i)                       = position_1(peak_index_1(i))-3.5/histogram_bins;                                                                                              % Sets the lower center relative to the peak

            t_hist_upper_a(i)                       = (1+amplitude_fract)*peak_height_1(i);                                                                                             % Sets the upper amplitude bound for the fit
            t_hist_upper_s(i)                       = max_spacing_stdev(j);                                                                                                             % Sets the upper standard deviation
            t_hist_upper_x(i)                       = position_1(peak_index_1(i))+3.5/histogram_bins;                                                                                              % Sets the upper center relative to the peak

            t_hist_start_a(i)                       = 1.00*peak_height_1(i);                                                                                                            % Sets the starting amplitude bound for the fit
            t_hist_start_s(i)                       = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;                                                                                    % Sets the starting standard deviation
            t_hist_start_x(i)                       = position_1(peak_index_1(i));                                                                                                      % Sets the starting center relative at the peak
            func_form                               = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];                                                % Generates the functional form to fit
            if i < min(3,length(peak_height_1))
                func_form                           = [func_form,'+'];                                                                                                                  % Generates the functional form to fit
            end
        end

        %func_form2 = ['gauss', num2str(min(3,length(peak_height_1)))];

        % Combines all the bouds into a form for the fits
        t_hist_lower                                = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];                                                                                   % Combines the lower bounds
        t_hist_upper                                = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];                                                                                   % Combines the upper bounds
        t_hist_start                                = [t_hist_start_a,t_hist_start_s,t_hist_start_x];                                                                                   % Combines the start point

        % Fits the data with the speficied conditions
        t_hist = fitoptions('Method','NonlinearLeastSquares',...                                                                                                                        % Sets the conditions for fitting
                   'Lower',      t_hist_lower,...
                   'Upper',      t_hist_upper,...
                   'Startpoint', t_hist_start,...
                   'MaxIter',    1e2,...
                   'MaxFunEvals',3e2);
        f_Number_hist = fittype(func_form,  'options' ,t_hist);
        %f_Number_hist = fittype(func_form2);
        [peak_1_gaussian,peak_1_gaussian_gof]       = fit( position_1',...                                                                                                              % Fits the data
                                                            height_1,...
                                                            f_Number_hist);

        peaks_nb = min(3,length(peak_height_1));
         for i                                   = 1:min(3,length(peak_height_1))       % Detects peaks smaller than min_peak_height; peaks found by findpeaks are bigger than this value, but if gaussian fit not good, peaks smaller can 'apppear'
             eval(['peak_ampl = peak_1_gaussian.a', num2str(i), ';']);
             if peak_ampl < min_peak_height
                 peaks_nb = peaks_nb - 1;
             end
         end

        % Combines all the amplitudes, the positions and the standard
        % deviations into one vector and sets the lower and upper threshold
        % for all peaks
         peak_1_gaussian_ampl = [];
         peak_1_gaussian_pos = [];
         peak_1_gaussian_std = [];
         for i = 1:peaks_nb
            eval(['peak_1_gaussian_ampl(', num2str(i), ')= peak_1_gaussian.a', num2str(i), ';']);
            eval(['peak_1_gaussian_pos(', num2str(i), ')= peak_1_gaussian.x', num2str(i), ';']);
            eval(['peak_1_gaussian_std(', num2str(i), ')= peak_1_gaussian.s', num2str(i), ';']);
         end

         % For the overlap checking, sort the peaks in ascending order (in
         % term of position)
         if peaks_nb == 1
             peak_gaussian_pos = peak_1_gaussian_pos;
             peak_gaussian_std = peak_1_gaussian_std;
             peak_gaussian_ampl = peak_1_gaussian_ampl;
         else
             [peak_gaussian_pos, ind] = sort(peak_1_gaussian_pos, 'ascend');
             for i = 1:peaks_nb
                peak_gaussian_std(i) = peak_1_gaussian_std(ind(i));
                peak_gaussian_ampl(i) = peak_1_gaussian_ampl(ind(i));
             end
         end

           % For each peak, calculates the lower and upper threshold
           for i = 1:peaks_nb
            eval(['thresh_low_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') - thresh_range*peak_gaussian_std(', num2str(i), ');']);
            eval(['thresh_high_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') + thresh_range*peak_gaussian_std(', num2str(i), ');']);
           end

%          disp(['thresh_low_1 = ', num2str(thresh_low_1), '; thresh_high_1 = ', num2str(thresh_high_1)]);
%          if peaks_nb  > 1
%          disp(['thresh_low_2 = ', num2str(thresh_low_2), '; thresh_high_2 = ', num2str(thresh_high_2)]);
%          end
%          if peaks_nb > 2
%             disp(['thresh_low_3 = ', num2str(thresh_low_3), '; thresh_high_3 = ', num2str(thresh_high_3)]);
%          end

         % For each pair of peaks, look if they overlap and what is their
         % intersection point
         if peaks_nb  == 1                                                                                                                                              % Only one peak, so no problem of overlap
             % do nothing
         else
             for i = 1:peaks_nb-1
                 % Does some plotting
%                  x_plot = linspace(peak_gaussian_pos(1)-5*peak_gaussian_std(1), peak_gaussian_pos(peaks_nb)+5*peak_gaussian_std(peaks_nb), 500);
%                  figure; plot(x_plot, peak_gaussian_ampl(i)*exp(-(x_plot-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2)), 'b', x_plot, peak_gaussian_ampl(i+1)*exp(-(x_plot-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)), 'c');
%                  hold on;
%                  curr_thresh_low = peak_gaussian_pos(i)-thresh_range*peak_gaussian_std(i);
%                  plot(curr_thresh_low, 1:peak_gaussian_ampl(i), 'b');
%                  curr_thresh_high = peak_gaussian_pos(i)+thresh_range*peak_gaussian_std(i);
%                  plot(curr_thresh_high, 1:peak_gaussian_ampl(i), 'b');
%                  curr_thresh_low = peak_gaussian_pos(i+1)-thresh_range*peak_gaussian_std(i+1);
%                  plot(curr_thresh_low,1:peak_gaussian_ampl(i+1), 'c');
%                  curr_thresh_high = peak_gaussian_pos(i+1)+thresh_range*peak_gaussian_std(i+1);
%                  plot(curr_thresh_high, 1:peak_gaussian_ampl(i+1), 'c');

                 fct = @(x)(peak_gaussian_ampl(i)*exp(-(x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2))- peak_gaussian_ampl(i+1)*exp(-(x-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)));
                 eval(['thresh_start_x = (thresh_low_', num2str(i), '+ thresh_high_', num2str(i+1), ')/2;']);
                [intersect_x, fval, exit_flag] = fzero(fct, thresh_start_x);                                                                                                         % Finds intersection of the two peaks
%                 plot(intersect_x, peak_gaussian_ampl(i)*exp(-(intersect_x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2)), 'rx', 'MarkerSize', 12);
%                 hold off;
                if intersect_x > peak_gaussian_pos(i)+ 5*peak_gaussian_std(i) || intersect_x < peak_gaussian_pos(i) - 5*peak_gaussian_std(i)        % Wrong zero was found; it can happen when the intersection point is where the gaussian curves tend to zero
                exit_flag = 0;                                                                                                                      % Sets the exit_flag to a value different of 1 so it doesn't do anything
                end
                    if exit_flag == 1                                                                                                                                                         % An intersection point between the two peaks was found
                    if i == 1
                         if thresh_high_1 > intersect_x
                            thresh_high_1 = intersect_x;                                                                                                            % Sets the upper threshold of the first peak to the intersection position
                         end
                         if thresh_low_2  < intersect_x
                            thresh_low_2 = intersect_x;                                                                                                              % Sets the lower threshold of the second peak to the intersection position
                         end
                    end
                    if i == 2       % (i can be equal to 1 or 2)
                         if thresh_high_2 > intersect_x
                            thresh_high_2 = intersect_x;                                                                                                            % Sets the upper threshold of the first peak to the intersection position
                         end
                         if thresh_low_3  < intersect_x
                            thresh_low_3 = intersect_x;                                                                                                              % Sets the lower threshold of the second peak to the intersection position
                         end
                    end
                end
             end
         end

         % Put the peaks in their originial order (largest amplitude first)
          [peak_gaussian_ampl, ind] = sort(peak_gaussian_ampl, 'ascend');
           for i = 1:peaks_nb
                eval(['thresh_low(', num2str(i), ') = thresh_low_', num2str(ind(i)), ';']);
                eval(['thresh_high(', num2str(i), ') = thresh_high_', num2str(ind(i)), ';']);
           end

         % Set the lower and upper threshold of the peak to be flattened
         if peaks_nb >= peak_to_flatten(j)
             eval(['thresh_low = thresh_low_', num2str(peak_to_flatten(j)), ';']);
             eval(['thresh_high = thresh_high_', num2str(peak_to_flatten(j)), ';']);
         else
             thresh_low = thresh_low_1;
             thresh_high = thresh_high_1;
         end

%           disp(['new_thresh_low_1 = ', num2str(thresh_low_1), '; thresh_high_1 = ', num2str(thresh_high_1)]);
%          if peaks_nb > 1
%          disp(['new_thresh_low_2 = ', num2str(thresh_low_2), '; thresh_high_2 = ', num2str(thresh_high_2)]);
%          end
%          if peaks_nb > 2
%             disp(['new_thresh_low_3 = ', num2str(thresh_low_3), '; thresh_high_3 = ', num2str(thresh_high_3)]);
%          end

%         thresh_low_1                                  = peak_1_gaussian.x1 - thresh_range*peak_1_gaussian.s1;                                                                                        % Sets the lower threshold for the next round of flattening
%         thresh_high_1                                = peak_1_gaussian.x1 + thresh_range*peak_1_gaussian.s1;                                                                                        % Sets the upper threshold for the next round of flattening
        image_test                                  = images(:,:,j);                                                                                                                    % Duplicates the images as a placeholder for the thresholding
        image_test((images(:,:,j) < thresh_low))    = NaN;                                                                                                                              % Thresholds out all of the data below
        image_test((images(:,:,j) > thresh_high))   = NaN;                                                                                                                              % Thresholds out all of the data above


        First_Thresh_File = [ImageDir,'Matlab',filesep,'first_threshold.txt'];

    if  index == 1 && j == 1
        fid = fopen(First_Thresh_File, 'a');
        fprintf(fid, [  'Image','\t',...
                        'Channel','\t',...
                        'Thresh Low','\t',...
                        'Thresh High','\n'],...
                        'char');

        fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                        channel_info(j).Name,'\t',...
                        num2str(thresh_low,'%0.0f'),'\t',...
                        num2str(thresh_high,'%0.3f'),'\n'],...
                        'char');
        fclose(fid);
    else
        fid = fopen(First_Thresh_File, 'a');
        fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                        channel_info(j).Name,'\t',...
                        num2str(thresh_low,'%0.0f'),'\t',...
                        num2str(thresh_high,'%0.3f'),'\n'],...
                        'char');
        fclose(fid);
    end

        % Fit the thresholded surface to an n x m polynomial
        % Loops for the number of images
        coef_length                                 = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);                                                                      % Determines the number of coeffecients for the user selected polynomial order
        [xInput_poly, yInput_poly, zOutput_poly]    = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);                                                 % Preps vectors for surface fitting
        poly_type                                   = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );                                                     % Sets the fit type to be an X-Y polynomial of user choice
        poly_opts                                   = fitoptions( poly_type );                                                                                                          % Sets the fit  construction

        if index == 1
            poly_opts.Lower                         = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                         =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                          = coeffvalues(Second_poly_coef{index-1,j})-0.5*abs(coeffvalues(Second_poly_coef{index-1,j}));
            poly_opt.Upper                          = coeffvalues(Second_poly_coef{index-1,j})+0.5*abs(coeffvalues(Second_poly_coef{index-1,j}));
        end                                                                                                       % Sets the upper bounds

        poly_result                                 = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );                                                            % Fits the polynomial to the actual data
        images(:,:,j)                               = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Second_poly_coef{index,j}                   = poly_result;

        images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));

        if S.GUI_flag                       % program was called from the GUI
            if S.processAll  ==  0          % process only one image
                poly_flat_1_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
            end
        end

        if plot_intermediate_steps == 1
            figure(summary_fig)
            subplot(2, 7, 4);
            plot_nice_image;
            title('After 1st polynomial flattening')
            subplot(2, 7, 11);
            image_inside_range                      = images(:,:,j);
            image_inside_range((images(:,:,j) < min_hist))   = NaN;
            image_inside_range((images(:,:,j) > max_hist))   = NaN;
            hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
            plot_peaks;
            title('Histogram')
            xlabel('This step improved the image');
        end

        % Checks if the polynomial flattening brought some improvement or
        % not
        disp('--- First Poly Flattening ----')
        [improv_param, peak_2_gaussian_std, peak_2_gaussian_ampl] = check_improvement(images(:,:,j), ...
                                                                    peak_1_gaussian_std, peak_1_gaussian_ampl,...
                                                                    min_spacing_stdev(j), max_spacing_stdev(j), ...
                                                                    amplitude_fract,  min_peak_distance, min_peak_height, peak_threshold, ...
                                                                    position_est, peak_height_est, peak_index_est, ...
                                                                    peak_pos_tolerance, histogram_bins, smoothing_level,...
                                                                    min_hist, max_hist ,ref_peak_pos, peaks_fit(j));
        if max(peak_gaussian_ampl) >  max_peak_ampl
            max_peak_ampl = max(peak_gaussian_ampl);
        end
        if improv_param < 0                         % If true, the polynomial flattening didn't bring an improvement
            % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' poly 1 fail.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
            images(:,:,j) = saved_image;            % Cancels the polynomial flattening
            peak_2_gaussian_std = peak_1_gaussian_std;
            peak_2_gaussian_ampl = peak_1_gaussian_ampl;
            disp('The first polymomial flattening did NOT improve the image')
            if plot_intermediate_steps == 1
                xlabel('This step did NOT improve the image');
            end
        else
             % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' poly 1 win.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
        end
         if improv_param == 0
            if plot_intermediate_steps == 1
               xlabel('No change')
            end
         end

          [height position]                           = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
%           min_peak_distance                           = floor(min(peak_2_gaussian_std) * histogram_bins / (max(position)-min(position)))              % Readjusts the value of min_peak_distance
%           if (min_peak_distance < 4)
%               min_peak_distance                       = 4;
%           end
%

%             images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));
%             subplot(2, 7, 4);
%             plot_nice_image;
%             title('After 1st polynomial flattening')
%

    % Refine reference peaks:
        image_inside_range                      = images(:,:,j);
        image_inside_range((images(:,:,j) < min_hist))   = NaN;
        image_inside_range((images(:,:,j) > max_hist))   = NaN;

        [height position]                    = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);                         % Calculates the histogram of the corrected height data (First image only)
        height                               = smooth(height,smoothing_level,'loess');
        [peak_height peak_index]             = findpeaks(height,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                  % Finds up to four peaks in the corrected image to identify different terraces

        peak_index                           = find_closest(ref_peak_pos, peak_index(1:min(7, length(peak_index))), position); % peak_index_1(1:...): to keep only the 5 higher peaks

        ref_peak_pos = position(peak_index);

After the first round of thresholded plane correction this does another median correction, respecting the above thresholds to correct scaring

        if Med_corr_2(j) == 1

            saved_image(:,:) = images(:,:,j);

            disp('--- Second median correction ----')
            [images(:,:,j), peak_gaussian_std, peak_gaussian_ampl, improv_param]       = medianCorrection2(images(:,:,j), 2.5, ...
                                                                                        peak_2_gaussian_std, peak_2_gaussian_ampl, ...
                                                                                        min_spacing_stdev(j), max_spacing_stdev(j), ...
                                                                                        amplitude_fract, min_peak_distance, min_peak_height, peak_threshold, ...
                                                                                        position_est, peak_height_est, peak_index_est, ...
                                                                                        peak_pos_tolerance, histogram_bins, smoothing_level, ...
                                                                                        min_hist, max_hist, ref_peak_pos, peaks_fit(j));  % Med_corr_2

            images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));

            if S.GUI_flag                       % program was called from the GUI
                if S.processAll  ==  0          % process only one image
                    med_corr_2_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
                end
            end


            if plot_intermediate_steps == 1
                figure(summary_fig)
                subplot(2, 7, 5);
                plot_nice_image;
                title('After 2nd median correction')
                subplot(2, 7, 12);
                image_inside_range                      = images(:,:,j);
                image_inside_range((images(:,:,j) < min_hist))   = NaN;
                image_inside_range((images(:,:,j) > max_hist))   = NaN;
                hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
                plot_peaks;
                title('Histogram')
            end

             if max(peak_gaussian_ampl) >  max_peak_ampl
            max_peak_ampl = max(peak_gaussian_ampl);
             end

            if improv_param < 0
                % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' 2 med fail.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
                images(:,:,j)         = saved_image(:,:);
                peak_gaussian_std = peak_2_gaussian_std;
                peak_gaussian_ampl = peak_2_gaussian_ampl;
                if plot_intermediate_steps == 1
                    xlabel('This step did NOT improve the image');
                end
            elseif improv_param == 0
                if plot_intermediate_steps == 1
                    xlabel('No change')
                end
            else
                if plot_intermediate_steps == 1
                    xlabel('This step improved the image');
                end

                 % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' 2 med improv.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);

            end

        end


          [height position]                           = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);

This does a second round of thresholding and polynomial flattening to improve on the results from before.

        saved_image                                 = images(:,:,j);

         image_inside_range                      = images(:,:,j);
         image_inside_range((images(:,:,j) < min_hist))   = NaN;
         image_inside_range((images(:,:,j) > max_hist))   = NaN;

        [height_final position_final]               =  hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);       % Calculates the histogram of the corrected height data (First image only)
        height_final                                = smooth(height_final,smoothing_level,'loess');
        [peak_height_final peak_index_final]        = findpeaks(height_final,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                % Finds peaks in the corrected image to identify different terraces

        peak_index_final                            = find_closest(ref_peak_pos, peak_index_final(1:min(peaks_fit(j)+1, length(peak_index_final))), position_final); % peak_index_final(1:...) : to keep only the 5 higher peaks

        ref_peak_pos = position_final(peak_index_final);        % Refine peaks reference positions

        if sort_peaks(j) == 1
            peak_index_final = sort(peak_index_final(1:min(peaks_fit(j), length(peak_index_final))),'ascend');
        else
            peak_index_final = peak_index_final(1:min(peaks_fit(j), length(peak_index_final)));
        end
        peak_height_final = height_final(peak_index_final);

        func_form                                   = [];
        t_hist_lower_a                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the lower bounds for fitting
        t_hist_lower_s                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the lower bounds for fitting
        t_hist_lower_x                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the lower bounds for fitting
        t_hist_upper_a                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the upper bounds for fitting
        t_hist_upper_s                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the upper bounds for fitting
        t_hist_upper_x                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the upper bounds for fitting
        t_hist_start_a                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the start point for fitting
        t_hist_start_s                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the start point for fitting
        t_hist_start_x                              = zeros(1,min(3,length(peak_height_final)));                                                                                        % Initializes the start point for fitting

        for i                                       = 1:min(3,length(peak_height_final))                                                                                                % Loops for the minimum of the number of peaks or 3
            t_hist_lower_a(i)                       = (1-amplitude_fract)*peak_height_final(i);                                                                                         % Sets the lower amplitude bound for the fit
            t_hist_lower_s(i)                       = min_spacing_stdev(j);                                                                                                             % Sets the lower standard deviation
            t_hist_lower_x(i)                       = position_final(peak_index_final(i))-3.5/histogram_bins;                                                                                      % Sets the lower center relative to the peak

            t_hist_upper_a(i)                       = (1+amplitude_fract)*peak_height_final(i);                                                                                         % Sets the upper amplitude bound for the fit
            t_hist_upper_s(i)                       = max_spacing_stdev(j);                                                                                                             % Sets the upper standard deviation
            t_hist_upper_x(i)                       = position_final(peak_index_final(i))+3.5/histogram_bins;                                                                                      % Sets the upper center relative to the peak

            t_hist_start_a(i)                       = 1.00*peak_height_final(i);                                                                                                        % Sets the starting amplitude bound for the fit
            t_hist_start_s(i)                       = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;                                                                                    % Sets the starting standard deviation
            t_hist_start_x(i)                       = position_final(peak_index_final(i));                                                                                              % Sets the starting center relative at the peak
            func_form                               = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];                                                % Generates the functional form to fit

            if i < min(3,length(peak_height_final))
                func_form                           = [func_form,'+'];                                                                                                                  % Generates the functional form to fit
            end
        end

        clear peak_gaussian_std peak_gaussian_ampl peak_gaussian_pos;

        % Combines all the bouds into a form for the fits
        t_hist_lower                                = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];                                                                                   % Combines the lower bounds
        t_hist_upper                                = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];                                                                                   % Combines the upper bounds
        t_hist_start                                = [t_hist_start_a,t_hist_start_s,t_hist_start_x];                                                                                   % Combines the start point

        % Fits the data with the speficied conditions
        t_hist = fitoptions('Method','NonlinearLeastSquares',...                                                                                                                        % Sets the conditions for fitting
                   'Lower',      t_hist_lower,...
                   'Upper',      t_hist_upper,...
                   'Startpoint', t_hist_start,...
                   'MaxIter',    1e2,...
                   'MaxFunEvals',3e2);
        f_Number_hist                               = fittype(func_form, 'options' ,t_hist);
        [peak_final_gaussian,peak_final_gaussian_gof] = fit( position_final',...
                                                             height_final,...
                                                             f_Number_hist);

         peaks_nb = min(3,length(peak_height_final));
         for i                                   = 1:min(3,length(peak_height_final))       % Detects peaks smaller than min_peak_height; peaks found by findpeaks are bigger than this value, but if gaussian fit not good, peaks smaller can 'apppear'
             eval(['peak_ampl = peak_final_gaussian.a', num2str(i), ';']);
             if peak_ampl < min_peak_height
                 peaks_nb = peaks_nb - 1;
             end
         end

        % Combines all the amplitudes, the positions and the standard
        % deviations into one vector and sets the lower and upper threshold
        % for all peaks
         peak_final_gaussian_ampl = [];
         peak_final_gaussian_pos = [];
         peak_final_gaussian_std = [];
         for i = 1:peaks_nb
            eval(['peak_final_gaussian_ampl(', num2str(i), ')= peak_final_gaussian.a', num2str(i), ';']);
            eval(['peak_final_gaussian_pos(', num2str(i), ')= peak_final_gaussian.x', num2str(i), ';']);
            eval(['peak_final_gaussian_std(', num2str(i), ')= peak_final_gaussian.s', num2str(i), ';']);
         end

         % For the overlap checking, sort the peaks in ascending order (in
         % term of position)
         if peaks_nb == 1
             peak_gaussian_pos = peak_final_gaussian_pos;
             peak_gaussian_std = peak_final_gaussian_std;
             peak_gaussian_ampl = peak_final_gaussian_ampl;
         else
             [peak_gaussian_pos, ind] = sort(peak_final_gaussian_pos, 'ascend');
             for i = 1:peaks_nb
                peak_gaussian_std(i) = peak_final_gaussian_std(ind(i));
                peak_gaussian_ampl(i) = peak_final_gaussian_ampl(ind(i));
             end
         end

           % For each peak, calculates the lower and upper threshold
           for i = 1:peaks_nb
            eval(['thresh_low_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') - thresh_range*peak_gaussian_std(', num2str(i), ');']);
            eval(['thresh_high_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') + thresh_range*peak_gaussian_std(', num2str(i), ');']);
           end

%          disp(['thresh_low_1 = ', num2str(thresh_low_1), '; thresh_high_1 = ', num2str(thresh_high_1)]);
%          if peaks_nb  > 1
%          disp(['thresh_low_2 = ', num2str(thresh_low_2), '; thresh_high_2 = ', num2str(thresh_high_2)]);
%          end
%          if peaks_nb > 2
%             disp(['thresh_low_3 = ', num2str(thresh_low_3), '; thresh_high_3 = ', num2str(thresh_high_3)]);
%          end

         % For each pair of peaks, look if they overlap and what is their
         % intersection point
         if peaks_nb  == 1                                                                                                                                              % Only one peak, so no problem of overlap
             % do nothing
         else
             for i = 1:peaks_nb-1
                 % Does some plotting
%                  x_plot = linspace(peak_gaussian_pos(1)-5*peak_gaussian_std(1), peak_gaussian_pos(peaks_nb)+5*peak_gaussian_std(peaks_nb), 500);
%                  figure; plot(x_plot, peak_gaussian_ampl(i)*exp(-(x_plot-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2)), 'b', x_plot, peak_gaussian_ampl(i+1)*exp(-(x_plot-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)), 'c');
%                  hold on;
%                  curr_thresh_low = peak_gaussian_pos(i)-thresh_range*peak_gaussian_std(i);
%                  plot(curr_thresh_low, 1:peak_gaussian_ampl(i), 'b');
%                  curr_thresh_high = peak_gaussian_pos(i)+thresh_range*peak_gaussian_std(i);
%                  plot(curr_thresh_high, 1:peak_gaussian_ampl(i), 'b');
%                  curr_thresh_low = peak_gaussian_pos(i+1)-thresh_range*peak_gaussian_std(i+1);
%                  plot(curr_thresh_low,1:peak_gaussian_ampl(i+1), 'c');
%                  curr_thresh_high = peak_gaussian_pos(i+1)+thresh_range*peak_gaussian_std(i+1);
%                  plot(curr_thresh_high, 1:peak_gaussian_ampl(i+1), 'c');

                 fct = @(x)(peak_gaussian_ampl(i)*exp(-(x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2))- peak_gaussian_ampl(i+1)*exp(-(x-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)));
                 eval(['thresh_start_x = (thresh_low_', num2str(i), '+ thresh_high_', num2str(i+1), ')/2;']);
                [intersect_x, fval, exit_flag] = fzero(fct, thresh_start_x);                                                                                                         % Finds intersection of the two peaks
%                 plot(intersect_x, peak_gaussian_ampl(i)*exp(-(intersect_x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2)), 'rx', 'MarkerSize', 12);
%                 hold off;
                if intersect_x > peak_gaussian_pos(i)+ 5*peak_gaussian_std(i) || intersect_x < peak_gaussian_pos(i) - 5*peak_gaussian_std(i)        % Wrong zero was found; it can happen when the intersection point is where the gaussian curves tend to zero
                exit_flag = 0;                                                                                                                      % Sets the exit_flag to a value different of 1 so it doesn't do anything
                end
                    if exit_flag == 1                                                                                                                                                         % An intersection point between the two peaks was found
                        if i == 1
                             if thresh_high_1 > intersect_x
                                thresh_high_1 = intersect_x;                                                                                                            % Sets the upper threshold of the first peak to the intersection position
                             end
                             if thresh_low_2  < intersect_x
                                thresh_low_2 = intersect_x;                                                                                                              % Sets the lower threshold of the second peak to the intersection position
                             end
                        end
                         if i == 2       % (i can be equal to 1 or 2)
                             if thresh_high_2 > intersect_x
                                thresh_high_2 = intersect_x;                                                                                                            % Sets the upper threshold of the first peak to the intersection position
                             end
                             if thresh_low_3  < intersect_x
                                thresh_low_3 = intersect_x;                                                                                                              % Sets the lower threshold of the second peak to the intersection position
                             end
                         end
                    end
             end
         end

         % Put the peaks in their originial order (largest amplitude first)
          [peak_gaussian_ampl, ind] = sort(peak_gaussian_ampl, 'ascend');
           for i = 1:peaks_nb
                eval(['thresh_low(', num2str(i), ') = thresh_low_', num2str(ind(i)), ';']);
                eval(['thresh_high(', num2str(i), ') = thresh_high_', num2str(ind(i)), ';']);
           end

         % Set the lower and upper threshold of the peak to be flattened
         if peaks_nb >= peak_to_flatten(j)
             eval(['thresh_low = thresh_low_', num2str(peak_to_flatten(j)), ';']);
             eval(['thresh_high = thresh_high_', num2str(peak_to_flatten(j)), ';']);
         else
             thresh_low = thresh_low_1;
             thresh_high = thresh_high_1;
         end

%           disp(['new_thresh_low_1 = ', num2str(thresh_low_1), '; thresh_high_1 = ', num2str(thresh_high_1)]);
%          if peaks_nb > 1
%          disp(['new_thresh_low_2 = ', num2str(thresh_low_2), '; thresh_high_2 = ', num2str(thresh_high_2)]);
%          end
%          if peaks_nb > 2
%             disp(['new_thresh_low_3 = ', num2str(thresh_low_3), '; thresh_high_3 = ', num2str(thresh_high_3)]);
%          end

%         thresh_low                                  = peak_final_gaussian.x1 - thresh_range*peak_final_gaussian.s1;                                                                     % Sets the lower threshold for the next round of flattening
%         thresh_high                                 = peak_final_gaussian.x1 + thresh_range*peak_final_gaussian.s1;
        image_test                                  = images(:,:,j);                                                                                                                    % Makes a duplicate of the image
        image_test((images(:,:,j) < thresh_low))    = NaN;                                                                                                                              % Thresholds data below
        image_test((images(:,:,j) > thresh_high))   = NaN;                                                                                                                              % Thresholds data above

        Second_Thresh_File = [ImageDir,'Matlab',filesep,'second_threshold.txt'];
        if  index == 1 && j == 1
            fid = fopen(Second_Thresh_File, 'a');
            fprintf(fid, [  'Image','\t',...
                            'Channel','\t',...
                            'Thresh Low','\t',...
                            'Thresh High','\n'],...
                            'char');

            fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                            channel_info(j).Name,'\t',...
                            num2str(thresh_low,'%0.0f'),'\t',...
                            num2str(thresh_high,'%0.3f'),'\n'],...
                            'char');
            fclose(fid);
        else
            fid = fopen(Second_Thresh_File, 'a');
            fprintf(fid, [  num2str(file_number,'%0.0f'),'\t',...
                            channel_info(j).Name,'\t',...
                            num2str(thresh_low,'%0.0f'),'\t',...
                            num2str(thresh_high,'%0.3f'),'\n'],...
                            'char');
            fclose(fid);
        end


        % Fit the thresholded surface to an n x m polynomial
        coef_length                                     = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);                                                                      % Calculates the number of coeffecients
        [xInput_poly, yInput_poly, zOutput_poly]        = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);                                                 % Preps vectors for surface fitting
        poly_type                                       = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );                                                     % Sets the fit type to be an X-Y polynomial
        poly_opts                                       = fitoptions( poly_type );                                                                                                          % Tells it to fit a polynomial
        if index == 1
            poly_opts.Lower                             = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
            poly_opts.Upper                             =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
        else
            poly_opt.Lower                              = coeffvalues(Third_poly_coef{index-1,j})-0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
            poly_opt.Upper                              = coeffvalues(Third_poly_coef{index-1,j})+0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
        end                                                                                                       % Sets the upper bounds

        poly_result                                     = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );                                                            % Fits the polynomial to the actual data
        images(:,:,j)                                   = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
        Third_poly_coef{index,j}                        = poly_result;


        images(:,:,j) = images(:,:,j) - median(median(images(:,:,j)));

        if S.GUI_flag                       % program was called from the GUI
            if S.processAll  ==  0          % process only one image
                poly_flat_2_images(:,:,channels_indices(j)) = images(:,:,j); % Saves the raw image
            end
        end

        if plot_intermediate_steps == 1
            figure(summary_fig)
            subplot(2, 7, 6);
            plot_nice_image;
            title('After 2nd polynomial flattening')
            subplot(2, 7, 13);
            image_inside_range                      = images(:,:,j);
            image_inside_range((images(:,:,j) < min_hist))   = NaN;
            image_inside_range((images(:,:,j) > max_hist))   = NaN;
            hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
            plot_peaks;
            title('Histogram')
            xlabel('This step improved the image');
        end


        % Checks if the polynomial flattening brought some improvement or
        % not
        disp('--- Second Poly Flattening ----')
        [improv_param, peak_gaussian_std, peak_gaussian_ampl] = check_improvement(images(:,:,j), ...
                                                                    peak_final_gaussian_std, peak_final_gaussian_ampl,...
                                                                    min_spacing_stdev(j), max_spacing_stdev(j), ...
                                                                    amplitude_fract,  min_peak_distance, min_peak_height, peak_threshold, ...
                                                                    position_est, peak_height_est, peak_index_est, ...
                                                                    peak_pos_tolerance, histogram_bins, smoothing_level, ...
                                                                    min_hist, max_hist, ref_peak_pos, peaks_fit(j));

        if max(peak_gaussian_ampl) >  max_peak_ampl
            max_peak_ampl = max(peak_gaussian_ampl);
        end
        if improv_param < 0                         % If true, the polynomial flattening didn't bring an improvement
            % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' 2nd poly fail.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
            images(:,:,j) = saved_image;            % Cancels the polynomial flattening
            peak_gaussian_std = peak_final_gaussian_std;
            peak_gaussian_ampl = peak_final_gaussian_ampl;
            disp('The second polymomial flattening did NOT improve the image')
            if plot_intermediate_steps == 1
                figure(summary_fig)
                xlabel('This step did NOT improve the image')
            end
        else
             % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' 2nd poly win.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
        end
         if improv_param == 0 && plot_intermediate_steps == 1
            figure(summary_fig)
            xlabel('No change')
         end

%         images(:,:,j) = images(:,:,j) - median(median(images(:,:,j)));
%         subplot(2, 7, 6);
%         plot_nice_image;
%         title('After 2nd polynomial flattening')

%         % Another median correction (Med_corr_2)
%         if strcmp(Med_corr_3(j), 'Y')
%
%             [images(:,:,j), peak_gaussian_std, peak_gaussian_ampl, improv_param]       = medianCorrection2(images(:,:,j), 2.5, peak_gaussian_std, peak_gaussian_ampl, min_spacing_stdev(j), max_spacing_stdev(j), amplitude_fract, min_peak_distance, min_peak_height, peak_threshold, histogram_bins);  % Med_corr_3
%
%             figure(summary_fig)
%             subplot(2, 7, 7);
%             plot_nice_image;
%             title('After 3rd median correction')
%             subplot(2, 7, 14);
%             hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
%             plot_peaks;
%             axis([-10 10 0 max(peak_height_1)*2])
%             title('Final histogram')
%
%             if improv_param < 0
%                 xlabel('This step did NOT improve the image');
%             elseif improv_param == 0
%                 xlabel('No change');
%             else
%                 xlabel('This step improved the image');
%             end
%         end

      %  saveas(gcf, [ImageDir,'Matlab',filesep,'Figure Figs/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_summary_image_processing.fig']);

Computes the final mask for direct one-step correction of the initial image (image after first median correction)

        if sequential_flattening == 1
            % already done
        else % Applies the final mask to the initial image one_step_corrected_image
            % Applies the final mask to the initial image
            % final_flattening_image and flatten
            image_test                                  = final_flattening_image(:,:);                                                                                                  % Makes a duplicate of the image
            image_test((images(:,:,j) < thresh_low))    = NaN;                                                                                                                              % Thresholds data below
            image_test((images(:,:,j) > thresh_high))   = NaN;                                                                                                                              % Thresholds data above

            % Fit the thresholded surface to an n x m polynomial
            coef_length                                     = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);                                                                      % Calculates the number of coeffecients
            [xInput_poly, yInput_poly, zOutput_poly]        = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);                                                 % Preps vectors for surface fitting
            poly_type                                       = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );                                                     % Sets the fit type to be an X-Y polynomial
            poly_opts                                       = fitoptions( poly_type );                                                                                                          % Tells it to fit a polynomial
            if index == 1
                poly_opts.Lower                             = -Inf*ones(coef_length,1);                                                                                                         % Sets the lower bounds
                poly_opts.Upper                             =  Inf*ones(coef_length,1);                                                                                                         % Sets the upper bounds
            else
                poly_opt.Lower                              = coeffvalues(Third_poly_coef{index-1,j})-0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
                poly_opt.Upper                              = coeffvalues(Third_poly_coef{index-1,j})+0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
            end                                                                                                       % Sets the upper bounds

            poly_result                                     = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );                                                            % Fits the polynomial to the actual data
            images(:,:,j)                                   = final_flattening_image(:,:) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));

            images(:,:,j)                                   = images(:,:,j) - median(median(images(:,:,j)));

            % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/text/','image ',Filename{index},' global fit alternate.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),'delimiter', '\t','precision', 5);
        end

Applies an offset to the whole image to put the reference peak at a given height reference_height

        if Offset_data(j) == 1

            image_inside_range                      = images(:,:,j);
            image_inside_range((images(:,:,j) < min_hist))   = NaN;
            image_inside_range((images(:,:,j) > max_hist))   = NaN;

            [height_offset_calc position_offset_calc]   = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);          % Calculates the histogram of the corrected height data (First image only)
            [peak_height_hist peak_index_hist]          = findpeaks(smooth(height_offset_calc,smoothing_level,'loess'),'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                         % Finds peaks in the corrected image to identify different terrace

            peak_index_hist                                = find_closest(ref_peak_pos, peak_index_hist(1:min(peaks_fit(j)+1, length(peak_index_hist))), position_offset_calc); % peak_index_hist: only keep the 5 higher peaks

            peak_index_hist = sort(peak_index_hist(1:min([3, peaks_fit(j), length(peak_index_hist)])),'ascend');       % Sorts peaks from smallest index (position) to highest one
            peak_height_hist = height_offset_calc(peak_index_hist);

            func_form                                   = [];
            t_hist_lower_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the lower bounds for fitting
            t_hist_lower_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the lower bounds for fitting
            t_hist_lower_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the lower bounds for fitting
            t_hist_upper_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the upper bounds for fitting
            t_hist_upper_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the upper bounds for fitting
            t_hist_upper_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                              % Initializes the upper bounds for fitting
            t_hist_start_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
            t_hist_start_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
            t_hist_start_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting

            for i                                       = 1:min(peaks_fit(j),length(peak_height_hist))                                                                                                % Loops for the minimum of the number of peaks or 3
                t_hist_lower_a(i)                       = (1-amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the lower amplitude bound for the fit
                t_hist_lower_s(i)                       = min_spacing_stdev(j);                                                                                                                     % Sets the lower standard deviation
                t_hist_lower_x(i)                       = position_offset_calc(max([1,peak_index_hist(i)-4]));                                                                                                   % Sets the lower center relative to the peak

                t_hist_upper_a(i)                       = (1+amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the upper amplitude bound for the fit
                t_hist_upper_s(i)                       = max_spacing_stdev(j);                                                                                                                     % Sets the upper standard deviation
                t_hist_upper_x(i)                       = position_offset_calc(min([peak_index_hist(i)+4,length(position_offset_calc)]));                                                                                                   % Sets the upper center relative to the peak

                t_hist_start_a(i)                       = 1.00*peak_height_hist(i);                                                                                                                % Sets the starting amplitude bound for the fit
                t_hist_start_s(i)                       = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;                                                                                            % Sets the starting standard deviation
                t_hist_start_x(i)                       = position_offset_calc(peak_index_hist(i));                                                                                                              % Sets the starting center relative at the peak
                func_form                               = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];                                                        % Generates the functional form to fit
                if i                                    < min(peaks_fit(j),length(peak_height_hist))
                    func_form                           = [func_form,'+'];                                                                                                                          % Generates the functional form to fit
                end
            end

            % Combines all the bouds into a form for the fits
            t_hist_lower                                = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];                                                                                           % Combines the lower bounds
            t_hist_upper                                = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];                                                                                           % Combines the upper bounds
            t_hist_start                                = [t_hist_start_a,t_hist_start_s,t_hist_start_x];                                                                                           % Combines the start point

            % Fits the data with the speficied conditions
            t_hist                                      = fitoptions(   'Method','NonlinearLeastSquares',...                                                                                        % Sets the conditions for fitting
                                                                        'Lower',      t_hist_lower,...
                                                                        'Upper',      t_hist_upper,...
                                                                        'Startpoint', t_hist_start,...
                                                                        'MaxIter',    1e2,...
                                                                        'MaxFunEvals',3e2);
            f_Number_hist                               = fittype(func_form,  'options' ,t_hist);
            [peak_hist_gaussian,peak_hist_gaussian_gof] = fit( position_offset_calc',...
                                                                 smooth(height_offset_calc,smoothing_level,'loess'),...
                                                                 f_Number_hist);

        % Combines all the amplitudes, the positions and the standard
        % deviations into one vector and sets the lower and upper threshold
        % for all peaks
         peak_hist_gaussian_ampl = [];
         peak_hist_gaussian_pos = [];
         peak_hist_gaussian_std = [];
         for i = 1:length(peak_height_hist)
            eval(['peak_hist_gaussian_ampl(', num2str(i), ')= peak_hist_gaussian.a', num2str(i), ';']);
            eval(['peak_hist_gaussian_pos(', num2str(i), ')= peak_hist_gaussian.x', num2str(i), ';']);
            eval(['peak_hist_gaussian_std(', num2str(i), ')= peak_hist_gaussian.s', num2str(i), ';']);
         end

         [peak_hist_pos, ind] = sort(peak_hist_gaussian_pos, 'ascend');
         for i = 1:length(peak_hist_pos)
                peak_hist_ampl(i) = peak_hist_gaussian_ampl(ind(i));
         end


         if length(peak_hist_pos) >= peak_reference_height(j)
                eval(['peak_hist_gaussian_x_ref = peak_hist_pos(', num2str(peak_reference_height(j)), ');']);
         elseif peak_reference_height(j) > 1 && length(peak_height_hist) >= peak_reference_height(j)-1
                eval(['peak_hist_gaussian_x_ref = peak_hist_pos(', num2str(peak_reference_height(j)-1), ');']);
         else
                peak_hist_gaussian_x_ref = peak_hist_pos(1);
         end

         images(:,:,j)                                     = images(:,:,j) + z_height_reference(j) - peak_hist_gaussian_x_ref;                           % Sets first peak to reference height                                                                                         % Adds the user selected Z-offset

        else
        end

        image_inside_range                      = images(:,:,j);
        image_inside_range((images(:,:,j) < min_hist))   = NaN;
        image_inside_range((images(:,:,j) > max_hist))   = NaN;
        [height_final position_final]               = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);          % Calculates the histogram of the corrected height data (First image only)
        [peak_height_final peak_index_final]          = findpeaks(smooth(height_final,smoothing_level,'loess'),'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');                                         % Finds peaks in the corrected image to identify different terrace
        % peak_index_final                                = find_closest(ref_peak_pos+z_height_reference(j)-peak_hist_gaussian_x_ref, peak_index_final, position_final);
        peak_index_final = sort(peak_index_final(1:min([3, peaks_fit(j), length(peak_index_final)])),'ascend');       % Sorts peaks from smallest index (position) to highest one
        peak_height_final = height_final(peak_index_final);
        peak_position_final = position_final(peak_index_final);

        while length(peak_height_final) < max(peaks_fit)
            peak_height_final = [peak_height_final NaN];
            peak_position_final = [peak_position_final NaN];
        end


        if plot_intermediate_steps == 1
            figure(summary_fig)
            % Plots all the histograms with the same axis
            subplot(2, 7, 10)
            axis([min_hist max_hist 0 max_peak_ampl*1.5])
            subplot(2, 7, 11)
            axis([min_hist max_hist 0 max_peak_ampl*1.5])
            subplot(2, 7, 12)
            axis([min_hist max_hist 0 max_peak_ampl*1.5])
            subplot(2, 7, 13)
            axis([min_hist max_hist 0 max_peak_ampl*1.5])
        end

        all_images(index, :, :, j) = images(:,:,j);
       % all_histo(index, :, :, j)  = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);

Computes a final histogram and CDF of the corrected image and fits Gaussians (like before) Displays the final histogram and fits and saves plots

        if S.GUI_flag == 0  % Program was not called from the GUI
             if index == 1 && j ==  channels_flatten(1)
                figure_hand_hist = figure('Name','Hist');
            else
                figure(figure_hand_hist);
            end
        else    % Program was called from the GUI
                axes(S.ax2)
        end

        [tok remain]                                    = strtok(Filename{index},'.');
        file_number                                     = str2double(remain(2:4));
        if index == 1 && ~overwrite_xlim(j) == 1
            [h_lim x_lim]                               = ecdf(reshape(images(:,:,j),size(images,1)*size(images,2),1));
            min_z(j)                                    = floor(x_lim(find(h_lim >= 0.01,1)));
            max_z(j)                                    =  ceil(x_lim(find(h_lim >= 0.99,1)));
            span_z                                      = max_z(j) - min_z(j);
            min_z(j)                                    =  ceil(min_z(j)-0.05*span_z);
            max_z(j)                                    = floor(max_z(j)+0.05*span_z);

        else
        end
        [y1_data x_data]                                = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),min_z(j):(max_z(j)-min_z(j))/(histogram_bins-1):max_z(j));
        y2_data                                         = zeros(size(y1_data,1),1);

        x_data                                          = x_data';
        y1_data                                         = y1_data';
        y1_data                                         = 100*y1_data/sum(y1_data);
        for i                                           = 1:size(y1_data,1)
            y2_data(i,1)                                = sum(y1_data(1:i));
        end

        [AX,H1,H2]                                      = plotyy(x_data,y1_data,x_data,y2_data,'plot');
        set(AX(1),'YColor','k');
        set(AX(2),'YColor','r');
        axis auto;

        set(AX(1),'YLim',[-0.002*ceil(max(y1_data)) 1.002*ceil(max(y1_data))],'YTick',0:ceil(max(y1_data))/10:ceil(max(y1_data)),'YTickLabel',0:ceil(max(y1_data))/10:ceil(max(y1_data)))
        set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
            axis_lim = [min_z(j),max_z(j)];
            set(AX(1),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
            set(AX(2),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));

        set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
        set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
        set(get(AX(1),'Ylabel'),'String','% Area','FontSize', FontSize, 'FontName', FontName,'Color','k');
        set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');

        set(H1,'Color','k','LineWidth',2);
        set(H2,'Color','r','LineWidth',2);

        xlabel([channel_info(j).Name,' (',channel_info(j).Unit,')'],'FontSize', FontSize, 'FontName', FontName );
        %title([channel_info(j).Name,'(',num2str(j),')',' ',num2str(file_number)],'FontSize', FontSize*1.5, 'FontName', FontName );
        set(AX(1),'FontSize', FontSize, 'FontName', FontName);
        set(AX(2),'FontSize', FontSize, 'FontName', FontName);

        if disp_percent == 1
        percent_val                                     = zeros(length(percent_search),1);
        text_disp                                       =  cell(length(percent_search),1);
        for i                                           = 1:length(percent_search)
            index_height                                = find(y2_data >= percent_search(i),1);
            percent_val(i)                              = x_data(index_height);
            text_disp{i,1}                              = [num2str(percent_search(i)),'% below ',num2str(percent_val(i),'%02.2f'),' ',channel_info(j).Unit];
        end
        else
        end

        x_range                                         = get(AX(1),'XLim');
        y_range                                         = get(AX(1),'YLim');
        text(x_range(2) -.02*(x_range(2)-x_range(1)),y_range(2)-.02*(y_range(2)-y_range(1)),text_disp,'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','top','HorizontalAlignment','right' );

        if fit_hist(j) == 1
            [peak_height_hist peak_index_hist]          = findpeaks(smooth(y1_data,smoothing_level,'loess'),'MinpeakDistance', 8, 'Sortstr','descend'); %min_peak_distance, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                       % Finds peaks in the corrected image to identify different terrace
            func_form                                   = [];
            t_hist_lower_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
            t_hist_lower_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
            t_hist_lower_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
            t_hist_upper_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
            t_hist_upper_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
            t_hist_upper_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
            t_hist_start_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
            t_hist_start_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
            t_hist_start_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting

            for i                                       = 1:min(peaks_fit(j),length(peak_height_hist))                                                                                                % Loops for the minimum of the number of peaks or 3
                t_hist_lower_a(i)                       = (1-amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the lower amplitude bound for the fit
                t_hist_lower_s(i)                       = min_spacing_stdev(j);                                                                                                                     % Sets the lower standard deviation
                t_hist_lower_x(i)                       = x_data(max([1,peak_index_hist(i)-3]));                                                                                                   % Sets the lower center relative to the peak

                t_hist_upper_a(i)                       = (1+amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the upper amplitude bound for the fit
                t_hist_upper_s(i)                       = max_spacing_stdev(j);                                                                                                                     % Sets the upper standard deviation
                t_hist_upper_x(i)                       = x_data(min([peak_index_hist(i)+3,length(x_data)]));                                                                                                   % Sets the upper center relative to the peak

                t_hist_start_a(i)                       = 1.00*peak_height_hist(i);                                                                                                                % Sets the starting amplitude bound for the fit
                t_hist_start_s(i)                       = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;                                                                                            % Sets the starting standard deviation
                t_hist_start_x(i)                       = x_data(peak_index_hist(i));                                                                                                              % Sets the starting center relative at the peak
                func_form                               = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];                                                        % Generates the functional form to fit
                if i                                    < min(peaks_fit(j),length(peak_height_hist))
                    func_form                           = [func_form,'+'];                                                                                                                          % Generates the functional form to fit
                end
            end

            % Combines all the bouds into a form for the fits
            t_hist_lower                                = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];                                                                                           % Combines the lower bounds
            t_hist_upper                                = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];                                                                                           % Combines the upper bounds
            t_hist_start                                = [t_hist_start_a,t_hist_start_s,t_hist_start_x];                                                                                           % Combines the start point

            % Fits the data with the speficied conditions
            t_hist                                      = fitoptions(   'Method','NonlinearLeastSquares',...                                                                                        % Sets the conditions for fitting
                                                                        'Lower',      t_hist_lower,...
                                                                        'Upper',      t_hist_upper,...
                                                                        'Startpoint', t_hist_start,...
                                                                        'MaxIter',    1e2,...
                                                                        'MaxFunEvals',3e2);
            f_Number_hist                               = fittype(func_form,  'options' ,t_hist);
            [peak_hist_gaussian,peak_hist_gaussian_gof] = fit( x_data,...
                                                                 smooth(y1_data,smoothing_level,'loess'),...
                                                                 f_Number_hist);

        hold on
        plot(x_data,peak_hist_gaussian(x_data),'.b','LineWidth',2);
        for i                                           = 1:min(peaks_fit(j),length(peak_height_hist))
        text(1.05*eval(['peak_hist_gaussian.x',num2str(i)]),...
         0.85*eval(['peak_hist_gaussian.a',num2str(i)]),...
         { [num2str(sum(eval(['peak_hist_gaussian.a',num2str(i),'*exp(-(x_data-peak_hist_gaussian.x',num2str(i),').^2/(2*(peak_hist_gaussian.s',num2str(i),')^2))'])),'%02.2f'),'% @'];...
           [num2str(eval(['peak_hist_gaussian.x',num2str(i)]),'%02.2f'),' +- ',num2str(eval(['peak_hist_gaussian.s',num2str(i)]),'%02.2f'),' ',channel_info(j).Unit]},...
         'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','bottom','HorizontalAlignment','left','color','b' );
        end

        else
        end

        if S.GUI_flag == 0
            print('-loose',  '-djpeg' ,  [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.jpg']);
            saveas(gcf,                                             [ImageDir,'Matlab',filesep,'Figure Figs/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.fig']);
            hold off;
            %close gcf;
        else
            if S.processAll == 1
%                 print('-loose',  '-djpeg' ,  [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.jpg']);
%                 saveas(gcf,                                             [ImageDir,'Matlab',filesep,'Figure Figs/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.fig']);
%                 hold off;
            end
        end

        hold off;

Plots the final corrected image with user set color and scale bars and saves the output

        ScanZ                                       = max_z(j)-min_z(j);                                                                                                                    % Sets the user defined data scale for each channel

%         if index == 1 && j ==  channels_flatten(1)
%             figure_hand_image = figure('Name','Image');
%         else
%             figure(figure_hand_image);
%         end


    if S.GUI_flag == 0  % Program was not called from the GUI
        if index == 1 && j ==  channels_flatten(1)
            figure_hand_image = figure('Name','Image');
        else
            figure(figure_hand_image);
        end
    else  % Program was called from the GUI
        axes(S.ax1)
    end


        imshow(imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),[min_z(j),max_z(j)])                                      % Generates a new Figure with the name of the correct channel,  Resizes the image to a square the size of the largest dimension.  This is useful when using nonsquare pixels.  Sets the display range from the min to the max

        if S.GUI_flag == 0
            print('-dtiff' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.tiff']);                                   % Saves the images.  I'm not sure if i dig the naming just yet.  I'll deal with it later
        else
%             if S.processAll == 1
%                 print('-dtiff' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.tiff']);
%             end
        end

        % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.txt'],imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]));

        % Loads the Georg approved colormap
        colormap(mycmap_sky);                                                                                                                                                               % Sets the colormap
            hcb= colorbar('YTick',min_z(j):ScanZ/10:max_z(j),...
                          'YTickLabel',{            [num2str(min_z(j)+0*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Inserts a color bar.  Sets the 1st tick with approriate labels
                                                    [num2str(min_z(j)+1*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 2nd tick with approriate labels
                                                    [num2str(min_z(j)+2*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 3rd tick with approriate labels
                                                    [num2str(min_z(j)+3*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 4th tick with approriate labels
                                                    [num2str(min_z(j)+4*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 5th tick with approriate labels
                                                    [num2str(min_z(j)+5*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 6th tick with approriate labels
                                                    [num2str(min_z(j)+6*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 7th tick with approriate labels
                                                    [num2str(min_z(j)+7*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
                                                    [num2str(min_z(j)+8*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
                                                    [num2str(min_z(j)+9*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 8th tick with approriate labels
                                                    [num2str(min_z(j)+10*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit]},...                                                                  % Sets the 9th tick with approriate labels
                                                    'FontSize', FontSize, 'FontName', FontName);                                                                                            % Sets the font and size for the color bar
            set(hcb,'YTickMode','manual','FontSize', FontSize, 'FontName', FontName);                                                                                                       % Apply the font to the colorbar


        %This plots the integral scale bar.  It works...  It's not all my code.
        DistFactor                                  = 0.05;                                                                                                                                 % Scalebar height
        Length                                      = ScaleBar;                                                                                                                             % User set scalebar length
        Pos                                         = [0.5*max(size(images(:,:,j))),0.95*max(size(images(:,:,j)))];                                                                         % positions scale bar
        Scale                                       = ScanSize/(size(images,2)*(1+crop_fraction));                                                                                                              % Image scale
        UnitsName                                   = ' nm';                                                                                                                                % Units for the scale bar
        NextPlot                                    = get(gca,'NextPlot');
        hold on;

        XLim                                        = get(gca,'XLim');
        YLim                                        = get(gca,'YLim');
        Xdiff                                       = abs(diff(XLim));
        Ydiff                                       = abs(diff(YLim));

        LineX                                       = Pos(1) + 0.5.*Length./Scale.*[-1;+1];
        LineY                                       = Pos(2).*[+1;+1];
        DistXdir                                    = 0;
        DistYdir                                    = -1;


        %--- plot line ---
        plot(LineX,LineY,'w-','LineWidth',2);


        %--- plot text ---
        DistX                                       = DistXdir.*DistFactor.*Xdiff;
        DistY                                       = DistYdir.*DistFactor.*Ydiff;

        Htext                                       = text(Pos(1)+DistX,Pos(2)+DistY,sprintf('%5.1f %s',Length,UnitsName));
        set(Htext,'HorizontalAlignment','center','Color', 'w', 'FontSize', FontSize, 'FontName', FontName);
        set(gca,'XLim',XLim);
        set(gca,'YLim',YLim);
        set(gca,'NextPlot',NextPlot);
        set(gca, 'XTickMode', 'manual');

       if S.GUI_flag == 0
            %print( '-r300', '-loose',  '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'.jpg']);
            print('-loose',  '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'.jpg']);                                   % Saves the images.  I'm not sure if i dig the naming just yet.  I'll deal with it later
            hold off;
            %close gcf;
        else
        end


      hold off;

      if S.GUI_flag       % program was called from the GUI
         if S.processAll  == 0     % If only process one image: save the final image and raw histogram for display
            final_images(:,:,channels_indices(j)) = images(:,:,j);
         end
      end

   pause(0.1) % just to get a good display...
     catch exception

        msgbox([{'An error occured while processing image:'}, {file_name}, {'Channel: ', channel_info(j).Name}], 'Error', 'error');
         if S.GUI_flag       % program was called from the GUI
            raw_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            med_corr_1_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            poly_flat_1_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            med_corr_2_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            poly_flat_2_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            final_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
            all_images(index, :, :, j) = zeros(size(images, 1), size(images, 2));
            if S.processAll  == 0     % If only process one image: save the variables before exiting the program
                if exist('min_spacing_stdev', 'var') == 0
                    min_spacing_stdev = min_spacing_stdev_user;
                    max_spacing_stdev = max_spacing_stdev_user;
                end
                if exist('min_hist', 'var') == 0
                    min_hist = 0;
                    max_hist = 0;
                end
                if exist('std_est', 'var') == 0
                    std_est = 0;
                end
                fit_hist(j) = 0;
                pause(2)
            end
         end

     end  % end of try


    end   % end of "for j" loop

      if S.GUI_flag       % program was called from the GUI
           if S.processAll  == 0     % If only process one image: save the final image and raw histogram for display
%                  % Get names of variables in caller's workspace
%                 w = evalin('caller','who');
%                 % Loop through variables and put them in base workspace.
%                 for i = 1:length(w)
%                 assignin('base',w{i},evalin('caller',w{i}));
%                 end
                save('Nano6_variables', 'channel_info', 'Filename', ...            % Saves only variables needed for plotting to be faster
                    'final_images', 'raw_images', 'med_corr_1_images', 'poly_flat_1_images', 'med_corr_2_images', 'poly_flat_2_images', ...
                    'min_spacing_stdev', 'max_spacing_stdev', 'mycmap_sky', 'min_z', 'max_z', 'histogram_bins', 'ScanSize', ...
                    'smoothing_level', 'min_hist', 'max_hist', 'std_est');
                break;
           end
      end
end

In case of "Process all", prints all the images and all the histograms in jpg and matlab figures

    if S.GUI_flag       % program was called from the GUI
       if S.processAll  ==  1 % process all images
           FontSize                    = 16;                                                                                              % Font Size for images
           FontName                    = 'Times';                                                                                          % Font for images

           for index = 1:length(Filename)                                                     % Loops for hte number of images
               for j                                           = 1:length(find(channels_flatten == 1))             % Loops for the number of different channels in the image
                   if max(max(all_images(index, :, :, j))) == 0 && min(min(all_images(index, :, :, j))) == 0 % => the program failed at this point and the image was set to zero
                       break;
                   end

--- Prints all histograms --- %%%

                    if index == 1 && j ==  1
                            figure_hand_hist = figure('Name','Saving all histograms...');
                     else
                             figure(figure_hand_hist);
                     end

                    [tok remain]                                    = strtok(Filename{index},'.');
                    file_number                                     = str2double(remain(2:4));

                    if index == 1 && ~overwrite_xlim(j) == 1
                        [h_lim x_lim]                               = ecdf(reshape(all_images(index, :, :, j),size(all_images,2)*size(all_images,3),1));
                        min_z(j)                                    = floor(x_lim(find(h_lim >= 0.01,1)));
                        max_z(j)                                    =  ceil(x_lim(find(h_lim >= 0.99,1)));
                        span_z                                      = max_z(j) - min_z(j);
                        min_z(j)                                    =  ceil(min_z(j)-0.05*span_z);
                        max_z(j)                                    = floor(max_z(j)+0.05*span_z);

                    else
                    end

                    [y1_data x_data]                                = hist(reshape(all_images(index, :, :, j),size(all_images,2)*size(all_images,3),1),min_z(j):(max_z(j)-min_z(j))/(histogram_bins-1):max_z(j));
                    y2_data                                         = zeros(size(y1_data,1),1);

                    x_data                                          = x_data';
                    y1_data                                         = y1_data';
                    y1_data                                         = 100*y1_data/sum(y1_data);
                    for i                                           = 1:size(y1_data,1)
                        y2_data(i,1)                                = sum(y1_data(1:i));
                    end

                    [AX,H1,H2]                                      = plotyy(x_data,y1_data,x_data,y2_data,'plot');
                    set(AX(1),'YColor','k');
                    set(AX(2),'YColor','r');
                    axis auto;

                    set(AX(1),'YLim',[-0.002*ceil(max(y1_data)) 1.002*ceil(max(y1_data))],'YTick',0:ceil(max(y1_data))/10:ceil(max(y1_data)),'YTickLabel',0:ceil(max(y1_data))/10:ceil(max(y1_data)))
                    set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
                        axis_lim = [min_z(j),max_z(j)];
                        set(AX(1),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
                        set(AX(2),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));

                    set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
                    set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
                    set(get(AX(1),'Ylabel'),'String','% Area','FontSize', FontSize, 'FontName', FontName,'Color','k');
                    set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');

                    set(H1,'Color','k','LineWidth',2);
                    set(H2,'Color','r','LineWidth',2);

                    xlabel([channel_info(j).Name,' (',channel_info(j).Unit,')'],'FontSize', FontSize, 'FontName', FontName );
                    title([channel_info(j).Name,'(',num2str(j),')',' ',num2str(file_number)],'FontSize', FontSize*1.5, 'FontName', FontName );
                    set(AX(1),'FontSize', FontSize, 'FontName', FontName);
                    set(AX(2),'FontSize', FontSize, 'FontName', FontName);

                    if disp_percent == 1
                    percent_val                                     = zeros(length(percent_search),1);
                    text_disp                                       =  cell(length(percent_search),1);
                    for i                                           = 1:length(percent_search)
                        index_height                                = find(y2_data >= percent_search(i),1);
                        percent_val(i)                              = x_data(index_height);
                        text_disp{i,1}                              = [num2str(percent_search(i)),'% below ',num2str(percent_val(i),'%02.2f'),' ',channel_info(j).Unit];
                    end
                    else
                    end

                    x_range                                         = get(AX(1),'XLim');
                    y_range                                         = get(AX(1),'YLim');
                    text(x_range(2) -.02*(x_range(2)-x_range(1)),y_range(2)-.02*(y_range(2)-y_range(1)),text_disp,'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','top','HorizontalAlignment','right' );

                    if fit_hist(j) == 1
                        [peak_height_hist peak_index_hist]          = findpeaks(smooth(y1_data,smoothing_level,'loess'),'MinpeakDistance', 8, 'Sortstr','descend'); %min_peak_distance, 'Threshold', peak_threshold, 'Sortstr','descend');                                                                       % Finds peaks in the corrected image to identify different terrace
                        func_form                                   = [];
                        t_hist_lower_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
                        t_hist_lower_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
                        t_hist_lower_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the lower bounds for fitting
                        t_hist_upper_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
                        t_hist_upper_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
                        t_hist_upper_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the upper bounds for fitting
                        t_hist_start_a                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
                        t_hist_start_s                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting
                        t_hist_start_x                              = zeros(1,min(peaks_fit(j),length(peak_height_hist)));                                                                                        % Initializes the start point for fitting

                        for i                                       = 1:min(peaks_fit(j),length(peak_height_hist))                                                                                                % Loops for the minimum of the number of peaks or 3
                            t_hist_lower_a(i)                       = (1-amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the lower amplitude bound for the fit
                            t_hist_lower_s(i)                       = min_spacing_stdev(j);                                                                                                                     % Sets the lower standard deviation
                            t_hist_lower_x(i)                       = x_data(max([1,peak_index_hist(i)-3]));                                                                                                   % Sets the lower center relative to the peak

                            t_hist_upper_a(i)                       = (1+amplitude_fract)*peak_height_hist(i);                                                                                                 % Sets the upper amplitude bound for the fit
                            t_hist_upper_s(i)                       = max_spacing_stdev(j);                                                                                                                     % Sets the upper standard deviation
                            t_hist_upper_x(i)                       = x_data(min([peak_index_hist(i)+3,length(x_data)]));                                                                                                   % Sets the upper center relative to the peak

                            t_hist_start_a(i)                       = 1.00*peak_height_hist(i);                                                                                                                % Sets the starting amplitude bound for the fit
                            t_hist_start_s(i)                       = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;                                                                                            % Sets the starting standard deviation
                            t_hist_start_x(i)                       = x_data(peak_index_hist(i));                                                                                                              % Sets the starting center relative at the peak
                            func_form                               = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];                                                        % Generates the functional form to fit
                            if i                                    < min(peaks_fit(j),length(peak_height_hist))
                                func_form                           = [func_form,'+'];                                                                                                                          % Generates the functional form to fit
                            end
                        end

                        % Combines all the bouds into a form for the fits
                        t_hist_lower                                = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];                                                                                           % Combines the lower bounds
                        t_hist_upper                                = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];                                                                                           % Combines the upper bounds
                        t_hist_start                                = [t_hist_start_a,t_hist_start_s,t_hist_start_x];                                                                                           % Combines the start point

                        % Fits the data with the speficied conditions
                        t_hist                                      = fitoptions(   'Method','NonlinearLeastSquares',...                                                                                        % Sets the conditions for fitting
                                                                                    'Lower',      t_hist_lower,...
                                                                                    'Upper',      t_hist_upper,...
                                                                                    'Startpoint', t_hist_start,...
                                                                                    'MaxIter',    1e2,...
                                                                                    'MaxFunEvals',3e2);
                        f_Number_hist                               = fittype(func_form,  'options' ,t_hist);
                        [peak_hist_gaussian,peak_hist_gaussian_gof] = fit( x_data,...
                                                                             smooth(y1_data,smoothing_level,'loess'),...
                                                                             f_Number_hist);

                    hold on
                    plot(x_data,peak_hist_gaussian(x_data),'.b','LineWidth',2);
                    for i                                           = 1:min(peaks_fit(j),length(peak_height_hist))
                    text(1.05*eval(['peak_hist_gaussian.x',num2str(i)]),...
                     0.85*eval(['peak_hist_gaussian.a',num2str(i)]),...
                     { [num2str(sum(eval(['peak_hist_gaussian.a',num2str(i),'*exp(-(x_data-peak_hist_gaussian.x',num2str(i),').^2/(2*(peak_hist_gaussian.s',num2str(i),')^2))'])),'%02.2f'),'% @'];...
                       [num2str(eval(['peak_hist_gaussian.x',num2str(i)]),'%02.2f'),' +- ',num2str(eval(['peak_hist_gaussian.s',num2str(i)]),'%02.2f'),' ',channel_info(j).Unit]},...
                     'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','bottom','HorizontalAlignment','left','color','b' );
                    end

                    end


                    print(gcf, '-loose',  '-djpeg' ,  [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.jpg']);
                    saveas(gcf, [ImageDir,'Matlab',filesep,'Figure Fig/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.fig']);

                    hold off;

Does Kumar's Grain Analysis if turned on

                grain_analysis      = get(S.disp_grain_analysis, 'value');
                image_test          = [];
                image_particles     = [];
                if grain_analysis   == 1
                    image_test      = all_images(index, :, :, j);
                    image_particles(1:size(image_test,2),1:size(image_test,3)) = image_test(1,1:size(image_test,2),1:size(image_test,3));
                    image_particles = wiener2(image_particles,[5,5]);
                    image_test_min  = (min(min(image_particles)));
                    image_test_span = (max(max(image_particles))-min(min(image_particles)));
                    part_thresh     = (peak_hist_gaussian.a1+peak_hist_gaussian.s1 - image_test_min)/image_test_span;
                    image_particles = (image_particles - image_test_min)./image_test_span;
                    image_particles = im2bw(image_particles, part_thresh);
                    B               = bwconncomp(image_particles);
                    disp(B);



                    stats = regionprops(B, 'Area');
                    idx = find([stats.Area] > 25);
                    BW2 = ismember(labelmatrix(B), idx);
                    graindata       = regionprops(BW2,'basic');
                    BW2 = bwconncomp(BW2);
                    disp(BW2);
                    grain_areas     = [graindata.Area]*(ScanSize/size(image_test,2))^2;
                    grain_centroids = [graindata.Centroid];

                    [y1_area_data x_area_data]       = hist(grain_areas,[min(grain_areas):(max(grain_areas)-min(grain_areas))/(length(grain_areas)*5):max(grain_areas)]);
                    y2_area_data                     = zeros(size(y1_area_data,1),1);
                    x_area_data                      = x_area_data';
                    y1_area_data                     = y1_area_data';
                    y1_area_data                     = 100*y1_area_data/sum(y1_area_data);
                    for i                            = 1:size(y1_area_data,1)
                        y2_area_data(i,1)            = sum(y1_area_data(1:i));
                    end

                    [AX,H1,H2]                       = plotyy(x_area_data,y1_area_data,x_area_data,y2_area_data,'plot');
                    set(AX(1),'YColor','k');
                    set(AX(2),'YColor','r');
                    axis auto;

                    set(AX(1),'YLim',[-0.002*ceil(max(y1_area_data)) 1.002*ceil(max(y1_area_data))],'YTick',0:ceil(max(y1_area_data))/10:ceil(max(y1_area_data)),'YTickLabel',0:ceil(max(y1_area_data))/10:ceil(max(y1_area_data)))
                    set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
                    set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
                    set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
                    set(get(AX(1),'Ylabel'),'String','% of total Particles','FontSize', FontSize, 'FontName', FontName,'Color','k');
                    set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');

                    set(H1,'Color','k','LineWidth',2);
                    set(H2,'Color','r','LineWidth',2);

                    xlabel('Particle area [nm^2]','FontSize', FontSize, 'FontName', FontName );
                    title( 'Particle area [nm^2]','FontSize', FontSize*1.5, 'FontName', FontName );
                    set(AX(1),'FontSize', FontSize, 'FontName', FontName);
                    set(AX(2),'FontSize', FontSize, 'FontName', FontName);

                    x_range                                         = get(AX(1),'XLim');
                    y_range                                         = get(AX(1),'YLim');

                    print(gcf, '-loose',  '-djpeg' ,  [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_part_area_hist_cdf.jpg']);



                    part_vol    = zeros(length(grain_areas),1);
                    part_height = zeros(length(grain_areas),1);
                    image_test      = all_images(index, :, :, j);
                    image_particles(1:size(image_test,2),1:size(image_test,3)) = image_test(1,1:size(image_test,2),1:size(image_test,3));
                    image_particles = wiener2(image_particles,[5,5]);
                    for i       = 1:length(grain_areas)
                        grain   = false(size(BW2));
                        grain(BW2.PixelIdxList{i}) = true;
                        disp(size(grain));
                        disp(size(image_particles));
                        %part_height(i)  = max(max(immultiply(image_particles,grain)));
                        %part_vol(i)     = sum(sum(immultiply(image_particles,grain)*(ScanSize/size(image_test,2))^2));
                    end

                    disp(part_height);
                else
                end
               end
           end

           close(figure_hand_hist)

--- Prints all images --- %%%

            for index = 1:length(Filename)                                                     % Loops for hte number of images
               for j                                           = 1:length(find(channels_flatten == 1))             % Loops for the number of different channels in the image

                   [tok remain]                                    = strtok(Filename{index},'.');
                   file_number                                     = str2double(remain(2:4));

                    % Plots the final corrected image with user set color and scale bars and saves the output

                    ScanZ                                       = max_z(j)-min_z(j);            %Sets the user defined data scale for each channel

                     if index == 1 && j ==  1
                            figure_hand_image = figure('Name','Saving all images...');
                     else
                             figure(figure_hand_image);
                     end

                     if max(max(all_images(index, :, :, j))) == 0 && min(min(all_images(index, :, :, j))) == 0 % => the program failed at this point and the image was set to zero
                       break;
                     end


                    images(:,:, j) = all_images(index, :, :, j);
                    imshow(imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),[min_z(j),max_z(j)])                                      % Generates a new Figure with the name of the correct channel,  Resizes the image to a square the size of the largest dimension.  This is useful when using nonsquare pixels.  Sets the display range from the min to the max

                     print(gcf, '-dtiff' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.tiff']);                                   % Saves the images.  I'm not sure if i dig the naming just yet.  I'll deal with it later

                    % dlmwrite([ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.txt'],imresize(all_images(index, :, :, j),[max(size(all_images(index, :, :, j))),max(size(all_images(index, :, :, j)))]));

                    % Loads the Georg approved colormap
        color_string    = get(S.ColorBarSelect, 'string');
        color_value     = get(S.ColorBarSelect, 'value');
        switch color_string{color_value}
            case 'Autumn'
                colormap(autumn);
            case 'Bone'
                colormap(bone);
            case 'ColorCube'
                colormap(colorcube);
            case 'Cool'
                colormap(cool);
            case 'Copper'
                colormap(copper);
            case 'Flag'
                colormap(flag);
            case 'Gray'
                colormap(gray);
            case 'Hot'
                colormap(hot);
            case 'HSV'
                colormap(hsv);
            case 'Jet'
                colormap(jet);
            case 'Lines'
                colormap(lines);
            case 'Pink'
                colormap(pink);
            case 'Prism'
                colormap(prism);
            case 'Spring'
                colormap(spring);
            case 'Summer'
                colormap(summer);
            case 'Winter'
                colormap(winter);
            case 'Sky'
                load('Sky','mycmap_sky')
                colormap(mycmap_sky);
        end

                    % Sets the colormap
                        hcb= colorbar('YTick',min_z(j):ScanZ/10:max_z(j),...
                                      'YTickLabel',{            [num2str(min_z(j)+0*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Inserts a color bar.  Sets the 1st tick with approriate labels
                                                                [num2str(min_z(j)+1*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 2nd tick with approriate labels
                                                                [num2str(min_z(j)+2*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 3rd tick with approriate labels
                                                                [num2str(min_z(j)+3*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 4th tick with approriate labels
                                                                [num2str(min_z(j)+4*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 5th tick with approriate labels
                                                                [num2str(min_z(j)+5*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 6th tick with approriate labels
                                                                [num2str(min_z(j)+6*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 7th tick with approriate labels
                                                                [num2str(min_z(j)+7*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
                                                                [num2str(min_z(j)+8*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
                                                                [num2str(min_z(j)+9*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...                                                                    % Sets the 8th tick with approriate labels
                                                                [num2str(min_z(j)+10*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit]},...                                                                  % Sets the 9th tick with approriate labels
                                                                'FontSize', FontSize, 'FontName', FontName);                                                                                            % Sets the font and size for the color bar
                        set(hcb,'YTickMode','manual','FontSize', FontSize, 'FontName', FontName);                                                                                                       % Apply the font to the colorbar

disp_scale_bar      = get(S.disp_scale_bar, 'value');
if disp_scale_bar == 1
                    %This plots the integral scale bar.  It works...  It's not all my code.
                    DistFactor                                  = 0.05;                                                                                                                                 % Scalebar height
                    Length                                      = ScaleBar;                                                                                                                             % User set scalebar length
                    Pos                                         = [0.5*max(size(all_images(index, :, :, j))),0.95*max(size(all_images(index, :, :, j)))];                                                                         % positions scale bar
                    Scale                                       = ScanSize/(size(all_images,3)*(1+crop_fraction));                                                                                                              % Image scale
                    UnitsName                                   = ' nm';                                                                                                                                % Units for the scale bar
                    NextPlot                                    = get(gca,'NextPlot');
                    hold on;

                    XLim                                        = get(gca,'XLim');
                    YLim                                        = get(gca,'YLim');
                    Xdiff                                       = abs(diff(XLim));
                    Ydiff                                       = abs(diff(YLim));

                    LineX                                       = Pos(1) + 0.5.*Length./Scale.*[-1;+1];
                    LineY                                       = Pos(2).*[+1;+1];
                    DistXdir                                    = 0;
                    DistYdir                                    = -1;


                    %--- plot line ---
                    plot(LineX,LineY,'w-','LineWidth',2);


                    %--- plot text ---
                    DistX                                       = DistXdir.*DistFactor.*Xdiff;
                    DistY                                       = DistYdir.*DistFactor.*Ydiff;

                    Htext                                       = text(Pos(1)+DistX,Pos(2)+DistY,sprintf('%5.1f %s',Length,UnitsName));
                    set(Htext,'HorizontalAlignment','center','Color', 'w', 'FontSize', FontSize, 'FontName', FontName);
                    set(gca,'XLim',XLim);
                    set(gca,'YLim',YLim);
                    set(gca,'NextPlot',NextPlot);
                    set(gca, 'XTickMode', 'manual');
else
end
                    print(gcf, '-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'.jpg']);                                   % Saves the images.  I'm not sure if i dig the naming just yet.  I'll deal with it later

                  hold off;


               end
            end

            close(figure_hand_image)
       end
    end