数字图像处理-基础

一、闲谈

还是按照以往的阅读方式，先分析整章的框架，再逐点分析。本章从人类视觉系统剖析；光、电磁波的成像特点；数字图像的生成；图像基本操作；像素间基本关系以及本书的数学工具几方面进行展开。即图像的来源、获取方式与处理方式的介绍。

二、人类视觉系统

根据光线沿直线传播的原理以及凸透镜折射原理，物体对光的反射通过晶状体在视网膜上呈现倒立缩小的图像。人眼的光感受器有两类：锥状体和杆状体。锥状体对颜色敏感；杆状体没有色彩感觉，对低照明度敏感。

三、光和电磁波谱

彩色光源的质量有三种：发光强度，光通量和亮度。发光强度是从光源流出的能量总量。光通量用流明数（lm）来度量，用于表示观察者从光源感受到的能量。亮度是光感知的主观描绘方式，无法用于度量。通常场景的明暗程度多采用光通量来评估。

没有颜色的光称为单色光或无色光。单色光的唯一属性是它的强度或大小。因为感知单色光的强度从黑色到灰色变化，最后到白色，灰度级一词通常用来表示单色光的强度。

同样，电磁波图像的获取也是通过评估其能量强度实现的。

四、图像感知与获取

其本质是通过传感器获取能量值将其转化为电压，从而实现量化效果，为图像构建提供数据。

传感器获取方式有：

1. 单个传感器获取图像
2. 条带传感器获取图像
3. 传感器阵列获取图像

五、图像取样与量化

采样可以理解为将图像沿某一与坐标轴平行的数据等比例分割成若干份，而将每一份样本的幅值映射成数值的形式可以理解为量化。采样与量化后的数据更便于存储，

下图为提取图片第800行和第500列的像素值的变化情况。不难发现，以列提取，天空的渐亮而像素值也逐渐变高，在灯带处到达顶峰；地面突然变漆黑，像素值呈现断崖式下降。

贴一波代码：

image = Image.open("img/image.jpg")  # 替换为你的图像文件路径
image_gray = image.convert("L")  # 转换为灰度图像
row_index = 800 # 要获取的行索引
column_index = 500  # 要获取的列索引
image_array = np.array(image)
image_array[row_index, :] = (255,0,0) #将行高亮为红色
image_array[: ,column_index] = (255,0,0) #将列高亮为红色
new_image = Image.fromarray(image_array)
row_pixels = list(image_gray.getdata())[row_index * image_gray.width : (row_index + 1) * image_gray.width]
column_pixels = [image_gray.getpixel((column_index, y)) for y in range(image_gray.height)]

5.1. 数字图像的表示

令$f(s,t)$表示一幅具有两个连续变量$s$和$t$的连续图像函数，通过取样和量化，可以把该函数转换成数字图像。假设把这幅连续图像取样为一个二维阵列$f(x,y)$，该阵列包含有$M$行和$N$列，其中$(x,y)$是离散坐标。通常图像在任何坐标$(x,y)$处的值记为$f(x,y)$，其中$x$和$y$都是整数。由一幅图像的坐标张成的实平面部分称为空间域，$x$和$y$称为空间变量或空间坐标。

数值阵列通常用于处理和算法开发。以公式形式，可将一个$M\times N$的数值阵列表示为:
$$f(x,y)=\left[\begin{array}{lclc}f(0,0)& f(0,1)& \dots & f(0,N-1)\\ f(1,0)& f(1,1)& \dots & f(1,N-1)\\ & \dots \dots & & \\ f(M-1,0)& f(M-1,1)& \dots & f(M-1,N-1)\end{array}\right]$$
根据公式中的坐标可以很清楚的了解到数字图像的原点位于左上角，其中正x轴向下延伸，正y轴向右延伸。这样表示的好处基于：许多图像显示（譬如电视显示屏）扫描都是从左上角开始的，然后依次向下移动一行。更重要的事实是矩阵的第一个元素按照惯例应在阵列的左上角，这种表示方式使用的是我们熟悉的标准右手笛卡尔坐标系统。

同样采用上图，将其每个坐标的像素值大小定义为z轴的值（此处采用red通道的值），进而构建三维图像如下：

贴个代码：

image = o3d.io.read_image("img/image.jpg")
pixels = np.asarray(image)
point_cloud = o3d.geometry.PointCloud()
height, width, _ = pixels.shape
points = []
for i in range(height):
    for j in range(width):
        points.append([i, j, pixels[i, j, 0]])  # 使用R通道的像素值作为高度
