
% Samantha Sutton 2008
%
% This document is a concatenation of 3 scripts: main_final.m, CNR_final.m,
% and BF_mask_final.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script main_final.m
%
% This is the main script to run when analyzing cells. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


YFP  = YFP.tif       % YFP image of cells expressing the Phospholocator. Can be of any bit size (8-bit, 12-bit, 16-bit, etc.)
BF   = BF-8bit.tif   % BF image of cells expressing the Phospholocator. Must be 8-bit.
RFP  = RFP-8bit.tif  % RFP image of cells expressing the Phospholocator. Must be 8-bit.

YFP_back  = YFP.tif       % YFP image of cells not expressing the Phospholocator. Can be of any bit size (8-bit, 12-bit, 16-bit, etc.)
BF_back   = BF-8bit.tif   % BF image of cells not expressing the Phospholocator. Must be 8-bit.
RFP_back  = RFP-8bit.tif  % RFP image of cells not expressing the Phospholocator. Must be 8-bit.


[number_cells, YFP_total_average, cytoplasmic_nuclear_vector, geometric_cytosolic_nuclear,...
    std_error_average_total, cytosolic_YFP_average, nuclear_YFP_average, cell_size]...
    = CNR_final(BF, RFP, YFP, BF_back, RFP_back, YFP_back)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Measurements output by CNR_final:
%
% -- number_cells reports the number of cells identified and analyzed.
% -- YFP_total_average is a vector containing the average YFP value (per pixel) of each cell analyzed.
% -- cytoplasmic_nuclear_vector is a vector containing the average cytoplasmic:average nuclear fluorescence ratios
%    of each cell analyzed.
% -- geometric_cytosolic_nuclear is the geometic mean of cytoplasmic_nuclear_vector.
% -- std_error_average_total is the standard error of cytoplasmic_nuclear_vector.
% -- cytosolic_YFP_average is a vector containing the average YFP value (per pixel) of the cytosol of each cell analyzed.
% -- nuclear_YFP_average is a vector containing the average YFP value (per pixel) of the nucleus of each cell analyzed.
% -- cell_size is a vector containing the total pixel size of each cell analyzed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script CNR_final.m
%
% This script is called by main_final.m. The script measures nuclear boundaries, and computes final  fluorescence levels. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



function [number_cells, YFP_total_average, cytoplasmic_nuclear_vector, geometric_cytosolic_nuclear,...
    std_error_average_total, cytosolic_YFP_average, nuclear_YFP_average, cell_size]...
    = CNR_final(BF, RFP, YFP, BF_back, RFP_back, YFP_back)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For our microscope, RFP images were slightly shifted from the BF and YFP
% images (5 pixels down and 2 pixels to the left). Thus, we must shift the 
% RFP image to align with the others. Here, we define the region
% encompassing the upper 5 rows and rightmost 2 columns of RFP images.
% This is the "danger zone."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

picture_size = size(RFP);
danger_zone_counter = 1;
shift_up = 5;
shift_right = 2;

% cordon off the upper "shift_up" rows

for i = 1:picture_size(2)
        danger_zone(danger_zone_counter) = (i - 1) * picture_size(1) + 1;
        danger_zone(danger_zone_counter + 1) = danger_zone(danger_zone_counter) + 1;
        danger_zone(danger_zone_counter + 2) = danger_zone(danger_zone_counter) + 2;
        danger_zone_counter = danger_zone_counter + shift_up; 
end

% cordon off the right-most "shift_right rows", now called the danger_zone_annex

danger_zone_annex = [(picture_size(1) * (picture_size(2) - shift_right) + 1): (picture_size(1) * picture_size(2))];

% concatenate the two danger zones together

danger_zone = cat(2,danger_zone, danger_zone_annex);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert RFP image of cell expressing the Phospholocator to Binary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
level_RFP = graythresh(RFP);
bw_RFP = im2bw(RFP,level_RFP);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Label cells
%
% Notes:
% "labeled_..." is a matrix showing each object and where it is located
% "numObjects_..." gives you the number of objects in labeled.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[labeled_RFP,numObjects_RFP] = bwlabel(bw_RFP,8); 


