%%Cell Segmentation v1.0%%

%Step 1: Load all images (Replace with respective filenames accordingly)
image1 = double(imread(strcat('Brightfield.tif')));
image2 = double(imread(strcat('DAPI.tif')));
image3 = double(imread(strcat('QUAD.tif')));
image4 = double(imread(strcat('HER2.tif')));
image5 = double(imread(strcat('ERPR.tif')));

%Step 2: Identify potential cell candidates in brightfield
image1_m = mat2gray(image1);

%%%%Search for edges
[~,threshold] = edge(image1_m,'Prewitt');
fudgeFactor = 0.22; %Adjust according to results
BWs = edge(image1_m,'Prewitt',threshold * fudgeFactor,'nothinning');

%%%%Dilate the selected regions
se90 = strel('line',3,90);
se0 = strel('line',3,0);
BWsdil = imdilate(BWs,[se90 se0]);

%%%%Fill interrior gaps
BWdfill = imfill(BWsdil,'holes');

%%%%Remove border cells
BWfinal = imclearborder(BWdfill,4);

%%%%Smooth the cells
seD = strel('diamond',1);
BWfinal = imerode(BWfinal,seD);

%%%%Define the perimeter of all cells
bw3 = imopen(BWfinal, ones(5,5));
bw4 = bwareaopen(bw3, 100);
bw4_perim = bwperim(bw4);

%Step 3: Identify cells with nucleus using DAPI image
image2_m = mat2gray(image2);

%%%%Identify nucleus
mask_em = imextendedmax(image2_m, 0.1);
mask_em = imclose(mask_em, ones(7,7));
mask_em = imfill(mask_em, 'holes');
mask_em = bwareaopen(mask_em, 20);

%%%%Correlate segmented cells with nucleus
images1_c = imcomplement(image1_m);
images1_mod = imimposemin(images1_c, ~bw4 | mask_em);
L = watershed(images1_mod);
LabeledCells = bwlabel(L);
mask = imbinarize(LabeledCells, 1);
mask_perim = bwperim(mask);
mask_perim = imdilate(mask_perim, strel('disk',1));
mask_final = imdilate(mask_perim, strel('disk',1));
mask_final = imfill(mask_final,'holes');

%%%%Show the end results of the cell segmentation
overlay_cell = imoverlay(image1, mask_final, [1 .3 .3]);
figure; imshow(overlay_cell);

%Step 4: Retrieve the mean intensity of all segmented cells in each fluorescence images
LabeledCells = logical(L);
CellMeasurements_DAPI = regionprops(LabeledCells, image2, 'MeanIntensity');
CellMeasurements_QUAD = regionprops(LabeledCells, image3, 'MeanIntensity');
CellMeasurements_HER2 = regionprops(LabeledCells, image4, 'MeanIntensity');
CellMeasurements_ERPR = regionprops(LabeledCells, image5, 'MeanIntensity');
Cell_No = size(CellMeasurements_DAPI, 1);

CellMeasurements_DAPI = CellMeasurements_DAPI(2:Cell_No);
CellMeasurements_QUAD = CellMeasurements_QUAD(2:Cell_No);
CellMeasurements_HER2 = CellMeasurements_HER2(2:Cell_No);
CellMeasurements_ERPR = CellMeasurements_ERPR(2:Cell_No);
Cell_No = Cell_No - 1;