point_cloud.points = o3d.utility.Vector3dVector(points)
o3d.visualization.draw_geometries([point_cloud])

量化时每个像素的值会定义为灰度级数跨越范围$[0,L-1]$中的整数。而灰度跨越的值域有一个非正式的别名：动态范围。动态范围定义为系统中最大可度量灰度与最小可检测灰度之比。上限取决于饱和度，下限取决于噪声。一幅图像中最高和最低灰度级间的灰度差定义为对比度。

出于存储和量化硬件的考虑，灰度级数通常会采用2的整数幂，即$L=2^k$，这也意味着每个像素需要占用$k$bit的存储。一张图像需要占用$M\times N\times k$bit的容量。拥有$L=2^k$灰度级的图像也被称为“K比特图像”。

5.2. 空间和灰度分辨率

空间分辨率与灰度分辨率，个人理解是采样与量化的产物。

空间分辨率是图像中可辨别的最小细节的度量。而图像分辨率是单位距离内可分辨的最大线对数量。线对数量是指在单位距离内能够分辨的最大特征线对数目。每个线对由一条明线和一条暗线组成，它们之间有明显的对比。线对数量越高，表示图像能够显示更多的细节和较小的特征。这与图像分辨率有关，因为较高的图像分辨率意味着图像中有更多的像素，可以更好地捕捉和表示细节。

灰度分辨率是指在灰度级中可分辨的最小变化。灰度中可分辨的真实变化不仅受噪声和饱和度的影响，也受人类感知能力的影响。

5.3. 图像内插

初始化了一张图像，在这个图像上显示不同内插方法的差异。

5.3.1. 最近邻内插

最近邻算法通过计算新图坐标点$(X_{dst},Y_{dst})$在原图中对应坐标$(X_{src},Y_{src})$，从而获取新坐标点对应的像素值。
$$\begin{gathered} X_{src}=X_{dst}\times (\frac{Widt{h}_{src}}{Widt{h}_{dst}}), \\ {Y}_{src}=Y_{dst}\times (\frac{Heigh{t}_{src}}{Heigh{t}_{dst}}) \end{gathered}$$
将图像放大1.5倍效果如下：

5.3.2. 双线性内插

新图中目标位置为$(x,y)$，四个邻近像素为$(x_1,y_1)$，$(x_1,y_2)$，$(x_2,y_1)$和$(x_2,y_2)$，对应的像素值为$p_{11}$，$p_{12}$，$p_{21}$和$p_{22}$。目标位置与邻近像素的关系如下：

• $(x_1,y_1)$：位于目标位置的左上方
• $(x_1,y_2)$：位于目标位置的右上方
• $(x_2,y_1)$：位于目标位置的左下方
• $(x_2,y_2)$：位于目标位置的右下方

(x1, y1) ------Q1------ (x1, y2)
  |            |           |
  |         (x, y)         |
  |            |           |
(x2, y1) ------Q2------ (x2, y2)

首先设定$(x,y)$在该区域水平和垂直方向上的相对位置$\alpha$与$\beta$：
$$\begin{gathered} \alpha =\frac{y-y_1}{y_2-y_1}, \\ \beta =\frac{x-x_1}{x_2-x_1}, \end{gathered}$$
然后计算$Q_1$与$Q_2$的像素值：
$$\begin{gathered} Q_1=\alpha p_{11}+(1-\alpha )p_{12}, \\ {Q}_2=\alpha p_{21}+(1-\alpha )p_{22} \end{gathered}$$
有了$Q_1$与$Q_2$的像素值后就可以计算目标点的像素值：
$$f(x,y)=\beta Q_1+(1-\beta )Q_2$$
整理得计算公式如下：
$$f(x,y)=\alpha \beta p_{11}+\alpha (1-\beta )p_{12}+(1-\alpha )\beta p_{21}+(1-\alpha )(1-\beta )p_{22}$$

当邻近像素没有值的时候，例如在边缘情况下，通常会采用边界处理策略来确定缺失像素的值。常见的边界处理策略有以下几种：

