图像处理：基于灰度信息的直方图均衡化

图像处理课老师介绍的方法（Region-Based Weighted Histogram and Its Application to Image Contrast Enhancemen)，自己在matlab下实现了一下：

• 在对大背景的图像直接直方图均衡化，会由于背景信息权重过大，影响结果：

I=imread('Xray.jpg');
J=histeq(I);
subplot(2,2,1);imshow(I);title('原图');
subplot(2,2,2);imshow(J);title('均衡化');
subplot(2,2,3);imhist(I,256);
subplot(2,2,4);imhist(J,256);

结果：

文章对图片进行边缘检测，采用灰度梯度大的像素点作为均衡化依据，这样便可排除背景干扰。
sobel梯度边缘检测滤波器：

[R,C] = size(img);   % 获得图像的高度和宽度
U = double(img);
map = img;%uint8格式，不必结尾再转换;存储跳跃值图
for i = 2:R - 1   %类似sobel边缘检测
    for j = 2:C - 1
        Gx = (U(i+1,j-1) + 2*U(i+1,j) + U(i+1,j+1)) - (U(i-1,j-1) + 2*U(i-1,j) + U(i-1,j+1));
        Gy = (U(i-1,j+1) + 2*U(i,j+1) + U(i+1,j+1)) - (U(i-1,j-1) + 2*U(i,j-1) + U(i+1,j-1));
        map(i,j) = sqrt(Gx^2 + Gy^2); %四周一圈没变，沿用img的值了
    end
end 

跳跃值的图像map：

然后将所得到的灰度梯度排序，从小到大分为5部分。取灰度梯度最大的五分之一在原图上对应的像素数据，依据此做直方图均衡化：

A=reshape(map,1,R*C);%矩阵排成行向量
A=sort(A);%排序
a=A(round(R*C*4/5));%取最高五分之一阈值

最大的1/5的值是50，所以，统计跳跃值>50的点对应的像素值的直方图

%% 直方图
cnt = zeros(1, 256);
sum=0;
for i = 1 : R
    for j = 1 : C
        if map(i,j)>a %统计1/5
            cnt(1, I(i, j) + 1) = cnt(1, I(i, j) + 1) + 1;
            sum=sum+1;
        end
    end
end

统计得到的直方图cnt：

比原始直方图要占灰度值均匀许多，用该直方图去均衡化

%% 作均衡化
f = zeros(1, 256);
f = double(f); cnt = double(cnt);
% 统计每个像素值出现的概率， 得到概率直方图
for i = 1 : 256
    f(1, i) = cnt(1, i) / sum;
end
% 求累计概率，得到累计直方图
for i = 2 : 256
    f(1, i) = f(1, i - 1) + f(1, i);
end
% 用f数组实现像素值[0, 255]的映射。 
for i = 1 : 256
    f(1, i) = f(1, i) * 255;
end

直方图均衡化函数f：

% 完成每个像素点的映射
I = double(I);
for i = 1 : R
    for j = 1 : C
        I(i, j) = f(1, I(i, j) + 1);
    end
end
img = uint8(I);%格式转换

结果：

由于这个自己构造的边缘检测方法比较粗糙，有些边缘点没滤出来，掉到了灰度梯度小的四部分中去了。
所以，对小的那部分按照其存在边缘点（采用matlab自带的sobel或canny滤波器获得）的个数比重，进行加权。区别不大有空再码（懒了）。边缘检测map = edge(I,'canny');

• 不过好像直接采用现成的边缘检测滤波器也没什么问题，效果几乎完全一致，就不放图了，代码如下
  比如采用canny边缘检测

I = imread('...');  
% img = rgb2gray(img);        %转化为灰度图
map = edge(I,'canny');

%% 得到灰度直方图
[R,C] = size(I);
cnt = zeros(1, 256);
sum=0;
for i = 1 : R
    for j = 1 : C
        if map(i,j)==1 %统计1/5
            cnt(1, I(i, j) + 1) = cnt(1, I(i, j) + 1) + 1;
            sum=sum+1;
        end
    end
end

原文连接Improving histogram-based image contrast enhancement using gray-level information histogram with application to X-ray images