灰度图像直方图均衡化公式及实现

图像的直方图：直方图是图像中像素强度分布的图形表达方式。它统计了每一个强度值所具有的像素个数。

直方图均衡化：是通过拉伸像素强度分布范围来增强图像对比度的一种方法。是图像处理领域中利用图像直方图对对比度进行调整的方法。

均衡化指的是把一个分布(给定的直方图)映射到另一个分布(一个更宽更统一的强度值分布)，所以强度值分布会在整个范围内展开。映射函数应该是一个累积分布函数(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;
}

#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;
}