图像处理——直方图匹配

1 原理

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

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

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

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

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

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->z
		std::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中搜索sk
				int 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;
	}