数字图像处理——章节总复习

数字图像基础

对应梁老师课件第二章，着重掌握基本概念。

图像的获取

略

图像的采样与量化

  • 大多数传感器的输出是连续电压波形
  • 为了产生一幅数字图像，需要把连续的感知数据转化为数字形式
  • 这包括两种处理：采样和量化

采样

  • 图像空间坐标的数字化（离散化）
  • 决定了图像的空间分辨率
  • 行采样像素\(M\)个，列采样\(N\)个构成\(M \times N\)的实数矩阵
  • 为编程方便采用矩阵坐标系
  • 1

量化

  • 图像像素灰度值的数字化（离散化）
  • 从模拟量到数字量（离散量）
  • 决定了图像的幅度（灰度级）分辨率

采样点数和量化级数的关系

对一幅图像，当量化级数一定时，采样点数对图像质量有着显著的影响。采样点数越多，图像质量越好； 当采样点数减少时，图上的块状效应就逐渐明显。

同理，当图像的采样点数一定时，采用不同量化级数的图像质量也不一样。量化级数越多，图像质量越好，当量化级数越少时，图像质量越差，量化级数最小的极端情况就是二值图像，图像出现假轮廓。

像素间基本关系

通路

2
2

距离

3
3

主要区分\(D_e\)欧式距离，\(D_4\)距离和\(D_8\)距离。

4
4
5
5
6
6

图像类型

  • 二值图像(Binary images)

    • 二值图像也叫黑白图像，编程实现中就指是图像象素只存在 0,1 两个值的逻辑数组。

  • 亮度图像(Intensity images)

    • 亮度图像是包含亮度级的图像，如 64 级，256 级等。
    • 如当亮度图像像素用 unit8 或 unit16 表示时，每个像素的整数取值范围分别是[0,255]和[0,65535]。
    • 如当亮度图像像素用 double 表示时，则像素的取值为浮点数，规定双精度归一化亮度图像的取值范围为[0,1]。

  • 索引图像(Indexed images)

    • 索引图像把像素值直接作为索引颜色的序号。
    • 根据索引颜色的序号就可以找到该像素的实际颜色。
    • 当把索引图像读入计算机时，索引颜色将被存储到调色板中。
  • RGB 图像(RGB images)
    • 一副 RGB 图像就是彩色像素的一个 M×N×3 数组，其中每一个彩色像素点都是在特定空间位置的彩色图像对应的红、绿、蓝三个分量。也可视为由三幅灰度图像形成的“堆”。
    • 这类图像不使用单独的调色板，每一个像素的颜色由存储在相应位置的红，绿，蓝颜色分量共同决定。
    • RGB 图像是 24 位图像，红绿蓝分量分别占用 8 位，理论上可以包含 16M 种不同的颜色。

空间域图像处理

对应梁老师课件第三章，王老师第二、三、四单元。着重掌握二值化，gamma 矫正，直方图均衡/匹配。

点运算

空间域增强：\(g(x,y) = T [ f(x,y )]\)

若将邻域限制为\(1 \times 1\)，则可简化为\(s=T(r)\)即点运算。

  • 点运算将输入图象映射为输出图象，输出图象每个象素点的灰度值仅由对应的输入象素点的值决定。
    • 常用于改变图象的灰度范围及分布
    • 也称为对比度增强、对比度拉伸或灰度变换；
  • 点运算可以是线性的，也可以是平方的,对数的,或其它任意单调函数的灰度变换；
  • 点运算可以利用一个 LUT（Look-up table）容易实现（或在彩色至少 R、G、B 三个 LUT）。

基本的变换

  • 图像反转
    • \(s = L-1-r\)
  • 对数运算
    • \(s = c \log(1+r)\)
    • 有时原图的动态范围太大，超出某些显示设备的允许动态范围，如直接使用原图，则一部分细节可能丢失。
    • 解决办法是对原图进行灰度压缩，如对数变换
  • 幂次变换
    • \(s = cr^\gamma\)
    • \(\gamma=1\)等幂变换与原图像相同。
    • \(\gamma>1\)用于过亮图像变暗。
    • \(\gamma<1\)用于过暗图像变亮。
    • 伽马矫正
  • 对比拉伸
    • 单调增函数保证灰度值变化保序。
    • 将需要提升差值的部分拉大灰度值极差从而达到提高对比度的效果。
  • 灰度切割
    • 提高图像中特定灰度范围的亮度
  • 位图切割

gamma 矫正和二值化的代码实现

  • gamma 矫正

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
      # 建表以后LUT
    def GammaTable(gamma):
      invGamma = 1.0 / gamma
      table = np.array([((i / 255.0) ** invGamma) * 255
                        for i in np.arange(0, 256)]).astype("uint8")
      return table
    def LUT(image , lutTable):
      # image: grayscale or RGB color image
      # luTable: [255,] 1D numpy array mapping 0-255 values to ohter values
      lut = lambda x: lutTable[x]
      return lut(image)
      # 直接伽马函数映射
    def DirectGammaFunc(image,gamma):
      invGamma = 1.0/gamma
      fuc =  lambda x: ((x/255.0)**invGamma)*255
      return fuc(image).astype("uint8")
      # R为单通道灰度图
      # 调用
      R_crrt = LUT(R,GammaTable(1.5))
      R_crrt = DirectGammaFunc(R,1.5)

  • 二值化（三种写法，见注释）

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    # 朴素循环
    def ImThresh(im, minv, maxv):
        BinImg = np.zeros(im.shape, dtype=im.dtype)
        for i in range(im.shape[0]):
            for j in range(im.shape[1]):
                if im[i,j]>=minv and im[i,j]<=maxv:
                    BinImg[i,j]=1
                else:
                    BinImg[i,j]=0
        return BinImg

    # 向量运算
    def ImThreshv2(image, minv, maxv):
        assert(len(image.shape)==2)
        # 二值逻辑
        group1 = image >= minv
        group2 = image <= maxv
        # 与操作 即要求满足 像素值>=minv 又要 <=maxv
        return (group1*group2).astype(np.uint8)

    # 杰哥的写法（也是向量运算）
    def ImThreshv3(image, minv, maxv):
        assert(len(image.shape)==2)
        # 只是这里用的非逻辑值，数字化为0，1（uint8类型）
        newimg = np.copy(image)
        newimg[newimg > maxv] = 0
        newimg[newimg < minv] = 0
        newimg[newimg != 0] = 1
        return newimg.astype(np.uint8)

    # 调用
    Threshimg = Imthresh(R,120,256)
    Threshimg2  = Imthreshv2(R,120,256)

