数字图像处理

1. 数字图像概述

1.1 图像的基本概念
图像因其表现方式的不同分为连续图像和离散图像两大类。

• 连续图像：指在二维坐标系中连续变化的图像，即图像的像点是无限稠密的，同时具有灰度值。又称为模拟图像。
• 离散图像：指用一个数字序列表示的图像。该阵列中的最小单位称为像素，每个像素由数字表示图像上每点的信息，又称为数字图像。

1.2 图像的数字化采样
目前的计算机只能处理数字信号，因此必须将连续的图像信息转化为数字形式，即需要数字化的采样与量化。

• 图像采样就是按照图像空间的坐标测量该位置上像素的灰度值。一般采用等间隔采样：将$(x,y)$平面上的图像分成均匀的小网格，采样值对应每个位置上网格的灰度值。
• 将模拟图像在空间上离散化为像素后，各个像素点的灰度值仍是连续量，此时需要把像素的灰度值进行量化，分为若干个灰度等级。

一幅图像经过采样和量化后便得到一幅数字图像，通常用一个矩阵来表示。矩阵中的元素称为像素。Matlab图像工具箱中像素值的表示支持各种数据类型：整数、浮点数、字符、逻辑算符，图像类型之间的转换函数有ind2gray、gray2ind、rgb2ind、ind2rgb、ntsc2rgb、rgb2ntsc等，可使用imtool命令查看一个图像文件的信息。

1.3 图像类型
在计算机中，按颜色和灰度的多少可以将图像分为二值图像、灰度图像、索引图像、RGB图像四种基本类型。

• 二值图像：矩阵元素仅由0、1两个值构成，0代表黑色，1代表白色。二值图像在Matlab中是一个取值只有0和1的逻辑数组，一个取值只有0和1的整型数组Matlab并不认为是二值图像。

B=logical(A)	%将数值型数组A转为逻辑数组B。数值0转为逻辑0，其他数值转为逻辑1

• 灰度图像：矩阵元素的整数取值范围通常为[0,255]，此时其数据类型一般为8位无符号整数（uint8）。若灰度图像的像素值是uint16类，则其整数取值范围为[0,65535]。若图像是double类，则像素取值范围为[0,1]。
• RGB彩色图像：若一幅RGB图像的数据类型是double，则它的取值范围就是[0,1]。类似地，数据类型也可以是uint8、uint16等。
• 索引图像：索引图像有两个分类，数据矩阵X和彩色映射矩阵map。map是一个大小为m×3且范围在[0,1]间由浮点数构成的double类数组，map的长度m同它所定义的颜色数目相等，map的每一行定义红、绿、蓝三个分量。每个像素的颜色由对应矩阵X的值作为指向map的一个指针决定。

2. 亮度变换与空间滤波

2.1 亮度变换
Matlab中使用函数imadjust对图像进行亮度变换。

g=imadjust(f,[low_in,high_in],[low_out,high_out],gamma)	%将图像f中的亮度值映射到g中的新值
%即将low_in至high_in之间的值映射到low_out至high_out之间的值
%参数gamma为调节权重，若gamma小于1则被映射到更亮的输出值，大于1责备映射到更暗的输出值
g=imadjust(f,[0,1],[1,0])	%图像亮度反转

2.2 线性空间滤波器
拉普拉斯滤波器增强图像的基本公式为$$g(x,y)=f(x,y)+c{\nabla }^{2}f(x,y)$$式中：f为输入图像；g为输出图像；c取1或-1。
对离散的数字图像，拉普拉斯算子${\nabla }^2$对图像f的作用就相当于如下矩阵${\mathbf{T}}_1$与f相乘$${\mathbf{T}}_1=\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]$$称${\mathbf{T}}_1$为滤波器或掩膜。或使用如下的矩阵${\mathbf{T}}_2$近似拉普拉斯算子$${\mathbf{T}}_2=\left[\begin{array}{ccc}1& 1& 1\\ 1& -8& 1\\ 1& 1& 1\end{array}\right]$$该矩阵更逼近二阶导数，对图像的改善作用更好。这些近似矩阵中心原始为负数时c=-1，反之c=1。
也可以选取其他的拉普拉斯近似矩阵，Matlab中函数fspecial(‘laplacian’,α)构造一个如下的拉普拉斯算子掩膜$$\left[\begin{array}{ccc}\frac{\alpha }{1+\alpha }& \frac{1-\alpha }{1+\alpha }& \frac{\alpha }{1+\alpha }\\ \frac{1-\alpha }{1+\alpha }& \frac{-4}{1+\alpha }& \frac{1-\alpha }{1+\alpha }\\ \frac{\alpha }{1+\alpha }& \frac{1-\alpha }{1+\alpha }& \frac{\alpha }{1+\alpha }\end{array}\right]$$

