《数字图像处理》第三章—

《数字图像处理》第三章——图像增强

图像增强是突出图像中的某些部分，将图像变得易于人眼观察，通常用于去噪，平滑，锐化，阈值处理。

处理方式：

空间域处理
- 点处理
- 邻域处理
频率域处理

点像素处理法

阈值处理：最常用的是二分阈值，也就是像素值大于某个阈值全部变成255，小于该阈值全部变成0。
图像反转：常用于灰度图的黑白反转，使得想要显示的图像细节变成低像素值，其他变成高像素值，公式如下： $$ s = L - 1 -r $$
伽马校正：原理是将图像的像素点再有效范围内分布更加均匀。它通常用于常用显示器显示图像或者对整体过亮的图像进行延展。公式如下： $$ s = cr^{\gamma} $$ 注：这里通过调整上标伽马的参数来进行校正
像素拉伸：原理是在低像素和高像素区域的像素值都乘以一个小于1的系数，处于中间像素区域的像素值都乘以一个大于1的系数。

图像均衡

直方图均衡化

直方图的定义：横坐标为像素值，纵坐标为像素的个数。
直方图归一化：横坐标依然为像素值，纵坐标变成了该像素值出现的概率。

直方图的目的是通过一个函数来使原图像的直方图分布从不均匀到一个均匀分布，公式如下： $$ s = T(r) = (L - 1)\int_{0}^{r}P_{r}(w)dw $$

注：其中L为256（8bit的最大数加1），P为概率分布函数

该函数为一个递增函数，由旧的像素值（整数）变成新的像素值（浮点数），又一个取整的过程，所以直方图均衡化之后，直方图分布并不是完全平均

直方图匹配

直方图匹配是指给定一张图像，要将其直方图分布变换成特定直方图分布。

直方图匹配的实现：

1.先对原图像进行直方图均衡化处理，得到结果T
2.对目标直方图进行直方图均衡化处理，得到映射关系
3.对T进行步骤2的逆映射处理

邻域矩阵处理法

这也就是我们通常所说的滤波。通常使用一个n×n矩阵和图像中的所有n×n子块进行求内积的操作

均值滤波（局部）：使用全为1的矩阵进行求内积。注意矩阵乘法之前必须进行归一化，防止数据超过原有的范围。

$$\sum_{i}(a - a_{i})^2取最小时就是其均值（误差较大）$$

高斯滤波：使用高斯滤波进行局部平均，相当于离中间像素值越大的其权重分配的越小，差值越小的权重分配的越大。
中值滤波：对图像的n×n子块，排序取中间值来替代原像素的值。（通常用于去除椒盐噪声，效果很好）

$$\sum_{i}\lvert a - a_{i} \rvert取最小时就是其中值（误差鲁棒性更小）$$

非局部均值滤波：根据图像的自相似性，我们可以将相似的局部区域一起做一个平均。

1
2
3
4
5
6
7


-p0指所有相似区域p1～pn的相同的部分
-s1～sn指所有区域不同的部分

p1 = p0 + s1
p2 = p0 + s2
···
pn = p0 + sn

如何进行相似区域的寻找？

其实只要比较区域间的均方误差，保证其值很小就可以了！

锐化和反锐化掩膜

离散情况下，x方向和y方向的导数为： $$ \frac{\partial F}{\partial x} = f(x + 1) - f(x)$$

$$ \frac{\partial^2 F}{\partial x^2} = f(x + 1) + f(x - 1) - 2f(x)$$ 那么x和y方向的2阶导数之和为： $$ \frac{\partial^2 F}{\partial x^2} + \frac{\partial^2 F}{\partial y^2} = f(x + 1) + f(x - 1) - 2f(x) + f(y + 1) + f(y - 1) - 2f(y)$$

这其实就是拉普拉斯算子的推导公式，其掩膜demo：

0 1 0
1 -4 1
0 1 0

锐化应用流程：

输入图像->拉普拉斯算子处理->修正处理后的图像->与输入图像叠加->效果为对原图像做了锐化处理