直方图均衡

直方图

  • 定义：将所收集的测定值或数据之全距分为几个相等的区间作为横轴,并将各区间内之测定值所出现次数累积而成的面积,用柱子排起来的图形；
  • 表示图像中具有某种属性（如灰度、颜色等）的像素的个数,反映了图像中每种属性级出现的频率,是图像的基本统计特征之一。

均衡

  • 直方图均衡化是通过灰度变换将一幅图象转换为另一幅具有均衡直方图，即在每个灰度级上都具有相同（离散情况下是相近）的象素点数的过程。这样就增加了像素灰度值的动态范围，从而达到增强图像整体对比度的效果
  • 使用的方法是灰度级变换：\(s = T(r)\)

为什么概率累计曲线可以用来实现直方图均衡化的效果

概率论基础（问题基础）

根据课本 P74 页，\(p_r(r)\)和\(p_s(s)\)分别表示随机变量\(r\)和\(s\)的概率密度函数，概率论的一个基本结果是，如果\(p_r(r)\)和\(T(r)\)已知且\(s = T(r)\)具备单调性，在感兴趣的值域上是连续且可微的，则变换后的变量\(s\)的概率密度函数可由下面的简单公式得到：

\[p_s(s) = p_r(r) \left| \frac{\mathrm{d}r}{\mathrm{d}s} \right|\]

这个式子怎么推导的呢？实际上我们就是在求一个随机变量函数的概率分布，概率论是学过的，这里回忆一下给出推导。

设\(r\)为随机变量，其概率分布符合\(p_r(r)\)，另有一单调的函数（变换）\(s = T(r)\)，现在求\(s\)的概率分布，设\(r\)值域均为\([0,L-1]\)。

由分布函数公式

\[F(x) = \int_{0}^{x} p(m) \mathrm{d}m\]

得\(r\)得分布函数

\[F_r(r) = \int_{0}^{r} p_r(m) \mathrm{d}m \]

又因为\(s = T(r)\)为单调函数，存在反函数\(r = T^{-1}(s)\)

\[F_s(s) = F_r(T^{-1}(s)) = \int_{0}^{ T^{-1}(s) } p_r(m) \mathrm{d}m \]

概率密度函数为分布函数的导数

\[p_s(s) = \left| \frac{\mathrm{d} F_s(s)}{\mathrm{d}s} \right|\]

应用链式求导法则

\[p_s(s) = \left| \frac{\mathrm{d}F_s(s)}{\mathrm{d}s} \right|\]

\[p_s(s) = \frac{\mathrm{d}F_r(r)}{\mathrm{d}r} \left| \frac{ \mathrm{d} r}{ \mathrm{d}s} \right|\]

\[p_s(s) = p_r(r) \left| \frac{ \mathrm{d} r}{ \mathrm{d}s} \right|\]

则得到了书本上的公式。

均衡变换（回答问题）

问题即为，为何任何随机变量的概率分布函数，若变换为累积分布函数，则变换后的随机变量的概率分布为均匀分布？这实际上是累积分布函数的一个特殊的性质，下面给出证明。

  1. 若\(T(r)\)恰为\(r\)的分布函数，则 \[s = T(r) = \int_{0}^{r} p_r(m) \mathrm{d}m \] \[\frac{\mathrm{d}s}{\mathrm{d}r} = p_r(r)\] \[p_s(s) = p_r(r) \left| \frac{ \mathrm{d} r}{ \mathrm{d}s} \right| = p_r(r) \frac{1}{p_r(r)} = 1\]

    则\(s\)为均匀分布\(s \in [0,1]\) 且\(p_s(s)\)值恒为\(1\)

  2. 若\(T(r)\)为\(r\)分布函数映射到\([0,L-1]\)，则 \[s = T(r) = (L-1) \int_{0}^{r} p_r(m) \mathrm{d}m\] \[p_s(s) = p_r(r) \frac{1}{(L-1) p_r(r)} = \frac{1}{L-1}\] 则\(s\)为均匀分布，且\(p_s(s)=\frac{1}{L-1}\)

这便是我们想要的函数变换，无论\(p_r(r)\)分布如何，经过这种类分布函数的函数变换，得到的\(p_s(s)\)就是一个均匀分布。

直方图均衡化代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def CaculateHistogram(input_image):
    # 参数1. 图像
    # 输出1. 总像素值
    # 输出2. 灰度值分布计数
    # 输出3. 灰度值分布前缀和（用来做后续直方图平局奴化）
    # 区分单通道还是三通道图像
    if len(np.shape(input_image)) == 3:
           height,width,level = np.shape(input_image)
           summ = height*width*level
    else :
           height,width = np.shape(input_image)
           summ = height*width
    caculate_num,index_x = np.histogram(input_image,np.arange(0,256))
    caculate_num = np.append(caculate_num,1)
    # 前缀和
    sum_num = np.copy(caculate_num)
    for i in range(1,256):
        sum_num[i] = sum_num[i-1] + sum_num[i]
    return summ,caculate_num,sum_num

# 计算直方图数组
def cumsum(img, bins):
    histogram = np.zeros(bins)
    for pixel in np.arange(0, bins, 1):
        # 向量化计算像素值为piexel的像素数量
        histogram[pixel] += len(img[img==pixel])
    return histogram

# 直方图均衡接口函数，利用了LUT
def HistogramEqualizationLUT(input_image):
    size,data,data_sum = CaculateHistogram(input_image)
    fxy  = lambda x: (255*data_sum[x])//size
    table = np.array([fxy(i) for i in range(256)])
    lut = lambda x: table[x]
    return lut(input_image),table

直方图匹配

方法：求两图均衡化，有逆函数存在，可实现匹配。

证明

我们易知\(s\)服从均匀分布，我们寻找一个\(z\)随机变量符合另一特殊分布\(p_z(z)\)，易知存在一函数（变换）：

\[s = G(z) = (L-1) \int_{0}^{z} p_z(m) \mathrm{d}m\]

由于\(G(z)\)为累积分布函数，单调，故存在反函数\(G^{-1}(z)\)，则有新变换 \(N = G^{-1} \cdot T\)，使得\(r\)随机变量经\(N\)变换后符合\(z\)的分布特征。

\[z = G^{-1}(s) = G^{-1}(T(r)) = N(r)\]

由此，我们找到了某值域内求从某一分布随机变量\(r\)变换到指定分布随机变量\(z\)的方法。

空间域滤波基础

书本 P93 页有详细概念。