1. 边界复制（Edge replication）：使用最靠近缺失像素的边界像素的值来填充缺失像素。这意味着在边缘情况下，邻近像素的值会被复制到缺失的位置。
2. 常数填充（Constant filling）：将缺失像素的值设置为一个预定义的常数值，例如0或255。这种方法假设边缘外部的像素都具有相同的值。
3. 对称反射（Symmetric reflection）：通过对邻近像素的值进行对称反射，将缺失像素的值填充为反射后的像素值。这种方法假设图像在边缘外部是对称的。
4. 零填充（Zero filling）：将缺失像素的值设置为零。这种方法假设边缘外部的像素都具有零值。

将图像放大1.5倍效果如下：

5.3.3. 双三次内插

该方法通过计算16个最近邻点的加权值进而计算当前点的像素值，其公式如下：
$$f(x,y)=\sum _{i=0}^3\sum _{j=0}^3{\alpha }_{ij}{x}^i{y}^j$$
详细过程可以看这篇论文：Cubic convolution interpolation for digital image processing

将图像放大1.5倍效果如下：

简单观察会发现双三次内插与双线性内插效果似乎一致，实际上像素值会有一些不同，而且双线性插值其实是三次内插的特例形式。

六、像素间的关系

6.1. 相邻像素

符号：			  坐标：
[[a,b,c],		[[(x-1,y-1),(x-1,y),(x-1,y+1)],
 [d,e,f],		 [ (x,y-1) , (x,y) , (x,y+1) ],
 [g,h,i]]		 [(x+1,y-1),(x+1,y),(x+1,y+1)]]

在如上例子中，bh,df分别是e的垂直和水平的相邻像素，记作$N_4(p)$。acgi是e的对角相邻像素，记作$N_D(p)$。这八个点称为e的8邻域，记作$N_8(p)$。

6.2. 邻接性、连通性、区域和边界

6.2.1. 邻接性

将灰度值范围取一个子集$V$作为邻接像素。则邻接有以下三种定义：

4邻接：如果q在集合$N_4(p)$中，则具有$V$中数值的两个像素p和q是4邻接的。
8邻接：如果q在集合$N_8(p)$中，则具有V中数值的两个像素p和q是8邻接的。
m邻接（混合邻接）：以下情形可以认定p和q是m邻接的

1. q在$N_4(p)$中。
2. q在$N_D(p)$中，且集合$N_4(p)$∩$N_4(p)$中没有来自$V$中数值的像素。

6.2.2. 连通性

存在一条路线从一个像素到另一个像素，并且该路线上任意相邻两点均为邻接，则称这两个点是连通的。

若从一个像素出发能找到一条通路回到该点，则称这条通路是闭合通路。

令S是图像中的一个像素子集。如果 S的全部像素之间存在一个通路，则可以说两个像素p和q 在S中是连通的。对于S中的任何像素p，S中连通到该像素的像素集称为 S的连通分量。如果S仅有一个连通分量，则集合 S称为连通集。

6.2.3. 区域

令R是图像中的一个像素子集。如果R是连通集，则称R为一个区城。如果两个区域联合形成一个连通集，则区域尽$R_i$和$R_j$，称为邻接区域。不邻接的区域称为不连接区域。

6.2.4. 边界

假设一幅图有$K$个不连接的区域，且他们都不接触图像的边界。令$R_u$为$K$个区域的并集，令$R_u^c$为其补集，则称$R_u$中的所有点为图像的前景，$R_u^c$中的所有点为图像的背景。

区域$R$中与$R$的补集相邻的点称为称为$R$的边界（边缘或轮廓）。

6.3.距离度量

6.3.1. 欧式距离

像素$p(x,y)$与$q(s,t)$的欧式距离，即直角坐标系下两点间的距离：
$$D_e(p,q)=[(x-s{)}^2+(y-t{)}^2{]}^{0.5}$$

6.3.2. 城市街区距离（曼哈顿距离）

可以理解为通路步数的计算：
$$D_4(p,q)=|x-s|+|y-t|$$

6.3.3. 棋盘距离

可以理解为以$p$点为中心的矩形填充层数：
$$D_8(p,q)=max(|x-s|,|y-t|)$$

欧式距离给出的结果比较准确，但计算时要进行平方和开方运算，计算量大。 城市距离和棋盘距离均为非欧式距离， 计算量小，但有一定的误差。

七、工具介绍

7.1. 阵列与矩阵操作

阵列相乘：
$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\cdot \left[\begin{array}{cc}x& y\\ z& w\end{array}\right]=\left[\begin{array}{cc}ax& by\\ cz& dw\end{array}\right]$$
矩阵相乘：
$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]\cdot \left[\begin{array}{cc}x& y\\ z& w\end{array}\right]=\left[\begin{array}{cc}ax+bz& ay+bw\\ cx+dz& cy+dw\end{array}\right]$$