反锐化掩膜（目的是为了增强图像的轮廓）：

原图像-> 平滑处理（局部均值）-> 求图像差值（原图-平滑处理后的图，一般在边缘会有差别）->乘一个修正因子（加强边缘）-> 再与原图像进行想加（得到轮廓增强的效果）

边缘检测和彩色边缘检测

梯度: $$ \nabla f(x,y) = (\frac{\partial F}{\partial x}, \frac{\partial F}{\partial y}) $$ 梯度方向变化最快
变化值： $$ \lvert \nabla f(x,y) \rvert = \sqrt{(\frac{\partial F}{\partial x})^2 + (\frac{\partial F}{\partial y})^2} $$ 沿着梯度方向，有最大的变化值。沿着梯度法线的方向，有最小的变化值

如果是RGB图像，可以对其三个通道分别求其梯度方向，再加上这三个通道方向的矢量叠加，就是其变化最快的方向。

课后作业

task

实现

  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
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440


import cv2 as cv #OpenCV
import matplotlib.pyplot as plt #图表
import numpy as np #矩阵处理
import random



# 显示图像
def cv_show(name, img):
    cv.imshow(name, img)
    cv.waitKey(0)
    cv.destroyAllWindows()

# 随机噪声
def random_noise(img, n):
    '''
    添加随机噪点（实际上就是随机在图像上将像素点的灰度值变为255即白色）
    :param image: 需要加噪的图片
    :param noise_num: 添加的噪音点数目，一般是上千级别的
    :return: img_noise
    '''
    #
    # 参数image：，noise_num：
    img_noise = img
    noise_num = 3000 * n
    rows, cols, chn = img_noise.shape
    # 加噪声
    for i in range(noise_num):
        x = np.random.randint(0, rows)#随机生成指定范围的整数
        y = np.random.randint(0, cols)
        img_noise[x, y, :] = 255
    return img_noise

# 高斯噪声
def noise_Gauss(img, loc=0, scale=0.005):
    """
    添加高斯噪声
    :param img: 输入灰度图像
    :param loc: 高斯均值
    :param scale: 高斯方差
    :param scale**2:高斯标准差
    :return: img：加入高斯噪声后的图像
    """

    img = np.float32(img) / 255  # (0-255)->(0-1)
    Gauss_noise = np.random.normal(loc, scale ** 0.5, img.shape)
    img = img + Gauss_noise  # 加性噪声
    if img.min() < 0:
        low_clip = -1
    else:
        low_clip = 0
    img = np.clip(img, low_clip, 1.0)  # 约束结果在(0-1)
    img = np.uint8(img * 255)  # (0-1)->(0-255)

    return img

# 椒盐噪声
def pepper_and_salt(img,percentage):
    num=int(percentage*img.shape[0]*img.shape[1])#  椒盐噪声点数量
    random.randint(0, img.shape[0])
    img2=img.copy()
    for i in range(num):
        X=random.randint(0,img2.shape[0]-1)#从0到图像长度之间的一个随机整数,因为是闭区间所以-1
        Y=random.randint(0,img2.shape[1]-1)
        if random.randint(0,1) ==0: #黑白色概率55开
            img2[X,Y] = (255,255,255)#白色
        else:
            img2[X,Y] =(0,0,0)#黑色
    return img2

# 非局部均值的子函数
def arraycompare(array1, array2, height, width):
    resultarray = np.zeros((height, width))
    for row in range(0, height):
        for col in range(0, width):
            resultarray[row, col] = max(array1[row, col], array2[row, col])
    return resultarray

# 非局部均值算法的子函数
def integralImgSqDiff2(paddedimg_val, Ds_val, t1_val, t2_val):

    lengthrow = len(paddedimg_val[:, 0])
    lengthcol = len(paddedimg_val[0, :])
    Dist2 = (paddedimg_val[Ds_val:-Ds_val, Ds_val:-Ds_val] -
             paddedimg_val[Ds_val + t1_val:lengthrow - Ds_val + t1_val,
                           Ds_val + t2_val:lengthcol - Ds_val + t2_val]) ** 2
    Sd_val = Dist2.cumsum(0)
    Sd_val = Sd_val.cumsum(1)
    return Sd_val