7
7
8
8

平滑和锐化

公式太多就贴 PPT 了

9
9
10
10
11
11
12
12
13
13
14
14
15
15
16
16

梯度算子

17
17
18
18

添加噪声代码实现

注意椒盐

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# 参考资料上的
def noisy(noise_typ,image):
    if noise_typ == "gauss":
        row,col,ch= image.shape
        mean = 0
        var = 0.1
        sigma = var**0.5
        gauss = np.random.normal(mean,sigma,(row,col,ch))
        gauss = gauss.reshape(row,col,ch)
        noisy = image + gauss
        return noisy
    elif noise_typ == "s&p":
        row,col,ch = image.shape
        s_vs_p = 0.5
        amount = 0.004
        out = np.copy(image)
        # Salt mode
        num_salt = np.ceil(amount * image.size * s_vs_p)
        coords = [np.random.randint(0, i - 1, int(num_salt)) for i in image.shape]
        out[coords] = 1
        # Pepper mode
        num_pepper = np.ceil(amount* image.size * (1. - s_vs_p))
        coords = [np.random.randint(0, i - 1, int(num_pepper)) for i in image.shape]
        out[coords] = 0
        return out
    elif noise_typ == "poisson":
        vals = len(np.unique(image))
        vals = 2 ** np.ceil(np.log2(vals))
        noisy = np.random.poisson(image * vals) / float(vals)
        return noisy
    elif noise_typ =="speckle":
        row,col,ch = image.shape
        gauss = np.random.randn(row,col,ch)
        gauss = gauss.reshape(row,col,ch)
        noisy = image + image * gauss
        return noisy

# 自己写的
# 生成随机椒盐噪点
def MakeNoise(input_image,number_of_noise=1000):
    res_image = np.copy(input_image)
    row,column = np.shape(input_image)
#     还可以加对三通道的特判
    for i in range(number_of_noise):
        x = np.random.randint(0,row)
        y = np.random.randint(0,column)
        borw = np.random.randint(0,2)
        if borw == 0:
            res_image[x][y] = 0
        else:
            res_image[x][y] = 255
    return res_image

均值滤波代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def Smooth(input_image,kernal_size = 3):
    ans = input_image.copy()
    border = np.uint8(kernal_size/2)
    # 根据核大小扩充边界
    addBorder = cv2.copyMakeBorder(input_image,border,border,border,border,cv2.BORDER_REFLECT_101)
    filWin = np.ones((kernal_size,kernal_size),dtype=np.int64)
    row,column = np.shape(addBorder)
    for i in range(border,row-border):
        for j in range(border,column-border):
            # 以图像i,j为中心，大小为2border+1的区域
            temp = addBorder[i-border:i+border+1,j-border:j+border+1]
            temp_sum = np.sum(temp*filWin)
            ans[i-border][j-border] = temp_sum/(kernal_size**2)
    return ans

Laplace 锐化滤波代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def Laplace(input_image,lap_template):
    x,y = np.shape(lap_template)
    if(x!=y):
        print("Shape of Laplace template wrong!")
        pass

    print("The Laplace template is: ")
    print(lap_template)
    border = x//2
    ans = np.zeros(np.shape(input_image), dtype = np.int64)

    # 扩充边界
    addBorder = cv2.copyMakeBorder(input_image,border,border,border,border,cv2.BORDER_REFLECT_101)

    row,column = np.shape(addBorder)
    for i in range(border,row-border):
        for j in range(border,column-border):
            ans[i-border][j-border] = np.sum(lap_template*addBorder[i-border:i+border+1 , j-border: j+border+1])
    return ans

# Laplace 常用算子 [[-1,-1,-1],[-1,8,-1],[-1,-1,-1]]
laplace_template_4 = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]])
laplace_template_8 = np.array([[-1,-1,-1],[-1,8,-1],[-1,-1,-1]])


lap_img_4 = Laplace(im,laplace_template_4)

频率域图像处理

频率域分析基础

这次课，我们首次接触了“频率域分析”方法。我们从高数中学习过的“无穷级数”概念切入，讲到基于三角函数的傅里叶级数。其中的关键思路是：任意周期函数都可以用傅里叶级数展开，而傅里叶级数的基函数是频率各异的三角函数；也就是说任意周期函数都可以用一组频率各异的三角函数的线性组合来表示。而线性组合系数，即基函数的加权系数，就可以解读为原函数在基函数上的“投影系数”（或坐标，或两者的相似度）。

如果能够求出傅里叶级数中的加权系数，那么就可以实现对原函数的“频域分解”，即原函数主要可以由哪些频率的三角函数叠加而成，也可以解读为原函数主要包含哪些频率的分量。而“频率”对于图像而言是有直观意义的。一般而言，低频成分对应图像中的总体形状和总体色彩和明暗分布；而高频成分对应图像中的细节和边缘信息。如果我们把图像看作是二维函数，并能够将其进行频域分解，那么增强其低频分量就是相当于平滑操作，增强其高频分量就相当于锐化操作。这就是我们除了前几章学习的“空间域分析”方法之外，又多了一个全新的视角和全新的方法。

傅里叶变换就是求解傅里叶级数中的加权系数的方法。其思路很简单，由于加权系数对应原函数与基函数的相似程度，我们可以用两者的离散采样点组成的向量的点积来估计二者的相似程度。因为向量的点积与向量的夹角的余弦成正比，夹角越小代表两个向量越相似，其余弦值就越大。工程上有一种名为“快速傅里叶变换”（FFT）的数值运算方法，可以快速求解傅里叶变换。但对于图像而言，它要求图像的长宽为 2 的整数次幂。

傅里叶逆变换是运用傅里叶级数来重构原函数的方法。其表达式就是傅里叶级数的表达式。有了基函数和加权系数，进行傅里叶逆变换是很直接的。

在工程上，傅里叶变换（FT）一般是对离散数据进行处理的，故名为离散傅里叶变换(DFT)。又因为常用快速傅里叶变换算法来求解，所以常用快速傅里叶变换(FFT)表示，其逆变换为 iFFT。无论是 FFT 还是 iFFT 都有现成的高度优化过的软件包使用，一般不用你自己实现。在 Python 语言中的 Numpy 包中就有相关实现，在 C++语言中也有 OpenCV 包中的对应实现。

