分水岭算法 matlab实现

ゝ一世哀愁。 2022-09-11 07:19 132阅读 0赞

## 背景 ##

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

## 大概思想 ##

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

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

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

## 流程框图 ##

![在这里插入图片描述][0afe26705f3442e2be03dc5a277b2a17.png]

![在这里插入图片描述][ae926199bc64472ebf702389b8adeb4b.png]

## 代码 ##

`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);

[0afe26705f3442e2be03dc5a277b2a17.png]: /images/20220828/c92f0d193d5b4b8caa490f06800bd5d8.png
[ae926199bc64472ebf702389b8adeb4b.png]: /images/20220828/3da261346a964571859c75d8ec8117a6.png