前言

直方图匹配可以使得影像A的灰度分布变换为影像B的灰度分布。

代码

#coding=utf-8
import os
import matplotlib
import numpy as np
from numba import jit
from osgeo import gdal
from PIL import Image
from matplotlib import pyplot as plt
import cv2 as cv

@jit
def arrayToHist(grayArray,nums):
    if(len(grayArray.shape) != 2):
        print("length error")
        return None
    w,h = grayArray.shape
    hist = {}
    for k in range(nums):
        hist[k] = 0
    for i in range(w):
        for j in range(h):
            if(hist.get(grayArray[i][j]) is None):
                hist[grayArray[i][j]] = 0
            hist[grayArray[i][j]] += 1
    n = w*h
    for key in hist.keys():
        hist[key] = float(hist[key])/n
    return hist
@jit
def histMatch(grayArray,h_d):
    tmp = 0.0
    h_acc = h_d.copy()
    for i in range(256):
        tmp += h_d[i]
        h_acc[i] = tmp

    h1 = arrayToHist(grayArray,256)
    tmp = 0.0
    h1_acc = h1.copy()
    for i in range(256):
        tmp += h1[i]
        h1_acc[i] = tmp
    M = np.zeros(256)
    for i in range(256):
        idx = 0
        minv = 1
        for j in h_acc:
            if (np.fabs(h_acc[j] - h1_acc[i]) < minv):
                minv = np.fabs(h_acc[j] - h1_acc[i])
                idx = int(j)
        M[i] = idx
    des = M[grayArray]
    return des

@jit
def equalization(grayArray,h_s,nums):
    tmp = 0.0
    h_acc = h_s.copy()
    for i in range(256):
        tmp += h_s[i]
        h_acc[i] = tmp

    if(len(grayArray.shape) != 2):
        print("length error")
        return None
    w,h = grayArray.shape
    des = np.zeros((w,h),dtype = np.uint8)
    for i in range(w):
        for j in range(h):
            des[i][j] = int((nums - 1)* h_acc[grayArray[i][j] ] +0.5)
    return des

@jit
def drawHist(hist,name):
    keys = hist.keys()
    values = hist.values()
    x_size = len(hist)-1
    axis_params = []
    axis_params.append(0)
    axis_params.append(x_size)

    if name != None:
        plt.title(name)
    plt.bar(tuple(keys),tuple(values))


def read_img(filename):
    dataset=gdal.Open(filename)

    im_width = dataset.RasterXSize
    im_height = dataset.RasterYSize

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0,0,im_width,im_height)

    del dataset 
    return im_proj,im_geotrans,im_data


def mkdir(path):
    if not os.path.exists(path):
        os.mkdir(path)


if __name__ == '__main__':
    import matplotlib.pyplot as plt
    imdir_match = "./test/0000000005_image.tif"
    imdir = "./test5/aircraft_229.jpg"
    out_file = './test5_/'
    mkdir(out_file)

    data1 = cv.imread(imdir)
    data2 = cv.imread(imdir_match)
    # data1 = cv.resize(data1, dsize=(1000, 1000))
    com = np.zeros_like(data1)
    for i in range(3):
        f1 = data1[:,:, i]
        f2 = data2[:,:, i]

        hist_m = arrayToHist(f2, 256)
        im_f1 = histMatch(f1, hist_m)

        com[:, :, i] = im_f1
    
    out_full_path = os.path.join(out_file, '201612_match.tif')
    cv.imwrite(out_full_path, com)

图像处理之直方图匹配_直方图匹配