% STEP 4: Identification and Counting of Polonies

% The focus of Step 4 is to localize and count polonies bearing each of the 1023 possible signatures.  

% 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. Initialize object list matrix and final polony matrix
    
OBJECT_LIST = [];

FINAL_POLONY_LIST = [];

% 4. Set "view image" flag.  If this flag is set to one, then the object identification process can be viewed as it proceeds.  The signature image is displayed and black plus-marks 
% are added at predicted polony centroids as the algorithm is running.

VIEW_FLAG = 0;
    
% 5. 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),'/');
    
    % Load Signature Image
    
    load(strcat(DIR_NAME,'signature_image.mat'));
    
    % Apply Gel Mask
    
    SIG_IMAGE = imsubtract(SIG_IMAGE,GEL_MASK);
    
    % Load background-subtracted, median-filtered Exon 10 scans.  Signal intensity distributions in these images will be used to identify situations where multiple polonies with 
    % the same signature are overlapping, such that they can be appropriately counted as multiple objects rather than as a single joined object.  There are a handful of rare 
    % isoforms that do not contain Exon 10, but the odds of these overlapping are diminishingly low.
 
    SCAN_E10 = imread(strcat(DIR_NAME,'SL',num2str(D),'_SC01_E10_Cy3.Tif.Sub'),'Tif');

    % Initialize image

    if (VIEW_FLAG == 1) 
   
       hold off;
       close all;
       hold on;
       close all;
       hold on;
       imagesc(SIG_IMAGE);
        
    end
       
    % Iterate over signature values 
    
    for SIG_INDEX = 1:1023

        [D SIG_INDEX size(FINAL_POLONY_LIST)]
        
        % Grab signature to consider from the index array

        SIG = SIG_INDEX;
        
        % Define the set of pixels with signatures matching the current signature being analyzed
           
        OBJECT_IMAGE    = uint16(bwlabel(SIG_IMAGE==SIG));
        
        % Initialize temporary matricies for storing candidate object information
        
        PHASE_1_CANDIDATE_COMMON_OBJECT_LIST = [];
        PHASE_2_CANDIDATE_COMMON_OBJECT_LIST = [];
        
        % Proceed to object finding only if there are one or more pixels matching this signature
           
        if (find(OBJECT_IMAGE) ~= []) 
                              
            % Collect, in parallel, morphological information on the full set of objects (e.g. connected/adjacent groups of pixels) bearing this signature
                
            OBJECT_PROPERTIES = [];
                
            OBJECT_PROPERTIES(:,1) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'Area')))';
            OBJECT_PROPERTIES(:,2) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'Solidity')))';
            OBJECT_PROPERTIES(:,3) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'Eccentricity')))';
            OBJECT_PROPERTIES(:,4) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'Extent')))';
            OBJECT_PROPERTIES(:,5) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'MajorAxisLength')))';
            OBJECT_PROPERTIES(:,6) = cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'MinorAxisLength')))';
            OBJECT_PROPERTIES = horzcat(OBJECT_PROPERTIES, reshape(cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'Centroid'))),2,[])');
            OBJECT_PROPERTIES = horzcat(OBJECT_PROPERTIES, reshape(cell2mat(struct2cell(RegionProps(OBJECT_IMAGE,'BoundingBox'))),4,[])');

            OBJECT_COUNT = size(OBJECT_PROPERTIES);
                           
            % Phase 1: Iterate over the full set of objects bearing this signature
                
            for (j = 1:OBJECT_COUNT(1))
                    
                % Apply a set of morphological criteria to decide whether this object is polony-like
                    
                if ((OBJECT_PROPERTIES(j,1) > 100 & OBJECT_PROPERTIES(j,2) > 0.50 & (OBJECT_PROPERTIES(j,5)/OBJECT_PROPERTIES(j,6)) < 3) | (OBJECT_PROPERTIES(j,1) > 25 & OBJECT_PROPERTIES(j,2) > 0.80 & (OBJECT_PROPERTIES(j,5)/OBJECT_PROPERTIES(j,6)) < 1.5) | OBJECT_PROPERTIES(j,1) > 500)
          
                    % Examine fraction of pixels that are 'perimeter pixels' for this object
                    
                    PERIM_ARE = OBJECT_IMAGE((OBJECT_PROPERTIES(j,10)-5.5):((OBJECT_PROPERTIES(j,10)+OBJECT_PROPERTIES(j,12)+5.5)),(OBJECT_PROPERTIES(j,9)-5.5):((OBJECT_PROPERTIES(j,9)+OBJECT_PROPERTIES(j,11)+5.5)));                    
                    PERIM_BOX = medfilt2((PERIM_ARE==j),[3 3]);  
                    PERIM_PIX = uint16(bwperim(PERIM_BOX));
                    PERIM_TP = size(find((reshape(PERIM_BOX,[],1))>0));
                    PERIM_PP = size(find((reshape(PERIM_PIX,[],1))>0));
                    PERIM_FR = PERIM_PP(1)/PERIM_TP(1);
                    OBJECT_PROPERTIES(j,3) = PERIM_FR;
                    
                    % End perimeter section
                    
                    % What would you expect for a filled circle of this area?  (approximation based on perimeter of a box housing a circle of radius based on area of object under consideration)
                    
                    r = (2 * pi * sqrt(OBJECT_PROPERTIES(j,1)/pi));
                    OBJECT_PROPERTIES(j,4) = (r / PERIM_TP(1));
                                        
                    % Add an additional morphological criterion based on perimeter fractions
                    
                    if ((OBJECT_PROPERTIES(j,3) / OBJECT_PROPERTIES(j,2)) < 1)
                    
                        APPENDED_PROPERTIES = horzcat(j,SIG,OBJECT_PROPERTIES(j,:));
                        
                        PHASE_1_CANDIDATE_COMMON_OBJECT_LIST = vertcat(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST,APPENDED_PROPERTIES);
        
                    end
                    
                end
                    
            end
                       
            % Phase 2: Reiterate over the list and decide whether each given object actually consists of multiple polonies of the same signature
            
            PHASE_1_OBJECT_COUNT = size(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST);
                    
            if (PHASE_1_OBJECT_COUNT > 0 & SIG < 512)
                                       
                PHASE_2_CANDIDATE_COMMON_OBJECT_LIST = horzcat(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(:,:),[PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(:,9) PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(:,10)]);
                               
            elseif (PHASE_1_OBJECT_COUNT > 0)

                for (j = 1:PHASE_1_OBJECT_COUNT(1))
                
                    % Pull sub-region of image (Bounding Box for object) from Exon 10 scan and from object-image for this signature
                
                    Box1 = SCAN_E10     ((PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,12)+0.5):(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,12)-0.5+PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,14)),(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,11)+0.5):(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,11)-0.5+PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,13)));              
                    Box2 = OBJECT_IMAGE ((PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,12)+0.5):(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,12)-0.5+PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,14)),(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,11)+0.5):(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,11)-0.5+PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,13)));
               
                    % Determine overlap of polony object with Exon 10 scan
                
                    Box3 = uint16(double(Box1) .* double(Box2==PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,1)));
                
                    % Define regional maxima of object based on Exon 10 scan intensity values
                
                    Box4 = imextendedmax(Box3,3000);
                    Box5 = uint16(bwlabel(Box4));
                
                    % Define perimeter pixels of object
                
                    Box6 = bwperim(Box3,8);

                    % Evaluate regional maxima
                
                    if (find(Box5) ~= [])
                
                        % Collect information on regional maxima within this object
                    
                        RM_Area = cell2mat(struct2cell(RegionProps(Box5,'Area')))';
                       RM_Centroid = reshape(cell2mat(struct2cell(RegionProps(Box5,'Centroid'))),2,[])';

                        % Iterate over the regional maxima within this object
                    
                        RM_number = size(RM_Area);

                        Validated_Count = 0;
                    
                        for (k = 1:RM_number(1))
                        
                            % If the regional maxima pixels do not touch the border of this object OR 
                            % there is only one regional maxima within the polony OR 
                            % the number of pixels within the regional maximum focus area that are not touching the perimeter of the polony is greater than four...
                            % ...then add this regional maximum as a polony object location
                                   
                             BorderPixCount = size(find(double(Box6) .* double(Box5==k)));
                        
                            InternalPixCount = RM_Area(k)- BorderPixCount(1);
                        
                            if (BorderPixCount(1) == 0 | InternalPixCount > 4 | RM_number(1) == 1)
                        
                                MyCentroid_X = RM_Centroid(k,1) + PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,11);
                                MyCentroid_Y = RM_Centroid(k,2) + PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,12);
                            
                                NewLine = horzcat(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,:),[MyCentroid_X MyCentroid_Y]);
                                PHASE_2_CANDIDATE_COMMON_OBJECT_LIST = vertcat(PHASE_2_CANDIDATE_COMMON_OBJECT_LIST,NewLine);
                            
                                Validated_Count = Validated_Count + 1;

                                if (VIEW_FLAG == 1)
                                    
                                    % Mark polony locations with thick plus signs
                                       
                                    plot(NewLine(15),NewLine(16),'black+');
                                    pause(0.01);
                                    
                                end
    
                            end 
                                                
                        end
                    
                        % If none of the >1 regional maxima met the criteria, go ahead and add the centroid of the polony as a legitamate object
                            
                        if (Validated_Count == 0)
                            
                            NewLine = horzcat(PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,:),[PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,9) PHASE_1_CANDIDATE_COMMON_OBJECT_LIST(j,10)]);
                            PHASE_2_CANDIDATE_COMMON_OBJECT_LIST = vertcat(PHASE_2_CANDIDATE_COMMON_OBJECT_LIST,NewLine);
                        
                            if (VIEW_FLAG == 1)
                            
                                % Mark polony locations with thick plus signs
                                       
                                plot(NewLine(15),NewLine(16),'black+');
                                pause(0.01);
                                    
                           end
                            
                        end
                    
                    end
                
                end
                
            end
            
        end
        
        PHASE_2_OBJECT_COUNT = size(PHASE_2_CANDIDATE_COMMON_OBJECT_LIST);
        
        if (PHASE_2_OBJECT_COUNT(1) > 0) 
            
            PHASE_2_CANDIDATE_COMMON_OBJECT_LIST(:,1) = D;        

            FINAL_POLONY_LIST = vertcat(FINAL_POLONY_LIST,PHASE_2_CANDIDATE_COMMON_OBJECT_LIST);
    
        end
        
    end % End iteration over 1023 signatures
      
