clear;
close all;
clc;
%第一步 读取图像，转为灰度图
rgb=imread('gecko.bmp');
I = rgb2gray(rgb);
%第二步，转为梯度图
hy = fspecial('prewitt');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure
imshow(gradmag,[]);title('sobel梯度图');
% 第三步，标记前景对象
%使用两组模板，每组模板有两个大小不同的模板，使得大区域与小区域都标记为目标
se = strel('disk', 8);
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
se1 = strel('disk', 3);
Ie1 = imerode(I, se1);
Iobr1 = imreconstruct(Ie1, I);
Iobrd1 = imdilate(Iobr1, se1);
Iobrcbr1 = imreconstruct(imcomplement(Iobrd1), imcomplement(Iobr1));
Iobrcbr1 = imcomplement(Iobrcbr1);
Nse2 = strel('disk', 9);
Ie2 = imerode(I, Nse2);
Iobr2 = imreconstruct(Ie2, I);
Iobrd2 = imdilate(Iobr2, Nse2);
Iobrcbr2 = imreconstruct(imcomplement(Iobrd2), imcomplement(Iobr2));
Iobrcbr2 = imcomplement(Iobrcbr2);
figure;
imshow(Iobrcbr);title('111');
figure;
imshow(Iobrcbr1);title('222');
fgm = imregionalmax(Iobrcbr);
fgm1 = imregionalmax(Iobrcbr1);
Nfgm1 = imregionalmax(Iobrcbr2);
figure
imshow(fgm);title('aaaa');
I2 = I;
I2(Nfgm1|fgm1) = 255;
figure
imshow(I2), title('前景标定图像')
%去除斑点
se2 = strel(ones(2,2));
fgm2 = imclose(fgm|fgm1, se2);
fgm3 = imerode(fgm2, se2);
fgm4 = bwareaopen(fgm3, 5);
Nfgm2=imclose(fgm1|Nfgm1, se2);
Nfgm3 = imerode(Nfgm2, se2);
Nfgm4 = bwareaopen(Nfgm3, 5);
I3 = I;
I3(Nfgm4) = 255;
figure
imshow(I3);title('去除斑点前景标定图');
%第四步，标记背景区域
%bw = im2bw(Iobrcbr, 0.42);
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));
bw1 = im2bw(Iobrcbr1, graythresh(Iobrcbr1));
bw2=adaptivethresh(Iobrcbr1,23,9);
lvbo1=strel(ones(2,2));
bw2=imerode(bw2,lvbo1);
bw2=imfill(bw2,'hole');
figure
imshow(bw2);title('自适应二值化法')
D = bwdist(bw);
DL = watershed(D);
D1 = bwdist(bw1);
DL1 = watershed(D1);
D2 = bwdist(bw2);
DL2 = watershed(D2);
bgm = DL == 0;
bgm1 = DL1 == 0;
bgm2 = DL2 == 0;
figure;
imshow(bw2);title('阈值分割图')
figure;
imshow(bgm2, []);title('分水岭分割')
%第五步：计算分割函数的分水岭变换
gradmag2 = imimposemin(gradmag, bgm2 |fgm4);
gradmag3 = imimposemin(gradmag, bgm2 |Nfgm4);
Logicgradmag2=(gradmag2==-inf);
gradmag2(Logicgradmag2)=gradmag3(Logicgradmag2);
%gradmag2=gradmag2*0.5+gradmag3*0.5;
figure
imshow(gradmag2,[]);title('分割处理后的梯度图');
L = watershed(gradmag2);
Stats=regionprops(L,'Area');
Stats1=struct2cell (Stats);
Stats2=cell2mat(Stats1);
[m n]=size(Stats2);
[idx, Codebook] = kmeans(Stats2(1,2:n)',2);
[a b]=sort(idx,'descend');
[m1 n1]=size(L);
zero1=zeros(size(L));
zero2=zero1;
numberBigArea=sum(a==size(Codebook,1));
for i=1:numberBigArea
c(i)=b(i)+1;
zero1(L==c(i))=2;
end
zero2(L==1)=1;
L1=zero1+zero2;
%第六步：查看结果
I4 = I;
%I4(imdilate(L == 0, ones(3, 3)) | bgm2|bgm1 | fgm4) = 255;
% Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
Lrgb = label2rgb(L1, 'jet', 'w', 'shuffle');
figure
imshow(I4);
figure
imshow(Lrgb)
title('Colored watershed label matrix (Lrgb)')
figure
imshow(I)
hold on
himage = imshow(Lrgb);
himage.AlphaData = 0.3;
