# 直方图均衡算法(Histogram Equalized)

Lab1: Histogram Equalization

### 1. 实验环境(C++)

• 操作系统版本 MacOS Catalina 10.15
• OpenCV4.0 (imgcodecs | core | highgui | imgproc)
• Cmake-3.14
• Clang-1100.0.33.8

### 2. 实验步骤

1. Calculate the histogram H for src

2. Normalize the histogram.

	std::array<double, 256> calNormalizedHist(cv::Mat& source) {
std::array<double, 256> acc{0};
// Calculate the histogram H for src
for(int i = 0; i < source.rows; i++)
for (int j = 0; j < source.cols; j++)
acc[ source.ptr<uchar>(i)[j]] ++;
// Normalize the histogram.
for(int i = 0; i < acc.size(); i++)
acc[i] /= source.rows * source.cols;
return acc;
}

3. Compute the integral of the histogram: H'

4. Transform the image using H′ as a look-up table: $$𝚍𝚜𝚝(x,y)=H′(𝚜𝚛𝚌(x,y)) * 255$$

	void equalizeHist(cv::Mat& source, cv::Mat& result) {
source.copyTo(result);
auto hist = calNormalizedHist(source);
// Compute the integral of the histogram: H'
for(int i = 1; i < 256; i++) hist[i] += hist[i-1];
// Transform the image using H′ as a look-up table: 𝚍𝚜𝚝(x,y)=H′(𝚜𝚛𝚌(x,y)) * 255
for(int i = 0; i < source.rows; i++)
for (int j = 0; j < source.cols; j++)
result.ptr<uchar>(i)[j] = hist[source.ptr<uchar>(i)[j]] * 255;
}

5. Test Code:

	cv::Mat cv, lab, m = cv::imread("../data/lena.png", cv::IMREAD_GRAYSCALE);
cv::equalizeHist(m, cv);
equalizeHist(m, lab);

cv::imwrite("../out/m.png", m);
cv::imwrite("../out/cv.png", cv);
cv::imwrite("../out/lab.png", lab);

cv::imwrite("../out/cv-m.png", abs(cv-m));
cv::imwrite("../out/lab-m.png", abs(lab-m));
cv::imwrite("../out/lab-cv.png", abs(lab-cv));

cv::imwrite("../out/hist.m.png", drawHist(calNormalizedHist(m)));
cv::imwrite("../out/hist.cv.png", drawHist(calNormalizedHist(cv)));
cv::imwrite("../out/hist.lab.png", drawHist(calNormalizedHist(lab)));

cv::imwrite("../out/acchist.m.png", drawHist(calNormalizedHist(m), true));
cv::imwrite("../out/acchist.cv.png", drawHist(calNormalizedHist(cv), true));
cv::imwrite("../out/acchist.lab.png", drawHist(calNormalizedHist(lab),true));


### 3. 实验结果

image: Origin

image: OpenCV

image: Result

dif: Origin&OpenCV

dif: Result&Origin

dif: Result&OpenCV

Hist: Origin

Hist: OpenCV

Hist: Result

AccHist: Origin

AccHist: OpenCV

AccHist: Result
###4. 实验结果分析
• 明显本文实现的算法与OpenCV实现的高度一致(不考虑指令集优化: SIMD|SEE4 etc.)

• 由Histogram与Accumulate Histogram上看, 连续变量的均匀分布意味着概率分布函数满足线性分布, 与理论推导得出的性质一致

• 从直方图直观观察和算法分析, 易知基于直接灰度映射(点对点映射)的经典直方图均衡算法会导致灰阶的减少.因而加剧色带(Banding)现象:

• Lookup-Table $$H’$$在定义上不是单射

所以这个意义上, 直方图的均衡是一个可以把相邻柱子合并的挪动柱子的游戏~

### 5 Attachment

• 怎么得到均衡后的直方图 ?

• 比较常见的方法: 图像均衡后再统计一次直方图

• 比较少用的方法: 直接均衡原直方图

	std::array<double, 256> equalizeHist(std::array<double, 256>&& hist) {
std::array<double, 256> hist_eq{0}, hist_acc{0};
hist_acc[0] = hist[0];
// Compute the integral of the histogram: H'
for(int i = 1; i < 256; i++) hist_acc[i] = hist_acc[i-1] + hist[i];

for(int i = 0; i < 256; i++)
hist_eq[hist_acc[i]*255] += hist[i];

return hist_eq;
}

• 先计算均衡后的直方图再均衡图像对吗?

• 不对, 均衡是一种方法, 可以直接均衡图像然后计算新的直方图, 也可以直接均衡直方图, 两直方图一致
• 直方图均衡的均衡指什么?

• 均衡一种概率变换, 可以将线性的目标概率分布函数拓展到任意的分布函数:
​ Transformation between two distribution funtion.

#### 5.* 下面给出指定累计直方图的均衡结果的部分例子

• $y=x^2$

AccHist

Image

Hist:

Hist: Directly
- $$y=x^5$$

AccHist

Image

Hist:

Hist: Directly
- $$y=sin(x)$$

AccHist

Image

Hist:

Hist: Directly
- $$y=\sqrt{2x-x^2}$$

AccHist

Image

Hist:

Hist: Directly
- $$y=e^{x-1}$$

AccHist

Image

Hist:

Hist: Directly
#include <iostream>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <array>

cv::Mat drawHist(std::array<double, 256>&& acc, bool accumulate = false) {
if(accumulate) for(int i = 1; i < acc.size(); i++) acc[i] += acc[i-1];

double max_v = *std::max_element(acc.begin(), acc.end());

cv::Mat visual = cv::Mat::zeros(256, 256, CV_8UC3);
for(int i = 0; i < 256; i++)
cv::line(visual, cv::Point(i, 256.0*(1-acc[i]/max_v)), cv::Point(i, 256), 255);

return visual;
}

std::array<double, 256> calNormalizedHist(cv::Mat& source) {
std::array<double, 256> acc{0};

// Calculate the histogram H for src
for(int i = 0; i < source.rows; i++)
for (int j = 0; j < source.cols; j++)
acc[ source.ptr<uchar>(i)[j]] ++;

// Normalize the histogram.
for(int i = 0; i < acc.size(); i++)
acc[i] /= source.rows * source.cols;

return acc;
}

std::array<double, 256> equalizeHist(std::array<double, 256>&& hist,
std::function<double(double)> && f=[](double x){return x;}) {
std::array<double, 256> hist_eq{0};

std::array<int, 256> H{0};

double sum = 0;
for(int i = 0, j = 0; i < 256; i++)
for(;sum < f(i/255.0); j++) {
sum += hist[j]; H[j] = i;
}

for(int i = 0; i < 256; i++)
hist_eq[H[i]] += hist[i];

return hist_eq;
}

// Only for single channel
void equalizeHist(cv::Mat& source, cv::Mat& result,
std::function<double(double)> && f=[](double x){return x;}) {
source.copyTo(result);

auto hist = calNormalizedHist(source);

std::array<int, 256> H{0};

// Compute the integral of the histogram: H'
double sum = 0;
for(int i = 0, j = 0; i < 256; i++)
for(;sum < f(i/255.0); j++) {
sum += hist[j]; H[j] = i;
}

// directly gray mapping
// Transform the image using H′ as a look-up table: 𝚍𝚜𝚝(x,y)=H′(𝚜𝚛𝚌(x,y)) * 255
for(int i = 0; i < source.rows; i++)
for (int j = 0; j < source.cols; j++)
result.ptr<uchar>(i)[j] = H[source.ptr<uchar>(i)[j]];
}

int main() {

cv::Mat cv, lab;

cv::equalizeHist(m, cv);
equalizeHist(m, lab);

cv::imwrite("../out/m.bmp", m);
cv::imwrite("../out/cv.bmp", cv);
cv::imwrite("../out/lab.bmp", lab);

cv::imwrite("../out/cv-m.bmp", abs(cv-m));
cv::imwrite("../out/lab-m.bmp", abs(lab-m));
cv::imwrite("../out/lab-cv.bmp", abs(lab-cv));

cv::imwrite("../out/hist.m.bmp", drawHist(calNormalizedHist(m)));
cv::imwrite("../out/hist.cv.bmp", drawHist(calNormalizedHist(cv)));
cv::imwrite("../out/hist.lab.bmp", drawHist(calNormalizedHist(lab)));
cv::imwrite("../out/hist.lab_direct.bmp", drawHist(equalizeHist(calNormalizedHist(m))));

cv::imwrite("../out/acchist.m.bmp", drawHist(calNormalizedHist(m), true));
cv::imwrite("../out/acchist.cv.bmp", drawHist(calNormalizedHist(cv), true));
cv::imwrite("../out/acchist.lab.bmp", drawHist(calNormalizedHist(lab),true));

{
auto f = [&](std::string f_name, std::function<double(double)>&& f) {
cv::Mat m_f;
equalizeHist(m, m_f, std::move(f));
cv::imwrite("../out/lab."+f_name+".bmp", m_f);
cv::imwrite("../out/hist."+f_name+".bmp", drawHist(calNormalizedHist(m_f)));
cv::imwrite("../out/acchist."+f_name+".bmp", drawHist(calNormalizedHist(m_f), true));
cv::imwrite("../out/hist."+f_name+"_direct.bmp",drawHist(equalizeHist(calNormalizedHist(m), std::move(f))))    ;
};

f("x^2", [](double x){return x*x;});
f("sin(x)", [](double x){return sin(3.1415926 * 0.5 * x);});
f("sqrt(x*(2-x))", [](double x){return sqrt(x*(2-x));});
f("exp(x-1)", [](double x){return std::exp(x-1);});
f("x^5", [](double x){return std::pow(x,5);});
}

//	cv::imshow("lena", m);
//	cv::imshow("lena: equalizeHist", cv);
//	cv::imshow("lena: equalizeHist zlb", lab);
//
//
//	cv::imshow("origin", drawHist(calNormalizedHist(m)));
//	cv::imshow("cvHist", drawHist(calNormalizedHist(cv), true));
//	cv::imshow("zlbHist", drawHist(calNormalizedHist(lab), true));
//	cv::waitKey();
return 0;
}

posted @ 2019-10-31 14:14  xiconxi  阅读(2710)  评论(0编辑  收藏  举报