h1=fspecial('laplacian',0)	%构造一个T1滤波器
g1=f-imfilter(f,h1)		%相当于c=-1，即从原图像中减去拉普拉斯算子处理的结果

2.3 非线性空间滤波器
Matlab中非线性滤波器的函数为ordfilt2，它可以生成统计排序滤波器。该滤波器基于对图像邻域中所包含的像素进行排序，然后使用排序结果确定的值来替代邻域中的中心像素的值。

g=ordfilt2(f,order,domain)	%使用邻域的一组排序元素中的第order个元素来替代f中的每个元素，该邻域由domain中的非零元素指定
g=ordfilt2(f,1,ones(m,n))	%最小滤波器
g=ordfilt2(f,m*n,ones(m,n))	%最大滤波器
g=ordfilt2(f,median(1:m*n),ones(m,n))	%中值滤波器
%基于中值滤波器的实用性，工具箱提供了一个二维中值滤波函数medfilt2
g=medfilt2(f,[m,n])
g=medfilts(f)	$使用3*3的邻域来计算中值，并用0来填充输入图像的边界

3. 频域变换

为了快速有效地对图像进行处理和分析，如图像增强、图像分析、图像复原、图像压缩等，常需要将原定义在图像空间的图像以某种形式转换到频域空间，并利用频域空间的特有性质方便地进行一定的加工，最后再转换回图像空间，以得到所需要的效果。

3.1 傅里叶变换
傅里叶变换将图像从空域变换到频域。

3.1.1 二维连续傅里叶变换$$F(u,v)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }f(x,y)e^{-iux}{e}^{-iuy}dxdy$$式中，x和y为空间变量；u和v为频率变量；F(u,v)称为f(x,y)的频率表示。
二维傅里叶逆变换定义为$$f(x,y)=\frac{1}{4{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(u,v)e^{iux}{e}^{iuy}dudv$$

3.1.2 二维离散傅里叶变换（DFT）
令f(x,y)表示一幅大小为M×N的图像，其中$x=0,\dots ,M-1$和$y=0,\dots ,N-1$，二维离散傅里叶变换定义如下$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-i2\pi (ux/M+vy/N)}$$式中，$u=0,\dots ,M-1$和$v=0,\dots ,N-1$。由u和v定义的M×N矩形区域称为频率矩形，频率矩形的大小与输入图像的大小相同。
二维离散傅里叶逆变换为$$f(x,y)=\frac{1}{MN}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{i2\pi (ux/M+vy/N)}$$

3.1.3 基于离散傅里叶变换的频域滤波
在频域中滤波首先计算输入图像的傅里叶变换F(u,v)，然后用滤波器H(u,v)对F(u,v)作变换，最后对其得到的变换结果作傅里叶逆变换就得到频域滤波后的图像。对离散情况主要为以下5步：

1. 用$(-1{)}^{x+y}$乘以输入图像进行中心变换：$f_c(x,y)=(-1{)}^{x+y}f(x,y)$
2. 计算图像$f_c(x,y)$的离散傅里叶变换：$F(u,v)={\sum }_{x=0}^{M-1}{\sum }_{y=0}^{N-1}{f}_c(x,y)e^{-i2\pi (ux/M+vy/N)}$
3. 用滤波器H(u,v)作用F(u,v)：$G(u,v)=H(u,v)F(u,v)$
4. 计算G(u,v)的离散傅里叶逆变换，并取实部：$g_p(x,y)=\mathrm{Re}[F^{-1}[G(u,v)]]$
5. 用$(-1{)}^{x+y}$乘以$g_p(x,y)$得到中心还原滤波图像：$g(x,y)=(-1{)}^{x+y}{g}_p(x,y)$

理想低通滤波器具有传递函数$$H(u,v)=\{\begin{array}{l}1,D(u,v)⩽D_0\\ 0,D(u,v)>D_0\end{array}$$式中，$D_0$为指定非负数；D(u,v)为点(u,v)到滤波器中心的距离。理想低通滤波器相当于切断了圆外的所有分量，而圆上和圆内的分量不变。
n阶巴特沃兹低通滤波器的传递函数为$$H(u,v)=\frac{1}{1+[D(u,v)/D_0{]}^{2n}}$$该滤波器的传递函数并不是在$D_0$处突然不连续。对于平滑传递函数的滤波器，通常要定义一个截止频率。
高斯低通滤波器的传递函数为$$H(u,v)=\exp (-\frac{D^2(u,v)}{2D_0^2})$$

cf=fft2(cm)		%傅里叶变换
cm1=real(ifft2(cf1))	%傅里叶逆变换后取实部

3.2 离散余弦变换（DCT）
通过DCT变换可以将图像空间域上的信息变换到频率域上，较好地利用了人类视觉系统的特点。
一个M×N图像f(x,y)的二维DCT变换定义为$$F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)h(x,y,u,v)$$DCT逆变换公式为$$f(x,y)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}F(u,v)h(x,y,u,v)$$DCT变换在Matlab中有两种实现方法：

