数字图像及其基础处理|Image Processing

基本概念

图像

• 图：光的折射或者反射光的分布

• 像：人的视觉系统对图的感知

一幅图像可以定义为一个二元函数$f(x,y)$，其中$x$和$y$是空间坐标，幅值$f$是对应点的强度或灰度。当$x、y、z$为有限的离散值时，称为数字图像。

一幅图像可分解为许多个单元。每个基本单元叫做图像元素，简称像素。

数字图像处理|Image Processing 就是借助数字计算机对数字图像进行处理。

• 低级：图像处理 输入输出都是图像

• 中级：图像分析分割 输入图像 输出特征

• 高级：计算机视觉识别 输入图像 输出理解

应用领域

安防监控、图像搜索、工业视觉、人机交互、视觉导航、虚拟现实、生物医学、遥感测绘

图像类型

宇宙射线图像、伽马射线图像、X射线图像、紫外线图像、可见光图像、红外线图像、无线电波图像、交流电波图像

图像的表示

我们一般通过 单个传感器、条带状传感器或阵列式传感器 进行图像的获取。

而要产生一幅数字图像，就需要把连续感测的数据转为数字形式，包括两种处理：采样、量化

• 采样：对坐标值进行数字化
• 量化：对幅度值进行数字化

进而，我们利用上述的$2-D$ 函数对其进行表示，如下图所示：

特别注意，在图像的表示中，$x$ 和$y$ 的方向在不同的教材和资料下的定义有所不同。请以具体场景下的定义为准。

文件格式

表示图像常使用两种不同的形式：矢量表示、光栅表示（也称位图、像素图）。

• 矢量表达 中，图像是由一系列线段及其组合体表示的，矢量文件内有一系列命令和数据，计算机根据命令讲图像绘制出来。
• 光栅表达 与人对图像的直观理解一致，是由规则排列的图像单位组合而成。适用于色彩、阴影或形状变换比较复杂的图像，但是缺少了像素间关系的表示，并且限定了 空间分辨率。

常见的文件格式有：

1. BMP 格式：是 Windows 环境下的标准格式，BMP图像文件也称位图文件。
2. GIF 格式：一种公用图像文件格式，8 位存储，即最多只能存储 256 色图像。
3. TIFF 格式：一种独立于操作系统和文件系统的格式，便于在软件之间进行图像数据交换。
4. JPEG 格式：源自于对静止灰度或彩色图像的一种压缩标准，有损压缩时可节省巨大空间。

动态范围

出于存储和量化硬件的考虑，灰度级数通常取为 2 的整数次幂，即：$L=2^k$

灰度区间为$[0,L-1]$.

动态范围上限取决于饱和度，下限取决于噪声

• 一般高动态范围对应高对比度；

• 相反低动态范围图像沉闷，灰度被冲淡.

空间分辨率

空间分辨率是图像中最小可辩别细节的测度。常用的测度是 点数/英寸，单位为$dpi$。

灰度分辨率

灰度分辨率是指在灰度级中可分辨的最小灰度级变化，$bit$。(被量化为256级的图像，其灰度分辨率为$8\ bit$)

像素关系

邻接像素

在二值图像中，把具有 1 值的像素归诸于邻接像素，$V=\{1\}$。

在灰度图像中，$V$可以 是$[0,255]$之间灰度的子集。

对于一个像素点$p$，常见的相邻像素的定义有：

1. 4-邻域，记为$N_4(p)$。如果像素$q$ 在集合$N_4(p)$ 中，且有$f(p)\in V$ 和$f(q)\in V$，则像素$p$ 和$q$ 是 4 邻接的（下同）。
2. 对角邻域，记为$N_D(p)$。
3. 8-邻域，记为$N_8(p)$。

m邻接（混合邻接：可以消除8-邻接的二义性）

如果数值在$V$ 中的$q$ 满足：$q\in N_4(p)$ 或$q\in N_D(p)$ 且$\forall e\in N_4(p) \cap N_4 (q),f(e)\notin V$