d_RFP = size(labeled_RFP);
q = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating a mask that passes the "skinny" nucleii. We pass nucleii that:
%       - aren't found in the "danger zone" (the area that we will shift 
%       up and over to align images)
%       - size of nucleii is > 100 and < 900
% Once a nucleus has passed our test, we shift it up 5 pixels and over 2
% to the right.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:numObjects_RFP
        pix_temp= find(labeled_RFP == i);
        if (length(pix_temp)>100 & length(pix_temp) < 900)
            whammie = 1;
            for j = 1:length(pix_temp)
                if (find(danger_zone == pix_temp(j)))
                    whammie = 0;
                    break 
                end
            end
            if(whammie)  % shift all pixels up and over
                for k = 1:length(pix_temp)
                    pix_temp(k) = pix_temp(k) + shift_right * picture_size(1) - shift_up;
                end
                q = q + 1;
                pix_RFP(q) = {pix_temp};
                if (q==1)
                    mask_RFP_skinny = pix_RFP{q};
                else
                    mask_RFP_skinny = cat(1,mask_RFP_skinny, pix_RFP{q});
                end
                v = floor(pix_RFP{q}(1)/d_RFP(1));
                u = mod((pix_RFP{q}(1)-v*d_RFP(1)),d_RFP(1));
                keepers_RFP(q,1) = i;
                keepers_RFP(q,2) = pix_temp(1);

            end
        end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We dilate bw_RFP so that we can define a "no man's land" around the 
% nucleus that is counted as neither nucleus nor cytoplasm. This is to 
% prevent bleed through.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% dilate the image twice
se90 = strel('line', 3, 90);
se0 = strel('line', 3, 0);

RFP_dil_temp = imdilate(bw_RFP, [se90 se0]);

RFP_dil = imdilate(RFP_dil_temp, [se90 se0]);

[labeled_dil,numObjects_dil] = bwlabel(RFP_dil,8); 


d_dil = size(labeled_dil);
q = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generating a mask that passes the "dilated" nucleii. We pass nucleii that:
%       - aren't found in the "danger zone" (the area that we will shift 
%       up and over to align images)
%       - size of nucleii is > 100 and < 900
% Once a nucleus has passed our test, we shift it up 5 pixels and over 2
% to the right.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for i = 1:numObjects_dil
        pix_dil_temp = find(labeled_dil == i);
        if (length(pix_dil_temp)>100)
            whammie = 1;
            for j = 1:length(pix_dil_temp)
                if (find(danger_zone == pix_dil_temp(j)))
                    whammie = 0;
                    break 
                end
            end
            if(whammie)  % shift all pixels up and over
                for k = 1:length(pix_dil_temp)
                    pix_dil_temp(k) = pix_dil_temp(k) + shift_right * picture_size(1) - shift_up;
                end
                q = q + 1;
                pix_dil(q) = {pix_dil_temp};
                pix_intensity_YFP_dil{q} = YFP(pix_dil{q});   
                
                if (q==1)
                    mask_RFP_dil = pix_dil{q};
                else
                    mask_RFP_dil = cat(1,mask_RFP_dil, pix_dil{q});
                end

                v = floor(pix_dil{q}(1)/d_dil(1));
                u = mod((pix_dil{q}(1)-v*d_dil(1)),d_dil(1));
                keepers_dil(q,1) = i;
                keepers_dil(q,2) = u;  % this stores the x-coord of first pixel of the RPF nucleus
                keepers_dil(q,3) = v + 1; % this stores the y-coord of first pixel of the RPF nucleus
              
            end
        end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BACKGROUND CALCULATIONS: Now repeat the skinny and dil process with the
% cells not expressing the Phospholocator to find the nucleii of
% "background" cells.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Skinny nucleus calculations 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert RFP_back image to Binary
level_RFP_back = graythresh(RFP_back);
bw_RFP_back = im2bw(RFP_back,level_RFP_back);

[labeled_RFP_back,numObjects_RFP_back] = bwlabel(bw_RFP_back,8); 

d_RFP_back = size(labeled_RFP_back);
q = 0;

