def gaussian_noise(img, mean, sigma):
    """
    此函数用将产生的高斯噪声加到图片上
    传入: img   :  原图   mean  :  均值  sigma :  标准差
    返回: gaussian_out : 噪声处理后的图片
    """
    # 将图片灰度标准化
    img = img / 255
    # 产生高斯 noise
    noise = np.random.normal(mean, sigma, img.shape)
    # 将噪声和图片叠加
    gaussian_out = img + noise
    # 将超过 1 的置 1，低于 0 的置 0
    gaussian_out = np.clip(gaussian_out, 0, 1)
    # 将图片灰度范围的恢复为 0-255
    gaussian_out = np.uint8(gaussian_out * 255)
    # 将噪声范围搞为 0-255
    # noise = np.uint8(noise*255)
    return gaussian_out  # 这里也会返回噪声，注意返回值

图1的待实验图片加入高斯噪声（以mean = 0 ， sigma = 0.1为例）前后如图所示

2.椒盐噪声

椒盐噪声：椒盐噪声就是给图片添加黑白噪点，椒指的是黑色的噪点(0,0,0)盐指的是白色的噪点(255,255,255)，通过设置amount来控制添加噪声的比例，值越大添加的噪声越多，图像损坏的更加严重。

加入椒盐噪声的代码实现（以灰度图的加噪为例）如下：

def sp_noise(noise_img, proportion):
    """
    添加椒盐噪声
    proportion的值表示加入噪声的量，可根据需要自行调整
    return: img_noise
    """
    height, width = noise_img.shape[0], \\
                    noise_img.shape[1]  # 获取高度宽度像素值
    num = int(height * width * proportion)  # 一个准备加入多少噪声小点
    for i in range(num):
        w = random.randint(0, width - 1)
        h = random.randint(0, height - 1)
        if random.randint(0, 1) == 0:
            noise_img[h, w] = 0
        else:
            noise_img[h, w] = 255
    return noise_img

图1的待实验图片加入高斯噪声（以mean = 0 ， sigma = 0.1为例）与椒盐噪声（proportion = 0.1为例）前后如图所示。

三、图像的去/降噪

噪声的存在严重影响了遥感图像的质量，因此在图像增强处理和分类处理之前，必须使用图像的去/降噪方法予以纠正。本实验以中值滤波与均值滤波为例。

1.均值滤波

包含在滤波器领域内像素的平均值，本质上就是对目标像素及周边像素取平均值后再填会目标像素来实现滤波目的的方法，加权平均公式如图所示。

使用均值滤波器的代码实现（以灰度图的线性非加权降噪为例）如下：

# 均值降噪转换器
def convert_2d(r):
    n = 3
    # 3*3 滤波器, 每个系数都是 1/9
    window = np.ones((n, n)) / n ** 2
    # 使用滤波器卷积图像
    # mode = same 表示输出尺寸等于输入尺寸
    # boundary 表示采用对称边界条件处理图像边缘
    s = scipy.signal.convolve2d(r,
      window, mode='same', boundary='symm')
    return s.astype(np.uint8)

加入高斯噪声（以mean = 0 ， sigma = 0.1为例）与椒盐噪声（proportion = 0.1为例）的图1的待实验图片降噪前后如图所示。

分析：从图可见，线性非加权的均值滤波器降噪效果较差，且滤波后图像模糊，丢失较多信息，因此实际使用时一般使用加权均值滤波器或其他滤波器。

2.中值滤波

类似均值滤波的方法，用像素领域内的中间值代替该像素，公式与原理如0所示。

使用中值滤波器的代码实现（以灰度图的降噪为例）如下：

# 自定义中值降噪滤波器
def medium_filter(im, x, y, step):
    sum_s = []
    for k in range(-int(step / 2), int(step / 2) + 1):
        for m in range(-int(step / 2), int(step / 2) + 1):
            sum_s.append(im[x + k][y + m])
    sum_s.sort()
    return sum_s[(int(step * step / 2) + 1)]

加入高斯噪声（以mean = 0 ， sigma = 0.1为例）与椒盐噪声（proportion = 0.1为例）的图1的待实验图片降噪前后如图所示。

分析：从图可见，图片中噪点明显减少，证明中值滤波器降噪效果较好（相对均值滤波器），且滤波在去除噪音的同时，可以比较好地保留边的锐度和图像的细节，因此本实验去/降噪方法最终使用中值滤波器。

四、直方图与均衡化