我们在讲解了上述原理后，就尝试对一张图像（黑色背景的中央有个白色长方形）进行了傅里叶变换，对变换系数矩阵组成的图像进行了对数化（压缩值的分布，以便看清楚微弱的分布细节），观察了图像的频率域分解效果。然后我们分别演示了基于频率域分析的图像平滑（将 FFT 系数矩阵的边缘区域置零，再逆变换回空间域）和图像锐化（将 FFT 系数矩阵的中央区域置零，再逆变换回空间域，得到边缘图，再与原图进行叠加以便锐化图像）。这部分的演示，请看附件中的代码。

使用 OpenCV 实现 DFT

OpenCV 也提供的工具用于实现离散傅里叶变换，分别是cv2.dft()和cv2.idft()函数。返回有两个通道，第一个通道是结果的实部，第二个通道是结果的虚部。所以在计算其幅度的时候需要先手动转化一次。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
## Magnitude spectrum


## 使用Numpy实现DFT
## 进行二维DFT变换
f_img1 = np.fft.fft2(img1)
## 对换象限
fshift_img1 = np.fft.fftshift(f_img1)
## 幅度谱 对数化？
manitude_img1 = 20*np.log(np.abs(fshift_img1))
## img2图像也做相同的处理
f_img2 = np.fft.fft2(img2)
fshift_img2 = np.fft.fftshift(f_img2)
manitude_img2 = 20*np.log(np.abs(fshift_img2))

## 使用OpenCV实现DFT
f_img3 = cv2.dft(np.float32(img3), flags = cv2.DFT_COMPLEX_OUTPUT)
fshift_img3 = np.fft.fftshift(f_img3)
# manitude_img3 = 20*np.log(np.abs(fshift_img3))
# 手动转化计算magnitude
manitude_img3 = 20*np.log(cv2.magnitude(fshift_img3[:,:,0],fshift_img3[:,:,1])+1e-15)


## 使用OpenCV实现DFT
f_img4 = cv2.dft(np.float32(img4), flags = cv2.DFT_COMPLEX_OUTPUT)
fshift_img4 = np.fft.fftshift(f_img4)
manitude_img4 = 20*np.log(cv2.magnitude(fshift_img4[:,:,0],fshift_img4[:,:,1]))

理想滤波器

所谓"理想"是指无法通过硬件实现的硬截断

理想低通滤波器 ILPF

在圆外“阻断”所有频率，而在圆内无衰减的通过所有频率，这种二维低通滤波器称为理想低通滤波器（ILPF），由下面的函数确定

\[ H_{ILPF}(u,v) = \left \{ \begin{aligned} 1, & D(u,v) \le D_0 \\ 0, & D(u,b) > D_0 \end{aligned} \right. \]

其中\(D_0\)是一个正常数，\(D(u,v)\)表示频率域中的点\((u,v)\)距离频率域中心\((\frac{P}{2},\frac{Q}{2})\)的距离。

理想高通滤波器 IHPF

与低通类似，高通是将阈值的圆内“阻断”所有频率，而在圆外无衰减的通过所有频率，描述如下

\[ H_{IHPF}(u,v) = \left \{ \begin{aligned} 0, & D(u,v) \le D_0 \\ 1, & D(u,b) > D_0 \end{aligned} \right. \]

代码实现

代码实现低通滤波器并展示其频率域透视图、频率域图像显示、空间域图像显示和径向剖面图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
# 获得指定大小的理想滤波器
def getIdealMask(mask_shape, filter_d0,hl_type):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "lpf":
                if dis <= filter_d0:
                    mask[i,j] = 1
                else:
                    mask[i,j] = 0
            elif hl_type == "hpf":
                if dis <= filter_d0:
                    mask[i,j] = 0
                else:
                    mask[i,j] = 1
    return mask

# 测试ILPF
# 参数设置
mask_shape = (100,100)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getIdealMask(mask_shape,d,filter_type)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "ILPF Perspective Axes3D", cmap = "gray")
# 不想进行列表解析,需要调用frompyfunc构建np可以用的分段函数
ufunc1 = np.frompyfunc(lambda x: 0 if (x-d)>0 else 1, 1, 1)
drawCurv(ax2,[ufunc1],["ILPF"],d,title = "ILPF Curv")
drawPanel(ax3,myfilter,title = "ILPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "ILPF Spatial Panel Axes2D")
plt.show()

# 测试IHPF
d = 20
filter_type = "hpf"
myfilter = getIdealMask(mask_shape,d,filter_type)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "IHPF Perspective Axes3D", cmap = "gray")
# 不想进行列表解析,需要调用frompyfunc构建np可以用的分段函数
ufunc1 = np.frompyfunc(lambda x: 1 if (x-d)>0 else 0, 1, 1)
drawCurv(ax2,[ufunc1],["IHPF"],d,title = "Curv")
drawPanel(ax3,myfilter,title = "IHPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "IHPF Spatial Panel Axes2D")
plt.show()

理想低通 ILPF 结果四图

png
png

理想高通 IHPF 结果四图

png
png

布特沃斯滤波器

可通过硬件实现，可以通过阶数进行控制，一些资料中又称之为“巴特沃斯滤波器”。

布特沃斯低通滤波器 BLPF

截止频率位于距原点\(D_0\)的\(n\)阶布特沃斯滤波器（BLPF）的传递函数定义为：

\[H_{BLPF}(u,v) = \dfrac{1}{ 1 + {[ \dfrac {D(u,v)}{D_0} ]}^{2n} }\]

布特沃斯高通滤波器 BHPF

对应的传递函数定义为：

\[H_{BHPF}(u,v) = \dfrac{1}{ 1 + {[ \dfrac {D_0}{D_(u,v)} ]}^{2n} }\]

（分母分子颠倒）

两式中\(n\)对应了即阶参数，下面的代码给出巴特沃斯滤波器的实现，其频率域透视图、频率域图像显示、空间域图像显示和径向剖面图，曲线图绘制出不同阶下的取值。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
# 获得指定大小的布特沃斯滤波器
def getButterworthMask(mask_shape,filter_d0,hl_type,butter_n = 1):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "lpf":
                mask[i,j] = 1.0/(1+(dis/filter_d0)**(2*butter_n))
            elif hl_type == "hpf":
                # 除以0情况特判一下
                if np.abs(dis)<eps:
                    mask[i,j] = 0
                else:
                    mask[i,j] = 1.0/(1+(filter_d0/dis)**(2*butter_n))

    return mask

# 测试BLPF
# 参数设置
mask_shape = (100,100)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getButterworthMask(mask_shape,d,filter_type, butter_n=2)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "BLPF(n=2) Perspective Axes3D", cmap = "gray")
funcs = []
labels = []
for i in range(1,6):
    labels.append("BLPF "+"n="+str(i))
funcs.append(lambda x:1.0/(1+(x/d)**(2*1)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*2)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*3)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*4)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*5)))

drawCurv(ax2,funcs,labels,d,title = "BLPF(n=1,2,3,4,5) Curv")
drawPanel(ax3,myfilter,title = "BLPF(n=2) Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "BLPF(n=2) Spatial Panel Axes2D")
plt.show()

# 测试BHPF
d = 20
filter_type = "hpf"
myfilter = getButterworthMask(mask_shape,d,filter_type, butter_n=1)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "BHPF(n=2) Perspective Axes3D", cmap = "gray")
funcs = []
labels = []
for i in range(1,6):
#     funcs.append(lambda x:1.0/(1+(d/x)**(2*i)))
#     ufunc = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*i)), 1, 1)
#     funcs.append(ufunc)
    labels.append("BHPF "+"n="+str(i))
ufunc1 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*1)), 1, 1)
ufunc2 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*2)), 1, 1)
ufunc3 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*3)), 1, 1)
ufunc4 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*4)), 1, 1)
ufunc5 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*5)), 1, 1)
funcs.append(ufunc1)
funcs.append(ufunc2)
funcs.append(ufunc3)
funcs.append(ufunc4)
funcs.append(ufunc5)


drawCurv(ax2,funcs,labels,d,title = "BHPF(n=1,2,3,4,5) Curv")
drawPanel(ax3,myfilter,title = "BHPF(n=2) Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "BHPF(n=2) Spatial Panel Axes2D")
plt.show()

布特沃斯低通 BLPF 结果四图

png
png

布特沃斯高通 BHPF 结果四图

png
png

高斯低通滤波器 GLPF

高斯低通滤波器二维形式由下式给处：

\[H_{GLPF}(u,v) = e^{\dfrac{-D^2(u,v)}{2 \sigma^2}}\]

\(\sigma\)描述了中心的扩散速度，和其他滤波器描述式统一，通过令\(\sigma = D_0\)，可以用表示其他滤波器的方法表示高斯滤波器。

\[H_{GLPF}(u,v) = e^{\dfrac{-D^2(u,v)}{2 D_0^2}}\]

高斯高通滤波器 GHPF

如下：

\[H_{GHPF}(u,v) =1 - e^{\dfrac{-D^2(u,v)}{2 D_0^2}}\]

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# 获得指定大小的高斯滤波器
def getGaussianMask(mask_shape,filter_d0,hl_type):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "hpf":
                mask[i,j] = 1-np.exp(-(dis**2) / (2*(filter_d0**2)))
            elif hl_type == "lpf":
                mask[i,j] = np.exp(-(dis**2)/(2*(filter_d0**2)))
    return mask

# 测试GLPF
# 参数设置
mask_shape = (200,200)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getGaussianMask(mask_shape,d,filter_type)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "GLPF Perspective Axes3D", cmap = "gray")
drawCurv(ax2,[lambda x:np.exp(-(x**2)/(2*(d**2)))],["GLPF"],d,title = "GLPF Curv")
drawPanel(ax3,myfilter,title = "GLPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "GLPF Spatial Panel Axes2D")
plt.show()

# 测试GHPF
d = 20
filter_type = "hpf"
myfilter = getGaussianMask(mask_shape,d,filter_type)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "GLPF Perspective Axes3D", cmap = "gray")
drawCurv(ax2,[lambda x:1-np.exp(-(x**2)/(2*(d**2)))],["GHPF"],d,title = "GLPF Curv")
drawPanel(ax3,myfilter,title = "GLPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "GLPF Spatial Panel Axes2D")
plt.show()

高斯低通 GLPF 结果四图

png
png

高斯高通 GHPF 结果四图

png
png

滤波器对比总结

在规定滤波器为 100x100，阈值为 20 时，可以明显观察到，理想滤波器->高阶布特沃斯滤波器->低阶布特沃斯滤波器->高斯滤波器，可以由函数 Curv 看出对应的过渡。

我们同时也发现理想滤波器确实会存在振铃特性，这个将在后面的文章中再做分析学习。

振铃特性

理想低通滤波器实验中，我们也可以很明显的发现振铃现象，间隙处原本统一的纹理由于模糊变得有明暗起伏。而随着被滤去的高频内容的数量的减少，图像的纹理变得越来越好，甚至我们仔细看第三幅图，也能发现振铃现象的纹理，课本是这么评价振铃现象和 ILPF 的。

这种振铃现象是理想滤波器的一种特性，从这个例子我们可以清楚地看到，理想低通滤波器并不是非常实用。然而，作为滤波概念发展的一部分，研究这种滤波器的特性非常有用。

振铃现象的一些见解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 绘制空间域表示图和水平线灰度剖面图

for d in d_list:
    fre_mask=getIdealMask((688,688),d,"lpf")
    spa_mask=frequencyToSpatial(fre_mask)
    X = [i for i in range(spa_mask.shape[0])]
    Y = spa_mask[spa_mask.shape[0]//2]
    plt.figure(figsize=(8,4))
    ax1 = plt.subplot(121)
    ax2 = plt.subplot(122)
    ax1.set_title(f"Spatial Img with $D_0 = {d}$")
    ax2.set_title(f"Grayscale with $D_0 = {d}$")
    ax1.imshow(spa_mask,cmap = "gray")
    ax2.plot(X,Y)
    ax2.spines['left'].set_color('none')
    ax2.spines['top'].set_color('none')
    ax2.spines['right'].set_color('none')
    ax2.set_yticks([])
    ax2.set_yticklabels([])
png
png
png
png
png
png
png
png
png
png

分别对应了上面的阈值，观察其滤波器的空间与表示，在\(D_0\)较小的时候有很明显的波动函数形状。

ILPF 的模糊和振铃特性可用卷积定理来解释。由于 ILPF 在频率域的剖面图类似于盒状滤波器，因此可以预料相应空间滤波器具有 sinc 函数形状。空间域滤波可通过\(h(x,y)\)于图像卷积来实现。讲图像中的每个像素想象为一个离散冲击，它的强度与所在位置的灰度成正比。一个 sinc 函数与一个冲激卷积就是在冲激处复制这个 sinc 函数。sinc 函数的中心波瓣是引起模糊的主因，而外侧较小的波瓣是造成振铃的主要原因。sinc 函数“展开度”与\(H(u,v)\)半径成反比，所以\(D_0\)越大，空间 sinc 函数就趋近于一个卷积时不会导致模糊但也不会产生振铃的冲激。

图像去模糊（图像复原）

引言

前面所接触的图像增强是一个，已知原图像，经处理后，得到增强图像的一个过程，而这一章将要深入的图像复原（重建），则是希望从（被污染）过的图像，经处理后，得到原图像的过程，是以预先确定的目标来改善图像。图像复原试图利用退化现象的某种先验知识来复原被退化的图像。因而，复原技术时面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。虽然图像增强和图像复原两者在覆盖的领域和使用的技术栈有所重叠，其中还是有几点区别我们还是要提起注意的。

  • 形象化的描述

    图像增强主要是一个主观的过程，而图像复原大部分是一个客观的过程。

  • 已知与未知的区别

    图像增强已知原始图像与变换（卷积核或者其频率域的谱），对于增强效果是未知且非预先确定的（只有一个大概方向，比如模糊还是锐化），而图像复原则已知污染图像，且对原始图像是预先确定的（测试条件下甚至是有标准比对的原始图像，是已知的），对于复原变换（污染变换的逆）常常是未知的，这就要求我们在做图像处理时常常需要“估计”我们的复原变换。

  • 期望与探索的区别 由于两者已知和未知上的差距，这就导致图像复原通常会涉及设立一个最佳准则来产生期望结果的最佳估计。相比之下，图像增强技术基本上是一个探索性过程，即根据人类视觉系统的生理特点来设计改善图像的方法。

背景知识

线性系统

可加性，\(x_1(t)+x_2(t) = y_1(t)+y_2(t)\)，从而有\(a \times x_1(t) = a \times y_1(t)\)

平移不变性

卷积

卷积、离散二维卷积

从而可以利用卷积这个工具。

冲激响应：输入为一个脉冲信号，输出是一个冲激响应\(h(x)\)，实际上可以就是之前接触过的卷积核。

图形复原

复原

  • 试图利用退化过程的先验知识使已退化的图像恢复本来面目
  • 根据退化的原因，分析引起退化的环境因素
  • 建立相应的数学模型
  • 沿着使图像降质的逆过程回复图像

原因

  • 由于对焦不准导致的图像模糊是很常见的
  • 此时，物体的外观信息并未丢失，而是被分散叠加显示了
  • 如果能够把分散叠加的信息彼此分离，就可能去除模糊

目的

  • 在于消除或减轻退化的影响

方法

  • 由于退化系统是黑盒的，盲复原往往很困难，噪声干扰也为复原过程带来了困难和不确定性
  • 图像复原是寻求在一定优化准则下的原始图象的最有估计。因此，不同的优化准则会获得不同的图像复原。评价指标的选择目前也是研究的方向之一，如峰值信号比等。

图像退化/复原过程模型

退化过程（污染过程）的描述：建模为一个退化函数和一个加性噪声项，对于输入图像（原图像）\(f(x,y)\)进行处理，产生一副退化后的图像\(g(x,y)\)，图像复原目的就是已知\(g(x,y)\)的前提下，希望得到原图像的一个估计，这个估计越接近原始输入图像越好。空间域中的退化图像可由下式给出：

\[g(x,y) = h(x,y) \star f(x,y) + \eta(x,y)\]

上式中\(h(x,y)\)是退化函数的空间表示，由第四章内容，我们可以将上式的模型写成等价的频率域表示：

\[g(u,v) = H(u,v) F(u,v) + N(u,v)\]

这两个式子是本章后面大部分复原内容的基础。

采用频率域滤波图像复原的原因

在频率域滤波进行图像复原主要在两个方面效果较好，其一是利用频率域滤波消除周期噪声，另一个是利用频率域做退化函数的逆滤波。

为什么要在频率域做逆滤波？观察之前退化模型的两个式子，我们不难发现：

\[g(x,y) = h(x,y) \star f(x,y) + \eta(x,y) \tag 1\]

\[g(u,v) = H(u,v) F(u,v) + N(u,v) \tag 2\]

对(1)式空间域来说，想要从\(g(x,y)\)恢复\(f(x,y)\)，避不开的是“卷积”的逆运算，这在定义和实现的复杂上都比较困难，而转化到频率域，从(2)式我们或许可以通过一个“除法”来实现逆滤波，结合之前噪声模型相关内容，我们尝试在频率域上对仅退化函数影响的图像，和更复杂一些的，退化函数和加性噪声双重影响的图像进行复原。

逆滤波

退化函数已给出，或者由上面退化函数的估计方法获得后，最简单的复原方法是直接做逆滤波，即：

\[\hat{F}(u,v) = \dfrac{G(u,v)}{H(u,v)}\]

然而根据前述我们知道，在噪声的影响下，\(\hat{F}(u,v)\)和\(F(u,v)\)仍有差别，即

\[\hat{F}(u,v) = F(u,v) \dfrac{N(u,v)}{H(u,v)}\]

这个式子两点启发：

  1. 知道退化函数也不能完全复原未退化图像，因为噪声函数未知。
  2. 如果退化函数是零或是非常小的值，那么噪声影响会被放大

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

实际上维纳滤波是在这里是相对逆滤波来说的，而并非指特别的滤波函数，且不仅应用在运动模糊滤波中。

\[\hat{F}(u,v) =[\dfrac{1}{H(u,v)} \dfrac{ {|H(u,v)|}^2 }{ {|H(u,v)|}^2+ S_{\eta}(u,v)/S_f(u,v) }] G(u,v)\]

\(S_{\eta}(u,v)\)为噪声的功率谱而\(S_f(u,v)\)是未退化图像的功率谱，比值为噪信比。而由于谱\({|N(u,v)|}^2\)是一个常数，这大大简化了处理。我们常用下面的表达式来近似。

\[\hat{F}(u,v) =[\dfrac{1}{H(u,v)} \dfrac{ {|H(u,v)|}^2 }{ {|H(u,v)|}^2+ K }] G(u,v)\]

理解与启发：

  • 当频域某处没有噪声时， 𝑁/𝑆 的值为 0，方括号中表达式的值为 1，表示此时不需要对 𝐺/𝐻 修正。（回到逆滤波）
  • 当频域某处噪声大时， 𝑁/𝑆 的值比较大，方括号表达式的值小于 1，相当于对此处 𝐺/𝐻 的值进行抑制，以降低噪声的影响。

实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def wienerFiltering(input_img, h, NSR ,htype = "frequency"):
    assert htype in ("frequency","spatial")
    # 输入图像的傅里叶变换
    input_img_fft = np.fft.fft2(input_img)
    input_img_fft = np.fft.ifftshift(input_img_fft)

    if(htype == "spatial"):
        # 模糊化函数的傅里叶变换
        h_fft = np.fft.fft2(h)
    else :
        h_fft = h
    # 退化函数模值的平方
    h_abs_square = np.abs(h_fft)**2
    # 维纳滤波
    # 使用的是 共轭/模方+NSR的形式，并非 9 10 式
    output_image_fft = np.conj(h_fft) / (h_abs_square + NSR)

    # 输出图像傅里叶反变换
    output_image = np.fft.ifft2(output_image_fft * input_img_fft)
    output_image = np.abs(np.fft.fftshift(output_image))
    return output_image

多重曝光融合

  • 曝光程度对照片质量有影响
  • 不同曝光程度的照片中都包含有用场景信息

成像质量计算

19
19

为什么不直接用直方图均衡化来解决过曝光和欠曝光问题

过曝光和曝光不足是完全丧失了信息，而多重曝光得到的多幅图像则保留了这部分信息。举曝光不足的例子。比如黑色背景和灰色城堡，在曝光不足的情况下灰度值均为 0，显示为相同的黑色，实际上已经丢失了城堡和背景区分的信息，即便使用直方图均衡，可预知的结果是城堡和背景灰度值增加但是仍保持相同，无法区分两者。所以直方图均衡无法解决此类信息彻底丢失的问题。多重曝光则有多幅原始图像，可以对局部选择信息保留的部分进行融合。过曝光同理

计算问题

直接按下式进行加权融合会出现问题

\(𝒑=𝒘^′∗𝒑^′+𝒘^{′′}∗𝒑^{′′}+𝒘^{′′′}∗𝒑^{′′′}\)

解决办法

  • 基于原图边缘信息的图像融合

    • 用照片的边缘图把质量系数矩阵中的边缘过滤掉
    • Laplace 边缘图*成像系数矩阵

  • 基于多尺度边缘信息的图像融合

    • 单个 Laplacian 边缘图只包含图像中某个尺度的边缘信息
    • 所以将图像缩放，计算多尺度的边缘图，进行融合

图像压缩

三种冗余

P335，编码冗余、空间和时间冗余、不相关信息

一些基本的压缩方法

哈夫曼编码、算术编码、LZW 编码（课本 P348-P350）

LZW 编码步骤

  • 在词典中搜索与输入字符串的当前位置形成最大匹配的词汇 w；
  • 输出 w 的索引值；
  • 把 wa 输入词典，其中 a 是最大匹配后面的字符
  • 可用字典树加速查询最大匹配

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
func (c *Encoder) Circle() {
	fmt.Println("Start Encode Circle!")
	for key, char := range c.befComp {
		outIndex, addflag, _ := c.Digest(char)
		if addflag == true {
			c.preString = []byte{}
			c.preString = append(c.preString, char)
			c.aftComp = append(c.aftComp, outIndex)
			if key == len(c.befComp)-1 {
				c.aftComp = append(c.aftComp, char)
			}
		} else {
			// 在末尾的时候特判
			if key == len(c.befComp)-1 {
				// c.preString = []byte{}
				c.aftComp = append(c.aftComp, outIndex)
			}
		}
	}
	// fmt.Println("Test print:")
	// fmt.Println("before compression: ", c.befComp)
	// fmt.Println("after compression: ", c.aftComp)
}

func (c *Encoder) Digest(scanChar byte) (outIndex byte, flag bool, err error) {
	// 将Pre字符串和扫描的字符组合成为checkString
	checkString := append(c.preString, scanChar)
	// 这列可以换成 lastroot 有待改进
	nowPtr, isFind := c.FindTreeptr(checkString, c.root)
	if isFind == true {
		// 变更pre continue
		// 不做输出
		c.preString = checkString
		return byte(nowPtr.index), false, nil
	} else {
		// 加入新结点
		outIndex, err := c.AddTreeptr(nowPtr, scanChar)
		// 先输出
		return outIndex, true, err
	}

}

func (c *Encoder) FindTreeptr(leftString []byte, treeRoot *TreeNode) (addptr *TreeNode, flag bool) {
	tempChar := leftString[0]
	if treeRoot.childNode[tempChar] != nil {
		// 已经有节点
		if len(leftString) == 1 {
			// 返回本结点
			return treeRoot.childNode[tempChar], true
		} else {
			// 递归
			return c.FindTreeptr(leftString[1:], treeRoot.childNode[tempChar])
		}
	} else {
		// 未有节点
		// 返回父亲结点（输出过程需要知道index）
		return treeRoot, false
	}
}

func (c *Encoder) AddTreeptr(faptr *TreeNode, addchar byte) (outIndex byte, err error) {
	outIndex = byte(faptr.index)
	if c.capacity >= 256 {
		return outIndex, errors.New("Dictionary Full!")
	}
	// 由于获取父结点index需要所以这么操作
	faptr.childNode[addchar] = new(TreeNode)
	faptr.childNode[addchar].index = int(c.capacity)
	c.capacity++
	// c.dict[newIndex] =
	return outIndex, nil
}

LZW 解码步骤

  • 解码步骤
  • 初始化词典；
  • 翻译首个索引值为 w，输出 w；
  • 把 w？放入词典，其中？为待定字符；

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
func (c *Decoder) Circle() {
	fmt.Println("Start Decoder Circle!")
	for _, value := range c.befComp {
		if len(c.storeString) == 0 {
			c.aftComp = append(c.aftComp, value)
			// 第一个 特判 变更storeString
			c.storeString = append(c.storeString, value)
			continue
		}
		// 其他情况
		// (错误)写： aftComp是byte切片   storeString 也是byte切片
		// 先从字典获取现在值对应的字符串
		lookUpString := c.dict[value]
		// 加入输出
		c.aftComp = append(c.aftComp, []byte(lookUpString)...)

		// XX? 变成XXX
		// 把查到的字典加进来
		if flag := c.CheckDicFull(); flag == true {
			continue
		} else {

			c.StoreDic(c.storeString, lookUpString)
			// 将storeString 转化为当前 XX？
			c.storeString = []byte(lookUpString)
		}

	}
}

func (c *Decoder) CheckDicFull() (flag bool) {
	if len(c.dict) == 256 {
		return true
	}
	return false
}

func (c *Decoder) StoreDic(original []byte, prefix string) {
	// 原始+这次扫描的第一个byte（）
	newcontent := append(original, prefix[0])
	length := len(c.dict)
	fmt.Println(length, string(newcontent))
	c.dict[byte(length)] = string(newcontent)
}

阶段总结

  • 行程编码是利用数据的稀疏分布特点进行编码
  • 词典编码是根据词典将输入字符转换为对应词条的索引输出
  • LZW 方法在编码时根据输入的字符串动态生成词典
  • LZW 方法在解码时根据编码值和待定字符法进行词典的复原

LZ77 编码图解

20
20
  • 搜索缓存区远大于编码缓存区(32KB v.s. 258 Bytes)
  • 成块地读取字符以减轻 I/O 负担
  • 编码完成后采用哈夫曼法对三元组进行再编码
  • 编码质量一般优于 LZW 法
    • LZW 法的词典中存在没有使用的词汇
    • 在 LZW 法中已经扫描过的字符不能重组为词汇
  • 应用非常广泛(zip, rar, gzip, png, arj)

阶段总结

  • LZ77 是隐式词典法
  • LZ77 存储隐式词典增大后，编码效率降低的问题
  • 窗口 LZ77 编码法同时限制隐式词典的大小和搜索缓冲区的大小，进而保证编码效率
  • LZ77 法的编码值是三元组，具有不随着词典增大而增大的优点
  • 窗口 LZ77 编码的三元组定义具有特殊性

其他课件内容

DCT 相关

21
21
22
22
23
23
24
24
25
25
26
26
27
27
28
28

常见压缩思路

  • 二值图像压缩思路：行程编码
  • 灰度值图像压缩思路：
    • 设置阈值转换为 n 幅二值图像，分别用游程编码压缩
    • 用预测编码法，降低图像的信息熵，如简单的 yt+1=yt
    • 变换编码法（如离散余弦变换 DCT）
  • 真彩色图像的压缩思路：RGB 分为三幅灰度图像转灰度图像压缩。

常见图片压缩格式实现思路

29
29
30
30
31
31
32
32
33
33
34
34
35
35

图像形态学

概述

  • 形态学一般指生物学中研究动物和植物结构的一个分支
  • 用数学形态学（也称图像代数）表示以形态为基础对图像进行分析的数学工具
  • 基本思想是用具有一定形态的结构元素去度量和提取图像中的对应形状以达到对图像分析和识别的目的
  • 形态学图像处理的数学基础和所用语言是集合论
  • 形态学图像处理的应用可以简化图像数据，保持它们基本的形状特性，并除去不相干的结构
  • 形态学图像处理的主要运算有 4 个：膨胀、腐蚀、开操作和闭操作

四种运算

36
36
37
37
38
38
39
39
40
40
  • 开操作的 3 条性质
    • AoB 是 A 的子集合
    • 如果 C 是 D 的子集，则 CoB 是 DoB 的子集
    • (AoB)oB= AoB
  • 闭操作的 3 条性质
    • A 是 A•B 的子集合
    • 如果 C 是 D 的子集，则 C•B 是 D•B 的子集
    • (A•B)•B= A•B

形态学主要应用

边界提取：B 先对 A 腐蚀，A 减去腐蚀得到边界

区域填充：

联通分量提取

图像分割

概述

分割的目的：将图像划分为不同区域。

三大类方法：

  • 根据区域间灰度不连续搜寻区域之间的边界，在间断检测、边缘连接和边界检测介绍
  • 以像素性质的分布进行阈值处理，在阈值处理介绍
  • 直接搜寻区域进行分割，在基于区域的分割中介绍

Q&A：多尺度的边缘检测仍然需要选择确定的尺度，不能智能辨识。

间断检测类型

点检测、线检测、边缘检测

差分、数值微分：\(\partial\) 一阶对应线检测、二阶对应点检测

间断检测：拉普拉斯算子

\[∇^2f ＝4z_5 - (z_2 + z_4 + z_6 + z_8)\] \[∇^2f ＝8z_5 - (z_1 + z_2 + z_3 + z_4 + z_5 + z_6 + z_7 + z_8 + z_9)\]

  • 缺点：
    • 拉普拉斯算子对噪声具有敏感性
    • 拉普拉斯算子的幅值产生双边缘
    • 拉普拉斯算子不能检测边缘的方向
  • 优点：
    • 可以利用零交叉的性质进行边缘定位
    • 可以确定一个像素是在边缘暗的一边还是亮的一边

Derivative of Gaussian Filter

根据交换律，先进性高斯模糊核与边缘提取核先计算，从而减少整体的计算量。

课后探索：边缘提取：Canny 算法 Mask-RCNN

彩色图像

色彩的本质

人对色彩的感知：视杆细胞，视锥细胞

色彩空间（英语：Color space）是对色彩的组织方式。借助色彩空间和针对物理设备的测试，可以得到色彩的固定模拟和数字表示。色彩空间可以只通过任意挑选一些颜色来定义，比如像彩通系统就只是把一组特定的颜色作为样本，然后给每个颜色定义名字和代码；也可以是基于严谨的数学定义，比如 Adobe RGB、sRGB。

CIE 1931 xyz 是第一次人眼对色彩感知度量建立色彩空间的尝试，几乎所有其它色彩空间的基础。有一些变体。

  • CIELUV
  • CIE L*A*B* 基于人眼

RGB 加法混色法，从暗处叠加产生颜色（sRGB，Adobe RGB）

CMYK 减法混色法

HSV 艺术家们常用，使用色相饱和度和明度 HSL 用亮度（Lightness）代替明度（Value/Brightness）

其他讨论区问题

调研 C/C++语言中 main 函数的两个常见参数的意义和用法，将调研结果写在下面的讨论区

主函数程序应当含有一个名为 main 的全局函数，它被指定为程序的启动点。它应当有下列形式之一：

1
2
3
int main () { body } (1)
int main (int argc, char _argv[]) { body } (2)
// 其他由实现定义的形式，返回类型为 int (3)
  • argc - 非负数，表示从程序运行的环境传递给程序的实参个数。

  • argv - 指针，指向包含 argc + 1 个指针的数组的首元素。数组末元素为空指针，若其前面有任何元素，则它们指向空终止多字节字符串，表示从执行环境传递给程序的若干参数。若 argv[0] 不是空指针，或等价地 argc > 0 ，则它指向表示用于调用程序的名称的字符串，或空字符串。

  • body - 主函数的函数体 名字 argc 和 argv 是任取的，而且形参的类型的表示方式也是：int main(int ac, char** av) 同样合法。

main() 的一种常见的由实现定义的形式还有（除 argc 和 argv 之外的）第三个参数，类型为 char*[] ，指向指向执行环境变量的指针数组。