则 像素$p$ 和$q$ 是 m邻接 的。

连通性

从坐标$(x_0,y_0)$ 的点$p$ 到坐标$(x_n,y_n)$ 的点$q$ 的数字通路由一系列坐标序列表示：

$$(x_0,y_0),(x_1,y_1),\dots,(x_i,y_i),\dots,(x_n,y_n).$$

其中，$\forall i$，有点$(x_i,y_i)$ 和点$(x_{i+1},y_{i+1})$ 都是邻接的。（$0\leq i\leq n-1$）

---

令$S$ 为图像中的一个像素子集。如果$S$ 的全部像素之间存在一个通路，则称两个像素$p$ 和$q$ 在$S$ 中是连通的。

对于$S$ 中任何元素$p$，在$S$ 中找到的能和$p$ 形成通路的像素的构成的像素集称为$S$ 的连通分量。

如果$S$ 仅有一个连通分量，即$S$ 中所有像素都互相连通，则集合$S$ 称为连通集。

像素距离

给定三个像素$p,q,r$ ，定义距离函数$D(p,q)$ 使得：

1. $D(p,q)\geq 0$，且 当且仅当$p=q$ 时，$D(p,q)=0$；
2. $D(p,q)=D(q,p)$；
3. $D(p,r)\leq D(p,q)+D(q,r)$.

在离散化数字图像中，对于给定像素点$p(x,y),\;q(s,t)$ 常用的距离函数有：

$$\begin{aligned} 1.&D_E(p,q)=\sqrt{(x-s)^2+(y-t)^2}\\ 2.&D_4(p,q)=|x-s|+|y-t|\\ 3.&D_8(p,q)=\max(|x-s|,|y-t|) \end{aligned}$$

上面三者分别是像素点$p,q$ 的欧氏距离（欧拉距离）、城区距离和棋盘距离。
如果把$p,q$ 的坐标视为列向量，那么上面三者还分别表示$p,q$ 的 二范数、一范数以及无穷范数。

更一般的，有 闵可夫斯基距离 定义如下：

$$D_w(p,q)=\left(|x-s|^w+|y-t|^w\right)^{1/w}$$

距离变换

对图像区域中的点进行距离变换后，变换得到的结果应正比于该点到区域边界的最近距离。
距离变换（Distance Transformation，DT）是一种特殊的变换，它将一副二值图像变换为了灰度图像。

给定一副$W\times H$ 的二值图像，设其边界集合$B$ 已知。则距离图是一个尺寸为$W\times H$ 的矩阵，矩阵中每一个元素（像素）$p(i,j)$ 的值为$DT(p)$。

下面给出一种基于迭代方法的距离图计算：

首先初始化距离图，此时迭代指数$t=0$：

$$DT^{(0)}(p)=\begin{cases}0,&p\in B\\\infty,&p\notin B\end{cases}$$

区域与边界

待更

图像读取与显示|_MATLAB

img = imread('文件路径') % 0~255 uint8型
img = double(imread('文件路径')) % 0~255 双精度double型
img = im2double(imread(‘文件路径’)) % 0~1 规定化 双精度double型

imshow(img) % double型图像会映射到[0,1] 因此不适用
imshow(img,[])
imagesc(img),colormap(gray(256)) % 上面两个语句等效

% axis off 或 axis(‘off’) :关闭所有坐标轴线、刻度标记和标签。
% axis tight:设置坐标轴的范围为数据的范围。
% axis equal 或axis(‘equal’) :表示x轴和y轴的单位长度相同。

其他

人类视觉系统

人眼的光感知器分为 锥状体(Cones) 和 杆状体(Rods).

锥状体对颜色高度敏感，锥状体视觉又被称为 明视觉 或 亮视觉；

杆状体对低光照度敏感，例如白天色彩鲜艳的物体在夜晚却没有颜色这种现象，称为 暗视觉 或 微光视觉。

