遥感影像16位转8位(python)

遥感影像16位转8位(python)

不同于一般的8位数据，数据库中的遥感图像都是16位的无符号整形数据，而且大部分集中在[0,500]的范围内，对于8位的RGB图像来说根本无法显示正确的图像。因此，必须对灰度级进行压缩，从65536个灰度级压缩到256个灰度级范围内。

但是，仅仅使用线性压缩是无法解决问题的，因为前面也说过，数据主要集中在[0,500]范围内，即它只占用了前面的几百个灰度级，后面的根本没用到。为此，我们必须丢掉出现概率小的灰度级，保留出现概率大的灰度级范围。

****1.读16位图像数据

def read_mat(input_file):
    """
    读取影像
    Inputs:
    input_file:16位图像数据的路径
    Outputs:
    array_data：保存图像数据的numpy.array数组
    rows:高度
    cols:宽度
    bands:深度
    """
    result_data = scio.loadmat(input_file)
    dataset= result_data['imgPS']

    cols = dataset.shape[0]
    rows = dataset.shape[1]
    bands = dataset.shape[2]

    array_data = np.zeros((rows, cols,bands))

    for i in range(bands):
        band = dataset[:,:,i]
        array_data[:, :,i] = band

    return array_data,  rows, cols, bands

****2.直方图统计

def cumulativehistogram(array_data, band_min, band_max):
    """
    累计直方图统计
    Inputs:
    array_data:16位图像数据
    band_min:最小像素值
    band_max:最大像素值
    Outputs：
    cutmax:累积分布函数（CDF）==0.98对应的灰度级
    cutmin:累积分布函数（CDF）==0.02对应的灰度级
    """

    # 逐波段统计最值
    rows=array_data.shape[0]
    cols=array_data.shape[1]

    gray_level = int(band_max - band_min + 1)
    gray_array = np.zeros(gray_level)

    for row in range(rows):
        for col in range(cols):
            gray_array[int(array_data[row, col] - band_min)] += 1

    count_percent2 = rows*cols * 0.02
    count_percent98 = rows*cols* 0.98

    cutmax = 0
    cutmin = 0

    for i in range(1, gray_level):
        gray_array[i] += gray_array[i - 1]
        if (gray_array[i] >= count_percent2 and gray_array[i - 1] <= count_percent2):
            cutmin = i + band_min

        if (gray_array[i] >= count_percent98 and gray_array[i - 1] <= count_percent98):
            cutmax = i + band_min

    return cutmin, cutmax

****3.灰度压缩

def compress(origin_16, output_8):
    """
    Input:
    origin_16:16位图像路径
    Output:
    output:8位图像路径
    """
    array_data, rows, cols, bands = read_mat(origin_16)

    compress_data = np.zeros(( rows, cols,bands))

    for i in range(bands):
        band_max = np.max(array_data[:, :,i])
        band_min = np.min(array_data[:, :,i])

        cutmin, cutmax = cumulativehistogram(array_data[:, :,i], band_min, band_max)

        compress_scale = (cutmax - cutmin) / 255

        for j in range(rows):
            for k in range(cols):
                if (array_data[j, k,i] < cutmin):
                    array_data[j, k,i] = cutmin

                if (array_data[j, k,i] > cutmax):
                    array_data[j, k,i] = cutmax

                compress_data[j, k,i] = (array_data[j, k,i] - cutmin) / compress_scale

    write_mat(output_8, compress_data)

****5.生成8位图像

def write_mat(output_file,array_data):
    output_file_order=["B","G","R","NIR"]
    rows=array_data.shape[0]
    cols=array_data.shape[1]
    bands=array_data.shape[2]
    for band in range(bands):
        cv2.imwrite(output_file+output_file_order[band]+".tif", array_data[:,:,band])

链接:
link.
link.