% Generating a mask that passes nucleus
for i = 1:numObjects_RFP_back
        pix_temp_back= find(labeled_RFP_back == i);
        if (length(pix_temp_back)>100 & length(pix_temp_back) < 900)
            whammie = 1;
            for j = 1:length(pix_temp_back)
                if (find(danger_zone == pix_temp_back(j)))
                    whammie = 0;
                    break 
                end
            end
            if(whammie)  % shift all pixels up and over
                for k = 1:length(pix_temp_back)
                    pix_temp_back(k) = pix_temp_back(k) + shift_right * picture_size(1) - shift_up;
                end
                q = q + 1;
                pix_RFP_back(q) = {pix_temp_back};
                if (q==1)
                    mask_RFP_back_skinny = pix_RFP_back{q};
                else
                    mask_RFP_back_skinny = cat(1,mask_RFP_back_skinny, pix_RFP_back{q});
                end
                v = floor(pix_RFP_back{q}(1)/d_RFP_back(1));
                u = mod((pix_RFP_back{q}(1)-v*d_RFP_back(1)),d_RFP_back(1));
                keepers_RFP_back(q,1) = i;
                keepers_RFP_back(q,2) = pix_temp_back(1);
            end
        end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We dilate bw_RFP so that we can define a "no man's land" around the 
% nucleus that is counted as neither nucleus nor cytoplasm. This is to 
% prevent bleed through.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
% dilate the image twice
se90 = strel('line', 3, 90);
se0 = strel('line', 3, 0);

RFP_back_dil_temp = imdilate(bw_RFP_back, [se90 se0]);

RFP_back_dil = imdilate(RFP_back_dil_temp, [se90 se0]);

[labeled_dil_back,numObjects_dil_back] = bwlabel(RFP_back_dil,8); 

d_dil_back = size(labeled_dil_back);
q = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Finding the dilated nucleii that past our criteria, including avoiding
% the danger zone.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:numObjects_dil_back
        pix_dil_back_temp = find(labeled_dil_back == i);
        if (length(pix_dil_back_temp)>100)
            whammie = 1;
            for j = 1:length(pix_dil_back_temp)
                if (find(danger_zone == pix_dil_back_temp(j)))
                    whammie = 0;
                    break 
                end
            end
            if(whammie)  % shift all pixels up and over
                for k = 1:length(pix_dil_back_temp)
                    pix_dil_back_temp(k) = pix_dil_back_temp(k) + shift_right * picture_size(1) - shift_up;
                end
                q = q + 1;
                pix_dil_back(q) = {pix_dil_back_temp};
                pix_intensity_YFP_dil_back{q} = YFP(pix_dil_back{q});   
                
                if (q==1)
                    mask_RFP_back_dil = pix_dil_back{q};
                else
                    mask_RFP_back_dil = cat(1,mask_RFP_back_dil, pix_dil_back{q});
                end

                v = floor(pix_dil_back{q}(1)/d_dil_back(1));
                u = mod((pix_dil_back{q}(1)-v*d_dil_back(1)),d_dil_back(1));
                keepers_dil_back(q,1) = i;
                keepers_dil_back(q,2) = u;  % this stores the x-coord of first pixel of the RPF nucleus
                keepers_dil_back(q,3) = v + 1; % this stores the y-coord of first pixel of the RPF nucleus
              
            end
        end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we call BF_mask_final to find the BF outlines of the 
% background cells whose nucleii we just identified.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[pix_BF_back, correlate_RFP_BF_back, cells_back, pix_skinny_back_expanded, pix_dil_back_expanded] = BF_mask_final(BF_back, ...
    YFP_back, RFP_back, pix_RFP_back, pix_dil_back, 0, mask_RFP_back_skinny);