7.2. 线性操作与非线性操作

将算子定义为$H$,若图像$f(x,y)$经过算子$H$,得到$g(x,y)$:
$$H[f(x,y)]=g(x,y)$$
且满足：
$$H[a_i{f}_i(x,y)+a_j{f}_j(x,y)]=a_{i}H[f_i(x,y)]+a_{j}H[f_j(x,y)]=a_i{g}_i(x,y)+a_j{g}_j(x,y)$$
则称之为线性操作，否则为非线性操作。

7.3. 算术操作

图像间的算术操作属于阵列操作，
$$\begin{gathered} s(x,y)=f(x,y)+g(x,y) \\ d(x,y)=f(x,y)-g(x,y) \\ p(x,y)=f(x,y)\times g(x,y) \\ v(x,y)=f(x,y)÷g(x,y) \end{gathered}$$
先搞两张图片：

试试加法：

试试减法：

试试乘法：（这个效果不是很好，只能在英文部分略微看到一点黄色）

试试除法：

乘除法的效果都不是很好，究其原因是由于数值乘除后处于图像可表示范围的边界，对比度低，可视化效果较差。

7.4. 集合和逻辑操作

集合以及逻辑操作其实在高中就已经频繁接触，这里不多做赘述，直接上图片：

result_and = cv2.bitwise_and(image1, image2)
result_or = cv2.bitwise_or(image1, image2)
result_not1 = cv2.bitwise_not(image1)
result_not2 = cv2.bitwise_not(image2)
result_xor = cv2.bitwise_xor(image1, image2)
result_xnor = cv2.bitwise_not(result_xor)

7.5. 空间操作

空间操作可以分为三大类：1、单像素操作；2、邻域操作；3、几何空间变换。
关于单像素操作，文中提及的方案为：改变像素的灰度值，但没有具体展开，会在下一章提及，敬请期待吧。

关于邻域操作，可以理解为对某一点的邻域内的所有点执行某一函数得到的值，作为该点的新值。卷积和池化操作就属于邻域操作。

关于几何空间变换，书中公式有部分错误，此处参见wiki仿射变换

1.恒等变换

仿射矩阵$T$:$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$，对应的坐标公式为：$\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]$

2.尺度变换，通俗来讲就是缩放

仿射矩阵$T$:$\left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& 1\end{array}\right]$，对应的坐标公式为：$\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}{c}_x& 0& 0\\ 0& c_y& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{c}{s}_{x}v\\ s_{y}w\\ 1\end{array}\right]$

3.旋转变换

仿射矩阵$T$:$\left[\begin{array}{ccc}\cos (\theta )& -\sin (\theta )& 0\\ \sin (\theta )& \cos (\theta )& 0\\ 0& 0& 1\end{array}\right]$，

对应的坐标公式为：$\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}\cos (\theta )& -\sin (\theta )& 0\\ \sin (\theta )& \cos (\theta )& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{c}v\cos (\theta )-w\sin (\theta )\\ v\sin (\theta )+w\cos (\theta )\\ 1\end{array}\right]$

4：平移变换

仿射矩阵$T$:$\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\end{array}\right]$，对应的坐标公式为：$\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{c}v+t_x\\ w+t_y\\ 1\end{array}\right]$

5.偏移变换

仿射矩阵$T$:$\left[\begin{array}{ccc}1& \tan ({\theta }_x)& 0\\ \tan ({\theta }_y)& 1& 0\\ 0& 0& 1\end{array}\right]$，

对应的坐标公式为：$\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{ccc}1& \tan ({\theta }_x)& 0\\ \tan ({\theta }_y)& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}v\\ w\\ 1\end{array}\right]=\left[\begin{array}{c}v+tan({\theta }_x)w\\ tan({\theta }_y)v+w\\ 1\end{array}\right]$

来个中心旋转变换：

h, w = original_img.shape[:2]
#用于实现围绕中心点旋转
trans_mat = np.array([[1,0,int(h/2)],
                      [0,1,int(w/2)],
                      [0,0,1]])
trans_mat_inv = np.linalg.inv(trans_mat)

angle = 45  # 旋转角度为90度
angle_rad = np.deg2rad(angle)  # 将角度转换为弧度
#旋转矩阵
Tr = np.array([[np.cos(angle_rad), -np.sin(angle_rad), 0],
               [np.sin(angle_rad), np.cos(angle_rad), 0],
               [0, 0, 1]])

