clc;

clear all;

close all;

[fna,pna]=uigetfile({'.jpg;.bmp;*.png'},'Select file');

fanme=[pna fna];

im = double(imread(fanme));

im=imresize(im,[200 300]);

figure;imshow(uint8(im));title('Input image');

tic

%  Lab conversion and color thresholding

im=im/255;

R=im(:,:,1);

G=im(:,:,2);

B=im(:,:,3);

% Set a threshold

T = 0.008856;

[M, N,plane] = size(R);

s = M * N;

RGB = [reshape(R,1,s); reshape(G,1,s); reshape(B,1,s)];

% RGB to XYZ

MAT = [0.412453 0.357580 0.180423;

       0.212671 0.715160 0.072169;

       0.019334 0.119193 0.950227];

XYZ = MAT * RGB;

% Normalize for D65 white point

X = XYZ(1,:) / 0.950456;

Y = XYZ(2,:);

Z = XYZ(3,:) / 1.088754; 

XT = X > T;

YT = Y > T;

ZT = Z > T;

Y3 = Y.^(1/3); 

fX = XT .* X.^(1/3) + (~XT) .* (7.787 .* X + 16/116);

fY = YT .* Y3 + (~YT) .* (7.787 .* Y + 16/116);

fZ = ZT .* Z.^(1/3) + (~ZT) .* (7.787 .* Z + 16/116);

L = reshape(YT .* (116 * Y3 - 16.0) + (~YT) .* (903.3 * Y), M, N);

a = reshape(500 * (fX - fY), M, N);

b = reshape(200 * (fY - fZ), M, N);

  L = cat(3,L,a,b);

figure;imshow(L);title('L*a*b converted image'); 

% color thresholding

      out2 = colorthreshold(lab2uint8(L),'otsu');

      figure, imshow(uint8(out2));

      title('Color Thresholded image');

     OUCT=uint8(out2); 

     ReLAb=zeros(size(L));

     ReRGB=zeros(size(L));

Th=100;

for i=1:M

    for j=1:N

        if OUCT(i,j,1)<200 && OUCT(i,j,1)>90

            ReLAb(i,j,:)=L(i,j,:);

            ReRGB(i,j,:)=im(i,j,:);

        end

    end

end

   figure, imshow(ReLAb);title('L*a*b Mask image');

   figure, imshow(ReRGB);title('RGB Mask image');

%----RGB 2 Gray conversion

Grayim=rgb2gray(ReRGB);

figure, imshow(Grayim);title('Gray Mask image');

L=imadjust(Grayim);

H = histeq(Grayim);

%% Brighten most of the details in the image except the nucleus

R1=imadd(L,H);

%% Highlight all the objects and its borders in the image including the cell nucleus

R2 = imsubtract(R1,H);

%% Remove almost all the other blood components while retaining the nucleus with minimum affect of distortion on the nucleus part of the white blood cell.

% R3=imadd(R1,R2);

figure;

subplot(221);

imshow(Grayim);title('Gray image')

subplot(222);

imshow(H);title('Enhanced image')

subplot(223);

imshow(R1);title('Image except the nucleus')

subplot(224);

imshow(R2);title('Highlighted all objects')

toc

%check histogram

figure;

hi=imhist(R2);

hi1=hi(1:2:256);

horz=1:2:256;

bar(horz,hi1);

axis([0 255 0 1400]);

set(gca, 'xtick', 0:50:255);

set(gca, 'ytick', 0:2000:15000);

 xlabel('Gray level' );

ylabel('No of pixels' );

title('Histogram before opening the image');

%=====================

% Otsu thresholding   

   rk=otsu(R2);

  imthO=(Grayim>rk); 

figure, imshow(imthO);title('Otsu Binary Image');

%implements a 3-by-3 average filter

R3 = average_filter(imthO);

R31=double(im2bw(R3));

figure; imshow(R31);title('Average Filtered image');

hy = fspecial('sobel');

hx = hy';

Iy = imfilter(double(R3), hy, 'replicate');

