
% STEP 2: Threshold Determination

% This script determines approriate thresholding values for each image.  As images are scanned at varying laser powers and PMTs, and as different single-base-extensions may yield
% varying amounts of signal for a given quantity of extended DNA, the slides essentially need to be normalized to one another.  In moving forward we will be 'thresholding' the images.
% In other words, we will be setting a cutoff, intensities above which are assigned a '1', and below which are assigned a '0'.  Due to the stated reasons we cannot use the same cutoff
% for every image.  Here we seek to determine the approximetely correct cutoff to use for each image.

% The algorithm proceeds by first defining the set of pixels that are above a minimum value on every image.  These pixels generally will correspond to polonies that include all 10 exons.
% How this set of 'focus pixels' varies across the ten images will be used to define appropriate thresholds.

% On manual inspection of results it was noticed that the algorithm did not perform well on slides for which the background (e.g. post-stripping) images did not accurately reflect
% the background, such that the background-subtracted image still had moderately high background.  For these four images (out of forty images in total), thresholds decided on via
% manual inspection are forcibly assigned.

% 1. Reset variables & close image windows to free up memory.

clear, close all;

% 2. Load Gel Mask.  The Gel Mask is an image that defines the borders of the gel, such that pixels outside of the gel area can be set to zero.

load('GEL_MASK_100');

GEL_MASK = uint16(NEW_GEL_MASK .* 65535);

clear NEW_GEL_MASK;

% 3. Iterate over the eight slides (as we are processing all images in parallel).

for D = 1:8

    % Define the directory to be used
    
    DIR_NAME = strcat('./IMAGES/SL',num2str(D),'/');
    
    % Define the list of 10 files in that directory (e.g. files produced by Step 1 that bear the ".sub" extension)
    
    FILELIST = dir(strcat(DIR_NAME,'SL',num2str(D),'*.sub'));
    
    % Define the set of pixels that have a value greater than a set minimum across the full set of 10 images
        
    MIN_SIGNAL = 4000;                                          % Define the minimum intensity
              
    FOCUS_PIXELS = ones(size(GEL_MASK));                        % Initially everyone is a focus pixel!
        
    FOCUS_PIXELS = imsubtract(uint16(FOCUS_PIXELS),GEL_MASK);   % Block out pixels outside of the gel area
        
    for i = 1:10                                                % Iterate over the ten images
              
        IMAGE = imread(strcat(DIR_NAME,FILELIST(i).name));      % Load a given image
            
        IMAGE = IMAGE < MIN_SIGNAL;                             % Define the subset of pixels less than the minimum intensity
        
        FOCUS_PIXELS = imsubtract(FOCUS_PIXELS,uint16(IMAGE));  % Block out these pixels from the focus set
            
    end                                                         % Loop to next image
 
    % Now that the 'focus pixels' have been defined, go back through the ten images and grab the intensity values at these pixel positions
    
    FOCUS_LIST = find(FOCUS_PIXELS);
    
    IMAGE = imread(strcat(DIR_NAME,FILELIST(1).name));

    clear SLIST;
    
    for i = 1:10
        
        IMAGE = imread(strcat(DIR_NAME,FILELIST(i).name));
   
        FILELIST(i).name;
        
        SLIST(:,i) = double(IMAGE(FOCUS_LIST));
                 
    end

    % Sort the array containing the intensity values at focus pixels according to the intensities in the first image (column 1)
    
    SLIST = sortrows(SLIST,1);

    % Define the thresholds as the median value in the first 100 rows of the sorted focus-pixel intensity array
    
    thresholds(D,:) = median(SLIST(1:100,:))
    
end

% 5. Force specific thresholds on the four slides where the algorithm appears to have failed on manual inspection, as discussed above.

thresholds(5,10)= 6000;
thresholds(6,1) = 7000;
thresholds(6,2) = 11000;
thresholds(6,3) = 11000;
thresholds(6,4) = 11000;

% 6. Save the information!!

save('thresholds.mat','thresholds');

% 7. Reset variables & close image windows to free up memory.

clear, close all