Tr_inv = np.linalg.inv(Tr)

out_img = np.zeros_like(original_img)
for x in range(h):
    for y in range(w):
        pos = np.dot(np.dot(np.dot(trans_mat,Tr),trans_mat_inv),np.array([[x],[y],[1]]))
        x1, y1 = int(pos[0]), int(pos[1])
        if 0 <= x1 < h and 0 <= y1 < w:
            out_img[x, y] = original_img[x1, y1]

而当放大的时候，由仿射矩阵计算出来的坐标并不能完全覆盖整张图片。这个时候前面学习的内插就派上了用途：

#仅放大
out_img = np.zeros_like(original_img)
Ts = np.array([[2, 0, 0], 
               [0, 1.5, 0], 
               [0, 0, 1]])
trans_mat = np.array([[1,0,int(h/2)],
                      [0,1,int(w/2)],
                      [0,0,1]])
trans_mat_inv = np.linalg.inv(trans_mat)

#找原图像在新图像上的对应点，将像素填充进去
for x in range(original_img.shape[0]):
    for y in range(original_img.shape[1]):
        pos = np.dot(np.dot(np.dot(trans_mat,Ts),trans_mat_inv),np.array([x,y,1]))
        x1, y1 = int(pos[0]), int(pos[1])  
        if 0 <= x1 < out_img.shape[0] and 0 <= y1 < out_img.shape[1]:
            out_img[x1, y1] = original_img[x, y] 

#最近邻插值
out_img = np.zeros_like(original_img)
Ts = np.array([[2, 0, 0], 
               [0, 1.5, 0], 
               [0, 0, 1]])
trans_mat = np.array([[1,0,int(h/2)],
                      [0,1,int(w/2)],
                      [0,0,1]])
trans_mat_inv = np.linalg.inv(trans_mat)
Ts_inv = np.linalg.inv(Ts)

#找新图像在原图像上的近似坐标，用其填充当前位置
for x in range(out_img.shape[0]):
    for y in range(out_img.shape[1]):
        pos = np.dot(np.dot(np.dot(trans_mat,Ts_inv),trans_mat_inv),np.array([x,y,1]))
        x1, y1 = int(pos[0]), int(pos[1])  
        if 0 <= x1 < original_img.shape[0] and 0 <= y1 < original_img.shape[1]:
            out_img[x, y] = original_img[x1, y1] 

如果你还记得线性代数，又恰好能看懂我的代码，那你很容易发现这两个过程其实就是等式两边同乘以仿射矩阵的逆。

7.6. 向量和矩阵操作

通常RGB图像在每个像素上包含红绿蓝三个分量，而多光谱图片包含n个分量。这就构成了三维向量（n维向量）的情况。那么两个像素点的像素向量的欧式距离的计算方式如下：
$$D(z,a)=[(z-a{)}^{⊺}(z-a){]}^{0.5}=[(z_1-a_1{)}^2+(z_2-a_2{)}^2+\dots \dots +(z_n-a_n{)}^2{]}^{0.5}$$
注意，这里是指像素向量的欧氏距离，与之前提到的坐标距离有所差异。

参见7.5的代码，虽然我们用的是彩色图像，但是写法方面只需考虑坐标，而省略了向量的计算。也表明了向量计算的支持使得算法的开发与实现更加便捷。关于向量的计算后续还会详细提及。

7.7. 图像变换

二维线性变换$T$的通式可表示为：
$$T(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)r(x,y,u,v)$$
其中$r(x,y,u,v)$是正变换核，$u$和$v$是变换变量。而反变换的通式为：
$$f(x,y)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}T(u,v)s(x,y,u,v)$$
那么$s(x,y,u,v)$就成为反变换核。

关于这部分的用法在第四章会有详细推导过程。

7.8. 概率方法

令$z_i,i=1,2,\dots \dots L-1$表示图像中所有可能出现的灰度值。那么通过计算$z_k$在图像中出现的次数$n_k$就能够估算$z_k$出现的概率：
$$p(z_k)=\frac{n_k}{MN}$$
那么每个灰度值的概率和必然为1。均值、方差均与概率论中计算公式类似，此处不过多叙述，主要是基本概念的介绍。

八、小结

本章介绍了很多图像处理的基础概念与基础计算方式，其功效如何，且看下回分解。