直方图匹配（规定化）

1. 直方图匹配

在上一篇中，我们使用matlab完成了直方图的均衡化，在这一篇中，我们要实现图像的直方图匹配。与直方图均衡不同的一点是，直方图匹配要得到的是具有特定形状的直方图的图像。

当图像中接近0的像素过多，进行直方图均衡的时候，会导致统计概率变大，直接映射到高灰度，直接导致图像整体变量，达不到原来想要的结果，此时就需要使用直方图匹配，让变化后的图像具有特定的直方图。有时候，我们为了使两幅图像的色调保持一致，也可以采用该方法。

2. 直方图匹配的计算过程：

1. 根据以下公式（直方图均衡时用到的公式）计算输入图像的s = T(r)

2. 根据相同的公式计算模板直方图的 G(z)，p(z)为模板直方图的pdf。

3. 令G(z) = T(r)，求出r→z，当不止一个z满足时，通常选择最小的值。

注意：a. 对于灰度级为L，取值为[0, L-1]的灰度直方图，在计算T(r)和G(z)的过程中要四舍五入。 b. 由于对G(z)和T(r)进行了四舍五入，因此求r对应的z时，有可能有多个z满足要求，一般选择z最小的值。

r→z的求解：r→z表示对于一个r(对于L = 256的灰度直方图，r取0,1,2,...,254,255)，找出满足T(r) = G(z)的z，若不存在，则找出最接近的z（即找出min(abs( T(r) - G(z) ) 对应的z）。前面已经说了，当z有多个时，通常取最小的。（若不考虑只取整数的情况，则当p(r)和p(z)处处不为0时，仅有一个z与r对应）

3. matlab代码如下：

1. 直方图匹配函数如下：

该匹配函数为histmatching, 输入需要进行直方图匹配的图像，模板直方图，输出为直方图匹配之后的图像，变换向量。

当给出模板图像时，需要在调用该函数之前使用imhist(I)作为函数的模板直方图。该函数使用灰度级别为256的灰度图像。

function [image_out, T3] = histmatching(image_in, hist)
% 直方图匹配（规定化）函数
% 输入为需要进行直方图匹配的灰度图像，模板直方图
% 输出为直方图匹配后的灰度图像，进行变换的向量。

% Level为灰度级别
% T1, T2分别为输入图像，模板直方图的均衡化用到的变换向量
% T3为输入图像匹配模板直方图用到的变换向量
Level = 256;
[m, n] = size(image_in);
image_hist = imhist(image_in);
image_out = image_in;
% 求解T1
ac1 = zeros(Level, 1);
T1 = zeros(Level, 1, 'uint8');
ac1(1) = image_hist(1);
for i = 2 : Level
    ac1(i) = ac1(i - 1) + image_hist(i);
end
ac1 = ac1 * (Level - 1);
for i = 1 : 256
    T1(i) = uint8(round((ac1(i)) / (m * n)));
end

% 求解T2
ac2 = zeros(Level, 1);
T2 = zeros(Level, 1, 'uint8');
ac2(1) = hist(1);
for i = 2 : Level
    ac2(i) = ac2(i - 1) + hist(i);
end
ac2 = ac2 * (Level - 1);
hist_sum = sum(hist);
for i = 1 : 256
    T2(i) = uint8(round((ac2(i)) / hist_sum));
end

% 求解T3
% T1映射到T2^(-1)时，若有多个值，选取最小的那个值。
% 产生0 到 255 之间的256个点，即产生0,1,2,...,255的大小为256的数组
temp = zeros(Level, 1, 'uint8');
T3 = T1;
for i = 1 : 256
    for j = 1 : 256
        temp(j) = abs(T1(i) - T2(j));
    end
    [~, B] = min(temp);
    T3(i) = B - 1;
end

% 根据T3转换输入图像的值
for i = 1 : m
    for j = 1 : n
        image_out(i, j) = T3(uint32(image_in(i, j)) + 1);
    end
end
end

2. 测试代码如下：

close all;
clear all;
clc;

image2 = imread('2.png');
image3 = imread('3.png');
hist3 = imhist(image3);
match1 = histeq(image2, hist3);
[match2, T] = histmatching(image2, hist3);
figure;
subplot(2, 4, 1), imshow(image2), title('原图像');
subplot(2, 4, 2), imshow(image3), title('模板图像');
subplot(2, 4, 3), imshow(match2), title('匹配后得到的图像');
subplot(2, 4, 4), imshow(match1), title('调用histeq的结果');
subplot(2, 4, 5), imhist(image2), title('原图像的直方图');
subplot(2, 4, 6), imhist(image3), title('模板图像的直方图');
subplot(2, 4, 7), imhist(match2), title('匹配后的直方图');
subplot(2, 4, 8), imhist(match1), title('调用histeq得到图像对应的直方图');

3. 关于用到的matlab函数的一些说明：

1. [A, B] = min(V);  V为一个向量，函数返回值A表示V中最小的值，B表示该最小值的位置，当有多个最小值时，返回第一个出现的位置。

2. matlab的函数调用中，对于有多个返回值的函数，例如有2个返回值[A, B]，若在之后的过程中不需要使用到所有的值，可以在接受时将相应位置的参数用~代替。例如代码中对于min(V)函数，我们不需要用到其最小值，只需要用到最小值所在的位置，可以用[~,B] = min(V)。若只需要用到前几个返回值，可以省略后面的。例如：[A] = min(V); 返回的是V的最小值（A = min(V) 与 [A] = min(V)等价）。

3. histeq函数(可以在matlab 通过 help histeq 查看该函数的相关说明)

J = histeq(I, hist);  I为输入图像，hist为匹配直方图，可以是用不同方式表示的灰度图像的直方图，例如doule/single类型（横坐标范围0~1），uint8类型（横坐标范围0~255），uint16和int16

J = histeq(I, n); 直方图均衡函数，I为输入图像，n为均衡后的灰度级别，缺省值为64。

[J, T] = histeq(I); T为将灰度图像I中灰度级映射到J中灰度级的灰度变换，值范围为[0, 1]。

另外，还有3个针对索引图像调色板的函数如下：

newmap = histeq(X,map,hgram)
newmap = histeq(X,map)

[newmap,T] = histeq(X,...)

4. 代码运行结果如下：