背景
做图像分割的时候用到了,就学习了一下
大概思想
把图像中的像素大小理解成山地的海拔,向山地灌水,海拔低的地方会积水,这些地方称之为谷底。随着水位上升,不同谷底的水会相遇,相遇的地方就是分水岭。
总体上是按照给每个点贴标签的方法实现的对每个在谷底的点贴上从1开始的不同的标签,同一个谷底标签相同。给分水岭贴上0的标签。未贴标签的点标为-1,其实就是初始化的时候全部-1。
- 1、将图像中的每个像素按像素值大小排好顺序。
- 2、从较小的像素值开始进行处理。
- 3、搜索统计当前点周围像素(8邻域,边界的话就少几个像素点)中标签的种类,
- a、如果只有一种标签,为该点贴上相同的标签;
- b、 如果有2种及以上的标签,认为该点为分水岭。
- c、如果没有标签(-1没处理过)或只存在0(分水岭),搜索当前点的连通区域,,这儿可能会出现两种情况
- (1)、当前点的连通域和之前处理过的接通了,如果是这样的话呢,把当前点的连通区域标记为和他连通的区域中较小的编号
- (2)、当前点的连通域没有连通,那么当前点就把认为该点为新的谷底,谷底数加一,并把当前点的连通域都标记为新的谷底编号
流程框图
代码
data
是一个图片的灰度值矩阵
clear all
clcload data
picture = data;
[row,col] = size(picture); % 获得图片尺寸
watershed_result = -ones(row,col); % 将结果矩阵全部赋值为-1表示所有点都未处理过valley_number = 0; % 谷底数为0
tic
[picture_value, picture_index] = sort(picture(:)); % 将图片中所有元素按像素值大小排序
toc
total_pixel_number = row*col; % 总元素个数tic
for now_index = 1:total_pixel_number % 对每个像素都要处理if watershed_result(picture_index(now_index)) ~= -1 % 如果标记为处理过跳过该像素continue;endnow_picture_index = picture_index(now_index); % 正在处理的像素的位置now_picture_value = picture(now_picture_index); % 正在处理的像素的像素值vector_now_pixel_neighbor = neighbor_find_C(row,col,now_picture_index); % 获取当前像素点的周围8个像素点的位置temp_vector = sort(watershed_result(vector_now_pixel_neighbor));%获取周围像素点的标签temp_vector = unique(temp_vector);temp_vector = temp_vector(temp_vector>0); %除了-1,0的标签种类temp_vector_length = length(temp_vector); if temp_vector_length == 0 %种类为0is_same_area_index = zeros(2*col+2*row,1,'double');%和处理点像素值相同的连通区域is_same_area_index(1) = now_picture_index;area_num = 1;%连通区域的像素个数while(1) %获得和处理点像素值相同的连通区域%获取联通区域的周围像素点的位置 is_same_area_neighbor_index = area_neighbor_find_C(row,col,is_same_area_index(1:area_num));%获取联通区域的周围像素点像素值大小和当前点像素值大小相同的点的像素值temp_vector = is_same_area_neighbor_index(picture(is_same_area_neighbor_index) == now_picture_value);temp_vector_length = length(temp_vector);if temp_vector_length == 0break;endis_same_area_index(area_num+1:area_num+temp_vector_length) = temp_vector;area_num = area_num + temp_vector_length;is_same_area_index = is_same_area_index(is_same_area_index>0);endtemp_vector = sort(watershed_result(is_same_area_neighbor_index));%获取连通区域周围像素点的标签temp_vector = unique(temp_vector);temp_vector = temp_vector(temp_vector>0);temp_vector_length = length(temp_vector);%除了-1,0的标签种类if temp_vector_length == 0%种类为0valley_number = valley_number+1;%谷底数加1watershed_result(is_same_area_index(is_same_area_index>0)) = valley_number;%认为该连通区域为一个新的谷底else watershed_result(now_picture_index) = temp_vector(1);%否则该点贴与之相邻的标签中除-1,0外较小的标签endelse if temp_vector_length == 1 % 相邻8个像素中标签种类除-1,0还有1个watershed_result(now_picture_index) = temp_vector(1); % 该像素点贴相同标签else % 相邻8个像素中标签种类除-1,0为1个以上watershed_result(now_picture_index) = 0; % 认为该点为分水岭endend
end
toc
figure,imshow(watershed_result);%输出分水岭
用到的两个函数
neighbor_find_C.m
function vector_neighbor_C = neighbor_find_C(row,col,nn)
total_number=row*col;
now_col= floor(nn/row)+1;
temp=rem(nn,row);
if temp==0now_row=row;
elsenow_row=temp;
end
if nn==1vector_neighbor_C=[2;row+1;row+2];
else if nn==rowvector_neighbor_C=[row;row+row;row+row-1];else if nn==total_number-row+1vector_neighbor_C=[total_number-row+2;total_number-row-row;total_number-row-row+1];else if nn==total_numbervector_neighbor_C=[total_number-1;total_number-row;total_number-row-1];else if now_row==1vector_neighbor_C=[nn+row;nn-row;nn+1;nn+row+1;nn-row+1];else if now_row==rowvector_neighbor_C=[nn+row;nn-row;nn-1;nn+row-1;nn-row-1];else if now_col==1vector_neighbor_C=[nn+1;nn-1;nn+row;nn+row+1;nn+row-1];else if now_col==colvector_neighbor_C=[nn+1;nn-1;nn-row;nn-row+1;nn-row-1];elsevector_neighbor_C=[nn+1;nn-1;nn-row;nn+row;nn-row+1;nn-row-1;nn+row+1;nn+row-1];endendendendendendend
end
area_neighbor_find_C.m
function vector_neighbor = area_neighbor_find_C(row,col,area_index)
I = zeros(row,col);
se=ones(3);
I(area_index) = 1;
I2 = imfilter(I,se); % 滤波
I2=I2-I.*9;
vector_neighbor = find(I2>0);