1.直方图

直方图：在图像处理中，也经常需要分析图像的亮度（即像素级的分布情况），这就需要用到直方图了，如颜色直方图、灰度直方图等。

直方图通过对灰度级在范围[0，L-1]的数字图像（彩图分解RGB）使用如图的离散函数+归一化后构造的。

直方图构造（以灰度图为例）核心代码如下：

def histogram(Image):  # 灰度图像直方图统计
    im = Image
    a = [0] * 256  # 构造储存像素数的一维数组
    w = im.shape[0]  # 得到图像的宽高
    h = im.shape[1]
    for i in range(w):  # 遍历所有像素点统计
        for j in range(h):
            gray = im[i, j]
            a[gray] += 1
    y = a  # y表示灰度值的像素个数
    x = [i for i in range(256)]  # x表示0-256的灰度

图1（原图）与图2（原图）的直方图如图所示。

2.直方图均衡化

直方图均衡化：把一个已知灰度概率密度分布的图像经过一种灰度变换（变换函数的离散形式如图所示），使之演变为一幅具有均匀灰度概率密度分布（且尽量占满全部灰度）的新图像，以此增加了像素灰度值的动态范围从而达到增强图像整体对比度的效果。

直方图均衡化的实现步骤如图所示

五、灰度图均衡化

1. 非调库均衡化

根据直方图与均衡化的定义对降噪后的图1进行直方图构建与均衡化，非调库灰度图均衡化代码实现如下：

def homogenisation(a, Image):  # 均衡化函数：a--灰度像素数，Image--图片
    im = array(Image)
    b = [0] * 256  # 储存灰度像素所占的比例
    c = [0] * 256  # 储存累计分布概率
    w = im.shape[0]  # 图像的宽和高
    h = im.shape[1]
    mn = w * h  # mn表示像素点总数
    for i in range(len(a)):  # 计算灰度分布密度
        b[i] = a[i] / mn
    for i in range(len(c)):  # 计算累计分布密度c和累计直方图数据a
        if i == 1:
            c[i] = b[i]
        else:
            c[i] = c[i - 1] + b[i]
            a[i] = int(255 * c[i])
    for i in range(w):  # 更新每一个像素点的灰度值
        for j in range(h):
            im[i, j] = a[Image[i, j]]
    return im

均衡化前后直方图如图所示：

均衡化前后图片如图所示：

分析：根据图可见，均衡化后直方图灰度分布均匀且基本占满所有可能灰度，满足相关定义，观察图，均衡化后图片对比度增强，可见细节增加，满足预期，均衡化成功。

2. 调库均衡化

Python中Opencv（cv2）库中灰度均衡化函数为cv2.equalizeHist(img)，返回的即为均衡化后的灰度图像。

六、彩图均衡化

彩图本质是为RGB三通道图，均衡化时可根据不同考虑使用不同方法均衡化，本实验以拆分合并（复现库函数）、HSI模型<->RGB模型互换、库函数三种方式为例实现彩图均衡化。

1.拆分合并均衡化

核心原理（类似opencv库函数）：先将彩图RGB三通道拆分，采取类似灰度均衡化的方法对每一通道进行均衡化，最后进行通道合并，得到均衡化后的彩图，核心代码如下：

def my_equalizehist(single):
    hist = cv.calcHist([single], [0], None, [256], [0, 256])
    num_of_pixels = single.size
    ratio = np.zeros(256)
    transf_map = np.zeros(256)
    result = single.copy()
    j = 0
    for i in hist:
        if j > 0:
            ratio[j] = i / num_of_pixels + ratio[j - 1]
        else:
            ratio[j] = i / num_of_pixels
        transf_map[j] = round(ratio[j] * 255)
        j = j + 1
    for i in range(single.shape[0]):
        for j in range(single.shape[1]):
            result[i][j] = transf_map[single[i][j]]
        # result[result==j]=k[j]
    return result

均衡化前后图片如图：

2.库函数均衡化

Python中Opencv（cv2）库中彩图均衡化的函数为cv.equalizeHist(img)，可用于均衡化各个通道并使用cv.merge合并，具体用法如下：

hist_res=cv.calcHist([res],[0],None,[256],[0,256])
plt.plot(hist_res);
plt.show()
EB = cv.equalizeHist(B)
EG = cv.equalizeHist(G)
ER = cv.equalizeHist(R)
equal_test = cv.merge((EB, EG, ER))  # merge it back