• 一种是基于快速傅里叶变换的算法，通过dct2函数实现，逆变换通过idct2函数实现。
• 另一种是DCT变换矩阵方法，通过dctmtx函数来计算变换矩阵，对于矩阵$\mathbf{A}$，若其DCT变换矩阵为$\mathbf{T}$，则A的二维DCT变换为$\mathbf{B}=\mathbf{T}\mathbf{A}{\mathbf{T}}^{\mathbf{T}}$，$\mathbf{B}$的二维逆DCT变换为${\mathbf{T}}^{\mathbf{T}}\mathbf{B}\mathbf{T}$。

T=dctmtc(8)		%生成8×8的DCT变换矩阵

对于通常的图像来说，大多数的DCT系数的值非常接近于0，如果舍弃这些接近于0的值，则在重构图像时并不会带来图像画面质量的显著下降。所以利用DCT进行图像压缩可以节省大量的存储空间。

3.3 图像保真度和质量
在图像压缩中为增加压缩率有时会放弃一些图像细节，这种情况下需要有对信息损失的测度以描述解码图像相对于原始图像的偏离程度，这些测度称为保真度准则。常用的保真度准则可分为两大类。

3.3.1 客观保真度准则
最常用的一个准则是输入图和输出图之间的均方根误差$$e_{rms}={\left\{\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}[\hat{f}(x,y)-f(x,y){]}^2\right\}}^{\frac{1}{2}}$$另一个客观保真度准则是均方信噪比（SNR）。

3.3.2 主观保真度准则：通过人主观观测的方法来评定图像的质量。

4. 数字图像的水印防伪

4.1 水印防伪概述
4.1.1 水印防伪特点
数字水印技术，指在数字化的数据内容中嵌入不明显的记号，从而达到版权保护或认证的目的。被嵌入的记号通常是不可见或不可察觉的，但通过一些计算操作可以被检测或提取。数字图像的内嵌水印必须具有下列特点：

• 透明性：水印后图像不能有视觉质量的下降，与原始图像对比，很难发现二者的差别。
• 鲁棒性：加入图像中的水印必须能够承受施加于图像的变换操作（如加入噪声、滤波、有损压缩等），不会因变换处理而丢失，水印信息经检验提取后应清晰可辩。
• 安全性：数字水印应能抵抗各种蓄意的攻击，必须能够唯一地标识原始图像的相关信息，任何第三方都不能伪造他人的水印图像。

4.1.2 水印防伪分类
针对图像水印技术的研究，主要体现在空间域和频率域两个层面。

• 空间域水印：将水印信息嵌入到载体图像的空间域特性上，如图像像素的最低有效位。
• 频率域水印：将水印信息嵌入到载体图像的变换域系数等特性上，如在图像的DFT、DCT系数上嵌入水印信息。