end % End iteration over eight slides

% 6. Screen out candidate polonies whose centroids are in excluded regions.  These rectangular blocks encompass regions were identified as poor-quality on visual inspection of
% the original and signature images.  Polony objects whose centroids fall within these blocks will be excluded from consideration.

C = FINAL_POLONY_LIST;
C2 = [];
s = size(C);

for (i = 1:s(1))

    flag = 0;
    if (C(i,1) == 1 & C(i,15) > 400 & C(i,15) < 800 & C(i,16) > 250 & C(i,16) < 500) flag = 1; end  % SL1,  250: 500, 400: 800
    if (C(i,1) == 2 & C(i,15) > 150 & C(i,15) < 500 & C(i,16) >1400 & C(i,16) <1800) flag = 1; end  % SL2, 1400:1800, 150: 500
    if (C(i,1) == 3 & C(i,15) > 850 & C(i,15) < 950 & C(i,16) > 350 & C(i,16) < 450) flag = 1; end  % SL3,  350: 450, 850: 950
    if (C(i,1) == 4 & C(i,15) > 800 & C(i,15) <1200 & C(i,16) >1900 & C(i,16) <2300) flag = 1; end  % SL4, 1900:2300, 800:1200
    if (C(i,1) == 4 & C(i,15) > 950 & C(i,15) <1050 & C(i,16) >2400 & C(i,16) <2600) flag = 1; end  % SL4, 2400:2600, 950:1050
    if (C(i,1) == 7 & C(i,15) > 800 & C(i,15) <1200 & C(i,16) > 450 & C(i,16) < 850) flag = 1; end  % SL7,  450: 850, 800:1200
    if (C(i,1) == 8 & C(i,15) >   0 & C(i,15) < 400 & C(i,16) >   0 & C(i,16) <2600) flag = 1; end  % SL8,    0:2600,   0: 400
    if (C(i,1) == 8 & C(i,15) >1200 & C(i,15) <1800 & C(i,16) > 700 & C(i,16) <1550) flag = 1; end  % SL8,  700:1550,1200:1800
    
    if (flag == 0)
        C2 = vertcat(C2,C(i,:));    
    end
    
end

FINAL_POLONY_LIST = C2;

% 7. Save final list of locations of identified polonies

save('final_polony_list.mat','FINAL_POLONY_LIST');

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

clear, close all
 
 