3.HSI模型<->RGB模型互换均衡化

HSI颜色模型是一个满足计算机数字化颜色管理需要的高度抽象模拟的数学模型。HSI模型是从人的视觉系统出发，直接使用颜色三要素–色调（Hue）、饱和度（Saturation）和亮度（Intensity，有时也翻译作密度或灰度）来描述颜色，RBG转换到HIS空间公式如图所示，HSI转化成RGB的公式如图所示。

代码如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt
import collections


def rgb_hsi(rgb_image):
    # 保存原始图像的行列数
    rows = int(rgb_image.shape[0])
    cols = int(rgb_image.shape[1])
    # 图像复制
    hsi_image = rgb_image.copy()
    # 通道拆分
    b = rgb_image[:, :, 0]
    g = rgb_image[:, :, 1]
    r = rgb_image[:, :, 2]
    # 归一化到[0,1]
    b = b / 255.0
    g = g / 255.0
    r = r / 255.0
    for i in range(rows):
        for j in range(cols):
            num = 0.5 * ((r[i, j]-g[i, j])+(r[i, j]-b[i, j]))
            den = np.sqrt((r[i, j]-g[i, j])**2+(r[i, j]-b[i, j])*(g[i, j]-b[i, j]))
            if den == 0:
                hsi_h = 0
            else:
                theta = float(np.arccos(num / den))
                if b[i, j] <= g[i, j]:
                    hsi_h = theta
                else:
                    hsi_h = 2*np.pi - theta

            min_RGB = min(min(b[i, j], g[i, j]), r[i, j])
            sum = b[i, j]+g[i, j]+r[i, j]
            if sum == 0:
                hsi_s = 0
            else:
                hsi_s = 1 - 3*min_RGB/sum

            hsi_h = hsi_h/(2*np.pi)
            hsi_i = sum/3.0
            # 输出HSI图像，扩充到255以方便显示，一般H分量在[0,2pi]之间，S和I在[0,1]之间
            hsi_image[i, j, 0] = hsi_h*255
            hsi_image[i, j, 1] = hsi_s*255
            hsi_image[i, j, 2] = hsi_i*255
    return hsi_image


def hsi_rgb(hsi_image):
    # 保存原始图像的行列数
    rows = np.shape(hsi_image)[0]
    cols = np.shape(hsi_image)[1]
    # 对原始图像进行复制
    rgb_image = hsi_image.copy()
    # 对图像进行通道拆分
    hsi_h = hsi_image[:, :, 0]
    hsi_s = hsi_image[:, :, 1]
    hsi_i = hsi_image[:, :, 2]
    # 把通道归一化到[0,1]
    hsi_h = hsi_h / 255.0
    hsi_s = hsi_s / 255.0
    hsi_i = hsi_i / 255.0
    B, G, R = hsi_h, hsi_s, hsi_i
    for i in range(rows):
        for j in range(cols):
            hsi_h[i, j] *= 360
            if 0 <= hsi_h[i, j] < 120:
                B = hsi_i[i, j] * (1 - hsi_s[i, j])
                R = hsi_i[i, j] * (1 + (hsi_s[i, j] * np.cos(hsi_h[i, j] * np.pi / 180)) / np.cos(
                    (60 - hsi_h[i, j]) * np.pi / 180))
                G = 3 * hsi_i[i, j] - (R + B)
            elif 120 <= hsi_h[i, j] < 240:
                hsi_h[i, j] = hsi_h[i, j] - 120
                R = hsi_i[i, j] * (1 - hsi_s[i, j])
                G = hsi_i[i, j] * (1 + (hsi_s[i, j] * np.cos(hsi_h[i, j] * np.pi / 180)) / np.cos(
                    (60 - hsi_h[i, j]) * np.pi / 180))
                B = 3 * hsi_i[i, j] - (R + G)
            elif 240 <= hsi_h[i, j] <= 300:
                hsi_h[i, j] = hsi_h[i, j] - 240
                G = hsi_i[i, j] * (1 - hsi_s[i, j])
                B = hsi_i[i, j] * (1 + (hsi_s[i, j] * np.cos(hsi_h[i, j] * np.pi / 180)) / np.cos(
                    (60 - hsi_h[i, j]) * np.pi / 180))
                R = 3 * hsi_i[i, j] - (G + B)
            rgb_image[i, j, 0] = B * 255
            rgb_image[i, j, 1] = G * 255
            rgb_image[i, j, 2] = R * 255
    return rgb_image