4.2 基于矩阵奇异值分解的数字水印算法
4.2.1 矩阵的奇异值分解（SVD）与图像矩阵的能量
矩阵的F范数定义为$$||\mathbf{A}|{|}_F^2=\mathrm{tr}({\mathbf{A}}^{\mathbf{T}}\mathbf{A})=\sum _{i=1}^r{\sigma }_i^2$$式中，${\sigma }_i$为$\mathbf{A}$的奇异值。矩阵的F范数的平方等于矩阵的所有奇异值的平方和。对于一幅图像，通常用图像矩阵的F范数来衡量图像的能量，所以图像的主要能量集中在矩阵那些数值较大的奇异值上。
Weyl定理指出，当图像被施加小的扰动时，图像矩阵的奇异值变化不会超过扰动矩阵的最大奇异值。即图像矩阵奇异值分解的稳定性非常好。想要将水印信息嵌入到图像矩阵的奇异值中，如果在嵌入过程中选择一个嵌入强度因子来控制水印信息嵌入的程度，那么当嵌入强度因子足够小时，图像在视觉上不会产生明显的变化。

4.2.2 水印嵌入
设一幅图像对应的矩阵$\mathbf{A}$大小为M×N，需要嵌入的水印对应的矩阵$\mathbf{W}$大小为m×n，在矩阵$\mathbf{A}$的左上角取一个大小为m×n的子块${\mathbf{A}}_0$，嵌入步骤如下：

1. 对${\mathbf{A}}_0$进行奇异值分解：${\mathbf{A}}_0={\mathbf{U}}_1{\mathbf{S}}_1{\mathbf{V}}_1^{\mathbf{T}}$，目标是将水印$\mathbf{W}$嵌入到矩阵${\mathbf{S}}_1$中。
2. 定义嵌入强度因子a，水印嵌入过程为：${\mathbf{A}}_1={\mathbf{S}}_1+\mathbf{a}\mathbf{W}$，矩阵${\mathbf{A}}_1$包含了所有的水印信息，水印信息的能量反映在${\mathbf{A}}_1$的奇异值中。
3. 对${\mathbf{A}}_1$进行奇异值分解：${\mathbf{A}}_1={\mathbf{U}}_2{\mathbf{S}}_2{\mathbf{V}}_2^{\mathbf{T}}$，矩阵${\mathbf{S}}_2$反映了嵌入水印的图像的全部信息。
4. 得到子块${\mathbf{A}}_0$嵌入水印后的图像子块：${\mathbf{A}}_2={\mathbf{U}}_1{\mathbf{S}}_2{\mathbf{V}}_1^{\mathbf{T}}$，完成了水印的嵌入。

嵌入强度因子a衡量了水印对原图像的扰动情况，a较小时有利于水印的透明性，但嵌入的水印信息易收到外界噪声的干扰；a较大时有利于增加算法的鲁棒性，但可能破坏水印的透明性，影响图像的质量。

4.2.3 水印提取
水印提取是水印嵌入的逆过程，设得到的是受扰动的图像矩阵${\mathbf{A}}_2^{\ast }$，提取步骤如下：

1. 对${\mathbf{A}}_2^{\ast }$进行奇异值分解：${\mathbf{A}}_2^{\ast }={\mathbf{U}}_3{\mathbf{S}}_2^{\ast }{\mathbf{V}}_3^{\mathbf{T}}$，得到包含有全部水印信息的奇异值矩阵${\mathbf{S}}_2^{\ast }$
2. 利用水印嵌入时得矩阵${\mathbf{U}}_2,{\mathbf{V}}_2$得到：${\mathbf{A}}_1^{\ast }={\mathbf{U}}_2{\mathbf{S}}_2^{\ast }{\mathbf{V}}_2^{\mathbf{T}}$
3. 由水印嵌入算法得：${\mathbf{W}}^{\ast }=\frac{1}{a}({\mathbf{A}}_1^{\ast }-{\mathbf{S}}_1)$，得到水印信息

4.3 基于DCT变换的水印算法
在图像的DCT系数上嵌入水印信息具有诸多优势：

• DCT变化是实数域变化，对实系数的处理更加方便，且不会使相位信息发生改变。
• DCT变换是有损图像压缩JPEG的核心，基于DCT变换的图像水印将兼容JPEG图像压缩。
• 图像的频域系数反映了能量分布，DCT变换后图像能量集中在图像的低频部分，即DCT图像中不为0的系数大部分集中在左上角，因此编码效率很高，将水印信息嵌入图像的中频系数上具有较好地鲁棒性。

4.3.1 水印嵌入算法
通过调整载体图像子块的中频DCT系数的大小来实现对水印信息的编码嵌入。算法描述如下：

1. 读取原始载体图像$\mathbf{A}$，对$\mathbf{A}$进行8×8分块，并对每块图像进行DCT变换。
2. 在8×8的子块中，中频系数的掩膜矩阵取$$\mathbf{H}=\left[\begin{array}{llllllll}0& 0& 0& 0& 1& 1& 0& 0\\ 0& 0& 0& 1& 1& 0& 0& 0\\ 0& 0& 1& 1& 0& 0& 0& 0\\ 0& 1& 1& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$每个子块的中频系数共有11个位置，每个子块的中频位置可以嵌入11个像素点的亮度值。对应地，将水印图像按照11个像素点一组进行分块，水印图像的最后一个分块如果不足11个像素点，则通过把亮度值置0进行扩充。
3. 对载体图像DCT系数进行修改，有$$g_i^{\prime }=g_i+\alpha f_i,i=1,\cdots ,11$$式中，$g_i$为载体图像中频系数的DCT值；$g_i^{\prime }$为变换后的DCT值；$\alpha$为水印嵌入的强度，取$\alpha =0.05$；$f_i$为对应位置的水印图像的灰度值（亮度值）。
4. 对修改后的DCT矩阵进行DCT逆变换，得到嵌入了水印的合成图像。

4.3.2 水印提取算法

1. 计算合成图像和原始载体图像的差图像$\Delta \mathbf{A}$
2. 对差图像$\Delta \mathbf{A}$进行8×8分块，并对每个小块作DCT变换
3. 从DCT小块中提取可能的水印序列$$f_i=(g_i^{\prime }-g_i)/\alpha ,i=1,\cdots ,11$$
4. 用下列函数计算可能的水印$\mathbf{W}$和原嵌入水印${\mathbf{W}}^{\ast }$的相关性$$C({\mathbf{W}}^{\ast },\mathbf{W})=\sum _{i=0}^{L-1}(f_i^{\ast }{f}_i)/\sqrt{\sum _{i=0}^{L-1}{f}_i^2}$$根据相似性的值就可以判断图像中是否含有水印，从而达到版权保护的目的。

5. 图像的加密和隐藏

图像隐藏技术是利用人类感觉器官的不敏感，以及多媒体数字信号本身存在冗余，将秘密信息隐藏于掩护体中，不被察觉到，且不影响载体的感觉效果，以此来达到隐秘传输秘密信息的目的。

5.1 加密算法
图像加密有很多方法，下面利用Henon混沌序列打乱图像矩阵的行序和列序，该序列为$$\begin{gathered} x_{n+1}=1-\alpha x_n^2+y_n \\ {y}_{n+1}=\beta x_n \end{gathered}$$式中，$\alpha =1,4$；$\beta =0.3$。
设保密图像$\mathbf{A}$的大小为$M_1\times N_1$，载体图像的大小为$M_2\times N_2$，为了方便处理，利用Matlab将载体图像处理成语保密图像同样的大小。记$L=\max (M_1,N_1)$，设$L=N_1$，使用密匙key=-0.400001作为混沌序列式的初始值$x_0,y_0$，生成两个长度为L的混沌序列$X_L,Y_L$，截取$X_L$的前$M_1$个分量，并把得到的子序列按照从小到大的次序排列，利用该子序列的排序地址打乱保密图像矩阵$\mathbf{A}$的行序；同样地，把混沌序列$Y_L$按照从小到大排序，利用该序列的排序地址打乱保密图像矩阵$\mathbf{A}$的列序。

5.2 图像的隐藏
上面所说的水印算法可以用于图像隐藏，除此之外，还有基于空域LSB的图像隐藏算法。
对于uint8格式的灰度图像，每个像素点有256个灰度级别，可以用8位二进制码来表示。由于人对视觉低位平面不敏感，所以可以利用低4位进行信息隐藏。LSB算法步骤如下：

1. 将载体图像的像素值转为二进制码，并将低4位置0
2. 将加密图像的像素值转为二进制码，并将低4位置0后右移4位
3. 将两图像处理后的像素值相加得到合成像素值，即把加密图像嵌入到了载体图像中