下面三个题目是我们老师上课的时候布置的课堂作业,要求手算,这里我用python实现一遍。
题目1:对于如图1所示的5x5的3比特图像,(1) 画出灰度直方图; (2)求直图均衡变换函数,给出计算过程,(3) 给出直方图均衡化以后的图像中每个像素的灰度级; (4) 画出均衡化后的直方图。
(1) 画出灰度直方图。不做解释,直接上代码:
'''
Description:
Author: Weijian Ma
Date: 2020-10-13 10:09:21
LastEditTime: 2020-10-13 10:19:32
LastEditors: Weijian Ma
'''
import numpy as np
import matplotlib.pyplot as plt
import cv2
plt.rcParams['font.sans-serif']=['SimHei'] # 用黑体显示中文
plt.rcParams['axes.unicode_minus']=False # 正常显示负号
img = np.array([[0, 1, 1, 0, 0],
[3, 4, 0, 2, 1],
[3, 2, 1, 0, 2],
[2, 1, 1, 0, 3],
[2, 3, 3, 4, 4]])
countNum = np.zeros(8)
for i in range(8):
countNum[i] = np.sum(img==i)
## 画出灰度直方图
fig = plt.figure(figsize=(10,4))
sub01 = plt.subplot(121)
sub02 = plt.subplot(122)
sub01.set_xticks(np.arange(8))
sub02.set_xticks(np.arange(8))
sub01.set_title("直方图")
sub02.set_title("归一化直方图")
sub01.set_xlabel("灰度")
sub01.set_ylabel("频数")
sub02.set_xlabel("灰度")
sub02.set_ylabel("频率")
sub01.bar(np.arange(8), countNum)
sub02.bar(np.arange(8), countNum/np.sum(countNum))
plt.show()
结果:
![](https://img-blog.csdnimg.cn/2020101310201589.png?watmark/2/text/aHR0cDovL2cuY3Nkbi5uZXQvV5bm1hbjIzMw==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
(2)求直图均衡变换函数,给出计算过程。
具体原理及公式可参考下面的链接,这里我直接写计算过程。
https://blog.csdn.net/weixin_45476502/article/details/108968328
对于原图像:
p
r
(
r
k
)
=
n
k
M
N
,
k
=
0
,
1
,
2
,
.
.
.
,
K
−
1
p_r(r_k)=\frac{n_k}{MN},k=0,1,2,...,K-1
pr(rk)=MNnk,k=0,1,2,...,K−1
r
0
=
0
,
r
1
=
1
,
r
2
=
2
,
r
3
=
3
,
r
4
=
4
,
r
5
=
5
,
r
6
=
6
,
r
7
=
7.
r_0=0,r_1=1,r_2=2,r_3=3,r_4=4,r_5=5,r_6=6,r_7=7.
r0=0,r1=1,r2=2,r3=3,r4=4,r5=5,r6=6,r7=7.
n 0 = 6 , n 1 = 6 , n 2 = 5 , n 3 = 5 , n 4 = 3 , n 5 = 0 , n 6 = 0 , n 7 = 0. n_0=6,n_1=6,n_2=5,n_3=5,n_4=3,n_5=0,n_6=0,n_7=0. n0=6,n1=6,n2=5,n3=5,n4=3,n5=0,n6=0,n7=0.
p r ( 0 ) = 6 25 , p r ( 1 ) = 6 25 , p r ( 2 ) = 1 5 , p r ( 3 ) = 1 5 , p r ( 4 ) = 3 25 , p r ( 5 ) = 0 , p r ( 6 ) = 0 , p r ( 7 ) = 0. p_r(0)=\frac{6}{25},p_r(1)=\frac{6}{25},p_r(2)=\frac{1}{5},p_r(3)=\frac{1}{5},p_r(4)=\frac{3}{25},p_r(5)=0,p_r(6)=0,p_r(7)=0. pr(0)=256,pr(1)=256,pr(2)=51,pr(3)=51,pr(4)=253,pr(5)=0,pr(6)=0,pr(7)=0.
s k = T ( r k ) = ( L − 1 ) ∑ j = 0 k p r ( r j ) = ( L − 1 ) M N ∑ j = 0 k n j , k = 0 , 1 , 2 , . . . , K − 1 s_k=T(r_k)=(L-1)\sum_{j=0}^kp_r(r_j)=\frac{(L-1)}{MN}\sum_{j=0}^kn_j,k=0,1,2,...,K-1 sk=T(rk)=(L−1)j=0∑kpr(rj)=MN(L−1)j=0∑knj,k=0,1,2,...,K−1
则,
∑ j = 0 0 n j = 6 , ∑ j = 0 1 n j = 12 , ∑ j = 0 2 n j = 17 , ∑ j = 0 3 n j = 22 , ∑ j = 0 4 n j = 25 , ∑ j = 0 5 n j = 25 , ∑ j = 0 6 n j = 25 , ∑ j = 0 7 n j = 25 \sum_{j=0}^0n_j=6,\sum_{j=0}^1n_j=12,\sum_{j=0}^2n_j=17,\sum_{j=0}^3n_j=22,\sum_{j=0}^4n_j=25,\sum_{j=0}^5n_j=25,\sum_{j=0}^6n_j=25,\sum_{j=0}^7n_j=25 ∑j=00nj=6,∑j=01nj=12,∑j=02nj=17,∑j=03nj=22,∑j=04nj=25,∑j=05nj=25,∑j=06nj=25,∑j=07nj=25
所以,
s 0 = 7 25 × 6 , s 1 = 7 25 × 12 , s 2 = 7 25 × 17 , s 3 = 7 25 × 22 , s 4 = 7 25 × 25 , s 5 = 7 25 × 25 , s 6 = 7 25 × 25 , s 7 = 7 25 × 25 s_0=\frac{7}{25}\times6,s_1=\frac{7}{25}\times12,s_2=\frac{7}{25}\times17,s_3=\frac{7}{25}\times22,\\s_4=\frac{7}{25}\times25,s_5=\frac{7}{25}\times25,s_6=\frac{7}{25}\times25,s_7=\frac{7}{25}\times25 s0=257×6,s1=257×12,s2=257×17,s3=257×22,s4=257×25,s5=257×25,s6=257×25,s7=257×25
则,直图均衡变换函数的图像如下:
![](https://img-blog.csdnimg.cn/20201013111454248.png?watmark/2/text/aHR0cDovL2cuY3Nkbi5uZXQvV5bm1hbjIzMw==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
代码:
'''
Description:
Author: Weijian Ma
Date: 2020-10-10 09:13:23
LastEditTime: 2020-10-13 11:12:43
LastEditors: Weijian Ma
'''
## 对图像进行直方图均衡
import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] # 用黑体显示中文
plt.rcParams['axes.unicode_minus']=False # 正常显示负号
## 定义一个函数,用于转换图像
def T(r_k, nkCum, height, width):
s_k = (7/(height*width))*nkCum[r_k]
return s_k
img = np.array([[0, 1, 1, 0, 0],
[3, 4, 0, 2, 1],
[3, 2, 1, 0, 2],
[2, 1, 1, 0, 3],
[2, 3, 3, 4, 4]])
height = img.shape[0]
width = img.shape[1]
## 计算每个灰度出现的次数
n_k = np.zeros(256,dtype=np.int32)
for i in range(height):
for j in range(width):
n_k[img[i][j]] = n_k[img[i][j]]+1
nkCum = np.cumsum(n_k)
plt.plot(np.arange(8),np.around(T(np.arange(8),nkCum,height,width)),'-o')
plt.xlabel('$r$')
plt.ylabel('$s$')
plt.grid()
plt.show()
(3) 给出直方图均衡化以后的图像中每个像素的灰度级.
应用上一小题中得出的直图均衡变换函数(函数图像):
![](https://img-blog.csdnimg.cn/20201013111454248.png?watmark/2/text/aHR0cDovL2cuY3Nkbi5uZXQvV5bm1hbjIzMw==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
直接上代码:
'''
Description:
Author: Weijian Ma
Date: 2020-10-10 09:13:23
LastEditTime: 2020-10-13 11:21:11
LastEditors: Weijian Ma
'''
## 对图像进行直方图均衡
import cv2
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] # 用黑体显示中文
plt.rcParams['axes.unicode_minus']=False # 正常显示负号
## 定义一个函数,用于转换图像
def T(r_k, nkCum, height, width):
s_k = (7/(height*width))*nkCum[r_k]
return s_k
img = np.array([[0, 1, 1, 0, 0],
[3, 4, 0, 2, 1],
[3, 2, 1, 0, 2],
[2, 1, 1, 0, 3],
[2, 3, 3, 4, 4]])
height = img.shape[0]
width = img.shape[1]
## 计算每个灰度出现的次数
n_k = np.zeros(256,dtype=np.int32)
for i in range(height):
for j in range(width):
n_k[img[i][j]] = n_k[img[i][j]]+1
nkCum = np.cumsum(n_k)
print(np.around(T(img,nkCum,height,width)))
结果:
[[2. 3. 3. 2. 2.]
[6. 7. 2. 5. 3.]
[6. 5. 3. 2. 5.]
[5. 3. 3. 2. 6.]
[5. 6. 6. 7. 7.]]
(4) 画出均衡化后的直方图。
也就是把第(1)题中的图像换成第(3)题中得到的结果。
代码如下,
'''
Description:
Author: Weijian Ma
Date: 2020-10-13 10:09:21
LastEditTime: 2020-10-13 11:26:00
LastEditors: Weijian Ma
'''
import numpy as np
import matplotlib.pyplot as plt
import cv2
plt.rcParams['font.sans-serif']=['SimHei'] # 用黑体显示中文
plt.rcParams['axes.unicode_minus']=False # 正常显示负号
img = np.array([[2, 3, 3, 2, 2],
[6, 7, 2, 5, 3],
[6, 5, 3, 2, 5],
[5, 3, 3, 2, 6],
[5, 6, 6, 7, 7]])
countNum = np.zeros(8)
for i in range(8):
countNum[i] = np.sum(img==i)
## 画出灰度直方图
fig = plt.figure(figsize=(10,4))
sub01 = plt.subplot(121)
sub02 = plt.subplot(122)
sub01.set_xticks(np.arange(8))
sub02.set_xticks(np.arange(8))
sub01.set_title("直方图")
sub02.set_title("归一化直方图")
sub01.set_xlabel("灰度")
sub01.set_ylabel("频数")
sub02.set_xlabel("灰度")
sub02.set_ylabel("频率")
sub01.bar(np.arange(8), countNum)
sub02.bar(np.arange(8), countNum/np.sum(countNum))
plt.show()
结果:
![](https://img-blog.csdnimg.cn/2020101311270422.png?watmark/2/text/aHR0cDovL2cuY3Nkbi5uZXQvV5bm1hbjIzMw==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
题目2:用大小为3x3的盒状滤波器对图2的图像进行滤波,图像上方下方各补一行0,左右侧各补一列0,给出滤波响应图。
3x3的盒装滤波:
![](https://img-blog.csdnimg.cn/20201013123707215.png?watmark/2/text/aHR0cDovL2cuY3Nkbi5uZXQvV5bm1hbjIzMw==/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/gravity/SouthEast)
图2 上下左右各补一行0后:
道理很简单,手算非常容易实现。
下面是python代码实现:
'''
Description:
Author: Weijian Ma
Date: 2020-10-13 12:44:16
LastEditTime: 2020-10-13 13:08:33
LastEditors: Weijian Ma
'''
import cv2
import numpy as np
import matplotlib.pyplot as plt
## 图像
img = np.array([[0, 1, 1],
[0, 4, 0],
[2, 2, 2]])
## 核
kernel = np.array([[1, 1, 1],
[1, 1, 1],
[1, 1, 1]]) / 9
## 图像上下左右补0
img = np.insert(img, 0, values=np.array([0, 0, 0]), axis=1)
img = np.insert(img, 4, values=np.array([0, 0, 0]), axis=1)
img = np.insert(img, 0, values=np.array([0, 0, 0, 0, 0]), axis=0)
img = np.insert(img, 4, values=np.array([0, 0, 0, 0, 0]), axis=0)
box = np.zeros((3,3))
## 滤波操作
for i in range(1,4):
for j in range(1,4):
box[i-1][j-1] = np.sum(img[i-1:i+2,j-1:j+2] * kernel)
print("滤波结果:")
print(box)
print("四舍五入之后:")
print(np.around(box))
结果:
滤波结果:
[[0.55555556 0.66666667 0.66666667]
[1. 1.33333333 1.11111111]
[0.88888889 1.11111111 0.88888889]]
四舍五入之后:
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
所以,滤波响应图(这里用矩阵表示):
[[1. 1. 1.]
[1. 1. 1.]
[1. 1. 1.]]
题目3:Roberts、Sobel、 Prewitt 算子对图3对应的图像进行边缘检测,给出检测的结果。
这道题的python代码不知道怎么实现,python opencv好像并不容易实现(总是出现一些不知名的BUG),如果有大佬能补上更好hhh
用matlab好像挺容易实现的,鉴于我并没有正版的matlab,下面的这段程序全部都是在Octave上运行的,当然,在matlab上同样能运行,效果是一样的。
Octave代码(貌似有点题文不符,我对不起大家):
img = [1 1 1 1 1 1;1 1 1 1 1 1;1 1 1 1 2 3;2 2 2 2 2 3;2 2 2 2 2 3;2 2 2 3 3 3;3 3 3 3 3 3;3 3 3 3 3 3];
h1 = [-1, 0; 0, 1];
h2 = [0, -1; 1, 0];
Roberts_Edge_result = abs(imfilter(img, h1)) + abs(imfilter(img, h2))
Prewitt = fspecial('prewitt');
Prewitt_Edge_result = abs(imfilter(img, Prewitt)) + abs(imfilter(img, Prewitt'))
Sobel = fspecial('sobel');
Sobel_Edge_result = abs(imfilter(img, Sobel)) + abs(imfilter(Img, Sobel'))
结果: