图像处理之直方图均衡和直方图匹配

1 直方图均衡

直方图均衡是一种将图像中的灰度分布转换成均匀分布，从而增强图像的对比度的图像处理方法。直方图均衡可以将原本偏白或者偏黑的图像转换成对比度符合人眼视觉的图像。

1.1 原理

连续空间
连续空间内的图像灰度r∈[0,L−1]，L表示灰度级r\in[0,L-1]，L表示灰度级r∈[0,L−1]，L表示灰度级，期望的均衡化转换函数为
s=T(r),0≤r≤L−1s=T(r),0\le r\le L-1 s=T(r),0≤r≤L−1
其中TTT为灰度变换函数，sss是经过变换后的灰度值，该变化函数满足两个条件：

T(r)T(r)T(r)在区间[0,L−1][0, L-1][0,L−1]内是严格单调递增函数。保证转换是一一对应的避免出现二义性；
当0≤r≤L−10 \le r\le L - 10≤r≤L−1时，0≤s≤L−10 \le s\le L - 10≤s≤L−1。保证输入输出的灰度范围相同。

我们期望的是将输入的灰度概率分布p(r)p(r)p(r)转换成输出的灰度概率分布p(s)=1L−1p(s)=\frac{1}{L - 1}p(s)=L−11。这里选择下面的变换函数来证明这个函数是可行的：
s=T(r)=(L−1)∫0rpr(w)dw，w为积分的假变量s=T(r)=(L-1)\int_{0}^{r}p_r(w)dw，w为积分的假变量 s=T(r)=(L−1)∫0rpr(w)dw，w为积分的假变量
首先，∫0rpr(w)dw\int_{0}^{r}p_r(w)dw∫0rpr(w)dw的值域为[0,1][0,1][0,1]因此sss的值域为[0,L−1][0,L-1][0,L−1]；另外该函数也是严格单调递增的。基本条件满足。下面我们证明p(s)=1L−1p(s)=\frac{1}{L-1}p(s)=L−11。因为
p(s)=pr(r)∣drds∣=pr(r)1dsdr=pr(r)1(L−1)ddr[∫0rpr(w)dw]=pr(r)1(L−1)pr(r)=1L−1\begin{equation} \begin{aligned} p(s)&=p_r(r)|\frac{dr}{ds}|\\ &=p_r(r)\frac{1}{\frac{ds}{dr}}\\ &=p_r(r)\frac{1}{(L-1)\frac{d}{dr}[\int_{0}^{r}p_r(w)dw]}\\ &=p_r(r)\frac{1}{(L-1)p_r(r)}\\ &=\frac{1}{L-1} \end{aligned} \end{equation} p(s)=pr(r)∣dsdr∣=pr(r)drds1=pr(r)(L−1)drd[∫0rpr(w)dw]1=pr(r)(L−1)pr(r)1=L−11
因此s=T(r)=(L−1)∫0rpr(w)dws=T(r)=(L-1)\int_{0}^{r}p_r(w)dws=T(r)=(L−1)∫0rpr(w)dw可以作为均衡化变换函数。
离散空间
而对于图像处理始终是离散空间，离散空间的灰度rk，k为灰度级，k∈[0,L−1]r_k，k为灰度级，k\in [0,L-1]rk，k为灰度级，k∈[0,L−1]的概率近似为
pr(rk)=nkMNp_r(r_k)=\frac{n_k}{MN} pr(rk)=MNnk
其中k为灰度级，k∈[0,L−1]k为灰度级，k\in [0,L-1]k为灰度级，k∈[0,L−1]；MMM和NNN分别为图像的宽高,nkn_knk为当前灰度级像素点的数量。则利用上面的变换函数得到的sks_ksk为
sk=T(rk)=(L−1)∑j=0kpr(rj)=L−1MN∑j=0knjs_k=T(r_k)=(L-1)\sum_{j=0}^{k}p_r(r_j)=\frac{L-1}{MN}\sum_{j=0}^{k}{n_j} sk=T(rk)=(L−1)j=0∑kpr(rj)=MNL−1j=0∑knj
另外利用上面的变化函数计算出来的灰度值可能是小数，因此需要最近取整。以及对于RGB图像分别对R、G、B三个通道的图像进行变换的话会出现色相问题，因此需要将图像转换成其他颜色空间比如YUV，只对Y做灰度变换。

1.2 实现

实现比较简单，首先逐像素计算每个灰度级的数量，根据改数量计算出每个灰度级的概率，最后利用上面的公式计算出最终的颜色映射关系，最终将根据映射关系将原图中颜色进行映射。