---

亮度适应现象： 视觉系统不能同时在一个范围内工作，其可以同时辨别的不同强度级别的总范围与整个适应范围相比很小。

**亮度适应水平：**对于给定条件，视觉系统当前的灵敏度水平。

---

• 定义韦伯比$\Delta I_c/I$ 来表达亮度辨别能力

  • 韦伯比越小，意味着可以辨别小亮度的百分比变化，即表示“较好”的亮度辨别能力，反之
  • $\Delta I_c$是背景亮度为$I$时恰可察觉亮度差的$50\%$.

• 马赫带效应 视觉系统往往会在不同强度区域的边界处出 现“下冲”或“上冲”现象（毛边）

图像内插

内插 通常在图像的放大、缩小、旋转和几何校正等任务中使用；内插是用已知数据来估计未知位置的值的过程

图像放大

最近邻内插

• 有严重的直边失真

双线性内插

$$v(x,y)=ax+by+cxy+d$$

双三次内插

$$v(x,y)=\sum_{i=0}^3\sum_{j=0}^3a_{ij}x^iy^i$$

是商业图像编辑程序，如Adobe Photoshop，的标准内插方法

图像坐标变换

待更

灰度变换

对于空间域的处理，定义有表达式：

$$g(x,y)=T[f(x,y)]$$

其中，$T(·)$是在领域内对$f$进行处理的算子，可以应用于单个像素，也可以是一组像素。

当$g$只依赖于单个像素时，$T$也被称作灰度级映射函数。

此时，可将上述表达式简化为：$s=T(r)$。$r、s$分别是处理前后的像素值。

局部处理

图像反转

$$s=L-1-r$$

• 对灰度级为$[0,L-1]$的图像进行反转（也叫 反色或负片）
• 可用于增强图像暗色区域的白色或灰色细节

线性变换

$$s=kr+b$$

• 拉伸：当斜率$|k|>1$，经过线性变换的图像的特定灰度段得到拉伸，分布空间更宽。
• 挤压：当斜率$|k|<1$ ，经过线性变换的图像的灰度分布被挤压在更窄的空间中。
• 反色：斜率$k=-1$ ，截距$b=255$时，达到反色的效果。

通常来说，我们可以通过线性变换进行图像的归一化（ min-max 法），即：

$$s=\frac{r-r_{min}}{r_{max}-r_{min}}$$

在 MATLAB 中也可以轻松实现：

B = im2double(I);
Bmax = max(max(B));
Bmin = min(min(B));
new_I = (B-Bmin)/(Bmax-Bmin);

对数变换

$$s=c\log(1+r)$$

式中，$c$是一个常数，且假定有$r\geq0$。

• 可将输入值范围较窄的低灰度值映射为范围较宽的灰度级，对亮区压缩，暗区扩展，方便观察暗区的细节

幂律变换(伽马变换)

$$s=cr^\gamma$$

考虑到偏移情况，有时也定义为：$s=c(r+\varepsilon)^\gamma$.

$\gamma \lt 1$提高灰度级，在正比函数上方，使图像变亮

$\gamma \gt 1$降低灰度级，在正比函数下方，使图像变暗

• 用分数值$\gamma$将窄范围的暗输入值映射为较宽范围的输出值，将高输入值映射为较窄范围的输出值

• 可用于对显示器的伽马矫正

在MATLAB中，imadjust() 函数可以帮助调整图像的强度值或颜色图，以实现灰度变换。
具体如下：

J = imadjust(I,[low_in high_in],[low_out high_out],gamma)

把原始图像中的小于 low_in 的像素值映射到 low_out。把大于 high_in 的像素值映射到 high_out。gamma 则是伽马变换中的参数。

而 [low_in, high_in] 除了手动限制以外，可以由 stretchlim(I,Tol) 根据 Tol 的最低与最大对比度返回。（默认值： Tol = [0.01, 0.99]）

此外，还有：

