小白学习图像处理——分水岭算法

主要参考《数字图像处理 原理与实践（matlab版）》一书

用到的函数介绍：

imreconstruct（）函数的功能是对图形形态修饰。

imreconstruct（）书写主要格式为

IM = imreconstruct(MARKER,MASK) 

MARKER——标记，标记和掩码可以是两个灰度图像或两个二进制图像具有相同的大小。标记必须是相同的大小。

MASK——面膜，作为面膜，其元素必须小于或等于面具的相应元素

https://blog.csdn.net/qq_26460507/article/details/55192961

功能：转换标记矩阵到RGB图像

用法：RGB = label2rgb(L) 

输入：L——标记矩阵(可由labelmatrix, bwlabel, bwlabeln, watershed返回) 

输出：RGB——彩色图像

注：根据L的数值对应，默认对应到colormap(jet)的色彩，返回RGB矩阵 

RGB = label2rgb(L, map)

输入：L——标记矩阵(可由labelmatrix, bwlabel, bwlabeln, watershed返回)、map——颜色映射表 

输出：RGB——彩色图像

注：map为n*3的矩阵，可以通过MATLAB的colormap函数来返回，比如colormap('jet')等。也可以根据要求自己定义。默认为colormap(jet)。 

RGB = label2rgb(L, map, zerocolor) defines the RGB color of the elements labeled 0 (zero) in the input label matrix L. As the value of zerocolor, specify an RGB triple or one of the strings listed in this table.

输入：L——标记矩阵(可由labelmatrix, bwlabel, bwlabeln, watershed返回)、map——颜色映射表、zerocolor——对应于标记0的颜色 

输出：RGB——彩色图像 b

注：zerocolor可以取值如表1.1，默认为[1 1 1]，即白色。 

表1.1 取值列表 

Value	Color
'b'	蓝色
'c'	蓝绿色
'g'	绿色
'k'	黑色
'm'	洋红色
'r'	红色
'w'	白色
'y'	黄色

RGB = label2rgb(L, map, zerocolor, order) 

输入：L——标记矩阵(可由labelmatrix, bwlabel, bwlabeln, watershed返回)、map——颜色映射表、zerocolor——对应于标记0的颜色、order——标记矩阵和颜色映射表对应方式
输出：RGB——彩色图像 

注：order默认为noshuffle，即根据L的数值来对应颜色。另外可以取值为shuffle，说明使用伪随机方式来对应。

https://blog.csdn.net/u012116229/article/details/44945205

可以使用MATLAB图形处理工具箱中的函数来确定图像的高亮度区域和低亮度区域。

imregionalmax函数和imregionalmin函数确定所有的极大值和极小值
imextendedmax函数和imextendedmin函数确定大于某阈值的极大值和小于某阈值的极小值
所以都是为了确定极大值和极小值。

这些函数把灰度图像作为输入参数，而把二值图像作为输出参数。在输出的二值图像中，局部极大值和局部极小值设为1，其他值设为0.

imextendedmax函数和imextendedmin函数，以imextendedmax为例，是首先将输入的阈值减去再用imregionalmax来做操作的。

https://blog.csdn.net/qing101hua/article/details/45559377

cat：用来联结数组

用法：C = cat(dim, A, B)       按dim来联结A和B两个数组。

      C = cat(dim, A1, A2, A3, ...)    按dim联结所有输入的数组。

      e.g.

      a=cat(3,A,B) 左括号后的3表示构造出的矩阵维数；在新的矩阵中第1、2维就是A和B这两个矩阵的行数和列数，第3维是A和B这两个矩阵的矩阵个数，即为2

cat(2, A, B)相当于[A, B];

cat(1, A, B)相当于[A; B].

https://www.cnblogs.com/saliency/archive/2014/03/11/3593308.html

matlab函数bwareaopen──删除小面积对象

格式：BW2 = bwareaopen(BW,P,conn)
作用：删除二值图像BW中面积小于P的对象，默认情况下使用8邻域。

https://blog.csdn.net/u013228046/article/details/40651663

bwdist函数用于计算元素之间的距离。

https://blog.csdn.net/hduxiejun/article/details/53011863

imimposemin函数

bwdist函数用于计算元素之间的距离。 
功能：施加极小值 
用法：I2 = imimposemin(I,BW) 
使用形态重构修改强度图像I，使得它在BW非零的地方只有区域极小值。BW是一个二值图像，和I的尺寸相同。 
默认情况下，imimposemin对二维图像使用8连通域，对三维图像使用26连通域。对于更高维数的情况，imimposemin使用conndef(ndims(I),’minimum’)。 
I2 = imimposemin(I,BW,conn) 

https://blog.csdn.net/jiejianmin2666/article/details/78148533

