分水岭算法 matlab实现

大家好，又见面了，我是你们的朋友全栈君。

背景

做图像分割的时候用到了，就学习了一下

大概思想

把图像中的像素大小理解成山地的海拔，向山地灌水，海拔低的地方会积水，这些地方称之为谷底。随着水位上升，不同谷底的水会相遇，相遇的地方就是分水岭。

总体上是按照给每个点贴标签的方法实现的对每个在谷底的点贴上从1开始的不同的标签，同一个谷底标签相同。给分水岭贴上0的标签。未贴标签的点标为-1，其实就是初始化的时候全部-1。

1、将图像中的每个像素按像素值大小排好顺序。
2、从较小的像素值开始进行处理。
3、搜索统计当前点周围像素（8邻域，边界的话就少几个像素点）中标签的种类，
a、如果只有一种标签，为该点贴上相同的标签；
b、如果有2种及以上的标签，认为该点为分水岭。
c、如果没有标签（-1没处理过）或只存在0（分水岭），搜索当前点的连通区域，，这儿可能会出现两种情况
（1）、当前点的连通域和之前处理过的接通了，如果是这样的话呢，把当前点的连通区域标记为和他连通的区域中较小的编号
（2）、当前点的连通域没有连通，那么当前点就把认为该点为新的谷底，谷底数加一，并把当前点的连通域都标记为新的谷底编号

流程框图

在这里插入图片描述

代码

data 是一个图片的灰度值矩阵

clear all 
clc

load 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;
    end
    
    now_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  %种类为0
        is_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 == 0
                break;
            end
            is_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);
        end

        temp_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%种类为0
            valley_number = valley_number+1;%谷底数加1
            watershed_result(is_same_area_index(is_same_area_index>0)) = valley_number;%认为该连通区域为一个新的谷底
        else 
            watershed_result(now_picture_index) = temp_vector(1);%否则该点贴与之相邻的标签中除-1，0外较小的标签
        end
    else 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;                            % 认为该点为分水岭
        end
    end
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==0
    now_row=row;
else
    now_row=temp;
end
if nn==1
    vector_neighbor_C=[2;row+1;row+2];
else if nn==row
    vector_neighbor_C=[row;row+row;row+row-1];
    else if nn==total_number-row+1
            vector_neighbor_C=[total_number-row+2;total_number-row-row;total_number-row-row+1];
        else if nn==total_number
                vector_neighbor_C=[total_number-1;total_number-row;total_number-row-1];
            else if now_row==1
                    vector_neighbor_C=[nn+row;nn-row;nn+1;nn+row+1;nn-row+1];
                else if now_row==row
                        vector_neighbor_C=[nn+row;nn-row;nn-1;nn+row-1;nn-row-1];
                    else if now_col==1
                            vector_neighbor_C=[nn+1;nn-1;nn+row;nn+row+1;nn+row-1];
                        else if now_col==col
                                vector_neighbor_C=[nn+1;nn-1;nn-row;nn-row+1;nn-row-1];
                            else
                                vector_neighbor_C=[nn+1;nn-1;nn-row;nn+row;nn-row+1;nn-row-1;nn+row+1;nn+row-1];
                            end
                        end
                    end
                end
            end
        end
    end
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);