imadjust(I) == imadjust(I, stretchlim(I,[0.01,0.99]), [0,1], 1)

分段线性阶梯量化

待更

全局处理

定义 灰度直方图 ：横坐标表示灰度值，纵坐标表示图像中具有该灰度值的像素的个数。

• 像素的空间信息丢失，只记录其灰度分布状况

直方图均衡化

在MATLAB中：
imhist() 函数可以用于显示一幅图的直方图；
histeq() 函数可以直接对一幅图进行直方图均衡化.

直方图均衡化是将原图的直方图通过变换函数修整为均匀的直方图，然后按均衡直方图修整原图像。

直方图均衡化是一种典型的通过对图像的直方图进行修正来获得图像增强效果的自动方法。

• 主要用于增强动态范围偏小的图像的反差

$$s_k=T(r_k)=(L-1)\sum_{j=0}^kp_r(r_j)=\frac{L-1}{MN}\sum_{j=0}^kn_j$$

其中，$n_k$表示灰度值为$r_k$的像素点个数，$p_k(r_k)$为归一化图像直方图，有$r_k\in[0,L-1]$，$k=0,1,...,L-1$.

一般地，我们再将结果进行处理，有$s_k\leftarrow[s_k]$，即四舍五入取整（保证在$[0,L-1]$内），进而再计算$p_s(s_k)$。

$p_s(s_k)$即是均衡化之后的直方图。

直方图规定化(直方图匹配)

直方图规定化：把直方图调整成尽可能满足我们期望的分布的方法

设原始图像的灰度概率密度函数为$p_r(r)$，目标图像的为$p_Z(Z)$，则有：

$$S=T(r)=\int_0^rp_r(r)dr\\ V=G(Z)=\int_0^Zp_Z(Z)dZ\\ Z=G^{-1}(V)$$

因此，其主要思想就是找到这样的变换函数$G$.

对于离散点，表达式如下：

$$s_k=T(r_k)=(L-1)\sum_{j=0}^kp_r(r_j)\\ s_k=G(z_q)=(L-1)\sum_{i=0}^qp_z(z_i)\\ z_q=G^{-1}(s_k)$$

一般地，我们不需要求反函数，只需将有限个点均计算出来后存放在一张表中，最后作映射即可。

其中，对于$G(z_q)$，我们同样可以要求有$G(z_q)\leftarrow[G(z_q)]$，即四舍五入取整（保证在$[0,L-1]$​内），进而得到规定化之后的直方图。

直方图统计量

待更

空间滤波

使用空间模板进行的图像处理，被称为空间滤波。模板本身被称为空间滤波器。

在数学上，有着相关和卷积的概念，特别地，关于一元函数的卷积，下面列出的本站文章也详尽阐述。

离散卷积算法与傅里叶

https://qslie.top/posts/edf5c9ad/

此处针对图像，我们引入二元函数的相关和卷积：

$$(f\star g)(x,y)=\sum_{s=-a}^a\sum_{t=-b}^bf(s,t)g(x+s,y+t)\\ (f* g)(x,y)=\sum_{s=-a}^a\sum_{t=-b}^bf(s,t)g(x-s,y-t)\\$$

$(f*G)$即卷积计算相当于将$G$旋转180度得到的$G'$与$f$进行互相关计算，即$(f\star G')$

其中，相关的计算可以简述为：对应加权求和并平滑，示意图如下：

通常，图像滤波则是指对一个图像$f(x,y)$应用特定的算子$G(s,t)$对其进行相关计算而得到.

但一般来说，我们还是把算子称作卷积核 ，这是因为卷积计算拥有分配律等等优良性质，所以我们再进行空域滤波时，我们先将算子旋转180度再与图像进行卷积，因为等效于算子直接与图像进行相关计算。

可以通过imfilter()进行滤波，而fspecial()函数提供有常见的滤波器（又称算子、卷积核、模板等）

平滑|低通滤波

作用：减弱或消除图像中的高频率分量，可用于消除图像中的噪声。

