目录
- 一：模拟计算图像直方图和累计直方图
- 二：计算图像的均值、标准差、相关系数和协方差
- 三：利用模板进行卷积运算
- 四：获取彩色图像的直方图
- 五：图像直方图均衡化
一：模拟计算图像直方图和累计直方图
① 调用的python工具包：
② 利用生成随机数的函数生成一个256*256的二位数组模拟图像的灰度值的排列：
③ 遍历影像中的每一个像元的像元值，将其存为一个一维数组并进行排序：
④ 由于空值通常用-3.4028235e+38表达，避免错误统计需要提前剔除（利用数组切片的方法：
⑤ 绘制灰度直方图：
import numpy as np
import matplotlib.pyplot as plt
im_data = np.random.randint(low=0, high=256, size=256 * 256)
im_data.shape = 256, 256
print(im_data.shape)
# ##########显示灰度直方图########
# 遍历影像中的每一个像元的像元值
data = []
for i in range(im_data.shape[0]):
for j in range(im_data.shape[1]):
# print(j)
data.append(im_data[i][j])
data.sort()
# 统计最大最小值
data = np.array(data)
print(data.min(), data.max())
# 由于空值通常用-3.4028235e+38表达，避免错误统计需要提前剔除
a = sum(data == 0)
print(a)
print(data.shape)
data = data[a:] # 切片，即a[1:]相当于a[1:len(data)]
# 根据影像中最大最小值设定坐标轴
bins = np.linspace(start=0, stop=256, num=256)
# 绘制直方图，设定直方图颜色
plt.hist(data, bins, facecolor="black")
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title('灰度分布直方图')
# 显示中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
# 展现绘制得到的图片
plt.show()
# 导出绘制得到的图片
# plt.savefig('C:/test2.jpg')
band1_data = im_data
运行结果：
⑥ 绘制累计灰度直方图（只需将hist函数的cumulative参数设置为True）：
运行结果：
（因为是随机生成数组，所以灰度值都分布比较均匀）
二：计算图像的均值、标准差、相关系数和协方差
① 在题目一的基础上利用相关函数可直接计算均值与标准差：
运行结果：
② 计算协方差矩阵和相关系数矩阵
a.先定义一个原始列表（模拟将图像灰度值矩阵转为数组后的数据形式），然后再计算其平均值：
b.再定义一个列表与原始列表结合，计算其协方差矩阵：
c. 再定义一个列表与原始列表结合，计算其相关系数矩阵：
origin_list = [[1, 0, 1, 0], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 0, 1]]
index_1 = 0
avg = [0, 0, 0, 0]
for li in origin_list:
for num in li:
avg[index_1] += num/len(li)
index_1 += 1
print('各个平均数: ', avg)
xiefangcha = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
i = -1
for li_1 in origin_list:
i += 1
j = -1
for li_2 in origin_list:
j += 1
for num_1, num_2 in zip(li_1, li_2):
xiefangcha[i][j] += (num_1 - avg[i])*(num_2 - avg[j])/3
print('协方差矩阵: ', xiefangcha)
xiangguanjuzheng = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
i = 0
for xie in xiefangcha:
j = 0
for xi in xie:
xiangguanjuzheng[i][j] = xi/((xiefangcha[i][i]**0.5)*(xiefangcha[j][j]**0.5))
j += 1
i += 1
print('相关矩阵: ', xiangguanjuzheng)
结果输出：
各个平均数: [0.5, 0.5, 0.5, 0.5]
协方差矩阵: [[0.3333333333333333, 0.3333333333333333, -0.3333333333333333, 0.0], [0.3333333333333333, 0.3333333333333333, -0.3333333333333333, 0.0], [-0.3333333333333333, -0.3333333333333333, 0.3333333333333333, 0.0], [0.0, 0.0, 0.0, 0.3333333333333333]]
相关矩阵: [[1.0, 1.0, -1.0, 0.0], [1.0, 1.0, -1.0, 0.0], [-1.0, -1.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]
三：利用模板进行卷积运算
① 定义图像与模板
② 进行卷积运算的方法
③ 参数的传入与结果输出
import numpy as np
# 定义图像与模板
input_data = [
[[1, 1, 7, 1, 8, 1, 7, 1, 1],
[1, 1, 1, 1, 5, 1, 1, 1, 1],
[1, 1, 1, 5, 5, 5, 1, 1, 7],
[7, 1, 1, 5, 5, 5, 1, 1, 1],
[1, 1, 1, 5, 5, 5, 1, 8, 1],
[1, 8, 1, 1, 5, 1, 1, 1, 1],
[1, 8, 1, 1, 5, 1, 1, 8, 1],
[1, 1, 1, 1, 5, 1, 1, 1, 1],
[1, 1, 7, 1, 8, 1, 7, 1, 1]],
]
weights_data = [
[[0, -1, 0],
[-1, 4, -1],
[0, -1, 0]],
]
def compute_conv(fm, kernel):
[h, w] = fm.shape
[k, _] = kernel.shape
r = int(k / 2)
# 保存计算结果
rs = np.zeros([h - 2, w - 2], np.float32)
# 将输入在指定该区域赋值，即除了4个边界后，剩下的区域
# 对每个点为中心的区域遍历
for i in range(1, h - 1):
for j in range(1, w - 1):
# 取出当前点为中心的k*k区域
roi = fm[i - 1:i + r + 1, j - 1:j + r + 1]
# 计算当前点的卷积,对k*k个点点乘后求和
rs[i - 1][j - 1] = np.sum(roi * kernel)
return rs
def my_conv2d(input, weights):
[c, h, w] = input.shape
[_, k, _] = weights.shape
outputs = np.zeros([h - 2, w - 2], np.float32)
# 对每个feature map遍历，从而对每个feature map进行卷积
for i in range(c):
# feature map==>[h,w]
f_map = input[i]
# kernel ==>[k,k]
w = weights[i]
rs = compute_conv(f_map, w)
outputs = outputs + rs
return outputs
def main():
# shape=[c,h,w]
input = np.asarray(input_data, np.float32)
# shape=[in_c,k,k]
weights = np.asarray(weights_data, np.float32)
rs = my_conv2d(input, weights)
print(rs)
if __name__ == '__main__':
main()
运行结果（未输出边缘的像素值）：
[[ 0. -6. -8. 5. -8. -6. 0.]
[ 0. -4. 8. 0. 8. -4. -6.]
[ -6. -4. 4. 0. 4. -4. -7.]
[ -7. -4. 8. 0. 8. -11. 28.]
[ 21. -7. -8. 8. -8. 0. -14.]
[ 21. -7. -4. 8. -4. -7. 28.]
[ -7. -6. -4. 5. -4. -6. -7.]]
四：获取彩色图像的直方图
这里用到的是PIL(Python Imaging Library)库，Python平台的图像处理标准库。
① 导入用到的相关的包：
②读入一张彩色tif图像，并分割RGB三种色彩：
③要对图像求直方图，就需要先把图像矩阵进行flatten操作，使之变为一维数组，然后再进行统计，分别提取R\G\B统计值，进行分区绘图：
④相关显示设置，显示原图像：
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
img = Image.open(r"D:\数据\影像test.tif")
r, g, b = img.split() # 分离R\G\B
# 要对图像求直方图，就需要先把图像矩阵进行flatten操作，使之变为一维数组，然后再进行统计
# 分别提取R\G\B统计值，叠加绘图
fig = plt.figure("ImageHist")
fig.add_subplot(2, 2, 1)
ar = np.array(r).flatten()
plt.hist(ar, facecolor='red', bins=256, edgecolor='red', )
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title('R 灰度分布直方图')
fig.add_subplot(2, 2, 2)
ag = np.array(g).flatten()
plt.hist(ag, bins=256, facecolor='green', edgecolor='green', )
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title('G 灰度分布直方图')
fig.add_subplot(2, 2, 3)
ab = np.array(b).flatten()
plt.hist(ab, bins=256, facecolor='blue', edgecolor='blue', )
# 横坐标轴名称
plt.xlabel('像元值')
# 纵坐标轴名称
plt.ylabel('频数')
# 图表头名称
plt.title('B 灰度分布直方图')
ax = fig.add_subplot(2, 2, 4)
plt.imshow(img)
# 显示中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
# 调整子图组合及子图之间的位置与间距
plt.subplots_adjust(left=0.1, right=0.9, hspace=0.5, wspace=0.3)
plt.show()
运行显示结果如下图：
五：图像直方图均衡化
直方图均衡化的基本原理是：对在图像中像素个数多的灰度值（即对画面起主要作用的灰度值）进行展宽，而对像素个数少的灰度值（即对画面不起主要作用的灰度值）进行归并，从而增大对比度，使图像清晰，达到增强的目的。
这里用到了openCV图像处理库，核心代码主要思路为：遍历图像每个像素的灰度，算出每个灰度的概率（n/MN-n是每个灰度的个数，MN是像素总数），用L-1乘以所得概率得到新的灰度：
import cv2
import numpy as np
img = cv2.imread("img.tif")
img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 直方图统计
def pix_gray(img_gray):
h = img_gray.shape[0]
w = img_gray.shape[1]
gray_level = np.zeros(256)
gray_level2 = np.zeros(256)
for i in range(1, h - 1):
for j in range(1, w - 1):
gray_level[img_gray[i, j]] += 1 # 统计灰度级为img_gray[i,j的个数
for i in range(1, 256):
gray_level2[i] = gray_level2[i - 1] + gray_level[i] # 统计灰度级小于img_gray[i,j]的个数
return gray_level2
# 直方图均衡化
def hist_gray(img_gary):
h, w = img_gary.shape
gray_level2 = pix_gray(img_gray)
lut = np.zeros(256)
for i in range(256):
lut[i] = 255.0 / (h * w) * gray_level2[i] # 得到新的灰度级
lut = np.uint8(lut + 0.5)
out = cv2.LUT(img_gray, lut)
return out
cv2.imshow('input', img_gray)
out_img = hist_gray(img_gray)
cv2.imshow('output', out_img)
cv2.waitKey(0)
# cv2.destroyWindow()
再利用前面同样的方法绘制原图与均衡化后的直方图，进行对比，如下图，可以看出均衡化后的图像的对比度得到了增强：
…
学习遥感数字图像处理中做的一些的编程练习。
借鉴参考：
【Python图像处理】图片读取/提取直方图
。。。。