# 计算灰度图的直方图
def draw_histogram(grayscale):
    # 对图像进行通道拆分
    hsi_i = grayscale[:, :, 2]
    color_key = []
    color_count = []
    color_result = []
    histogram_color = list(hsi_i.ravel())  # 将多维数组转换成一维数组
    color = dict(collections.Counter(histogram_color))  # 统计图像中每个亮度级出现的次数
    color = sorted(color.items(), key=lambda item: item[0])  # 根据亮度级大小排序
    for element in color:
        key = list(element)[0]
        count = list(element)[1]
        color_key.append(key)
        color_count.append(count)
    for i in range(0, 256):
        if i in color_key:
            num = color_key.index(i)
            color_result.append(color_count[num])
        else:
            color_result.append(0)
    color_result = np.array(color_result)
    return color_result


def histogram_equalization(histogram_e, lut_e, image_e):
    sum_temp = 0
    cf = []
    for i in histogram_e:
        sum_temp += i
        cf.append(sum_temp)
    for i, v in enumerate(lut_e):
        lut_e[i] = int(255.0 * (cf[i] / sum_temp) + 0.5)
    equalization_result = lut_e[image_e]
    return equalization_result


x = []
for i in range(0, 256):  # 横坐标
    x.append(i)

# 原图及其直方图
rgb_image = cv2.imread("lena.jpg")
cv2.imshow('rgb', rgb_image)
histogram = draw_histogram(rgb_image)
plt.bar(x, histogram)  # 绘制原图直方图
plt.savefig('before_histogram.png')
plt.show()

# rgb转hsi
hsi_image = rgb_hsi(rgb_image)
cv2.imshow('hsi_image', hsi_image)
cv2.imwrite('hsi_result.png', hsi_image)

# hsi在亮度分量上均衡化
histogram_1 = draw_histogram(hsi_image)
lut = np.zeros(256, dtype=hsi_image.dtype)  # 创建空的查找表
result = histogram_equalization(histogram_1, lut, hsi_image)  # 均衡化处理
cv2.imshow('his_color_image', result)
cv2.imwrite('his_color.png', result)  # 保存均衡化后图片

# hsi转rgb
image_equ = cv2.imread('his_color.png')  # 读取图像
rgb_result = hsi_rgb(image_equ)
cv2.imshow('rgb_image', rgb_result)
cv2.imwrite('gbr_result.png', rgb_result)

rgb = cv2.imread("gbr_result.png")
histogram_2 = draw_histogram(rgb)
plt.bar(x, histogram_2)
plt.savefig('after_histogram.png')
plt.show()

plt.plot(x, lut)  # 绘制灰度级变换曲线图
plt.savefig('Grayscale_transformation_curve.png')
plt.show()

cv2.waitKey(0)
cv2.destroyAllWindows()

七、图像增强

图像增强是使图像更适合于特定应用的图像处理技术，主要分为空间域图像增强与频率域图像增强，常见的图像增强方法框架如图所示。除直方图均衡化，几种常见的图像增强算法简介与原理如下：

（1）对数变换（非线性变换）：通过公式s= c log (1+r)对灰度变换，使窄带低灰度输入图像值映射为一宽带输出值，灰度变化前后如图35，图像增强样例（傅里叶频谱）如图36。

图35 对数变换灰度变化前后

图36 对数变换图像增强样例（图源ppt）

（2）减法运算（代数运算）：图像的减法运算属于代数运算，即两幅输入图像进行点对点的加、减、乘、除计算而得到输出图像的运算。减法运算则是通过公式C(x,y) = f(x,y) - h(x,y)得到图像，常用于除去背景噪声（缓慢变化的加性非随机噪声），图像增强样例如图37所示。

图37 减法图像增强样例（图源ppt）

（3）平滑滤波器：本实验用到的均值、中值滤波器均属于平滑滤波器，原理可见实验步骤与过程，核心为用像素领域内的某些值代替该像素。

（4）锐化滤波器：包含Laplace算子（微分算子，应用强调在图像中灰度的突变区域，并不强调灰度级缓慢变化的区域）与Sobel算子（主要用在工业检测图像的边缘检测），算子运行时采取类似卷积的方式，将模板在图像上移动并在每个位置上计算对应中心像素的梯度值，即对一幅灰度图像求梯度所得的结果是一幅梯度图，Sobel算子定义如图38所示，锐化前后如图39所示。