# 非局部均值算法
def nl_meansfilter(imagearray, h_=10, ds0=2, ds1=5):
    height, width = imagearray[:, :, 0].shape[:2]
    length0 = height + 2 * ds1
    length1 = width + 2 * ds1
    h = (h_ ** 2)
    d = (2 * ds0 + 1) ** 2
    imagearray_NL = np.zeros(imagearray.shape).astype('uint8')
    for i in range(0, 3):
        paddedimg = np.pad(imagearray[:, :, i], ds0 + ds1 + 1, 'symmetric')
        paddedimg = paddedimg.astype('float64')
        paddedv = np.pad(imagearray[:, :, i], ds1, 'symmetric')
        paddedv = paddedv.astype('float64')
        average = np.zeros((height, width))
        sweight = np.zeros((height, width))
        wmax = np.zeros((height, width))
        for t1 in range(-ds1, ds1 + 1):
            for t2 in range(-ds1, ds1 + 1):
                if t1 == 0 and t2 == 0:
                    continue
                Sd = integralImgSqDiff2(paddedimg, ds1, t1, t2)
                SqDist2 = Sd[2 * ds0 + 1:-1, 2 * ds0 + 1:-1] + Sd[0:-2 * ds0 - 2, 0:-2 * ds0 - 2] - \
                          Sd[2 * ds0 + 1:-1, 0:-2 * ds0 - 2] - Sd[0:-2 * ds0 - 2, 2 * ds0 + 1:-1]
                SqDist2 /= d * h
                w = np.exp(-SqDist2)
                v = paddedv[ds1 + t1:length0 - ds1 + t1, ds1 + t2:length1 - ds1 + t2]
                average += w * v
                wmax = arraycompare(wmax, w, height, width)
                sweight += w
        average += wmax * imagearray[:, :, i]
        average /= wmax + sweight
        average_uint8 = average.astype('uint8')
        imagearray_NL[:, :, i] = average_uint8
    return imagearray_NL


# 播放黑白原视频
def Vedio_play1(vc):
# 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False

    # 用黑白播放读取的彩色视频
    while open:
        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            gray = cv.cvtColor(frame, cv.COLOR_BGR2GRAY)  # 将彩色图像转化成灰度图
            cv.imshow("result", gray)  #播放每一帧
            if cv.waitKey(40) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
                break
    vc.release()
    cv.destroyAllWindows()

# 每帧进行直方图均衡化并播放视频（黑白）
def Vedio_play2(vc):
# 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False

    # 用黑白播放读取的彩色视频
    while open:

        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            gray = cv.cvtColor(frame, cv.COLOR_BGR2GRAY)  # 将彩色图像转化成灰度图
            equ = cv.equalizeHist(gray)
            cv.imshow("result", equ)  #播放每一帧
            if cv.waitKey(40) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
                break
    vc.release()
    cv.destroyAllWindows()


# 每一帧拼接成图像，并对整个图像进行直方图均衡化，再拆分并播放视频（黑白）
def Vedio_play3(vc):
    # 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False
    index = 0
    x, y = frame[:, :, 0].shape
    # 处理成大图均衡化后播放读取的彩色视频
    while open:
        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            gray = cv.cvtColor(frame, cv.COLOR_BGR2GRAY)  # 将彩色图像转化成灰度图
            if(index == 0):
                res = gray.copy()
            else:
                res = np.hstack((res, gray.copy()))
            index += 1
    print("黑白视视频共有" + str(index) + "帧")
    equ = cv.equalizeHist(res)
    for i in range(index):
        temp = equ[0 : x , i * y : i * y + y]
        cv.imshow("result", temp)
        if cv.waitKey(30) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
            break
    vc.release()
    cv.destroyAllWindows()

