灰度影像直方圖均衡化公式及實現

影像的直方圖：直方圖是影像中畫素強度分佈的圖形表達方式。它統計了每一個強度值所具有的畫素個數。

直方圖均衡化：是通過拉伸畫素強度分佈範圍來增強影像對比度的一種方法。是影像處理領域中利用影像直方圖對對比度進行調整的方法。

均衡化指的是把一個分佈(給定的直方圖)對映到另一個分佈(一個更寬更統一的強度值分佈)，所以強度值分佈會在整個範圍內展開。對映函式應該是一個累積分佈函式(cumulative distribution function(cdf))。

直方圖均衡化是通過調整影像的灰階分佈，使得在0~255灰階上的分佈更加均衡，提高了影像的對比度，達到改善影像主觀視覺效果的目的。對比度較低的影像適合使用直方圖均衡化方法來增強影像細節。

這種方法通常用來增加許多影像的全域性對比度，尤其是當影像的有用資料的對比度相當接近的時候。通過這種方法，亮度可以更好地在直方圖上分佈。這樣就可以用於增強區域性的對比度而不影響整體的對比度，直方圖均衡化通過有效地擴充套件常用的亮度來實現這種功能。這種方法對於背景和前景都太亮或者太暗的影像非常有用，這種方法尤其是可以帶來X光影像中更好的骨骼結構顯示以及曝光過度或者曝光不足照片中更好的細節。這種方法的一個主要優勢是它是一個相當直觀的技術並且是可逆操作，如果已知均衡化函式，那麼就可以恢復原始的直方圖，並且計算量也不大。這種方法的一個缺點是它對處理的資料不加選擇，它可能會增加背景噪聲的對比度並且降低有用訊號的對比度。

灰度直方圖均衡化演算法實現步驟：

(1)、統計原始影像各灰度級的畫素數目n_i,0≤i<L, L是影像中所有的灰度數(通常為256)；

(2)、影像中灰度為i的畫素的出現概率是:p_x(i)=p(x=i)=n_i/n,n是影像中所有的畫素數，p_x(i)實際上是畫素值為i的影像的直方圖，歸一化到[0, 1]；

(3)、p_x的累積分佈函式,是影像的累計歸一化直方圖：

(4)、直方圖均衡化計算公式, cdf_min為累積分佈函式最小值，M和N分別代表了影像的長寬畫素個數，而L則是灰度級數(如影像為8位深度，則灰度級別共有2^8=256級數，這也是最常見的灰度級數)，v為原始影像中為v的畫素值：

彩色影像直方圖均衡化：上面描述了灰度影像上使用直方圖均衡化的方法，但是通過將這種方法分別用於影像RGB顏色值的紅色、綠色和藍色分量，從而也可以對彩色影像進行處理。實際上，對彩色分量rgb分別做均衡化，會產生奇異的點，影像不和諧。一般採用的是用yuv空間進行亮度的均衡即可。

以上內容主要參考：維基百科

以下是分別採用C++和OpenCV實現的code，從執行結果可知，C++和OpenCV的結果完全一致：

histogram_equalization.cpp:

#include "funset.hpp"
#include <chrono>
#include <vector>
#include <algorithm>
#include "common.hpp"

int histogram_equalization_cpu(const unsigned char* src, int width, int height, unsigned char* dst, float* elapsed_time)
{
	//TIME_START_CPU

	const int hist_sz{ 256 };
	std::vector<int> hist(hist_sz, 0), lut(hist_sz, 0);
	for (int y = 0; y < height; ++y) {
		for (int x = 0; x < width; ++x) {
			++hist[src[y * width + x]];
		}
	}

	int i{ 0 };
	while (!hist[i]) ++i;

	int total{ width * height };
	if (hist[i] == total) {
		unsigned char* p = dst;
		std::for_each(p, p + total, [i](unsigned char& value) { value = i; });
		return 0;
	}

	float scale = (hist_sz - 1.f) / (total - hist[i]);
	int sum = 0;

	for (lut[i++] = 0; i < hist_sz; ++i) {
		sum += hist[i];
		lut[i] = static_cast<unsigned char>(sum * scale + 0.5f);
	}

	for (int y = 0; y < height; ++y) {
		for (int x = 0; x < width; ++x) {
			dst[y * width + x] = static_cast<unsigned char>(lut[src[y * width + x]]);
		}
	}

	//TIME_END_CPU

	return 0;
}

main.cpp:

#include "funset.hpp"
#include <random>
#include <iostream>
#include <vector>
#include <memory>
#include <string>
#include <algorithm>
#include "common.hpp"

int test_image_process_histogram_equalization()
{
	const std::string image_name{ "E:/GitCode/CUDA_Test/test_data/images/lena.png" };
	cv::Mat mat = cv::imread(image_name, 0);
	CHECK(mat.data);

	const int width{ mat.cols/*1513*/ }, height{ mat.rows/*1473*/ };
	cv::resize(mat, mat, cv::Size(width, height));

	std::unique_ptr<unsigned char[]> data1(new unsigned char[width * height]), data2(new unsigned char[width * height]);
	float elapsed_time1{ 0.f }, elapsed_time2{ 0.f }; // milliseconds

	CHECK(histogram_equalization_cpu(mat.data, width, height, data1.get(), &elapsed_time1) == 0);
	//CHECK(histogram_equalization_gpu(mat.data, width, height, data2.get(), &elapsed_time2) == 0);

	//fprintf(stdout, "image histogram equalization: cpu run time: %f ms, gpu run time: %f ms\n", elapsed_time1, elapsed_time2);

	cv::Mat dst;
	cv::equalizeHist(mat, dst);
	cv::imwrite("E:/GitCode/CUDA_Test/test_data/images/histogram_equalization.png", dst);

	CHECK(compare_result(data1.get(), dst.data, width*height) == 0);
	//CHECK(compare_result(data1.get(), data2.get(), width*height) == 0);

	save_image(mat, dst, width, height/2, "E:/GitCode/CUDA_Test/test_data/images/histogram_equalization_result.png");

	return 0;
}

執行結果：