盒式核

高斯核

$$G(s,t)=Ke^{-\frac{r^2}{2\sigma^2}}$$

其中，$r=(s^2+t^2)^{1/2}$.

统计排序滤波|非线性

中值滤波器 对椒盐噪声有着优秀的降噪能力

最大值滤波器

最小值滤波器

锐化|高通滤波_边缘检测

作用：减弱或消除图像中的低频率分量，可使图像反差增加，边缘明显。

二阶导数|Laplace

$$g(x,y)=f(x,y)+c\left[\nabla^2f(x,y) \right]$$

其中，有：

$$\begin{aligned} \nabla^2f&=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}\\ \frac{\partial^2 f}{\partial x^2}&=f(x+1,y)+f(x-1,y)-2f(x,y)\\ \frac{\partial^2 f}{\partial y^2}&=f(x,y+1)+f(x,y-1)-2f(x,y) \end{aligned}$$

由此得到拉普拉斯核：

MATLAB中的使用

%拉普拉斯
lap_H = [-1,-1,-1;-1,8,-1;-1,-1,-1]; %模板可选
img_lap = imfilter(img,lap_H);

%或 img_lap = imfilter(img, fspecial(’laplacian’,alpha));
% 参数alpha用于控制算子形状，取值范围为[0,1]

梯度|Sobel

$$\nabla f=\left[g_x,g_y\right]^T=\left[\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right]^T\\ ||\nabla f||=M(x,y)=\sqrt{g_x^2+g_y^2}\approx|g_x|+|g_y|$$

Sobel核

MATLAB中的使用