# 播放彩色视频
def Vedio_play4(vc):
# 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False

    # 播放读取的彩色视频
    while open:
        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            cv.imshow("result", frame)  #播放每一帧
            if cv.waitKey(40) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
                break
    vc.release()
    cv.destroyAllWindows()



# 彩色视频每一帧进行均衡化
def Vedio_play5(vc):
# 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False

    # 均衡化处理后播放读取的彩色视频
    while open:
        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            b, g, r = cv.split(frame)
            equ1 = cv.equalizeHist(b)
            equ2 = cv.equalizeHist(g)
            equ3 = cv.equalizeHist(r)
            equ = cv.merge((equ1, equ2, equ3))
            cv.imshow("result", equ)  #播放每一帧
            if cv.waitKey(40) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
                break
    vc.release()
    cv.destroyAllWindows()


# 每一帧拼接成图像，并对整个图像进行直方图均衡化，再拆分并播放视频（黑白）
def Vedio_play6(vc):
    # 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False
    index = 0
    x, y = frame[:, :, 0].shape
    # 处理成大图均衡化后播放读取的彩色视频
    while open:
        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            if(index == 0):
                res = frame.copy()
            else:
                res = np.hstack((res, frame.copy()))
            index += 1
    print("彩色视视频共有" + str(index) + "帧")
    b, g, r = cv.split(res)
    equ1 = cv.equalizeHist(b)
    equ2 = cv.equalizeHist(g)
    equ3 = cv.equalizeHist(r)
    equ = cv.merge((equ1, equ2, equ3))
    for i in range(index):
        temp = equ[0 : x , i * y : i * y + y]
        cv.imshow("result", temp)
        if cv.waitKey(30) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
            break
    vc.release()
    cv.destroyAllWindows()


# 每帧进行非局部均值去噪声并播放视频（彩色）
def Vedio_play7(vc):
# 检查是否正确打开
    if vc.isOpened():
        # open代表是否读取成功，为bool值，frame为每一帧的图像的三维数组
        open, frame = vc.read()  # vc.read()代表读取视频的每一帧（从第一帧开始读取）
    else:
        open = False

    # 用黑白播放读取的彩色视频
    while open:

        ret, frame = vc.read()
        if frame is None:  # 为空则break
            break
        if ret == True:
            nl = nl_meansfilter(frame)
            cv.imshow("result", nl)  #播放每一帧
            if cv.waitKey(40) & 0xFF == 27:  # waitKey的参数代表等待的时间（数值越大，播放的速度越慢）  0xFF代表退出键为27
                break
    vc.release()
    cv.destroyAllWindows()


