一、简介
分水岭算法是一种图像区域分割法，分割的过程中将图片转化为灰度图，然后我会将灰度值看作是海拔，然后向较低点注水，这种基于地形学的解释，我们着重考虑三种点：
极小值点，该点对应一个盆地的最低点，当我们在盆地里滴一滴水的时候，由于重力作用，水最终会汇聚到该点。注意：可能存在一个最小值面，该平面内的都是极小值点。
盆地的其它位置点，该位置滴的水滴会汇聚到局部最小点。
盆地的边缘点，是该盆地和其它盆地交接点，在该点滴一滴水，会等概率的流向任何一个盆地。
明白上述三种点之后，我们开始往盆地的极小值点注水，然后随着注水的深入，每一个极小值点慢慢的向外扩展，然后知道两个盆地的水汇合，汇合处就是我们需要的分水岭。
从下图可以直观理解一下，首先这三块区域都含有极小值点
然后逐渐填充就能获得分水岭（即分界线）
得到分界线就能完成图像分割：
二、源代码
-
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
三、运行结果
四、备注
完整代码添加QQ1575304183
转载：https://blog.csdn.net/weixin_50197058/article/details/116429315
查看评论