% sobel梯度
sobel_H = [-1,-2,-1;0,0,0;1,2,1];
img_sobel = abs(imfilter(img,sobel_H))+abs(imfilter(img,sobel_H'));

%或 img_sobel = abs(imfilter(img,fspecial('sobel')'))+...
            %   abs(imfilter(img,fspecial('sobel')));

其他

混合空间增强

此处以《数字图像处理》(第四版)中，人体骨骼处理(图3.57)为例子，根据书本中边缘增强的基本思想进行复现。

其基本思路如下图所示：

于是我们通过以下matlab代码进行复现：

% 图3.57复现
close all;clear;clc;
img = im2double(imread('bonescan.tif'));%读取原图并转为double

%拉普拉斯
lap_H = [-1,-1,-1;-1,8,-1;-1,-1,-1];
img_lap = imfilter(img,lap_H);
% 拉普拉斯的灰度显示/缩放
img_lap_gray= (img_lap-min(img_lap(:)))./max(img_lap(:))*255;

% sobel梯度
img_sobel = abs(imfilter(img,fspecial('sobel')'))+...
            abs(imfilter(img,fspecial('sobel')));

% 盒式滤波器
box_H = 1/25*ones(5,5);
img_box = imfilter(img_sobel,box_H);

% 边缘增强/Laplace与平滑sobel相乘
img_st = immultiply(img_lap,img_box);

% 幂律变换
img_gamma = imadjust(img+img_st,[0 1],[0 1],0.5);

% 图像显示
figure(1)
subplot(1,3,1),imshow(img),title('原图')
subplot(1,3,2),imagesc(img_lap_gray),colormap(gray),axis equal,axis tight,axis off,title('拉普拉斯核处理')
% subplot(1,3,2),imshow(img_lap_gray),title('拉普拉斯核处理')
subplot(1,3,3),imshow(img+img_lap),title('原图+拉普拉斯')
figure(2)
subplot(1,3,1),imshow(img_sobel),title('sobel梯度处理')
subplot(1,3,2),imshow(img_box),title('盒式滤波器平滑')
subplot(1,3,3),imshow(img+img_box),title('原图+sobel平滑')
figure(3)
subplot(1,3,1),imshow(img_st),title('边缘增强/Laplace与平滑sobel相乘')
subplot(1,3,2),imshow(img+img_st),title('原图+增强')
subplot(1,3,3),imshow(img_gamma),title('γ=0.5幂律变换')

频率滤波

为什么要在频率域研究图像滤波？

1. 可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务，在频率域中变得非常普通

2. 滤波在频率域更为直观，它可以解释空间域滤波的某些性质

3. 可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为空间域滤波器的指导

4. 一旦通过频率域试验选择了空间滤波，通常实施都在空间域进行

卷积定理：

$$f(x,y)\otimes g(x,y)\Leftrightarrow F(u,v)G(u,v)$$

注：此处列出的是连续函数的卷积定理。式中，$\otimes$是卷积符号

于是有基本的频域滤波公式：

$$g(x,y)=\ Real\{\mathfrak I^{-1}[H(u,v)F(u,v)]\}$$

其中，$\mathfrak I^{-1}$是$IDFT$.

基本滤波步骤

1. 已知一幅大小为$M×N$的输入图像$f(x,y)$

2. 使用零填充、镜像填充、复制填充，形成大小为$P×Q$的填充图像，$P= 2M,Q = 2N$

3. 将$f_P (x,y)$乘以$(-1)^{x+y}$，使傅里叶变换位于$P×Q$大小的频率矩形中心

4. 计算步骤 3 得到的图像的DFT，即$F(u,v)$

5. 构建一个实对称滤波器传递函数$H(u,v)$，其大小为$P×Q$ ，中心在$(P/2,Q/2)$

6. 采用对应像素相乘得到$G(u,v)=H(u,v)F(u,v)$，即$G(i,k)=H(i,k)F(i,k)$

7. 计算$G(u,v)$的IDFT得到滤波后的图像$G_P (x,y)=\{real|\mathfrak I^{-1} [G(u,v)]\}(-1)^{x+y}$

8. 从$G_P (x,y)$的左上象限提取一个大小为$M×N$的区域，得到与输入图像大小相同的滤波后结果$g(x,y)$

同态滤波

陷波滤波

$$H_{NR}(u,v)=\prod_{k=1}^Q H_k(u,v)H_{-k}(u,v)$$

周期干扰

% 图4.65 & 4.66 复现
close all;clear;clc;
img = double(imread("cassini-interference.tif"));%原图img
figure
subplot(2,2,1);imagesc(img),axis equal,axis tight,axis off,colormap(gray),title('原图')
IMG = fftshift(fft2(img));%频谱IMG
subplot(2,2,2);imagesc(log(1+abs(IMG))),axis equal,axis tight,axis off,colormap(gray),title('频谱')
[M,N] = size(img);
H = ones(M,N);%初始化滤波函数
H([1:M/2-15],[N/2-2:N/2+2]) = 0;
H([M/2+15:M],[N/2-2:N/2+2]) = 0;
subplot(2,2,3);imagesc(H),axis equal,axis tight,axis off,colormap(gray),title('滤波函数')
NEW_IMG = IMG.*H;
new_img = real(ifft2(fftshift(NEW_IMG)));
subplot(2,2,4);imagesc(new_img),axis equal,axis tight,axis off,colormap(gray),title('结果')
figure
model = real(ifft2(fftshift(IMG.*(1-H))));%干扰模式
subplot(1,2,1);imagesc(1-H),axis equal,axis tight,axis off,colormap(gray),title('提取的带通滤波函数')
subplot(1,2,2);imagesc(model),axis equal,axis tight,axis off,colormap(gray),title('干扰模式')

图像复原

图像复原技术的主要目的是以某种预定义的方式来改进图像。

图像复原利用退化现象的先验知识来复原已退化的图像。

复原技术主要是对退化建模并应用逆过程来恢复原图像。

估计退化函数的几种方法：观察法、试验法、数学建模法

实验证明，当退化模型的噪声较小时即轻度降质时，采用逆滤波复原方法可以获得较好的结果。逆滤波与半径范围有关。根据情况分析，不是越大和越小越好。

维纳滤波|最小均方误差滤波：通过在滤波过程中调节K值以得到准最佳结果

中值滤波器 降低随机噪声，模糊度小，对椒盐噪声效果显著

最大值滤波器 对胡椒噪声效果极好

最小值滤波器 对盐粒噪声效果极好

中点滤波器 对高斯和均匀噪声效果极好

彩色图像处理

光度学

发光强度 从光源流出的能量总量，通常用瓦特($W$)来度量

光通量 观察者从光源感受到的能量，单位为流明数($lm$)

亮度 光感知的主观描述子，体现了强度的无色概念， 描述彩色感觉的参数之一

照度 被光线照射的表面上的照度用照射在单位面积上的光通量来衡量，单位是$lx$（勒[克斯]，也有 用$lux$的）

• $1\ lx = 1\ lm/m^2$

彩色模型

RGB |红绿蓝

针对彩色显示器和彩色摄像机

CMY/CMYK |青深红黄+黑

运用在大多数在纸上沉积彩色颜料的设备，如彩色打印机和复印机

HSI 模型

人类描述和解释颜色的方式

色调H：描述的一种纯色（纯黄色、 纯橙色或纯红色）的颜色属性

饱和度S：度量纯色被白光稀释的程度

亮度I：主观描述子，实际不可度量，体现了无色的强度概念，并且 是描述颜色感觉的关键因子之一

伪彩色处理

伪彩色图像处理：按照规定的准则对灰度值赋予颜色的处理。

伪彩色这个词是用于区分对单色图像赋予彩色的处理和对真彩色图像赋予彩色的处理。

假彩色的主要应用是可视化和解释单幅图像或一序列图像中的灰度事件。

灰度分层与彩色编码

以 焊接的 X 射线图像 为例，设定两个分层进行彩色编码：

% 图6.19 复现
close all;clear;clc;
img = double(imread("weldXray.tif"));%原图img
[h,w] = size(img);
new_img = zeros(h,w,3);
for i=1:h
    for j=1:w
        if img(i,j)<255
            new_img(i,j,1)=0;
            new_img(i,j,2)=64/255;
            new_img(i,j,3)=115/255;
        else
            new_img(i,j,2)=1;
            new_img(i,j,1)=1;
        end
    end
end
figure
subplot(1,2,1),imagesc(img),axis equal,axis tight,axis off,colormap(gray),title('原图');
subplot(1,2,2),imagesc(new_img),axis equal,axis tight,axis off,title('彩色编码');

彩色边缘检测

% 图6.44 复现
close all;clear;clc;
img = im2double(imread("lenna-RGB.tif"));%原图img

h1 = [-1 -2 -1;
    0 0 0;
    1 2 1];
h2 = [-1 0 1;
    -2 0 2;
    -1 0 1];

R = img(:,:,1);G = img(:,:,2);B = img(:,:,3);
Rx = imfilter(R,h1,'replicate');
Ry = imfilter(R,h2,'replicate');
Gx = imfilter(G,h1,'replicate');
Gy = imfilter(G,h2,'replicate');
Bx = imfilter(B,h1,'replicate');
By = imfilter(B,h2,'replicate');
gxx = Rx.^2+Gx.^2+Bx.^2;
gyy = Ry.^2+Gy.^2+By.^2;
gxy = Rx.*Ry+Gx.*Gy+Bx.*Gy;
A = 0.5.*(atan(2.*gxy./(gxx-gyy+eps)));
G1 = 0.5.*((gxx+gyy)+(gxx-gyy).*cos(2*A)+2.*gxy.*sin(2*A));
A = A + pi/2;
G2 = 0.5.*((gxx+gyy)+(gxx-gyy).*cos(2*A)+2.*gxy.*sin(2*A));
G1 = G1.^0.5;
G2 = G2.^0.5;
RGB_transf = max(G1,G2);
R_sobel = sqrt(Rx.^2+Ry.^2);
G_sobel = sqrt(Gx.^2+Gy.^2);
B_sobel = sqrt(Bx.^2+By.^2);
RGB_sobel = R_sobel+G_sobel+B_sobel;

figure
subplot(2,2,1),imagesc(img),axis equal,axis tight,axis off,colormap(gray),title('原图');
subplot(2,2,2),imagesc(RGB_transf),axis equal,axis tight,axis off,colormap(gray),title('RGB空间直接梯度');
subplot(2,2,3),imagesc(RGB_sobel),axis equal,axis tight,axis off,colormap(gray),title('三分量单独求梯度');
subplot(2,2,4),imagesc(abs(RGB_transf-RGB_sobel)),axis equal,axis tight,axis off,colormap(gray),title('作差对比');

形态学处理

腐蚀

去掉小于结构元素的物体，选择不同大小的结构元素，可以去掉大小不同的物体
如果两物体之间有细小的连通，当结构元素足够大时，可以将物体分开

膨胀

把图像周围的背景点合并到物体中。如果两个物体距离比较近，通过膨胀可能连通在一起
对于填补图像分割后物体中的空洞十分有用。

开运算

先腐蚀再膨胀，平滑了轮廓，断开了细连线，删除了细凸起部分

闭运算

先膨胀再腐蚀，连接了狭窄的断点，填充了细长的沟槽和小于结构元的目标

SE = sterl('disk',20);

new_I = imerode(I,SE); % 腐蚀
new_I = imdilate(I,SE); % 膨胀
new_I = imopen(I,SE); % 开运算
new_I = imclose(I,SE); % 闭运算

图像分割

图像分割准则：不连续性（点/线/边缘检测）和相似性（阈值处理、图像分割）

不连续性

点、线、边缘检测的思想和方法

a) 点检测：一阶/二阶导数，孤立点：$Laplacian$核

b) 线检测：带有方向性的核函数

c) 边缘检测：步骤、常用的边缘检测算子（梯度算子（如Sobel、prewitt、Roberts）、$LoG$核、canny算子）

高斯拉普拉斯函数$LoG$：

$$\nabla^2G(x,y)=\left(\frac{x^2+y^2-2\sigma^2}{\sigma^4}\right)e^{-\frac{x^2+y^2}{2\sigma^2}}$$

相似性|阈值处理

阈值处理有：全局阈值、局部阈值（可变阈值：动态阈值或自适应阈值）、多阈值

全局阈值

1. 初始化一个阈值$T$
2. 用$T$分割出$G_1、G_2$
3. 计算$G_1、G_2$的平均灰度值$m_1、m_2$
4. 更新阈值$T\leftarrow\frac 1 2(m_1+m_2)$，并循环步骤2到4，直到阈值收敛到精度$\Delta T$内为止

Otsu算法

1. 计算输入图像的归一化直方图$p_i,i=0,1,..,L-1$
2. 计算全局灰度均值$m_G$
3. 在$[0,L-1]$灰度范围内选择阈值$k$，进行如下处理：
   • 根据阈值，将图像分为$C_1,C_2$两类($C_1$中的像素点的灰度值在0~k之间);
   • 计算$C_1$的概率$P_1$;
   • 计算$C_1$的灰度均值$m_1=m$;
   • 计算类间方差$\sigma_B=P_1(m_1-m_G)^2+P_2(m_2-m_G)^2=\frac{(m_GP_1-m)^2}{P_1(1-P_1)}$.
4. 最优化最大类间方差$\sigma_B$，找到最优参数$k^*$.

$$k^*=\arg \max_{0\le k \le L-1}\sigma_B(k)$$

其他分割算法

• 区域生长
• 区域分离和聚合
• $k$均值聚类分割
• $SLIC/Ncut$ 超像素分割
• 形态学分水岭分割

参考资料

1. 《数字图像处理》（第四版），冈萨雷斯[美] 著
2. Java基础图像处理之混合空间增强 ，勇者无敌，博客园