Ix = imfilter(double(R3), hx, 'replicate');

gradmag = sqrt(Ix.^2 + Iy.^2);

figure

imshow(gradmag,[]), title('Gradient magnitude')

% Finally we are ready to compute the watershed-based segmentation.

L = watershed(gradmag);

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');

figure

imshow(Lrgb)

title('Watershed segmentation')

%  Mark the Foreground Objects

se = strel('disk', 2);

Io = imopen(Grayim, se);

figure

imshow(Io), title('Opening (Io)')

Ie = imerode(Grayim, se);

Iobr = imreconstruct(Ie, Grayim);

figure

imshow(Iobr), title('Opening-by-reconstruction (Iobr)')

Ioc = imclose(Io, se);

figure

imshow(Ioc), title('Opening-closing (Ioc)')

Iobrd = imdilate(Iobr, se);

Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));

Iobrcbr = imcomplement(Iobrcbr);

figure

imshow(Iobrcbr), title('Opening-closing by reconstruction (Iobrcbr)')

fgm = Iobrcbr>0.1;

figure

imshow(fgm), title('Regional maxima of opening-closing by reconstruction (fgm)')

I2 = Grayim;

I2(fgm) = 255;

figure

imshow(I2), title('Regional maxima superimposed on original image (I2)')

% se2 = strel('disk', 1);

% % fgm2 = imclose(fgm, se2);

% fgm = imerode(fgm, se2);

% This procedure tends to leave some stray isolated pixels that must be removed. You can do this using bwareaopen, which removes all blobs that have fewer than a certain number of pixels.

% getin=inputdlg('Enter the Area opening dimension(small objects-100 & Large objects-250)');

% aREAIN=str2double(getin{1,1});

aREAIN=50;

fgm4 = bwareaopen(fgm,aREAIN);

I3 = Grayim;

I3(fgm4) = 255;

figure

imshow(I3)

title('Modified regional maxima superimposed on original image (fgm4)')

figure

imshow(fgm4)

title('After Eliminating Platelets')

% 

% Compute Background Markers

% Now you need to mark the background. In the cleaned-up image, Iobrcbr, the dark pixels belong to the background, so you could start with a thresholding operation.

D = bwdist(fgm4);

DL = watershed(D);

bgm = DL == 0;

figure

imshow(bgm), title('Broken Watershed Nucleus')

%% eliminate border touching segments

hy = fspecial('sobel');

hx = hy';

Iy = imfilter(double(fgm4), hy, 'replicate');

Ix = imfilter(double(fgm4), hx, 'replicate');

gradmag = sqrt(Ix.^2 + Iy.^2);

figure

imshow(gradmag,[]), title('Gradient magnitude')

% Finally we are ready to compute the watershed-based segmentation.

L = watershed(gradmag);

% Step 6: Visualize the Result

% One visualization technique is to superimpose the foreground markers, background markers, and segmented object boundaries on the original image. You can use dilation as needed to make certain aspects, such as the object boundaries, more visible. Object boundaries are located where L == 0.

I4 = Grayim;

I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;

figure

imshow(I4)

title('Markers and object boundaries superimposed on original image ')

Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');

figure

imshow(Lrgb)

title('Colored watershed label matrix (Lrgb)')

%% Final 

figure;

imshow(im);

hold on

boundaries = bwboundaries(fgm4);    

numberOfBoundaries = size(boundaries);

for k = 1 : numberOfBoundaries

    thisBoundary = boundaries{k};

    plot(thisBoundary(:,2), thisBoundary(:,1), 'r', 'LineWidth', 1.5);

end

hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=imthO;

aa=imresize(a,[256 256]);

u_bw_filename = aa;

b=im;

 bb=imresize(b,[256 256]);

u_GT_filename = im2bw(bb);

u_GT = [((u_GT_filename)) > 0 ];

u_bw = [((u_bw_filename)) > 0 ];

temp_obj_eval = objective_evaluation_core(u_bw, u_GT);

disp('PSNR value --');

disp(temp_obj_eval.PSNR);