(数字图像冈萨雷斯第二版教材)纯手打！给个赞吧！

一、基本原理

• 图像的读取、存储操作：

i= imread('filename') ; imwrite(i,’image.jpg’);  

• 图像显示的⽅法及区别：

imshow(i);    imshow(i,[]);%0~255映射      imshow(i,[min max])%指定映射区间min到max
解释：见博客：https://blog.csdn.net/Michael__Corleone/article/details/68483407
 问题：在使用imshow(A)显示一张灰度图片时，显示出的是一张纯白的图片。（A为double类型的图像矩阵）
   原因：在matlab中，为了保证精度，经过了运算的图像矩阵A其数据类型会从unit8型变成double型。
imshow()显示图像时对double型是认为在0~1范围内，即大于1时都是显示为白色。
imshow显示uint8型时是0~255范围。
   使用imshow(A,[])，即可把图像矩阵A显示为正常的灰度图像。
   原理：imshow(A,[])是将A的最大值（max(A)）和最小值（min(A)）分别作为纯白(255)和纯黑(0)，
中间的K值相应地映射为0到255之间的标准灰度值，这样就可以正常显示了。
相当于将double型的矩阵A拉伸成为了0-255的uint8型的矩阵，因此就可以正常显示。

• 图像的索引操作

见数字图像处理第二版P14  索引=下标  

二、图像灰度变换和空间滤波

• 图像常见灰度变换及适用场合

• gamma变换

教材P24 
g=imadjust(f,[low high],[new_low new_high],gamma)
1.将f中灰度值映射为图g中新值。
2.low high （输入的值）在0到1之间
3.gamma>1突出亮区； gamma<1突出暗区

• log(突出暗部细节)

G = intrans(F, 'log', C, CLASS)
g=c*log(1+f) 

• stretch变换（灰度拉伸变换）

G = intrans(F, 'stretch', M, E)
g=1./(1+(m./f).^E)   %m是阈值在0到1之间  E是斜率  
使图像亮的地方更亮，暗的地方更暗，从而增加图像的可视细节

• 图像直方图处理

  • 直方图含义及特点

    imhist（f）
    

    含义: 反映图片包含各个灰度所占百分比 （王美丽啥都会经常叨咕的）
    特点：与空间位置无关，对几何攻击不敏感

  • 直方图均衡化的思想、步骤、实现

  处理骨骼挺好的。
  思想：把输入图像否认灰度级扩展到较宽灰度级范围来实现图像增强。
  步骤：新的灰度是原始灰度的累加，累计概率密度映射

  实现代码：g=histeq(f,256)
  

• 直方图匹配（规定化）的思想、步骤、实现

  g=histeq(f,p)
  

  答：生成具有特定直方图的图像方法，称为直方图匹配或直方图规定化
  1.直方图规定化(直方图匹配)是将变换过程加以控制，
  能够修正直方图的形状，或得到具有指定直方图的输出图像。
  有选择地增强某个灰度范围内的对比度或使图像灰度值满足某种特定的分布。
  2.直方图规定化是在运用均衡化原理的基础上，
  通过建立原始图像和期望图像(待匹配直方图的图像)之间的关系，
  使原始图像的直方图匹配特定的形状，
  弥补直方图均衡化不具备交互作用的特性。

• 空间滤波

• 空间滤波的处理思想：利用模板对图像进行卷积

