直方图均衡化（histogram equalization）

一. 原理

直方图均衡化是想要将聚集在某一区间内分布的灰度值，变为均匀的在所有区间内分布。

为了达到这一目的，我们需要找出一个函数T，将r（原图像灰度）映射到s（新图像灰度）上。同时，由于不想将图像反转，我们需要保证函数单调不减（若需要逆运算，则要严格单调递增）
s = T ( r ) s=T(r) s=T(r)

设 p r ( r ) p_r(r) pr​(r)为r的概率分布函数， p s ( s ) p_s(s) ps​(s)为s的概率分布函数，则两者关系如下
p s ( s ) = p r ( r ) ∣ d r d s ∣ p_s(s)=p_r(r)|\frac{dr}{ds}| ps​(s)=pr​(r)∣dsdr​∣
也就是说， T ( r ) T(r) T(r)这个函数现在要满足条件：

1. 单调递增
2. 值域范围为[0,L-1]
3. p s ( s ) p_s(s) ps​(s)为一个常函数（表示每一个灰度值的概率都相等）

现在给出一个函数（至于为什么是这个函数，我也不知道啊，哪位看官可以帮忙推导一下QAQ）：
s = T ( r ) = ( L − 1 ) ∫ 0 r p r ( w ) d w s=T(r)=(L-1)\int_{0}^{r}{p_r(w)}dw s=T(r)=(L−1)∫0r​pr​(w)dw
我们将逐步证明，以上三点成立

1. 由于此函数为积分所得，因此函数单调递增

2. 由于 p r ( r ) p_r(r) pr​(r)为概率密度函数，因此，积分区间为[0,1]，乘以常数，区间为[0,L-1]

3. 推导过程：
   d s = d T ( r ) = d ( ( L − 1 ) ∫ 0 r p r ( w ) d w ) = ( L − 1 ) p r ( r ) d r ds=dT(r)=d((L-1)\int_{0}^{r}{p_r(w)}dw)=(L-1)p_r(r)dr ds=dT(r)=d((L−1)∫0r​pr​(w)dw)=(L−1)pr​(r)dr

   p s ( s ) = p r ( r ) ∣ d r d s ∣ = p r ( r ) ∣ 1 ( L − 1 ) p r ( r ) ∣ = 1 L − 1 p_s(s)=p_r(r)|\frac{dr}{ds}|=p_r(r)|\frac{1}{(L-1)p_r(r)}|=\frac{1}{L-1} ps​(s)=pr​(r)∣dsdr​∣=pr​(r)∣(L−1)pr​(r)1​∣=L−11​

   此时 p s ( s ) p_s(s) ps​(s)为一个均匀密度函数，符合

对于离散值，计数每个不大于r的概率和，就是这个r对应的s
s k = T ( r k ) = ( L − 1 ) ∑ j = 0 k p r ( r j ) s_k=T(r_k)=(L-1)\sum_{j=0}^{k}{p_r(r_j)} sk​=T(rk​)=(L−1)j=0∑k​pr​(rj​)

二. 代码

img = imread("14.jpeg");
r = img(:,:,1);
g = img(:,:,2);
b = img(:,:,3);

p = stat(r);
map = cal(p);
img_equa = my_hist_equation(r,map);

% 画图
%原图
subplot(331)
imshow(r);
title('原图');
subplot(332)
imhist(r);
title('原图直方图');

%自己的直方图均衡化图
subplot(334)
imshow(img_equa);
title('my直方图均衡化图');
subplot(335)
imhist(img_equa);
title('my直方图均衡化直方图');
subplot(336)
plot((0:255)/255,map/255);
title('my函数');

%标准的直方图均衡化图
[std_img,T] = histeq(r);
subplot(337)
imshow(std_img);
title('标准直方图均衡化图');
subplot(338)
imhist(std_img);
title('标准直方图均衡化图直方图');
subplot(339)
plot((0:255)/255,T);
title('std函数');

figure
plot((0:255)/255,map/255);
figure
plot((0:255)/255,T);

%计算图像的分布
function p=stat(img)
    [r,c] = size(img);
    p = linspace(0,0,256);
    for i=1:1:r
        for j=1:1:c
            p(img(i,j)+1)= p(img(i,j)+1)+1;
        end
    end
    p = p/(r*c);
end

%计算r与s的对应关系，矩阵坐标为r，值为s(动态规划)
function s=cal(probality)
    s=probality;
    l = length(probality);
    for i=2:l
        s(i) = s(i-1)+s(i);
    end
    s = round(s*255);
end

%将图片根据映射函数，映射
function img_equa=my_hist_equation(img,map)
    [r,c] = size(img);
    img_equa = zeros(r,c);
    for i=1:r
        for j=1:c
            img_equa(i,j) = map(img(i,j)+1);
        end
    end
    img_equa = uint8(img_equa);
end

放大我的映射函数图片为：

matlab自建的映射函数为：

对比两个映射函数可以看出来，自己实现的为比较连续的，而matlab自建的为阶段性的。这也解释了为什么两个直方图均衡化的后的直方图图像不同。因为对于matlab自建的的函数，将多个输入值映射为同一个输出值了。