static vector<float> countGrayProp(const Mat &img) {assert(img.channels() == 1);double pixelProp = 1.0 / (img.rows * img.cols);vector<float> ret(GRAPH_GRAY_LAYER_NUM, 0.0);for (int i = 0; i < img.rows; i++) {for (int j = 0; j < img.cols; j++) {ret[static_cast<int>(img.at<uchar>(i, j))] += pixelProp;}}return ret;
}
/** @brief 均衡化灰度图*/
static Mat avgGrayHistogram(const Mat &img) {assert(img.channels() == 1);vector<float> props = countGrayProp(img);vector<float> propSum(GRAPH_GRAY_LAYER_NUM, 0.0);for (int i = 0; i < GRAPH_GRAY_LAYER_NUM;i ++) {if (i == 0) {propSum[i] = props[i];}else {propSum[i] = (props[i] + propSum[i - 1]);}}Mat ret(img.rows, img.cols, CV_8UC1);for (int i = 0; i < img.rows; i++) {for (int j = 0; j < img.cols; j++) {int value = static_cast<int>(img.at<uchar>(i, j));ret.at<uchar>(i, j) = int((GRAPH_GRAY_LAYER_NUM - 1) * propSum[value]);}}return ret;
}Mat avgHistogram(const Mat &img) {if (img.channels() == 1) {return avgGrayHistogram(img);}else if (img.channels() == 3) {Mat yuvImg;cvtColor(img, yuvImg, COLOR_BGR2YUV);std::vector<Mat> yuvImgs;split(yuvImg, yuvImgs);yuvImgs[0] = avgGrayHistogram(yuvImgs[0]);Mat proYUV, ret;merge(yuvImgs, proYUV);cv::cvtColor(proYUV, ret, COLOR_YUV2BGR);return ret;}return img;
}

对于单通道的灰度图效果如下，能够看到均衡化后的灰度直方图中灰度更加均匀：

如果直接对rgb三个通道做均衡化会出现下面的情况（1,2,3,4四张图中1和3为处理前的图，34为处理后的）

首先将图像转换到其他颜色空间比如YUV，只对Y分量进行均衡化就不会有明显的的色差和奇异点。

2 直方图匹配

2.1 原理

利用直方图均衡对图像进行灰度级调整能够快速高效的将图像调整的符合人眼的感受。但是该方式过于简单粗暴，而且无法自定义目标的直方图，对于某些场景并不是很有效。直方图规定或者直方图匹配就是为了解决直方图均衡无法指定输出图像的直方图概率分布的方法，该方法能够将输入图像的概率分布映射到指定的输出直方图概率分布。
连续空间
任务目标：输入图像的概率密度分布pr(r)p_r(r)pr(r)，输出图像概率密度分布为pz(z)p_z(z)pz(z)。
在直方图均衡任务中，sss和rrr的映射关系为s=T(r)=(L−1)∫0rpr(w)dws=T(r)=(L-1)\int_{0}^{r}p_r(w)dws=T(r)=(L−1)∫0rpr(w)dw。同样，对于pz(z)p_z(z)pz(z)进行直方图均衡，则G(z)=(L−1)∫0zpz(w)dwG(z)=(L-1)\int_{0}^{z}p_z(w)dwG(z)=(L−1)∫0zpz(w)dw，其中G(z)=T(r)G(z)=T(r)G(z)=T(r)，即
T(r)=(L−1)∫0rpr(w)dw=G(z)=(L−1)∫0zpz(t)dt\begin{equation} \begin{aligned} T(r)=(L-1)\int_{0}^{r}p_r(w)dw=G(z)=(L-1)\int_{0}^{z}p_z(t)dt \end{aligned} \end{equation} T(r)=(L−1)∫0rpr(w)dw=G(z)=(L−1)∫0zpz(t)dt
则：
z=G−1(T(r))z=G^{-1}(T(r)) z=G−1(T(r))
因此，将一个概率密度为pr(r)p_r(r)pr(r)的图像规定化为pz(z)p_z(z)pz(z)的步骤为：

计算输入图像各个灰度级的概率分布pr(r)p_r(r)pr(r)；
计算输入图像的直方图均衡映射函数；
根据指定的pz(z)p_z(z)pz(z)计算映射函数G(z)G(z)G(z)；
求G(z)G(z)G(z)的反变换函数G−1(T(r))G^{-1}(T(r))G−1(T(r))，对图像中每一个像素进行处理得到目标灰度图像。