• 均匀滤波：

  h = fspecial(type,parameter）；
  g = imfilter(f,w,’replicate’）；

平均滤波器:  h = fspecial('average',hsize）
高斯滤波器:  h = fspecial('gaussian',hsize,sigma)
圆盘滤波器:  h = fspecial('disk',radius)
运动滤波器:  h = fspecial('motion',len,theta)

• 锐化滤波

教材P50
sobel算子    :  h = fspecial('sobel')
prewitte算子 :  h = fspecial('prewitte')
log算子      :  h = fspecial('log',hsize,sigma)
lapalase算子 ： h = fspecial('laplacian',alpha)

• 统计排序滤波器(非线性空间滤波器)

g=ordfilt2(f,order,domain)
最小滤波器：g=ordfilt2(f,1,ones(m,n));
最大滤波器：g=ordfilt2(f,m*n,ones(m,n));
最知名的中值滤波器：g=ordfilt2(f,(m*n+1)/2,ones(m,n));
中值滤波器专用实现：g=medfilt2(f,[m n],padopt);
%数组[m n]定义一个大小为m*n的邻域（在该邻域上计算中值）
%padopt指定三个可能的边界填充选项之一：
%‘zeros’(默认值)，‘symmetric’指出f按照镜像反射方式沿边界扩展
%‘indexed‘表示若f是double类的则用1填充，否则用0填充。
默认形式：g=medfilt2(f);使用一个大小为3*3的邻域并用0填充边界计算中值。
中值滤波器对椒盐噪声处理很好。

三、图像傅里叶变换和频率与处理

如果看了此文你还不懂傅里叶变换，那就过来掐死我吧【完整版】

• 傅里叶变换的性质：

  答：频谱共轭对称性，线性，周期性，旋转不变性，时域卷积对应频域乘积

• 傅立叶图像反映的信息有什么物理含义：

  答：四个角亮代表图像高频信息，中间大部分代表低频信息。

• 怎么对图像实现傅立叶变换，怎么对频谱中心进行平移？

  答：F=fft2(f); F=fftshift(F)

• 频域内图像处理的流程是什么？

  答：
  1.使用paddedsize获得填充参数：PQ = paddedsize（size（f））；%如果输入是彩色图像，必须要灰度化rgb2gray。
  2.得到使用填充的傅里叶变化：F = fft2(f, PQ(1), PQ(2));
  3.生成一个大小为PQ(1) X PQ(2) 的滤波函数H。如果该滤波函数已居中，使用前要令H = fftshift(H)。
  4.将变换乘以滤波函数：G = H.*F;
  5.获得G的傅里叶逆变换的实部：g = real(ifft2(G));
  6.将左上部分的矩形剪切为原来尺寸大小：
  g = g(1:size(f,1), 1:size(f, 2));

f=imread('C:\Users\liulang\Desktop\test.jpg');g=rgb2gray(f);% 利用imread读出一副图像，显示图像 f = imread (‘….’)；
f  = im2double(g); % 将图片类型转成double类型;
 h=fspecial('average');% 利用fspecial生成一个空间滤波器，滤波器类型自定义 h= fspecial(…)；
PQ=paddedsize(size(f))% 利用paddedsize设计频率滤波器的大小 PQ=paddedsize(size(f))；
H=freqz2(h,PQ(2),PQ(1))% 利用freqz2生成相应的频域滤波器H = freq2(h, PQ(2),PQ(1))；
F=fftshift(fft2(f,PQ(1),PQ(2)))% 利用函数fft2计算图像的傅里叶变换 F = fft2 (f, PQ(1),PQ(2)))；要中心平移，产生的滤波器原点在矩阵中心处
f=imfilter(f,h);% 利用imfilter和h进行空间滤波；
J=F.*H;J=ifft2(J);J=abs(J);% 将F和H相乘，在做ifft2变换，再取实部，再进行剪切，得到频域滤波结果  
J=J(1:size(f,1),1:size(f,2));%size里面两个参数
figure,subplot(131),imshow(g,[]),title('灰度图')% 对比频域滤波和空间滤波的区别
subplot(132),imshow(J,[]),title('频域滤波') 
subplot(133),imshow(f,[]),title('空间滤波')

频率域滤波的基本步骤.jpg

• 频域内的低通和高通滤波器有几种类型？每种有什么特点？适用于什么场合？

  答：（都是频域平滑）
  理想低通：理想低通滤波器的平滑作用非常明显，但由于变换有一个陡峭的波形，它的反变换h(x,y)有强烈的振铃特性，使滤波后图像产生模糊效果。因此这种理想低通滤波实用中不能采用。
  会有【振铃现象】。
  巴特沃兹低通：在通过频率与截止频率之间没有明显的不连续性，不会出现“振铃”现象，其效果好于理想低通滤波器。
  阶处没有振铃，2阶有但看不出来，随着阶数升高，振铃越明显。
  高斯低通：高斯低通滤波器从通过频率到截止频率之间没有明显的不连续性，而是存在一个平滑的过渡带，无明显的振铃现象。
  高斯则无振铃，因为其DFT也是高斯。但高斯与2阶巴特沃斯相比，效果差点。
  因此，在需要严格控制低频和高频截止频率的过渡的情况下，巴特沃斯低通是最好的选择。除此之外，高斯最佳。
  高频滤波器：除了上述的三种之外，还有高频强调滤波器，由于高通滤波器偏离了直流分分量，从而把图像的平均值降低到了0,。一种补偿方法是给高通滤波器加上一个偏移量。若偏移量与滤波器乘以一个大于1的常数结合，则这种方法就称为高频强调滤波，因为该常量乘数突出了高频部分。这个乘数同时增加了低频部分的幅度，但是只要偏移量与被乘数比较小，低频增强的影响就弱于高频增强的影响。

四、图像复原

image.png

• 图像降质模型的空间域和频率域表示

答：退化后的图像=退化函数*原图像+一个加性噪声
空间域：卷积

空间域卷积操作：h是退化函数，f是原图像，＋噪声函数,g是退化后图片

频率域：乘积

频域乘积操作：H是退化函数，F是原图像，+噪声，G是退化后的图像

• 几种噪声模型：

高斯噪声，Rayleigh噪声，指数噪声，均匀噪声，salt-and-pepper噪声，伽马噪声。

六种

• 图像中噪声的估计流程

答：
1.均匀区域法是评估噪声最简单的方法。均匀区域法从图像中选择均匀区域, 通过计算这些均匀区域的标准差获取图像的噪声评估结果。需要说明的是，均匀区域是指该区域是平坦的，且在区域内像素值之间的差异不大。
2.分块法是假设图像由大量均匀的小块构成, 利用各小块的方差进行噪声估计。

• 图像中降噪的方法

答：空间滤波：均值滤波器，统计排序滤波器，自适应滤波器。
频率域滤波消除周期噪声：带阻滤波器，带通滤波器，陷波滤波器，最佳陷波滤波

• 均值滤波器：

  • 算术均值（窗口越大模糊越大）
  • 几何均值（比算数均值好些）
  • 谐波均值滤波器
    • 处理盐噪声比较好
  • 逆谐波均值滤波器（处理单值）
    • Q为正：处理椒噪声
    • Q为负：处理盐噪声
    • Q为-1：谐波滤波器
    • Q为0：算数均值滤波器

• 统计排序滤波器

  • 中值滤波器
    • 处理椒盐噪声最好；
  • 最大值最小值滤波器
    • 最大值滤波器处理椒噪声；
    • 最小值滤波器处理盐噪声；
  • 中点滤波器
    • 简单计算最大值最小值中点
    • 处理高斯噪声或均匀噪声
  • 修正的阿尔法均值滤波器
    • d=0: 算术均值滤波器
    • d=(mn-1)/2 ：中值滤波器

• 自适应滤波器

  • 自适应局部降低噪声滤波器
  • 自适应中值滤波器
    • 去除椒盐噪声
    • 平滑非椒盐噪声
    • 并减少诸如物体边界细化或粗化等失真

在自适应中值滤波算法中，A步骤里面会先判断是否满足Zmin<Zmed<ZmaxZmin<Zmed<Zmax。
这一步骤实质是判断当前区域的中值点是否是噪声点，
通常来说是满足Zmin<Zmed<ZmaxZmin<Zmed<Zmax这个条件的，
此时中值点不是噪声点，跳转到B；
考虑一些特殊情况，如果Zmed=ZminZmed=Zmin或者Zmed=ZmaxZmed=Zmax，
则认为是噪声点，应该扩大窗口尺寸，在一个更大的范围内寻找一个合适的非噪声点，
随后再跳转到B，否则输出的中值点是噪声点； 
接下来考虑跳转到B之后的情况：判断中心点的像素值是否是噪声点，判断条件为Zmin<Zxy<ZmaxZmin<Zxy<Zmax，
原理同上，因为如果Zxy=ZminZxy=Zmin或者Zxy=ZmaxZxy=Zmax，
则认为是噪声点。
如果不是噪声点，我们可以保留当前像素点的灰度值；
如果是噪声点，则使用中值替代原始灰度值，滤去噪声。

空间噪声滤波器函数
f = spfilt(g,type,m,n,parameter)
----<线性和非线性空间滤波器>----
%   f = spfilt(g,'amean',m,n)   算术均值滤波器
%   f = spfilt(g,'gmean',m,n)   几何均值滤波器
%   f = spfilt(g,'hmean',m,n)   谐波均值滤波器
%   f = spfilt(g,'chmean',m,n,q) 逆谐波均值滤波器 ,Q 缺省值为 1.5

%   f = spfilt(g,'median',m,n)  中值滤波器
%   f = spfilt(g,'max',m,n) 最大值滤波器
%   f = spfilt(g,'min',m,n) 最小值滤波器
%   f = spfilt(g,'midpoint',m,n)    中点滤波器
%   f = spfilt(g,'atrimmed',m,n,d)  修正的阿尔法均值滤波器

• 周期噪声处理方法：

  • 带通滤波器
    允许在一定频率范围内的信号通过而阻止其他频率范围内信号通过的滤波器。
  • 带阻滤波器
    阻止一定频率范围内的信号通过而允许其他频率范围内信号通过的滤波器。
    如果利用带通滤波器把某个频率带中的频率分量提取出来，然后将其从
    图像中减去，也可获得消除或减弱图像中某个频率范围内分量的效果。
  • 陷波滤波器
    可以阻止或通过以某个频率为中心的一段范围里的频率，本质上是带阻或带通滤波器，且可分别称为陷波带阻滤波器和 陷波带通滤波器。

• 退化系统估计的方法有什么？

  答：
  1.图像观察估计法：根据经验。
  2.图像试验估计法：使用与被退化图像设备相似的装置，并得到一个脉冲
  的冲激响应，可以进行较准确的退化估计。

  3. 模型估计法：建立退化模型，考虑引起退化的环境因素。

• 图像复原的方法及实现

  答：

  • 逆滤波

    • 退化函数H（u,v）的推测。
    • 尽可能不让噪声项N（u,v）影响画质。

逆滤波公式

• 半径截止的逆滤波：

  • 带半径的逆滤波的半径要试，不容易确定。

• 最小均方误差（维纳）滤波

  • 未退化图像和噪声的功率谱必须是已知的；
  • 功率比(信噪比)常数K 的估计一般还是没有合适的解。
  • 当噪声信息比NSR等于0时，此时维娜滤波等同于逆滤波。因此可以直接使用matlab自带deconvwnr函数，将第三个参数NSR设置成0即可。

J = deconvwnr(I,PSF)
J = deconvwnr(I,PSF,NSR)
J = deconvwnr(I,PSF,NCORR,ICORR)

I是退化图像
          PSF系统函数（点扩散函数）
          NSR信噪比
          NCORR：噪声的自相关函数
          ICORR：退化图像的自相关函数
          J：复原图像
clc,clear,close all;
i=imread('C:\Users\liulang\Desktop\wind.jpg')
figure,subplot(321),imshow(i,[]),title('原图你懂的'),
L=100;
T=11;
PSF=fspecial('motion',L,T);
blu=imfilter(i,PSF,'circular','conv');
 subplot(322),imshow(blu,[]),title('跑得快')
 noise = 0.1*randn(size(i)); % 生成噪声信号
blurredNoisy = imadd(blu,im2uint8(noise)); % 加入图像
 subplot(323),
imshow (blurredNoisy,[]);
title('有噪声并且运动模糊')
i1=deconvwnr(blu,PSF);
subplot(324),
imshow(i1),title('逆滤波-运动模糊')

subplot(325),
i2=deconvwnr(blurredNoisy,PSF),imshow(i2),title('逆滤波-有噪声有运动模糊')

NSR=sum(noise(:).^2)/sum(im2double(i(:)).^2);
i3=deconvwnr(blurredNoisy,PSF,NSR); 

subplot(326),imshow(i3),title('空间信噪比-维纳滤波？？？')

NP=abs(fftn(noise)).^2;%对噪声傅里叶变换 fftn   能量谱
NCORR=fftshift(real(ifftn(NP)));%噪声自相关函数，平不平移都可以
IP=abs(fftn(im2double(i))).^2;%图像功率谱
ICORR=fftshift(real(ifftn(IP)));%图像自相关函数
figure, i4=deconvwnr(blurredNoisy,PSF,NCORR,ICORR);
imshow(i4),title('自相关-维纳滤波？？？？')

五、图像重建

• 图像的投影变换

  • 1、R = radon(I, theta)
    返回亮度图像在角度theta下的Radon变换R。
    Radon变换是一幅图像在一个特定的角度下的径向线方向的投影。如果theta是一个标量，R则是一个包含在theta的列向量。如果theta是一个向量，R则是一个矩阵，据真的每一列是对应其中一个theta的Radon变换。如果忽略掉theta，则其默认为0：179.
  • 2、[R,xp] = radon(...)
    应于R中的每一行，返回一个包含径向坐标的向量xp。xp中的径向坐标是沿着X’轴的数值，其为在theta下，X’轴逆时针方向映射来的。两个坐标系的原点为图像的中心点，且为floor((size(me)+1)/2)，例如在一个20-by-30的图像中，其中心点为(10,15)。
    算法讲解：一幅图像的Radon变换是每一个像素Radon变换的集合。此算法首先将一个像素分成四个子像素，并且四个子像素分别投影，如下图所示：
    通过投影的位置以及bin的中心，每一个子像素的作用是按比例分裂两个最相邻的bins。如果子像素的投影与一个bin的中心相接，则坐标轴上的bin获得子像素的值，或者原始大小像素的1/4的值。如果子像素投影时，投影到了量bins之间，则子像素的值按比例分给两个bins。

• 图像重建的方法及实现

  • 直接反投影重建：

    • 出现的问题：星状伪影 图像模糊 振铃现象

  • 傅里叶反投影重建:

降低振铃现象
中心切片定理指出：f(x,y)在某一方向上的投影函数gθ(R)的一维傅立叶变换函数Gθ(ρ)是原函数f(x,y)的二维傅立叶变换函数F(ρ, θ)在(ρ, θ)平面上沿同一方向且过原点的直线上的值。
1.计算每个投影的一维傅里叶变换
2.用一个滤波函数|w|乘以每个傅里叶变换，就是加窗。
3.得到每个滤波后的一维反傅里叶变换。
4.将3得到的求和

• 滤波反投影重建

滤波反投影法（加汉明窗）复原的图像很好的消除了“晕环现象”。

'Ram-Lak'   矩形窗
'Shepp-Logan'   正弦窗  平滑了图像，损失了部分高频信息
'Cosine'   余弦窗
'Hamming'    通用hanmming窗  降低了高频噪声, 可得到Hamming滤波器和Hanning滤波器 

'Hann'
'none'

clear;close all; clc;
 g = phantom ('Modified Shepp-Logan',600);
subplot(241);imshow(g,[]);title ('original image')
theta = 0:0.5:179.5;
[R,xp] = radon(g,theta);
f3 = iradon(R,theta,'Shepp-Logan');
subplot(245);imshow(f3,[]);title ('Reconstructed Image with Shepp-Logan filter')
f4 = iradon(R,theta,'Cosine');
subplot(246);imshow(f4,[]);title ('Reconstructed Image with Cosine filter')
f5 = iradon(R,theta,'Hamming');
subplot(247);imshow(f5,[]);title ('Reconstructed Image with Hamming filter')

六、图像分割

• 基于不连续点的分割

  • 点检测：如果一个孤立点与它周围的点不同，则可以使用上述模板进行检测。
  • 点检测的另外一种方法是在mn大小的邻域中，
    找到其最大像素值点和最小像素值点，其差值大于阈值的那些点则可认为是图像中的孤立点。
    g= ordfilter2(f，mn，ones(m,n))-ordfilter2(f,1,ones(m*n));
    g= g>=T
  • 线检测，4个模板 6个【-1】和3个【2】的排列组合；检测水平数值45度。
• Canny算子的检测的流程及实现
  Step1：利用高斯滤波对图像平滑 （去噪）
  Step2：利用梯度算子的计算图像的梯度和梯度方向（求导）
  Step3: 进行非极大值抑制（初步得到边缘点）
  Step4: 双阈值连接边缘

clc;clear
f = imread('Fig1006(a)(building).tif');
subplot(231);imshow(f);title('原始图像')

gv1 = edge(f,'sobel','vertical');
subplot(232);imshow(gv1);title('sobel边缘')

gv2 = edge(f,'prewitt','vertical');
subplot(233);imshow(gv2);title('prewitt边缘')

gv3 = edge(f,'roberts','vertical');
subplot(234);imshow(gv3);title('robert边缘')

gv4 = edge(f,'log');
subplot(235);imshow(gv4);title('log边缘')

gv5 = edge(f,'canny');
subplot(236);imshow(gv5);title('canny边缘')

• Hough变换
  • 点线的对偶性：图像空间中的一条线对应Hough空间中的一个点
    图像空间中的一个点(x0, y0)能映射为Hough空间中的一条直线
    Hough空间中两条线的交点用来表示过点(x0,y0)和点(x1,y1)的直线
    把在图象空间中的检测问题转换到参数空间 里，通过在参数空间里进行简单的累加统计完成检测任务
    在Hough空间中找某些点，通过这些点的线数最多。如左图所示的A点和B点，分别有三条线通过这两点。
    X-Y平面中的一点对应参数空间的一条正弦曲线
    把在图像空间中检测直线的问题转化为在极坐标参数空间中找通过点(r,θ)的最多正弦曲线数的问题。

[H, theta, rho] = hough(BW)
[H, theta, rho] = hough(BW, Parameter1, Parameter2)

BW是二值图像
Parameter1:指定了theta方向的分辨率
Parameter2:指定了lamda方向的分辨率
H：hough变换空间矩阵
Theta和rho：为返回的theta和lamda值

peaks = houghpeaks(H, numpeaks)
H是hough变换矩阵空间
Numpeaks指定峰值的数目
Peaks是容纳峰值的行列坐标的矩阵

lines = houghlines(BW, theta, rho, peaks)
lines = houghlines(..., param1, val1, param2, val2)
从hough变换空间选出候选峰值后，利用峰值所对应的theta和rho值找到直线的位置。
clc;clear;close all;
f = imread('building.tif');
subplot(231);imshow(f);title('原始图像')
[bw,tc] = edge(f,'canny',[0.04 0.10],1.5);
subplot(232);imshow(bw);title('边缘检测结果')

[H,theta,rho] = hough(bw,1); 
subplot(233);imshow(H,[]);title ('hough变换空间')
axis on, axis normal
[r,c] = houghpeaks(H,5);
subplot(234);plot(theta(c),rho(r),'linestyle','none',...
    'marker','s','color','k')
title('带有所选5个峰值的位置的 Hough 变换')
lines = houghlines(f,theta,rho,r,c);
subplot(235)
imshow(f);hold on;
for k = 1:length(lines)
    xy = [lines(k).point1;lines(k).point2];
    plot(xy(:,2),xy(:,1),'LineWidth',4,'Color','r');
    hold on;
end
title('Hough 变换峰值对应的线段')

• 基于阈值的分割

  • 全局阈值

    • 均值迭代阈值选择分割与实现
      • 选择一个T的初始估计值
        用T分割图像，生成两组像素：G1由所有灰度值大 于T的像素组成，而G2由所有灰度值小于或等于T的 像素组成
        对区域G1和G2中的所有像素计算平均灰度值µ1和µ2
        计算新的阈值 T = ( µ1+ µ2)/2
        重复步骤2到4，直到逐次迭代所得的T值之差小于事先定义的参数T0

clc;clear;close all;
f = imread('fingerprint.tif');
 
subplot(221);imshow(f,[]);title(‘原图');
subplot(222);imhist(f);title(‘直方图');
T = mean2(f);
done = false;
while ~done
    g = f>=T;
    Tnext = 0.5*(mean(f(g)) + mean(f(~g)));
    done = abs(T - Tnext) < 0.5;
    T = Tnext;
end
g = f<=T;
subplot(2,2,3:4);imshow(g,[])
title(‘分割结果')

• 最大类间方差方法选择阈值分割与实现。
  Step1：将初始阈值T设为0，类间方差初始值为0；
  Step2：用T分割图像，生成两组像素：G1由所有灰度值大于T的像素组成，而G2由所有灰度值小于或等于T的像素组成
  Step3：对区域G1和G2中的所有像素计算µ1，µ2 , ω1, ω2 ，进而计算类间方差，
  Step4：如果类间方差大于初始值，则将该类间方差赋给初始类间方差
  Step5：将T遍历整个图像灰度，找到最大类间方差所对应的图像灰度值即为最优阈值

clc;clear;close all;
f = imread('small-blob.tif');
subplot(221);imshow(f,[]);title(‘原始图像');
subplot(222);imhist(f);title('直方图');
T = mean2(f);
done = false;
while ~done
    g = f>=T;
    Tnext = 0.5*(mean(f(g)) + mean(f(~g)));
    done = abs(T - Tnext) < 0.5;
    T = Tnext;
end
g = im2bw(f,T/255);
subplot(2,2,3);imshow(g,[]); title('迭代阈值处理后的图像')
th = graythresh(f);
g1 = im2bw(f,th);
subplot(2,2,4);imshow(g1,[]);title(‘Otsu阈值处理后的图像')

• 局部阈值分割

  • 局部统计阈值分割的思想与实现
    Step1： 计算每个像素邻域内的均值mxy和标准差σxy
    Step2：计算可变阈值Txy ：
    Txy = aσxy+bmxy
    Step3：根据Txy分割图像

clc;clear;close all;
f = imread('cell.tif');
f= rgb2gray(f);
subplot(221);imshow(f,[]);title('original image');
 
zones = ones(3,3);
mean = uint8(imfilter(double(f),zones/9));
std = uint8(stdfilt(f,zones));
subplot(222);imshow(std,[]);title('standard image');

a=10;b=0.5;
g = (f>a*std & f>b*mean);
subplot(223);imshow(g,[]);title('local segmented image');
th = graythresh(f);
g1 = im2bw(f,th);
subplot(224);imshow(g1,[]);title('Otsu segmented image')

• 移动平均阈值分割
  局部阈值处理方法的一种特殊情况是：沿着一幅图像的扫描线计算移动平均。典型的扫描是以zigzag模式逐线执行，进而减少照明偏差。令Zk+1表示在扫描顺序中，在第k+1步遇到的一个点。

clc;clear;close all;
f = imread('word.tif'); [M,N]=size(f);
f0= double(f);
subplot(131);imshow(f,[]);title('original image');

f1 = f0(:);
maf = ones(1,20)/20;
ma =filter(maf,1,f1);

k =0.5;
g = f1 > k*ma;
g = reshape (g,M,N);
subplot(132);imshow(g,[]);title(‘moving thresh segmented image');
 
TH = graythresh(f);
g1 = im2bw(f,TH);
subplot(133);imshow(g1,[]);title('otsu segmented image');

• 基于区域的分割

  • 区域增长算法：流程，会计算
    • 定义：区域生长是把图像分割成特征相似的若干小区域，比较相邻小区域的特征，若相似则合并为同一区域，如此进行直到不能合并为止，最后生成特征不同的各区域。也称为区域扩张法。
    • 选择合适的种子点
      确定相似性准则(生长准则)
      确定生长停止条件
      采用的判断准则是：如果所考虑的象素与种子象素灰度值差的绝对值小于某个门限T，则将该象素包括进种子象素所在的区域
    • 生长准则：基于区域灰度差 基于区域灰度分布统计性质 基于区域形状

[G, NR, SI, TI] = regiongrow(F, S, T). 
F：被分割的图像
S：种子点图像或者一个标量（灰度值，所有具有这个灰度的像素都被记为种子点）
T：阈值数组或者标量，确定相似准则。
G：输出的分割图像
NR：输出的区域个数
SI：种子点图像
TI：阈值测试图像
clear; close all; clc;

I = imread('Fig1014(a)(defective_weld).tif');
subplot(221);imshow (I,[]);title('original image');
[G, NR, SI, TI] = regiongrow(I, 255, 0.26*255);
subplot(222);imshow (SI,[]);title('seed image');
subplot(223);imshow (TI,[]);title('threshold image');
subplot(224);imshow (G,[]);title('segmented image');

• 区域分裂与合并：会计算怎么分裂怎么合并
  • 另外一种分割的想法: 先从整幅图像开始通过不断分裂，得到各个区域（在实际中，先将图像分成任意大小且不重叠的区域，然后再合并或分裂这些区域，以满足分割的要求），在这类方法中，常根据图像的统计特性设定图像区域属性的一致性测度
  • 算法步骤：
    （1）对任一区域Ri，如果P(Ri)=FALSE, 就将其分裂成不重叠的四等分
    （2）对相邻的两个区域Ri和Rj（它们可以大小不同，即不在同一层），如果条件P(Ri∪Rj)=TRUE，就将它们合并
    （3）如果进一步的分裂或合并都不可能，则结束

[vals, r, c] = qtgetblk(I, S, dim)
[vals, idx] = qtgetblk(I, S, dim)
vals：四叉树分解中尺寸为dim*dim大小的块的值；
r，c为块的左上角的行列坐标;
S:为qtdecomp返回的稀疏矩阵；
I：输入图像。

• 分水岭分割算法：
  • 基于距离的分水岭分割算法：会计算二值图像距离图像
    • 二值图像的距离变换：每一个像素到最近非零像素的距离。
      值为1的像素距离变换为0。

clear; close all; clc;
 
I = imread('Fig1022(a)(gel-image).tif');
 subplot(231);imshow (I,[]);title('original image');
th = graythresh (I);
bw = im2bw(I,th);
subplot(232);imshow(bw,[]);title('bw image');
subplot(233);imshow(~bw,[]);title('negative of bw image');
d = bwdist(bw);
subplot(234);imshow (d,[]);title('bw distance image');
L = watershed (-d);
w = L==0;
subplot(235);imshow (w,[]);title('water shed lines');
final = I; final(w)=255;
subplot(236);imshow (final,[]);title('superpose image');

  - 局域图像梯度的分水岭分割算法：掌握梯度图像平滑的方法：开运算和闭运算。

腐蚀后膨胀的过程称为开运算。用来消除小物体、在纤细点处分离物体、平滑较大物体的边界的同时并不明显改变其面积。
先膨胀后腐蚀的过程称为闭运算。用来填充物体内细小空洞、连接邻近物体、平滑其边界的同时并不明显改变其面积。
- 梯度幅度图像在沿对象边缘处有较高的值，而在其它地方值较小。
理想情况下，在沿对象边缘处产生分水岭。

clear;close all; clc;
I = imread('Fig1022(a)(gel-image).tif');
subplot(231);imshow (I,[]);title('original image');
th = graythresh (I);
bw = im2bw(I,th);
subplot(232);imshow(bw,[]);title('bw image');
subplot(233);imshow(~bw,[]);title('negative of bw image');
 
h= fspecial ('sobel');
f= double (I);
fx = imfilter (f,h);
fy = imfilter (f,h');
fm = hypot (fx,fy); 
subplot(234);imshow (fm,[]);title('gradient image');
L = watershed (fm);
w = L==0;
subplot(235);imshow (w,[]);title('water shed lines');
final = f; final(w)=255;
subplot(236);imshow (final,[]);title('superpose image');

• 基于图像标记的分水岭分割算法
  上述方法通常会由于噪声和其它诸如梯度的局部不规则性的影响造成过度分割。而过度分割足以令结果变得毫无用处。
  解决方法：通过一系列的预处理步骤，加入附加知识，限制允许存在的区域数目。
  标记符属于图像的连通分量。内部标记符标记感兴趣对象的内部；外部标记符标记背景区域。

clear;close all; clc;
I = imread('Fig1022(a)(gel-image).tif');
subplot(231);imshow (I,[]);title('original image');
 
th = graythresh (I);
bw = im2bw(I,th);
h= fspecial ('sobel');
f= double (I);
fx = imfilter (f,h);
fy = imfilter (f,h');
fm = hypot (fx,fy); 
subplot(232);imshow (fm,[]);title('gradient image');
rm = imregionalmin (fm);
subplot(233);imshow (rm,[]);
title('localmin of gradient image');
im = imextendedmin (f,2);
fim = f;
fim(im)=175;
subplot(234);imshow (fim,[]);title('inner marker image');
Lim = watershed(bwdist(im));
w = Lim==0;
g = imimposemin(fm,im|w);
 
L = watershed(g);
w2 = L==0
 
subplot(235);imshow (w2,[]);title('water shed lines');
final = f; final(w2)=255;
subplot(236);imshow (final,[]);title('superpose image');

七、彩色图像处理

• Matlab中彩色图像的表示方法

  • 二值图像(Binary images)
    灰度图像(Intensity images)
    RGB图像(RGB images)
    索引图像(Indexed images)
  • 在MATLAB中一幅彩色图像要么被当作RGB图像，要么被当作索引图像进行处理。

• 灰度图像转成彩色图像的方法及实现

• 这一方法的基本概念是对任何输入像素的灰度级执行三个独立的变换。然后，三个变换结果分别送入彩色电视监视器的红、绿、蓝通道。这种方法产生一幅合成图像，其彩色内容受变换函数特性调制。

• 伪彩色（pseudocoloring，也称为假彩色）处理
  定义：指将灰度图像转化为彩色图像，或者将单色图像变换成给定彩色分布的图像。

• 基本原理：将灰度图像或者单色图像的各个灰度级匹配到彩色空间中的一点，从而使单色图像映射成彩色图像。

• 强度分层技术是伪彩色图像处理最简单的例子之一。

 I=imread('i_peppers_gray.bmp'); 
GS8=grayslice(I,8);
GS64=grayslice(I,64); 
subplot(1,3,1), imshow(I), title('原始灰度图像');
subplot(1,3,2), subimage(GS8,hot(8)), title('分成8层伪彩色');
subplot(1,3,3), subimage(GS64,hot(64)), title('分成64层伪彩色');  

• 彩色图像的处理方法

  • 针对每个通道进行处理
  • 把一个像素点的灰度看成一个向量处理

• 彩色图像的分割方法实现

图像分割是把图像分成各具特性的区域并提取出感兴趣目标的技术和过程。

• 基于RGB色彩空间
  假设目标是在RGB图像中分割特殊彩色区域的物体,给定一个感兴趣彩色的有代表性的彩色点样品集，可得到一个彩色“平均”估计，这种彩色是我们希望分割的彩色。
  对一幅RGB彩色图像，选择要分割的区域，计算该区域中的彩色点的平均向量a。
  盒子的中点在a，它的尺度沿每一个RGB轴以沿相应轴的数据标准差的1.25倍选择。
  例如，令R代表样点红分量的标准偏差，aR代表平均向量a的红分量：
  （aR-1.25(theta)R，aR+1.25(theta)R），
  这里在整个彩色图像中编码每一点的结果为：如果点位于盒子表面或内部为白色，否则为黑色。

rgb=imread('flower608.jpg');
rgb1=im2double(rgb);
r=rgb1(:,:,1);
g=rgb1(:,:,2);
b=rgb1(:,:,3);
r1=r(129:256,86:170);
r1_u=mean(mean(r1(:)));
[m,n]=size(r1);
sd1=0.0;
for i=1:m
        for j=1:n
            sd1=sd1+(r1(i,j)-r1_u)*(r1(i,j)-r1_u);
        end
end
r1_d=sqrt(sd1/(m*n));
r2=zeros(size(rgb1,1),size(rgb1,2));
ind=find((r>r1_u-1.25*r1_d)&(r<r1_u+1.25*r1_d));
r2(ind)=1;
result = cat(3,r2,g,b);
Imshow(result,[])

• 基于HSI色彩空间
  HSI模型（HSI Model）是面向颜色处理的，用色调(Hue)、饱和度(Saturation) 、亮度(Intensity)来描述颜色。
  用色调和饱和度描述色彩，用亮度描述光的强度。
  这个模型有二个特点：
  （1）I分量与图像的彩色信息无关；
  （2）H和S分量与人感受颜色的方式是紧密相连的。
  这些特点使得HSI模型非常适合于借助人的视觉系统来感知彩色特性的图像处理算法。

S1=(S>0.3*(max(max(S))));
F=S1.*H;