RS—|遥感数字图像处理编程练习(python)

一:模拟计算图像直方图和累计直方图

① 调用的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图像处理】图片读取/提取直方图
。。。。

  • 2
    点赞
  • 27
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值