概念：将嵌在一幅图像的恒定区域或亮度几乎不变的区域里的孤立点的检测，就是点检测。可以用点检测的模板来将孤立的点检测出来：这个模板的作用就是当模板中心是孤立点时，模板的相应最强，而在非模板中心时，相应为零。
− 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 \begin{array}{|c|c|c|} \hline-1 & -1 & -1 \\ \hline-1 & 8 & -1 \\ \hline-1 & -1 & -1 \\ \hline \end{array} −1−1−1−18−1−1−1−1
在MATLAB中进行点检测操作，基本语法为：
g = abs(imfilter(tofloat(f),w))>=T; %g是包含检测点的图像；f是输入图像；w是点检测模板运用上面的语法编写实验代码：
f = imread('point.png'); subplot(121),imshow(f); w = [-1 -1 -1;-1 8 -1;-1 -1 -1]; g = abs(imfilter(tofloat(f),w)); T = max(g(:)); %在滤波后图像中选择的最大值 g = g >=T; subplot(122),imshow(g);椭圆右上角的黑色点在滤波之后可以很明显的看出来
点检测的另一种方法是在大小为m×n的所有邻点中寻找一些点，最大值和最小值的差超出了T的值，可以用函数ordfilt2进行操作：
g = ordfilt2(f,m*n,ones(m,n))-ordfilt2(f,1,ones(m,n)); g = g>=T;通过比较典型模板的计算值，确定一个点是否在某个方向的线上。 也可以设计其他模板：模板系数之和为0，感兴趣的方向系数值较大。
举例
水平模板： − 1 − 1 − 1 2 2 2 − 1 − 1 − 1 \begin{array}{|c|c|c|} \hline-1 & -1 & -1 \\ \hline 2 & 2 & 2 \\ \hline-1 & -1 & -1 \\ \hline \end{array} −12−1−12−1−12−1 135°模板： 2 − 1 − 1 − 1 2 − 1 − 1 − 1 2 \begin{array}{|c|c|c|} \hline 2 & -1 & -1 \\ \hline-1 & 2 & -1 \\ \hline-1 & -1 & 2 \\ \hline \end{array} 2−1−1−12−1−1−12 用R1,R2,R3和R4分别代表水平、45°、垂直和 -45°方向线的模板响应，在图像中心的点，如果 |Ri|>|Rj| \text { |Ri|>|Rj| } |Ri|>|Rj| j不等于i，则此点被认为与在模板i方向上的线更相关。
检测指定方向的线 :
f = rgb2gray(imread('san.jpg')); w = [2 -1 -1;-1 2 -1;-1 -1 2]; % +45°方向检测线 g = imfilter(double(f),w); gtop = g(1:120,1:120); % 左上角区域 gtop = pixeldup(gtop,4); % 通过复制像素将图像扩大gtop*4倍 gbot = g(end-119:end,end-119:end); % 右下角区域 gbot = pixeldup(gbot,4); g1 = abs(g); % 检测图的绝对值 T = max(g1(:)); g2 = g1>=T; subplot(3,2,1);imshow(f);title('(a)连线模板图像'); subplot(3,2,2);imshow(g,[]);title('(b)+45°线处理后的结果'); subplot(3,2,3);imshow(gtop,[]);title('(c)(b)中左上角的放大效果'); subplot(3,2,4);imshow(gbot,[]);title('(d)(b)中右下角的放大效果'); subplot(3,2,5);imshow(g1,[]);title('(e)(b)的绝对值'); subplot(3,2,6);imshow(g2);title('(f)满足g>=T的所有点');实验分析： 在进行45°的边缘检测中，只能够检测到45°或者135°的边缘。绝对值的图像是将图像变成了二值图像（只有黑与白），比灰度图像更加清晰看出检测到的线条边缘。图二的45度线不是很明显，正说明了该图模板45方向上的线不怎么相关。
以上两种方法在图像分割中很重要，但是所常用的方法是边缘检测，这种方法是检测亮度的不连续性。这样的不连续是用一阶和二阶导数来检测的。
首先，我们先介绍二阶函数的梯度向量公式： ∇ f = [ g x g y ] = [ ∂ f ∂ x ∂ f ∂ y ] \nabla f=\left[\begin{array}{l} g_{x} \\ g_{y} \end{array}\right]=\left[\begin{array}{l} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{array}\right] ∇f=[gxgy]=[∂x∂f∂y∂f]
这个向量的幅值可以简化为如下形式：通常使用梯度的幅值或近似值来作为“梯度”。
∇ f = mag ( ∇ f ) = [ g x 2 + g y 2 ] 1 / 2 ≈ g x 2 + g y 2 ≈ ∣ g x ∣ + ∣ g y ∣ \nabla f=\operatorname{mag}(\nabla f)=\left[g_{x}^{2}+g_{y}^{2}\right]^{1 / 2} \approx g_{x}^{2}+g_{y}^{2} \approx\left|g_{x}\right|+\left|g_{y}\right| ∇f=mag(∇f)=[gx2+gy2]1/2≈gx2+gy2≈∣gx∣+∣gy∣
梯度向量的基本性质是：梯度向量指向(x,y)坐标处f的最大变化率方向。 最大变化率处发生的角度是：
α ( x , y ) = tan − 1 ( g y g x ) \alpha(x, y)=\tan ^{-1}\left(\frac{g_{y}}{g_{x}}\right) α(x,y)=tan−1(gxgy)
在图像处理工具箱中，可以使用函数edge来作为边缘估计器。
[g,t] = edge(f,'method',parameters) sobel边缘检测算子调用语法是：
[g,t] = edge(f,'sobel',T,dir) %T是阈值，dir是边缘检测首选方向：'horizontal'、'vertical'、‘both’ prewitt边缘检测算子 调用语法为： [g,t] = edge(f,'prewitt',T,dir) %T是阈值，dir是边缘检测首选方向：'horizontal'、'vertical'、‘both’3.Roberts边缘检测算子 调用语法为：
[g,t] = edge(f,'roberts',T,dir) %T是阈值，dir是边缘检测首选方向：'horizontal'、'vertical'、‘both’ LoG检测算子由高斯函数： G ( x , y ) = e − x 2 + y 2 2 σ 2 G(x, y)=e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}} G(x,y)=e−2σ2x2+y2
其中的σ 是标准差。这是平滑函数，如果和图像卷积，将会使图像变模糊。模糊的程度由σ 的值决定。这个函数的拉普拉斯算法是：
∇ 2 G ( x , y ) = ∂ 2 G ( x , x ) ∂ x 2 + ∂ 2 G ( x , x ) ∂ y 2 = [ x 2 + y 2 − 2 σ 2 σ 4 ] e − x 2 + y 2 2 σ 2 \nabla^{2} G(x, y)=\frac{\partial^{2} G(x, x)}{\partial x^{2}}+\frac{\partial^{2} G(x, x)}{\partial y^{2}}=\left[\frac{x^{2}+y^{2}-2 \sigma^{2}}{\sigma^{4}}\right] e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}} ∇2G(x,y)=∂x2∂2G(x,x)+∂y2∂2G(x,x)=[σ4x2+y2−2σ2]e−2σ2x2+y2
这个函数称为LoG算子。原理是：由于求二阶导数是线性操作，所以用这个函数卷积这幅图像与先用平滑函数对图像卷积，再对结果进行拉普拉斯计算的结果是一样的。
这样会得到两种效果：一种是平滑图像，减少了噪声；另一种是计算拉普拉斯，从而产生双边缘图像，然后在双边缘之间定位由发现的零交叉组成的边缘。
调用语法为：
[g,t] = edge(f,'log',T,sigma) %T是阈值，sigma是标准差，默认为2 零交叉检测算子 这种算子与LoG方法类似，不同之处在于卷积使用特殊的滤波函数H来完成的。调用语法为：
[g,t] = edge(f,'zerocross',T,H) %T是阈值 Canny检测算子 这种算子是edge函数中最强的边缘检测算子： （1）图像用指定了标准差σ的高斯滤波器来平滑，用来减少噪声； （2）局部梯度和边缘方向在每个点都进行计算。边缘点被定义为梯度方向上局部强度最大的点。 （3）在（2）中决定的边缘点在梯度幅度图像上给出脊。然后算法追踪所有脊的顶部，设置所有的不在脊的顶部的像素为零。因此在输出中给出一条细线，这是非最大值抑制处理。脊像素使用称为“滞后阈值”的技术进行阈值处理，这种技术以使用两个阈值为基础，即T1和T2，且T1<T2。值大于T2的脊像素称为强边缘像素，T1和T2之间的脊像素称为弱边缘像素。 （4）算法用合并8连接的弱像素点到强像素点的方法执行边缘连接。调用语法为：
[g,t] = edge(f,'canny',T,sigma) %T是向量=[T1,T2]；sigma是标准差，默认为1实验代码：
f = imread('house.jpg'); subplot(231),imshow(f),title('原始图像'); [gv,t] = edge(f,'sobel','vertical'); subplot(232),imshow(gv),title('默认阈值垂直的sobel图像'); t gv = edge(f,'sobel',0.15,'vertical'); subplot(233),imshow(gv),title('阈值0.15垂直的sobel图像') gboth = edge(f,'sobel',0.15); subplot(234),imshow(gboth),title('阈值为0.15垂直水平的sobel图像') wneg45 = [-2 -1 0;-1 0 1;0 1 2]; gneg45 = imfilter(tofloat(f),wneg45,'replicate'); T =0.3*max(abs(gneg45(:))); gneg45 = gneg45>=T; subplot(235),imshow(gneg45),title('阈值45°边缘图像') t = 0.0583实验分析： 当参数是默认垂直的sobel边缘检测时，得到边缘细节是垂直的，水平的边缘线是没有能够检测出来的。当加入参数0.15之后，之前能够检测到的细节却不能够检测出来了，是因为这已经超出了阈值的范围；当采用垂直水平的边缘检测，将之前没有能够检测出的水平线能够检测并显示出来了；当所检测的边缘线是45°的时候，该算法只能够检测到45°或135°的边缘线条，并显示出来。
将sobel、Log和canny边缘检测算子进行比较，编写实验代码：
% f = rgb2gray(imread('flower.jpg')) f = tofloat(f); subplot(221),imshow(f),title('原始图像') [gSobel_dafault,ts] = edge(f,'sobel'); [gLog_dafault,tlog] = edge(f,'log'); [gCanny_dafault,tc] = edge(f,'canny'); subplot(222),imshow(gSobel_dafault),title('默认的sobel图像'); subplot(223),imshow(gLog_dafault),title('默认的log图像'); subplot(224),imshow(gCanny_dafault),title('默认的canny图像'); gSobel_best = edge(f,'sobel',0.05); gLog_best = edge(f,'log',0.003,2.25); gCanny_best = edge(f,'canny',[0.04 0.1],1.5); subplot(222),imshow(gSobel_best),title('最好效果的sobel图像'); subplot(223),imshow(gLog_best),title('最好效果的log图像'); subplot(224),imshow(gCanny_best),title('最好效果的canny图像');分析：
感兴趣的主要边缘是建筑的角落、窗户、由亮的砖 结构形成的门口和门本身、屋顶线。由默认值开始，每个选项的参数分别随早些时候提到的显示主要特性的目标而变化，并尽可能减少不必要的细节。Sobel 得出的结果与我们想要的水泥带和门口的左边缘相差还较远。在 (d)中，LoG的结果相比 Sobel 的结果要好一些，比 LoG的默认值得出的结果要好得多。但是，门口的左边缘没有检测出来，建筑物周围水泥带的两个边也还是不够清晰。Canny得出的结果(图(f))到目前为止要远远好于前两种。同样的一幅图像，使用默认的sobel边缘检测得到的图像中线条是比较其他两种是比较少的，而且获取的边缘细节不是很完整，缺失了很多重要的细节。在理想情况下，只能产生边缘上的像素，为了更好的得到边缘特性，因此下面介绍霍夫变换的线检测方法。
图像处理工具箱提供了 3 个与霍夫变换有关的函数。
函数 hough 实现霍夫变换的概念。
函数 houghpeaks 寻找霍夫变换中的峰值（高计数累加器单元）。
函数 houghlines 基于前两个函数的结果，提取原始图像中的线段。
函数hough 语法：
[H, theta, rho] = hough(f) [H, theta, rho] = hough(f, ‘ThetaRes’, val1, ‘RhoRes’, val2)H 是霍夫变换矩阵，theta(以度计)和 rho 是 ρ 和 θ 值向量，在这些值上产生霍夫变换。val1是 0 到 90 的标量，指定了沿 θ 轴霍夫变换的间距(默认是 1)，val2 是 0<val2<hypot(size(I,1),size(I,2))的实标量，指定了沿 ρ 轴的霍夫变换的间隔(默认是 1)。
示例：
f = zeros(101,101); f(1,1) = 1;f(101,1) = 1;f(1,101) = 1; f(101,101) = 1;f(51,51) = 1; H = hough(f) imshow(H,[]) [H,theta,rho] = hough(f); imshow(H,[],'XData',theta,'YData',rho,'InitialMagnification','fit') axis on,axis normal xlabel('\theta'),ylabel('\rho')实验分析： 三条曲线（直线也可以考虑为曲线）在±45°处的交点指出：f中有两组三个共线的点。两条曲线在(ρ,θ)=(0,-90)、(-100,-90)、(0,0)和(100,0)处的交点指出：有4组位于垂直线和水平线上的共线点。（需要将图像放大以方便观察其坐标）
函数houghpeaks：寻找霍夫变换的峰值（累加单元的高计数）调用语法：
peaks = houghpeaks(H,NumPeaks) 或 peaks = houghpeaks(H,‘Threshold’,val1,'NHoodSize',val2)其中，val1可以从0到Inf变换，默认是0.5*max(H(：))，val2是指定量围绕峰值的邻域大小。
这个过程的基本思想是：通过把峰值的直接邻域中的霍夫变换单元置0来清理峰值。
函数houghlines：当峰值在霍夫变换中被识别出来，则可用此函数来决定线的起始点和终点。调用函数：
lines = houghlines(f,theta,rho,peaks) lines = houghlines(f,'FillGap',val1,'MinLength',val2)其中，theta和rho是来自函数hough的输出，peaks是来自houghpeaks的输出。val1指定了与相同的霍夫变换相关的两条线段的距离。当两条线段之间的鹅咀里小于指定的值时，则把线段合并为一条线段（默认20）；val2指定合并的线是保留还是去掉，如果合并的线比val2指定的值短，则丢掉（默认是40）。
示例：
f = imread('house.jpg'); subplot(131),imshow(f),title('原始图像'); [H,theta,rho] = hough(f,'ThetaRes',0.2) subplot(132),imshow(H,[ ],'XData',theta,'YData',rho,'InitialMagnification','fit') axis on,axis normal xlabel('\theta'),ylabel('\rho') peaks =houghpeaks(H,5); hold on plot(theta(peaks(:,2)),rho(peaks(:,1)),'linestyle','none','marker','s','color','w') lines = houghlines(f,theta,rho,peaks); subplot(133),imshow(f); for k = 1:length(lines) xy = [lines(k).point1;lines(k).point2]; plot(xy(:,1),xy(:,2),'LineWidth',4,'Color',[0.8 0.8 0.8]); end一般情况下，一张图片分为前景和背景，我们感兴趣的一般的是前景部分，所以我们一般使用阈值将前景和背景分割开来，使我们感兴趣的图像的像素值为1，不感兴趣的我0，有时一张图我们会有几个不同的感兴趣区域（不在同一个灰度区域），这时我们可以用多个阈值进行分割，这就是阈值处理。
单个阈值： 两个阈值： 示意图如下：
一般选取阈值就是图像直方图的视觉检测。将区分度大的两个灰度级部分之间进行划分，取T为阈值来分开它们。在此基础上学习一种自动地选择阈值的算法，方法如下：
针对全局阈值选择初始估计值T；用T分割图像，G1是所有灰度值大于T的像素组成，G2是所有灰度值小于等于T的像素组成；分别计算G1和G2区域内的平均灰度值m1和m2；计算出新的阈值，取他们的m1和m2平均值数；重复步骤2.-4.，直到在连续的重复中，T的差异比预先设定的参数小为止；使用函数im2bw分割图像：g = im2bw(f,T/den)示例：
f = imread('house.jpg'); f = imnoise(f,'gaussian'); count = 0; T = mean2(f); done = false; while ~done count = count+1; g = f>T; Tnext = 0.5*(mean(f(g))+mean(f(~g))); done = abs(T-Tnext)<0.5; T =Tnext; end count T g = im2bw(f,T/255); subplot(131),imshow(f),title('带噪声的图像'); subplot(132),imhist(f),title('图像的直方图'); subplot(133),imshow(g),title('全局阈值分割的结果')Otsu’s方法的最佳方法是选择阈值k,最大类间方差σ2(k)来定义的： 当设置的方差越大，则完全分割一幅图像的阈值就会越接近。公式中的k就是我们所要寻找的最佳阈值，当k不唯一时，则将所有的最佳阈值进行取平均值即可。
调用的语法是：
[T,SM] = graythresh(f) %SM是可分性度量,T是归一化的阈值示例：
f = imread('cell.jpg'); f = imnoise(f,'gaussian'); [T,SM] = graythresh(f) T*255（基本全局算法）一般来说，高SM值说明灰度分成两类的可能性比较高。
T = 0.4784 SM = 0.6949 ans = 122 f2 = imread('cell.jpg'); subplot(221),imshow(f2),title('原始图像'); subplot(222),imhist(f2),title('原始图像直方图'); count = 0; T = mean2(f2); done = false; while ~done count = count+1; g = f2>T; Tnext = 0.5*(mean(f2(g))+mean(f2(~g))); done = abs(T-Tnext)<0.5; T =Tnext; end g = im2bw(f2,T/255); subplot(223),imshow(g,[]),title('用基本全局算法分割结果') [T,SM] = graythresh(f2); SM T*255 g = im2bw(f2,T); subplot(224),imshow(g),title('用Otsu’s算法分割结果') SM = 0.7017 ans = 123实验分析： 可以看出，用基本全局阈值处理时细胞的分割效果并不是很好，没有能够分割清楚，是由于前景和背景的灰度级比较相近而不能够完全分离开，因此当使用了otsu’s算法进行分割时，可以完全将细胞（前景）和背景分割开来。尽管分离度的度量值比较低，但是还是可以准确地从背景中提取细胞。
噪声可以影响阈值的选择，当噪声不能够在源头减少，在阈值处理之前可以将图像进行平滑处理，这样可以更好地进行阈值处理。
编写实验代码：
f = imread('house.jpg'); subplot(241),imshow(f),title('原始图像'); fn = imnoise(f,'gaussian',0,0.038); subplot(242),imshow(fn),title('具有高斯噪声的图像'); subplot(243),imhist(fn),title('具有高斯噪声的直方图'); Tn = graythresh(fn); gn = im2bw(fn,Tn); subplot(244),imshow(gn),title('用otsu’s分割的图像'); w = fspecial('average',5); fa = imfilter(fn,w,'replicate'); subplot(245),imshow(fa),title('平滑后的图像'); subplot(246),imhist(fa),title('平滑后的直方图'); Ta = graythresh(fa); ga =im2bw(fa,Ta); subplot(247),imshow(ga),title('平滑后并用otsu’s分割的图像');边缘改进的阈值处理：主要是处理那些位于或接近物体和背景间边缘的像素，使得这些像素分离开的操作。
具体算法过程如下：
用一种边缘查找方式计算图像的模板的值。通过百分比指定阈值。由于计算的边缘模板值中有很多噪声，所以可以将计算值排序，并选择百分比相对高的值(大于百分下的值的阈值)作为阈值。根据指定的阈值，对第一步的图像边缘的值进行选择。使高于阈值的像素点值为1，低于的值为零，由此选择出部分边缘点的二值图像(模板)。用刚才计算出来的模板与原图像相乘，获得一幅新的图像。对新的图像使用otsu进行分割。用自定义的函数percentile2i来计算灰度值I的对应指定的百分比：
I = percentile2i(h,p)示例：
f = tofloat(imread('star.jpg')); subplot(231),imshow(f),title('原图像'); subplot(232),imhist(f),title('原图像直方图'); sx = fspecial('sobel'); sy = sx'; gx = imfilter(f,sx,'replicate'); gy = imfilter(f,sy,'replicate'); grad = sqrt(gx.*gx+gy.*gy); grad = grad/max(grad(:)); h =imhist(grad); Q = percentile2i(h,0.999); markerImage = grad>Q; subplot(233),imshow(markerImage),title('以99.9%进行阈值处理后的梯度幅值图像'); fp = f.*markerImage; subplot(234),imshow(fp),title('原图像与梯度幅值乘积的图像'); hp = imhist(fp); hp(1) = 0; subplot(235),bar(hp),title('原图像与梯度幅值乘积的直方图'); T = otsuthresh(hp); T*(numel(hp)-1) g = im2bw(f,T); subplot(236),imshow(g),title('改进边缘后的图像')实验分析： 拉普拉斯算法是用于图像的锐化，把图像的边缘提取出来。在这个部分中，用拉普拉斯边缘信息改进全局阈值处理也是一样的，首先将原始图像的边缘提取出来，然后将原始图像与已标记好的图像进行相乘，得到的就是实验中所需要的图像，最终图像就如边缘改进结果那样，修改边缘信息提取出来，再用阈值进行处理。
函数percentile2i :
function I=percentile2i(h,P) %PERCENTILE2I Computes an intensity value given a percentile. % I=PERCENTILE2I(H,P) Given a percentile ,p,and a histogram, % H, this function computes an intensity ,I representing % the Pth percentile and returns the value in I. P must be in the % range [0,1] and I is returned as a value in the range [0,1] also. % Example: % I=percentile2i(h,0.5) % Check value of P if P<0 || P>1 error('The percentile must be in the range [0,1].'); end % Normalized the histogram to unit area.if it is already normalized % the following computation has no effect. h=h/sum(h); % Camulative distribution C=cumsum(h); % Calculations. idx=find(C >= P,1,'first'); % Subtract 1 from idx because indexing starts at 1,but intensities % start at 0. Also ,normalize to the range [0,1]. I=(idx-1)/(numel(h)-1);当背景照明高度不均匀时，需要进行阈值处理的难度就增大，为了解决这个问题，运用局部统计的可变阈值处理的算法进行解决。我们用一幅图像中每个点的邻域中像素的标准差和均值，这是局部对比度和平均灰度的描述子。
为了计算局部标准差，使用函数stdfilt进行操作：
g = stdfilt(f,nhood);为了计算局部均值，可以用函数localmean进行操作：
function mean = localmean(f,nhood) if nargin ==1 nhood = ones(3)/9; else nhood = nhood/sum(nhood(:)); end mean = imfilter(tofloat(f),nhood,'replicate');用于执行局部阈值处理，可以用下面的函数来操作：
function g = localthresh(f,nhood,a,b,meantype) f = tofloat(f); SIG = stdfilt(f,nhood); if nargin == 5 && strcmp(meantype,'global') MEAN = mean2(f); else MEAN = localmean(f,nhood); end g = (f > a*SIG) & (f > b*MEAN);实验：
[TGlobal] = graythresh(f); gGlobal = im2bw(f,TGlobal); subplot(131),imshow(gGlobal),title('用otsu方法分割的图像'); g = localthresh(f,ones(3),30,1.5,'global'); SIG = stdfilt(f,ones(3)); subplot(132),imshow(SIG,[]),title('局部标准差图像'); subplot(133),imshow(g),title('局部阈值处理后的图像');实验分析： 观察原图像发现，图像可以分为三种灰度级，发现除了细胞核的部分，细胞核与细胞的灰度级差别并不是很大，因此当我们使用全局阈值处理时，不能将里面的细胞核与细胞分离开。
移动平均是一种局部阈值处理方法，该方法以一幅图像的扫描行计算移动平均为基础。移动平均分割能减少光照偏差，且计算简单。当感兴趣的物体与图像尺寸相比较小(或较细)时，基于移动平均的阈值处理会工作的很好。打印图像和手写文本图像满足这一条件。
移动平均需要计算图像中每一个点，因此分割用下面的表示式： 其中，K是[0,1]范围内的常数，mxy是输入图像的移动平均
可以用如下的函数movingthresh来进行移动平均的操作：
function g = movingthresh(f, n, K) f = tofloat(f); [M, N] = size(f); if (n < 1) || (rem(n, 1) ~= 0) error('n must be an integer >= 1.') end if K < 0 || K > 1 error('K must be a fraction in the range [0, 1].') end f(2:2:end, :) = fliplr(f(2:2:end, :)); f = f'; % Still a matrix. f = f(:)'; % Convert to row vector for use in function filter. maf = ones(1, n)/n; % The 1-D moving average filter. ma = filter(maf, 1, f); % Computation of moving average. g = f > K * ma; g = reshape(g, N, M)'; g(2:2:end, :) = fliplr(g(2:2:end, :));示例：
f = imread('book.jpg'); subplot(131),imshow(f),title('原始图像'); T =graythresh(f); g1 = im2bw(f,T); subplot(132),imshow(g1),title('用otsu全局阈值分割后的图像'); g2 = movingthresh(f,20,0.5); subplot(133),imshow(g2),title('用移动平均局部阈值分割后的图像');4.1 区域生长
区域生长：根据预先定义的生长准则将像素或子区域组合为更大的区域的过程。
基本方法是从一组 “种子” 点开始，将与种子性质相似的那些邻域像素附加到每个种子上来形成这些生长区域。
语法：
[g, NR, SI, TI] = regiongrow(f, S, T) S： 可以是一个数组（与 f 大小相同）或一个标量。若 S 是一个数组，则它在种子点的所有坐标处必须包含 1，而在其他地方包含 0。如果 S 是一个标量，则它定义一个灰度值，那么 f 中具有该值的所有点都会成为种子点。T： 可以是一个数组或一个标量。若 T 是一个数组，则它对 f 中的每个位置包含一个阈值。若 T 是一个标量，则它定义一个全局阈值。S 和 T 所有值必须标定到区间 [0, 1]，且与输入图像的类无关。g： 是分割后的图像。NR： 是所找到的区域的数量。SI： 是包含种子点的一幅图像。TI： 是包含处理连接性前就已通过阈值测试的像素的一幅图像。使用区域生长检测焊接空隙
f = rgb2gray(imread('hanjie.jpg')); subplot(2,2,1),imshow(f); title('(a)显示有焊接缺陷的图像'); %函数regiongrow返回的NR为是不同区域的数目，参数SI是一副含有种子点的图像 %TI是包含在经过连通前通过阈值测试的像素 [g,NR,SI,TI]=regiongrow(f,1,0.26);%种子的像素值为255，65为阈值 subplot(2,2,2),imshow(SI); title('(b)种子点'); subplot(2,2,3),imshow(TI); title('(c)通过了阈值测试的像素的二值图像(白色)'); subplot(2,2,4),imshow(g); title('(d)对种子点进行8连通分析后的结果');分析：如上图，焊接裂缝处并没有显示，是因为焊缝缺损区域中的一些像素不是趋向于有最大的数字，所以在第一步确定初始种子点选择的点并不是我们感兴趣的。
regiongrow函数：
function [g, NR, SI, TI] = regiongrow(f, S, T) f = tofloat(f); if numel(S) == 1 SI =f == S; S1 = S; else SI = bwmorph(S, 'shrink', Inf); S1 = f(SI); end TI = false(size(f)); for K = 1:length(S1) seedvalue = S1(K); S = abs(f - seedvalue) <= T; TI = TI | S; end [g, NR] = bwlabel(imreconstruct(SI,TI));4.2 区域分离和聚合 将一幅图像细分为一组任意的不相交区域，然后聚合或分离这些区域。过程如下：
分离任意区域 Ri 为 4 个不相连的象限，满足 P( R i R_i Ri)=FALSE。当无法进一步分离时，合并任何满足 P( R j R_j Rj∪ R k R_k Rk)=TURE 的区域 R j R_j Rj和$ R_k$。在无法进一步合并的时候停止。在工具箱中执行四叉树分解的函数是qtdecomp,语法如下：
Z = qtdecomp(f, @split_test, parameters)为了在四叉树分解中得到实际的四叉区域像素值，使用qtgetblk函数，语法是：
[vals, r, c] = qtgetblk(f, Z, m)如果其中的两个区域分别满足属性，它们将被合并。我们可以调用的 splitmerge函数拥有如下调用语法：
g = splitmerge(f, mindim, @predicate) g 是输出图像，输出图像每个连接区域都使用不同的整数标注。mindim 定义分解中允许的最小块尺寸，必须是 2 的正整数次幂。 flag = predicate(region)如果在 region 区域中的像素满足函数中由代码定义的属性，那么函数就必须写成返回 true 的形式，否则，flag 的值就必须是false。
f = imread('starbulk.jpg'); subplot(2,3,1),imshow(f); title('(a)区域分割原始图像'); g64=splitmerge(f,64,@predicate);d代表分割中允许最小的块 subplot(2,3,2),imshow(g64); title('(b)mindim为64时的分割图像'); g32=splitmerge(f,32,@predicate);2代表分割中允许最小的块 subplot(2,3,3),imshow(g32); title('(c)mindim为32时的分割图像'); g16=splitmerge(f,16,@predicate);代表分割中允许最小的块 subplot(2,3,4),imshow(g16); title('(d)mindim为16时的分割图像'); g8=splitmerge(f,8,@predicate);%8代表分割中允许最小的块 subplot(2,3,5),imshow(g8); title('(e)mindim为8时的分割图像'); g4=splitmerge(f,4,@predicate);%4代表分割中允许最小的块 subplot(2,3,6),imshow(g4); title('(f)mindim为4时的分割图像');分水岭变换将找到灰度图像中的集水盆地和脊线。集水盆地就是我们要识别的物体或区域。
二值图像的距离变换是指从每个像素到最接近零值的像素的距离。
使用距离变换和分水岭变换分割二值图像 ：
f = tofloat(imread('plantcell.jpg')); g=im2bw(f,graythresh(f)); %把图像转换为二值图像 subplot(2,3,1),imshow(f);title('(a)使用距离和分水岭分割原图'); subplot(2,3,2),imshow(g);title('(b)原图像阈值处理后的图像'); gc=~g; %对图像求补 subplot(2,3,3),imshow(gc),title('(c)阈值处理后取反图像'); D=bwdist(gc); %计算其距离变换 subplot(2,3,4),imshow(D),title('(d)使用距离变换后的图像'); L=watershed(-D); %计算负距离变换的分水岭变换 w=L==0; %L 的零值即分水岭的脊线像素 subplot(2,3,5),imshow(w),title('(e)距离变换后的负分水岭图像'); g2=g & ~w; %原始二值图像和图像 w 的 “补” 的逻辑 “与” 操作可完成分割 subplot(2,3,6),imshow(g2),title('(f)阈值图像与分水岭图像相与图像');在使用针对分割的分水岭变换之前，常常使用梯度幅度对图像进行预处理。梯度幅度图像沿着物体的边缘有较高的像素值，而在其他地方则有较低的像素值。在理想的情况下，分水岭变换可得到沿物体边缘的分水岭脊线。
使用梯度和分水岭变换分割灰度图像 ：
f = tofloat(rgb2gray(imread('starbulk.jpg'))); subplot(2,2,1),imshow(f); title('(a)使用梯度和分水岭变换分割灰度图像'); h=fspecial('sobel'); fd=double(f); g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2); %计算图像的幅度梯度 subplot(2,2,2),imshow(g,[]); title('(b)使用梯度和分水岭分割幅度图像'); L=watershed(g); wr=L==0; subplot(2,2,3),imshow(wr); title('(c)对梯度复制图像进行二值分水岭后图像'); g2=imclose(imopen(g,ones(3,3)),ones(3,3)); L2=watershed(g2); %计算梯度的分水岭变换，并找到分水岭脊线。 wr2=L2==0; f2=f; f2(wr2)=255; subplot(2,2,4),imshow(f2); title('(d)平滑梯度图像后的分水岭变换');分析：如图C所示，分割的效果不是很理想，有太多我们不感兴趣的分水岭边界线出现。这是过分分割的原因。针对这一问题的解决方法是在计算分水岭变换之前先平滑梯度图像。如图 D显示了叠加结果。虽然增强图C的目的达到了，但仍然存在较多的附加脊线，对于决定哪些汇水盆地真正与感兴趣物体有关还是很困难。
由于噪声和梯度的其他局部不规则性，常会导致过分分割。一种方法是基于标记符修改梯度图像。标记符是属于一幅图像的连通分量。
示例：
f = imread('cell.jpg'); h=fspecial('sobel'); fd=double(f); g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2); L=watershed(g); wr=L==0; figure,subplot(2,3,1),imshow(wr,[]); title('(a)控制标记符的分水岭分割幅度图像'); rm=imregionalmin(g);%梯度图像有很多较浅的坑，造成的原因是原图像不均匀背景中灰度细小的变化 subplot(2,3,2),imshow(rm,[]); title('(b)对梯度幅度图像的局部最小区域'); im=imextendedmin(f,2);%得到内部标记符 fim=f; fim(im)=175; subplot(2,3,3),imshow(f,[]); title('(c)内部标记符'); Lim=watershed(bwdist(im)); em=Lim==0; subplot(2,3,4),imshow(em,[]); title('(d)外部标记符'); g2=imimposemin(g,im | em); subplot(2,3,5),imshow(g2,[]); title('(e)修改后的梯度幅度值'); L2=watershed(g2); f2=f; f2(L2==0)=255; subplot(2,3,6),imshow(f2),title('(f)最后分割的结果');
实验分析： 当直接使用梯度幅度方法进行分割处理时，发现图像过度分割非常严重，而这样的分割是没有意义的，并不是我们所需要的。当采用控制标记符分水岭分割算法时，能调整过度分割的界限，设置内部标记符和外部标记符，把所需要的内外部标记符联合起来，再次进行分割，因此可以得到我们所需要的分割效果。
