一、簡介

分水嶺演算法是一種影象區域分割法，分割的過程中將圖片轉化為灰度圖，然後我會將灰度值看作是海拔，然後向較低點注水，這種基於地形學的解釋，我們著重考慮三種點：

極小值點，該點對應一個盆地的最低點，當我們在盆地裡滴一滴水的時候，由於重力作用，水最終會匯聚到該點。注意：可能存在一個最小值面，該平面內的都是極小值點。
盆地的其它位置點，該位置滴的水滴會匯聚到區域性最小點。
盆地的邊緣點，是該盆地和其它盆地交接點，在該點滴一滴水，會等概率的流向任何一個盆地。

明白上述三種點之後，我們開始往盆地的極小值點注水，然後隨著注水的深入，每一個極小值點慢慢的向外擴充套件，然後知道兩個盆地的水匯合，匯合處就是我們需要的分水嶺。

從下圖可以直觀理解一下，首先這三塊區域都含有極小值點

然後逐漸填充就能獲得分水嶺（即分界線）

得到分界線就能完成影象分割：

二、原始碼

function susanseg
clear all; close all; clc
image= imread('cell.jpg');
% 用SUSAN演算法進行邊緣檢測
image = susan(image,4);
figure, imshow(image,[]);
%imwrite(image, './susanout/susanout.jpg');
% 將image轉為二值影象儲存後，用影象處理工具
% 把其背景的所有連通區域處理為黑色，即只有細
% 胞體是白色，便於細胞數目的搜尋
BW = im2bw(image, graythresh(image));
bounder_area = length(find(BW==0));
%imwrite(BW, './susanout/bw.jpg');
figure, imshow(BW);


% 申明全域性變數
global B Dir m n;
B = imread('./blackbackground.jpg');
B = im2bw(B, graythresh(B));
[m,n] = size(B);
figure, imshow(B);

% 細胞的總面積，即細胞所佔的畫素數目，包括細胞的邊界
% 由於SUSAN提取出的邊界已被增寬，所以將邊界畫素數除以2
% 來作為細胞的邊界畫素數目
total_area = length(find(B==1)) + bounder_area/2;
NUM = 5; % 細胞面積閾值
count = 0; % 細胞總數
% 搜尋方向向量，4鄰域搜尋
Dir = [-1 0; 0 1; 1 0; 0 -1;];
% 搜尋方向向量，8鄰域搜尋
%Dir = [-1 0; -1 1; 0 1; 1 1; 1 0; 1 -1; 0 -1; -1 -1;];
for i = 1:m
    for j = 1:n
        if B(i,j)==1 % 是細胞畫素
            num = search(i,j,4) + 1; % 計算該細胞的畫素數目
            if num>NUM
                count = count  + 1;
            else
                total_area = total_area - num; % 減掉不是細胞的面積
            end
        end
    end
end
%fid = fopen('./susanout/results.txt', 'wt');
fprintf('影象尺寸: %d * %d, SUSAN閾值: 4, 細胞面積閾值: %d\n', ...
    n, m, NUM);
fprintf('細胞總數: %d, 細胞總面積: %.2f, 平均細胞面積: %.2f\n', ...
    count, total_area, total_area/count);
%fprintf(fid,'影象尺寸: %d * %d, SUSAN閾值: 4, 細胞面積閾值: %d\n', ...
%    n, m, NUM);
%fprintf(fid,'細胞總數: %d, 細胞總面積: %.2f, 平均細胞面積: %.2f\n', ...
%    count, total_area, total_area/count);
%fclose(fid);
end
% -----------------------------------------------------------------------
% 
% This function uses the SUSAN algorithm to find edges within an image
% 
%
% >>image_out = susan(image_in,threshold)
%
%
% Input parameters ... The gray scale image, and the threshold 
% image_out .. (class: double) image indicating found edges
% typical threshold values may be from 10 to 30
%
%
%The following steps are performed at each image pixel: 
% ( from the SUSAN webpage, http://www.fmrib.ox.ac.uk/~steve/susan/susan/node4.html )
% 
% Place a circular mask around the pixel in question. 
% Calculate the number of pixels within the circular mask which have similar brightness to 
% the nucleus. These define the USAN. 
% Subtract USAN size from geometric threshold to produce edge strength image. 
%
% Estimating moments to find the edge direction has not been implemented . 
% Non-maximal suppresion to remove weak edges has not been implemented yet.
%
% example:
%
% >> image_in=imread('test_pattern.tif');
% >> image = susan(image_in,27);
% >> imshow(image,[]) 
%
%
% Abhishek Ivaturi
% 
% -------------------------------------------------------------------------
 

function image_out = susan(im,threshold)

% check to see if the image is a color image...
%im= imread('test_pattern.tif')
%threshold=27;
d = length(size(im));
if d==3
    image=double(rgb2gray(im));
elseif d==2
    image=double(im);
end

% mask for selecting the pixels within the circular region (37 pixels, as
% used in the SUSAN algorithm

mask = ([ 0 0 1 1 1 0 0 ;0 1 1 1 1 1 0;1 1 1 1 1 1 1;1 1 1 1 1 1 1;1 1 1 1 1 1 1;0 1 1 1 1 1 0;0 0 1 1 1 0 0]);  


% the output image indicating found edges
R=zeros(size(image));


% define the USAN area
nmax = 3*37/4;

% padding the image
[a b]=size(image);
new=zeros(a+7,b+7);
[c d]=size(new);
new(4:c-4,4:d-4)=image;
  
for i=4:c-4
    
    for j=4:d-4
        
        current_image = new(i-3:i+3,j-3:j+3);
        current_masked_image = mask.*current_image;
   

%   Uncomment here to implement binary thresholding

%         current_masked_image(find(abs(current_masked_image-current_masked_image(4,4))>threshold))=0;         
%         current_masked_image(find(abs(current_masked_image-current_masked_image(4,4))<=threshold))=1;


%   This thresholding is more stable
                 
                   
        current_thresholded = susan_threshold(current_masked_image,threshold);
        g=sum(current_thresholded(:));
        
        if nmax<g
            R(i,j) = g-nmax;
        else
            R(i,j) = 0;
        end
    end
end
 

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423