def main():

    # 第一题实现
    img1 = cv.imread("huajianghu.jpg", 0) #彩色图像变灰度图像
    plt.hist(img1.ravel(), 256) # 原图的直方图
    plt.show()

    equ = cv.equalizeHist(img1)
    plt.hist(equ.ravel(), 256) #均衡化后的直方图
    plt.show()

    res1 = np.hstack((img1, equ))
    cv_show("原图与均衡化的图像", res1)

    # 第二题实现
    img2 = cv.imread("woman_gray.jpg")
    # 分别使用高斯噪声和椒盐噪声
    img2_zaosheng1 = noise_Gauss(img2) # 方差为0， 标准差为0.005
    img2_zaosheng2 = pepper_and_salt(img2, 0.04) #百分之4的椒盐噪音

    #cv_show("1", img2_zaosheng1)
    #cv_show("2", img2_zaosheng2)

    # 滤波器大小为3
    median1_3 = cv.medianBlur(img2_zaosheng1, 3)
    median2_3 = cv.medianBlur(img2_zaosheng2, 3)

    #滤波器大小为5
    median1_5 = cv.medianBlur(img2_zaosheng1, 5)
    median2_5 = cv.medianBlur(img2_zaosheng2, 5)

    #展示图像
    res2_1 = np.hstack((img2, img2_zaosheng1, median1_3, median1_5))  # 横着拼接
    res2_2 = np.hstack((img2, img2_zaosheng2, median2_3, median2_5)) # 横着拼接
    cv_show("高斯噪声中值滤波", res2_1)
    cv_show("椒盐噪声中值滤波", res2_2)

    # 第三题实现
    img3 = cv.imread("woman_color.jpg")
    cv_show("lena'", img3)
    img3_zaosheng1 = noise_Gauss(img3) # 方差为0， 标准差为0.005
    nl_means1 = nl_meansfilter(img3_zaosheng1)
    res3_1 = np.hstack((img3, img3_zaosheng1, nl_means1))
    cv_show("1", res3_1)


    # 第四题实现
    img4 = cv.imread("woman_color.jpg")
    # 分别用5，10，20 次添加随机噪声（其实就是扩大添加随机噪声的规模）
    img4_5 = random_noise(img4.copy(), 5)
    img4_10 = random_noise(img4.copy(), 10)
    img4_20 = random_noise(img4.copy(), 20)
    res4 = np.hstack((img4, img4_5, img4_10, img4_20))
    cv_show("1", res4)


    # 第五题实现
    img5 = cv.imread("shaizi.png")
    b, g, r = cv.split(img5)

    # 采用OpenCV自带的Canny函数进行边缘检测 （50和100代表双阈值检测中的minVal和maxVal）
    '''
    该方法的流程为：
    1.使用高斯滤波器平滑图像 
    2.计算每个像素点的梯度大小和方向（离散情况下有误差）
    3.使用非极大值抑制的方法消除杂散影响 
    4.应用双阈值检测的方法确定真实和潜在的边缘
    5.抑制孤立的弱边缘完成检测
    '''
    b1 = cv.Canny(b, 50, 100)
    g1 = cv.Canny(g, 50, 100)
    r1 = cv.Canny(r, 50, 100)
    temp1 = cv.merge([b1, g1, r1])
    temp2 = cv.merge([b1, b1, b1]) #假设三个通道的值都为blue通道的值
    res5_1 = np.hstack((img5, temp1))
    res5_2 = np.hstack((img5, temp2))
    cv_show("彩色边缘检测", res5_1)
    cv_show("三个通道相等时的情况", res5_2) #可见边缘全都变成了白色

    # 第六题实现
    #1.黑白视频
    vc = cv.VideoCapture("test.mp4")  # 读取视频
    Vedio_play1(vc)  # 播放原视频
    vc = cv.VideoCapture("test.mp4")  # 读取视频
    Vedio_play2(vc) #逐帧均衡化并播放视频

    #2.彩色视频
    vc = cv.VideoCapture("test3.mp4") # 读取彩色视频
    Vedio_play4(vc) # 播放彩色视频
    vc = cv.VideoCapture("test3.mp4") # 读取彩色视频
    Vedio_play5(vc) # 彩色视频逐帧均衡化后播放

    # 第七题实现
    #1.黑白视频
    vc = cv.VideoCapture("test.mp4")
    Vedio_play1(vc)
    vc = cv.VideoCapture("test.mp4")
    Vedio_play3(vc)

    #2.彩色视频
    vc = cv.VideoCapture("test3.mp4")
    Vedio_play4(vc)
    vc = cv.VideoCapture("test3.mp4")
    Vedio_play6(vc)
    '''
    结论就是感觉先变成大图均衡化后视频更加流畅了
    '''

    # 第八题实现
    # 第一小问
    vc = cv.VideoCapture("test3.mp4")
    Vedio_play4(vc) # 播放彩色视频
    vc = cv.VideoCapture("test3.mp4")
    Vedio_play7(vc) # 逐帧进行非局部去噪声（能实现但实在是太慢了，因为非局部去燥算法本身就很慢，视频需要30秒才能播放一帧）
    '''
    将其变成大图像也能实现，只要在第六题的基础上修改就ok了
    但是没有演示效果（因为非局部均值算法慢成屎了）
    '''
    #第二小问
    '''
    实在不会了，看不懂
    '''

if __name__ == '__main__':
    main()