Matlab中imshow（）函数的使用
    问题：在使用imshow(A)显示一张灰度图片时，显示出的是一张纯白的图片。。（A为double类型的图像矩阵）
    原因：在matlab中，为了保证精度，经过了运算的图像矩阵A其数据类型会从unit8型变成double型。imshow()显示图像时对double型是认为在0~1范围内，即大于1时都是显示为白色。imshow显示uint8型时是0~255范围。
    解决：使用imshow(A,[])，即可把图像矩阵A显示为正常的灰度图像。

    原理：imshow(A,[])是将A的最大值（max(A)）和最小值（min(A)）分别作为纯白(255)和纯黑(0)，中间的K值相应地映射为0到255之间的标准灰度值，这样就可以正常显示了。。。，相当于将double型的矩阵A拉伸成为了0-255的uint8型的矩阵，因此就可以正常显示.

matlab代码

clear
i=imread('D:\Documents\Desktop\egg.jpg');%读入彩色图像
ig=rgb2gray(i);%将图像转为灰度图像
hy=fspecial('sobel');%sobel 是梯度算子，用sobel算子构建计算模板
hx=hy';%获得垂直方向和水平方向两个方向的检测模板
Iy=imfilter(double(ig),hy,'replicate');%计算垂直方向上的梯度
Ix=imfilter(double(ig),hx,'replicate');%计算水平方向上的梯度
subplot(221),imshow(i),title('original image');
subplot(222),imshow(ig),title('gray image');
subplot(223),imshow(Ix),title('filter with hx');
subplot(224),imshow(Iy),title('filter with hy');

显示结果

gradmag = sqrt (Ix.^2+Iy.^2);%求得梯度向量的模值
L=watershed(gradmag);%直接对梯度幅值图使用matlab自带的分水岭算法
figure(2);
subplot(121),imshow(gradmag),title('gradient magnitude');
subplot(122),imshow(L),title('use watershed');

显示结果

se=strel('disk',12);%构建开运算的结构元素，圆-r=12
Ie=imerode(ig,se);%进行腐蚀
Iobr = imreconstruct(Ie,ig);%对腐蚀后的图像进行重建.即将腐蚀后的结果作为标记
Iobrc = imcomplement(Iobr);%对腐蚀后的图像进行求补

Iobrd = imdilate(Iobr,se);%进行膨胀
Iobrdc = imcomplement(Iobrd);%对膨胀后的图进行求补，函数imcomplent功能： 对图像数据进行取反运算（实现底片效果）

Iobrcbr = imreconstruct (Iobrdc,Iobrc);%对腐蚀&膨胀的求补的结果进行重建
Iobrcbrc = imcomplement (Iobrcbr);%对重建结果进行求补

figure(3);
subplot(331),imshow(Ie),title('直接进行腐蚀');
subplot(332),imshow(Iobr),title('腐蚀重建');
subplot(333),imshow(Iobrc),title('腐蚀重建后求补');
subplot(334),imshow(Iobrd),title('腐蚀重建后膨胀，即进行开运算');
subplot(335),imshow(Iobrdc),title('膨胀后进行求补');
subplot(336),imshow(Iobrcbr),title('对腐蚀和膨胀的求补结果进行重建');
subplot(337),imshow(Iobrcbrc),title('对腐蚀和膨胀的求补结果进行重建后求补');

fgm = imregionalmax(Iobrcbrc);
ir=i(:,:,1);%提取原彩色RGB图像中的红色分量
ig=i(:,:,2);%提取原彩色RGB图像中的绿色分量
ib=i(:,:,3);%提取原彩色RGB图像中的蓝色分量

ir(fgm)=255;
ib(fgm)=0;
ib(fgm)=0;

i2=cat(3,ir,ig,ib);%联结成3维数组

figure(4);
subplot(231),imshow(fgm,[]),title('检测局部最大值');%原图最亮的地方
subplot(232),imshow(ir),title('红色分量');
subplot(233),imshow(ig),title('绿色分量');
subplot(234),imshow(ib),title('蓝色分量');
subplot(235),imshow(i2),title('局部极大值叠加图像');
subplot(236),imshow(fgm),title('检测局部最大值,没[]');%结果表示和有[]的没啥区别

se2=strel(ones(15,15));%构建结构元素，15×15的1矩阵
fgm2=imclose(fgm,se2);%进行闭运算，即先膨胀再腐蚀
fgm3=imerode(fgm2,se2);%载进行腐蚀
fgm4=bwareaopen(fgm3,400);%删除小于400像素的连通块

figure(5);
subplot(221),imshow(fgm),title('检测局部最大值');
subplot(222),imshow(fgm2),title('对局部最大值进行闭运算');
subplot(223),imshow(fgm3),title('闭运算后进行腐蚀');
subplot(224),imshow(fgm4),title('删除小于400像素的连通块');

bw=im2bw(Iobrcbrc,graythresh(Iobrcbrc));%对Iobrcbrc进行二值化，二值化的阈值自动获取
D=bwdist(bw);
DL=watershed(D);
bgm=DL==0;
gradmag2=imimposemin(gradmag,bgm|fgm4);%bgm和fgm4进行或运算，使得梯度幅值图像gradmag上只在或运算结果为非零的位置上有局部最小值