离散空间
在离散空间，pr(r)p_r(r)pr(r)经过直方图均衡后的输出为sk=L−1MN∑j=0knjs_k=\frac{L-1}{MN}\sum_{j=0}^{k}{n_j}sk=MNL−1∑j=0knj，而G(z)=(L−1)∑i=0qpz(zj)G(z)=(L-1)\sum_{i=0}^{q}p_z(z_j)G(z)=(L−1)∑i=0qpz(zj)。根据连续空间的推论，输出的值为zq=G−1(sk)z_q=G^{-1}(s_k)zq=G−1(sk)。
从上面的推论中能够看出我们只需要找出逆变换G−1(sk)G^{-1}(s_k)G−1(sk)即可，但是实际上由于我们是用的灰度级都是有限的，比如256个灰度级，因此往往只需要进行穷举即可。

计算输入图像各个灰度级的概率分布pr(r)p_r(r)pr(r)；
对输入图像进行直方图均衡得到输入图像的r→skr\rightarrow s_kr→sk的映射；
计算目标概率密度分布的直方图均衡映射关系q→sqq\rightarrow s_qq→sq；
在步骤3中的映射关系中寻找最接近的G(z)G(z)G(z)，如果存在多个值时选择最小的，即得到sk→qs_k\rightarrow qsk→q的映射。

在实现中可能发现G(z)G(z)G(z)并不是严格单调递增的，不符合直方图均衡时提到的条件，但是这在离散空间并不是一个严重的问题。

2.2 实现

实现上比较简单，主要就是搜索内容。

 static Mat matchGrayHistogram(const cv::Mat &img, const std::vector<float>& targetProp) {assert(img.channels() == 1);vector<float> props = countGrayProp(img);vector<float> propSum(GRAPH_GRAY_LAYER_NUM, 0.0), targPropSum(GRAPH_GRAY_LAYER_NUM, 0.0);for (int i = 0; i < GRAPH_GRAY_LAYER_NUM; i++) {if (i == 0) {propSum[i] = props[i];targPropSum[i] = targetProp[i];}else {propSum[i] = (props[i] + propSum[i - 1]);targPropSum[i] = (targetProp[i] + targPropSum[i - 1]);}}//将输入和输出的标签映射整数化for (int i = 0; i < GRAPH_GRAY_LAYER_NUM; i++) {int layer = GRAPH_GRAY_LAYER_NUM - 1;propSum[i] = std::min<float>(GRAPH_GRAY_LAYER_NUM - 1, static_cast<int>(layer * propSum[i] + 0.5));targPropSum[i] = std::min<float>(GRAPH_GRAY_LAYER_NUM - 1, static_cast<int>(layer * targPropSum[i] + 0.5));}//映射关系搜索，propSum中存储r->sk，targPropSum中存储z->sk 目标搜索sk->zstd::vector<int> skzMap(GRAPH_GRAY_LAYER_NUM, -1);for (int i = 0; i < GRAPH_GRAY_LAYER_NUM; i++) {int sk = propSum[i];if (skzMap[sk] == -1) {//使用二分查找法在targetPropSum中搜索skint left = 0, right = targPropSum.size() - 1;int mid = 0;while (left < right) {mid = left + (right - left) / 2.0;if (targPropSum[mid] == sk) {while (mid > 1 && targPropSum[mid] == targPropSum[mid - 1]) { mid--; }skzMap[sk] = mid;break;}else if (targPropSum[mid] > sk) {right = mid - 1;}else {left = mid + 1;}}if (targPropSum[mid] != sk) {if (left != 0 && abs(sk - targPropSum[left - 1]) < abs(sk - targPropSum[left])) {skzMap[sk] = left - 1;}else {skzMap[sk] = left;}}else {printf("");}}}Mat ret(img.rows, img.cols, CV_8UC1);for (int i = 0; i < img.rows; i++) {for (int j = 0; j < img.cols; j++) {int value = static_cast<int>(img.at<uchar>(i, j));ret.at<uchar>(i, j) = skzMap[propSum[value]];}}return ret;}Mat GrayTransform::matchHistogram(const cv::Mat &img, const std::vector<float>& targetProp) {assert(targetProp.size() == 256);if (img.channels() == 1) {return matchGrayHistogram(img, targetProp);}else if (img.channels() == 3) {Mat yuvImg;cvtColor(img, yuvImg, COLOR_BGR2YUV);std::vector<Mat> yuvImgs;split(yuvImg, yuvImgs);yuvImgs[0] = matchGrayHistogram(yuvImgs[0], targetProp);Mat proYUV, ret;merge(yuvImgs, proYUV);cv::cvtColor(proYUV, ret, COLOR_YUV2BGR);return ret;}return img;}