r_back = size(correlate_RFP_BF_back); 
number_cells_back = r_back(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we tally all the YFP in the background cells and take nuclear and cytoplasmic averages.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:r_back(1)
    pix_BF_total_final_back{i} = pix_BF_back{correlate_RFP_BF_back(i,1)}; % defines total cell pixels.
    pix_RFP_skinny_final_back{i} = pix_skinny_back_expanded{correlate_RFP_BF_back(i,2)}; % defines total nuclear pixels (skinny).
    pix_RFP_dil_final_back{i} = pix_dil_back_expanded{correlate_RFP_BF_back(i,3)};
    
    pix_BF_cytoplasmic_final_back{i} = pix_BF_total_final_back{i};
    
    for k = 1:length(pix_RFP_dil_final_back{i})
       pix_BF_cytoplasmic_final_back{i} = pix_BF_cytoplasmic_final_back{i}(pix_BF_cytoplasmic_final_back{i} ~= pix_RFP_dil_final_back{i}(k));
    end
    
    pix_YFP_nuc_total_final_back{i} = YFP_back(pix_RFP_skinny_final_back{i});
    pix_YFP_cyto_total_final_back{i} = YFP_back(pix_BF_cytoplasmic_final_back{i});
    
    YFP_nuc_total_average_back(i) = mean(double(pix_YFP_nuc_total_final_back{i}));
    YFP_nuc_total_max_back(i) = max(double(pix_YFP_nuc_total_final_back{i}));
    YFP_cyto_total_average_back(i) = mean(double(pix_YFP_cyto_total_final_back{i}));
    YFP_cyto_total_max_back(i) = max(double(pix_YFP_cyto_total_final_back{i}));
    
end

% calculating average total fluorescence

'calculated YFP nuclear background'
YFP_nuc_back = mean(YFP_nuc_total_average_back);

'calculated YFP cytosolic background'
YFP_cyto_back = mean(YFP_cyto_total_average_back);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We found that for all cells tested, the background fluorescence levels
% of the nucleii were similar to those of the cytoplasm. Thus, we define 
% one background fluorescence value, YFP_back
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

YFP_back = YFP_cyto_back;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here, we call BF_mask_final to find the BF outlines of the 
% cells expressing the Phospholocator whose nucleii we identified at the
% beginning of this file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[pix_BF, correlate_RFP_BF, cells, pix_RFP_expanded, pix_dil_expanded] = BF_mask_final(BF, YFP, RFP, pix_RFP,...
    pix_dil, YFP_back, mask_RFP_skinny);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We make lists of the Cell and nuclear pixels of cell identified in
% BF_mask_final
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r = size(correlate_RFP_BF); 
number_cells = r(1);

for i = 1:r(1)
    pix_BF_total_final{i} = pix_BF{correlate_RFP_BF(i,1)};
    pix_RFP_skinny_final{i} = pix_RFP_expanded{correlate_RFP_BF(i,2)};
    pix_RFP_dil_final{i} = pix_dil_expanded{correlate_RFP_BF(i,3)};
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This loop calculates YFP image intensity per cell in various
% compartments.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i = 1:r(1)
    pix_BF_cytoplasmic_final{i} = pix_BF_total_final{i}; 
    
    % This loop subtracts out the nuclear pixels by going though each RFP pixel, and removing it from 
    % pix_BF_cytoplasmic_final
    
    for k = 1:length(pix_RFP_dil_final{i})
       pix_BF_cytoplasmic_final{i} = pix_BF_cytoplasmic_final{i}(pix_BF_cytoplasmic_final{i} ~= pix_RFP_dil_final{i}(k));
    end
    
    % calculating cytosolic YFP intensities
    YFP_cytoplasmic_final{i} = YFP(pix_BF_cytoplasmic_final{i});
    
    cytosolic_YFP_total(i) = sum(double(YFP_cytoplasmic_final{i}) - YFP_cyto_back);
    
    cytosolic_YFP_average(i) = mean(double(YFP_cytoplasmic_final{i}) - YFP_cyto_back);
    
    cytosolic_YFP_mode(i) = mode(double(YFP_cytoplasmic_final{i}) - YFP_cyto_back);
    

    % calculating nuclear YFP intensities
    YFP_nuclear_final{i} = YFP(pix_RFP_skinny_final{i});
    
    nuclear_YFP_total(i) = sum(cytosolic_YFP_total(i));
    
    nuclear_YFP_average(i) = mean(double(YFP_nuclear_final{i}) - YFP_nuc_back);
    
    nuclear_YFP_mode(i) = mode(double(YFP_nuclear_final{i}) - YFP_nuc_back);
    
    
    % calculating total YFP intensities
    
    YFP_total_final(i) = cytosolic_YFP_total(i) + nuclear_YFP_total(i);
    
    cell_size(i) = length(YFP_cytoplasmic_final{i}) + length(YFP_nuclear_final{i});
    
    YFP_total_average(i) = YFP_total_final(i)/cell_size(i);
    
    % Calculating masks
    
    if (i==1)
            mask_BF_final = pix_BF_total_final{i};
            mask_RFP_skinny_final = pix_RFP_skinny_final{i};
            mask_RFP_dil_final = pix_RFP_dil_final{i};
        else
            mask_BF_final = cat(1,mask_BF_final, pix_BF_total_final{i});
            mask_RFP_skinny_final = cat(1,mask_RFP_skinny_final, pix_RFP_skinny_final{i});
            mask_RFP_dil_final = cat(1,mask_RFP_dil_final, pix_RFP_dil_final{i});
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This loop calculates the fluorescence averages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% calculating average total fluorescence
YFP_average_total = mean(YFP_total_average);
YFP_std_error_total = 2 * std(YFP_total_average)/sqrt(length(YFP_total_average));

% calculating stats for the average
cytoplasmic_nuclear_vector_temp = cytosolic_YFP_average ./ nuclear_YFP_average;

% Ensuring that all cells have positive fluorescences. Discard those that
% don't.

real_counter = 1;
for i = 1:length(cytoplasmic_nuclear_vector_temp)
    if (cytoplasmic_nuclear_vector_temp(i) >0)
        cytoplasmic_nuclear_vector(real_counter) = cytoplasmic_nuclear_vector_temp(i);
        real_counter = real_counter + 1;
    end
end

geometric_cytosolic_nuclear = exp(mean(log(cytoplasmic_nuclear_vector)));

std_average_total = std(cytoplasmic_nuclear_vector);
std_error_average_total = 2 * std_average_total/sqrt(length(cytoplasmic_nuclear_vector));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script BF_mask_final.m
%
% This script is called by CNR_final.m. The script finds cell boundaries and matches them with nucleii identified in CNR_main.m.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [pix_BF, correlate_RFP_BF, m, pix_RFP_expanded, pix_dil_expanded] = BF_mask_final(BF, YFP, RFP, pix_RFP, pix_dil,...
    background_YFP, mask_RFP_skinny)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Convert BF image to Binary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

level_BF = graythresh(BF);
bw_BF = im2bw(BF,level_BF); 

% figure
% imshow(bw_BF),title('BF mask, graythresh')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Getting rid of the large "background objects"
%
% Note:
% "labeled_..." is a matrix showing each object and where it is located
% "numObjects_..." gives you the number of objects in labeled.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[labeled_BF_prelim,numObjects_BF_prelim] = bwlabel(bw_BF,8); 


BF_cells_only = bw_BF;

for i = 1:numObjects_BF_prelim
    pix_BF_prelim(i) = {find(labeled_BF_prelim == i)};
    if (length(pix_BF_prelim{i}) >= 6300)
        BF_cells_only(pix_BF_prelim{i}) = 0;
    end
end

% figure 
% imshow(BF_cells_only),title('BF mask, got rid of "big cells"')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Dilate BF image slightly to improve cell boundary connectivity.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dilate the image
se90 = strel('line', 3, 90);
se0 = strel('line', 3, 0);

BF_dil = imdilate(BF_cells_only, [se90 se0]);

% figure
% imshow(BF_dil),title('BF mask, dilated')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fill in the holes in the BW image with connectivity = 8
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bw_fill_BF_8 = imfill(BF_dil,8,'holes');

% figure
% imshow(bw_fill_BF_8),title('BF mask, pre-screen')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Label cells and gather info on them
%
% Notes:
% "labeled_..." is a matrix showing each object and where it is located
% "numObjects_..." gives you the number of objects in labeled.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[labeled_BF,numObjects_BF] = bwlabel(bw_fill_BF_8,8); 

bf_data = regionprops(labeled_BF,'all');

perimeter_measure = [bf_data.Perimeter];
solidity_measure = [bf_data.Solidity];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select BF cells to quantify. We choose BF cells that:
% -- are whole cells (using solidity measure)
% -- have average YFP fluorescence that is 2% above background
% -- are larger than 400 pixels
% -- match up with atleast one nucleus from RFP images.
% 
% Some of our mitotic cells have two nucleii (mother and daughter),
% and so we make sure to associate both nucleii with the same cell
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d_BF = size(labeled_BF);
d_RFP = size(pix_RFP);
d_dil = size(pix_dil);

pix_RFP_expanded = pix_RFP;
pix_dil_expanded = pix_dil;

m = 0;        

clear correlate_RFP_BF
clear pix_BF
counter_solid = 0;
counter_RFP = 0;
counter_dil = 0;
counter_all = 0;


for i = 1:numObjects_BF
    nuc_found = 0;
    pix_BF(i) = {find(labeled_BF == i)};
    if (solidity_measure(i) > 0.8 & mean(double(YFP(pix_BF{i}))) >= 1.02 * background_YFP & length(pix_BF{i})>400) % Testing solidity and YFP levels
        counter_solid = counter_solid + 1;
        for k = 1:d_RFP(2)
            clear g;
            g1 = find(pix_BF{i} == pix_RFP{k}(1));
            g2 = find(pix_BF{i} == pix_RFP{k}(50));
            g3 = find(pix_BF{i} == pix_RFP{k}(100));
            g4 = find(pix_BF{i} == pix_RFP{k}(25));
            g5 = find(pix_BF{i} == pix_RFP{k}(75));
            g6 = find(pix_BF{i} == pix_RFP{k}(85));

            if (g1|g2|g3|g4|g5|g6)
                counter_RFP = counter_RFP + 1;
                clear n;
                for n = 1:d_dil(2)
                    h1 = find(pix_dil{n} == pix_RFP{k}(50));                    % matching up RFP_skinny with RFP_dil
                    if(h1)
                        counter_dil = counter_dil + 1;
                        if (length(pix_BF{i}) > 1.1 * length(pix_dil{n}))
                            nuc_found = nuc_found + 1;
                            BF_skinny_dil_triplet(nuc_found) = {[i, k, n]};
                            counter_all = counter_all + 1;
                            break   % Don't need to search through anymore dil once found the matching one.
                        else
                            break   % Don't need to search through anymore dil once found the matching one.
                            
                        end
                    end
                end
            end
        end
        if(nuc_found == 1)
            'found a single nucleus cell'
            m = m +1;
            correlate_RFP_BF(m,1) = BF_skinny_dil_triplet{nuc_found}(1); % stores first element as the BF tag.
            correlate_RFP_BF(m,2) = BF_skinny_dil_triplet{nuc_found}(2); % stores second element as the RFP image from pix_RFP.
            correlate_RFP_BF(m,3) = BF_skinny_dil_triplet{nuc_found}(3); % stores third element as the dilated RFP image from pix_dil.
        elseif(nuc_found > 1)
            'found a multinucleate cell'
            total_nuc_skinny =  pix_RFP{BF_skinny_dil_triplet{1}(2)};
            total_nuc_dil    =  pix_dil{BF_skinny_dil_triplet{1}(3)};
            BF_skinny_dil_triplet{1};
            for q = 2:nuc_found
                BF_skinny_dil_triplet{q};
                skinny_nuc = BF_skinny_dil_triplet{q}(2) ;
                dil_nuc = BF_skinny_dil_triplet{q}(3); 
                total_nuc_skinny = cat(1,pix_RFP{skinny_nuc}, total_nuc_skinny);
                total_nuc_dil = cat(1,pix_dil{dil_nuc}, total_nuc_dil);
            end
            temp_length_skinny = size(pix_RFP_expanded);
            temp_length_dil = size(pix_dil_expanded);
            pix_RFP_expanded(temp_length_skinny(2) + 1) = {total_nuc_skinny};
            pix_dil_expanded(temp_length_dil(2) + 1) = {total_nuc_dil};
            m = m +1;
            correlate_RFP_BF(m,1) = BF_skinny_dil_triplet{nuc_found}(1); % stores first element as the BF tag.
            correlate_RFP_BF(m,2) = temp_length_skinny(2) + 1; % stores second element as the RFP image from pix_RFP_extended.
            correlate_RFP_BF(m,3) = temp_length_dil(2) + 1; % stores third element as the dilated RFP image from pix_dil_extended.
            
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If no BF cells pass the above tests, we run a remedial loop
% such that the program doesn't crash.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if (m == 0)
    'no cells detected: running BF consolation loop'
    for i = 1:numObjects_BF
        pix_BF(i) = {find(labeled_BF == i)};
        if (solidity_measure(i) > 0.8 | (length(pix_BF{i})> 700 & length(pix_BF{i}) < 1000))
            for k = 1:d_RFP(2)
                clear g;
                g = find(pix_BF{i} == pix_RFP{k}(100));
                if (g)
                    clear n;
                    for n = 1:d_dil(2)
                        h = find(pix_BF{i} == pix_dil{n}(100));
                        if (h & length(pix_BF{i}) > 1.1 * length(pix_dil{n}))
                                m = m +1;
                                correlate_RFP_BF(m,1) = i; % stores first element as the BF tag.
                                correlate_RFP_BF(m,2) = k; % stores second element as the RFP image from pix_RFP.
                                correlate_RFP_BF(m,3) = n; % stores third element as the dilated RFP image from pix_dil.
                                break
                        end
                    end
                end
            end
        end   
    end
end