L=watershed(gradmag2);
%第一种显示方法
Lrgb=label2rgb(L,'jet','w','shuffle');
figure(6);
subplot(121),imshow(Lrgb),title('分水岭显示结果1');
%第二种显示方法
subplot(122),imshow(i,[]),title('分水岭显示结果2');
hold on
himage=imshow(Lrgb);
set(himage,'AlphaData',0.3);

完整的程序

clear
i=imread('D:\Documents\Desktop\egg.jpg');%读入彩色图像
ig=rgb2gray(i);%将图像转为灰度图像
hy=fspecial('sobel');%sobel 是梯度算子，用sobel算子构建计算模板
hx=hy';%获得垂直方向和水平方向两个方向的检测模板
Iy=imfilter(double(ig),hy,'replicate');%计算垂直方向上的梯度
Ix=imfilter(double(ig),hx,'replicate');%计算水平方向上的梯度
subplot(221),imshow(i),title('original image');
subplot(222),imshow(ig),title('gray image');
subplot(223),imshow(Ix),title('filter with hx');
subplot(224),imshow(Iy),title('filter with hy');

%*********************************************************
gradmag = sqrt (Ix.^2+Iy.^2);%求得梯度向量的模值
L=watershed(gradmag);%直接对梯度幅值图使用matlab自带的分水岭算法
figure(2);
subplot(121),imshow(gradmag),title('gradient magnitude');
subplot(122),imshow(L),title('use watershed');

%*********************************************************
se=strel('disk',12);%构建开运算的结构元素，圆-r=12
Ie=imerode(ig,se);%进行腐蚀
Iobr = imreconstruct(Ie,ig);%对腐蚀后的图像进行重建.即将腐蚀后的结果作为标记
Iobrc = imcomplement(Iobr);%对腐蚀后的图像进行求补

Iobrd = imdilate(Iobr,se);%进行膨胀
Iobrdc = imcomplement(Iobrd);%对膨胀后的图进行求补，函数imcomplent功能： 对图像数据进行取反运算（实现底片效果）

Iobrcbr = imreconstruct (Iobrdc,Iobrc);%对腐蚀&膨胀的求补的结果进行重建
Iobrcbrc = imcomplement (Iobrcbr);%对重建结果进行求补

figure(3);
subplot(331),imshow(Ie),title('直接进行腐蚀');
subplot(332),imshow(Iobr),title('腐蚀重建');
subplot(333),imshow(Iobrc),title('腐蚀重建后求补');
subplot(334),imshow(Iobrd),title('腐蚀重建后膨胀，即进行开运算');
subplot(335),imshow(Iobrdc),title('膨胀后进行求补');
subplot(336),imshow(Iobrcbr),title('对腐蚀和膨胀的求补结果进行重建');
subplot(337),imshow(Iobrcbrc),title('对腐蚀和膨胀的求补结果进行重建后求补');

fgm = imregionalmax(Iobrcbrc);
ir=i(:,:,1);%提取原彩色RGB图像中的红色分量
ig=i(:,:,2);%提取原彩色RGB图像中的绿色分量
ib=i(:,:,3);%提取原彩色RGB图像中的蓝色分量

ir(fgm)=255;
ib(fgm)=0;
ib(fgm)=0;

i2=cat(3,ir,ig,ib);%联结成3维数组

figure(4);
subplot(231),imshow(fgm,[]),title('检测局部最大值');%原图最亮的地方
subplot(232),imshow(ir),title('红色分量');
subplot(233),imshow(ig),title('绿色分量');
subplot(234),imshow(ib),title('蓝色分量');
subplot(235),imshow(i2),title('局部极大值叠加图像');
subplot(236),imshow(fgm),title('检测局部最大值,没[]');%结果表示和有[]的没啥区别

se2=strel(ones(15,15));%构建结构元素，15×15的1矩阵
fgm2=imclose(fgm,se2);%进行闭运算，即先膨胀再腐蚀
fgm3=imerode(fgm2,se2);%载进行腐蚀
fgm4=bwareaopen(fgm3,400);%删除小于400像素的连通块

figure(5);
subplot(221),imshow(fgm),title('检测局部最大值');
subplot(222),imshow(fgm2),title('对局部最大值进行闭运算');
subplot(223),imshow(fgm3),title('闭运算后进行腐蚀');
subplot(224),imshow(fgm4),title('删除小于400像素的连通块');

bw=im2bw(Iobrcbrc,graythresh(Iobrcbrc));%对Iobrcbrc进行二值化，二值化的阈值自动获取
D=bwdist(bw);
DL=watershed(D);
bgm=DL==0;
gradmag2=imimposemin(gradmag,bgm|fgm4);%bgm和fgm4进行或运算，使得梯度幅值图像gradmag上只在或运算结果为非零的位置上有局部最小值

L=watershed(gradmag2);
%第一种显示方法
Lrgb=label2rgb(L,'jet','w','shuffle');
figure(6);
subplot(121),imshow(Lrgb),title('分水岭显示结果1');
%第二种显示方法
subplot(122),imshow(i,[]),title('分水岭显示结果2');
hold on
himage=imshow(Lrgb);
set(himage,'AlphaData',0.3);