%%obtain the gray scale distribution of the slice
II=imread('C:\Users\LYUKAI\Desktop\pretreated slice.tif');%%%input the pre-treated slice
I=II(1:1040,11:1270);%%choose the effective area
[H,D]=imhist(II);
for i=1:255
    G(i)=length(find(I<i));
end
G=G';
G=G./(length(find(I<255)))*100;
F=D(1:255);
plot(F,G,'k','LineWidth',3);
xlabel('Gray value','Fontname', 'Times New Roman','FontSize',20);
ylabel('Area segmented (%)','Fontname', 'Times New Roman','FontSize',20);
axis([0 255 0 100]);%%draw the gray scale distribution curve
title('Gray-scale distribution','Fontname', 'Times New Roman','FontSize',24);
hold on;
H=H./(length(find(I<255)))*100;
H=H(1:255);


%%linear fitting 
xx1=D(8:23);
yy1=G(8:23);
xx2=D(61:69);
yy2=G(61:69);%%%find two linear parts
K1=polyfit(xx1,yy1,1);
K2=polyfit(xx2,yy2,1);%%%fitting coefficients of the two parts

for i=1:100
    yyy1(i)=K1(1)*i+K1(2);
    yyy2(i)=K2(1)*i+K2(2);
end
x=1:100;
plot(x,yyy1,'r','LineWidth',2);
axis([0 255 0 100]);
hold on;
plot(x,yyy2,'r','LineWidth',2);
axis([0 255 0 100]);
hold;

%%calculation of the intersection point
syms x y;
s=solve(y==K1(1)*x+K1(2),y==K2(1)*x+K2(2),x,y);
X=double(s.x);
Y=double(s.y);%%intersection point and the X value will be used for final segmentation

%%change the grayscale image into binary image
figure(2);
a=find(I<X);
b=find(I>=X);
I(a)=0;
I(b)=255;
imshow(I);

