【数字图像处理实验】RGB与HSI互转、灰度直方图绘制、直方图均衡化、直方图规定化的MATLAB实现与Python实现

文中附上的代码可能有缺失，建议通过该 GitHub仓库访问代码。

1 原理

1.1 色彩模型转换原理

RGB转HSI利用以下公式：

HSI转RGB利用以下公式：
（1）H在[0, 2/3π]之间：

（2）H在[2/3π, 4/3π]之间：

（3）H在[4/3π,5/3π ]之间：

因此，遍历图像中每一个像素（或者直接使用矩阵运算），对每个像素点按公式进行计算即可。

1.2 灰度直方图绘制原理

• 灰度直方图是灰度级的函数，它表示图像中具有每种灰度级的像素的个数，反映图像中每种灰度出现的频率。
• 灰度直方图的横坐标是灰度级，纵坐标是该灰度级出现的频率，它是图像基本的统计特征。

因此，遍历图像中每一个像素，记录每个像素的灰度值（彩色图像记录转换为HSI色彩空间后的亮度通道值），累加统计。而后以灰度级为横坐标，灰度级出现频次或频率为纵坐标绘图即可。

1.3 直方图均衡化原理

直方图均衡化就是把原始图像的直方图变换为均匀分布的形式，以增加像素灰度值的动态范围，增加图像整体的对比度效果。过程可以归纳如下（参考自王科平《数字图像处理 MATLAB版》）：

1.4 直方图规定化原理

直方图规定化可以有选择性地增强某个灰度值范围内的对比度，或者使图像灰度值的分布满足特定的要求。主要有三个步骤：

可以认为是，将原始图像和需要的直方图均进行直方图均衡化后，两个均衡化后的结果对应上的部分，将其均衡化前的灰度值对应，形成原始图像和需要的直方图间的映射。用此映射去处理图像中的每一个像素即可。

2 实现源代码

2.1 MATLAB实现

2.1.1 RGB转HSI

function hsi=rgb2hsi(img)
% e.g.
%   img=imread('color.jpg');
%   hsi=rgb2hsi(img);
%   rgb=hsi2rgb(img);
rows=length(img(:,1));
cols=length(img(1,:)) / length(img(1,1,:));
img=im2double(img);
hsi=img;
for i=1:rows
    for j=1:cols
        % get r g b
        r=img(i,j,1);
        g=img(i,j,2);
        b=img(i,j,3);
        % cal h
        num=0.5*((r-g)+(r-b));
        den=sqrt(power(r-g,2)+(r-b)*(g-b));
        theta=acos(num/(den+eps));
        if b<=g
            hsi(i,j,1)=theta;
        else
            hsi(i,j,1)=2*pi-theta;
        end
        % cal s
        if r+g+b==0
            hsi(i,j,2)=1-3*min(min(r,g),b)/(r+g+b+eps);
        else
            hsi(i,j,2)=1-3*min(min(r,g),b)/(r+g+b);
        end
        if hsi(i,j,2)==0
            hsi(i,j,1)=0;
        end
        % cal i
        hsi(i,j,3)=1/3*(r+g+b);
    end
end
hsi=im2uint8(hsi);

2.1.2 HSI转RGB

function rgb=hsi2rgb(img)
% e.g.
%   img=imread('color.jpg');
%   hsi=rgb2hsi(img);
%   rgb=hsi2rgb(img);
rows=length(img(:,1));
cols=length(img(1,:)) / length(img(1,1,:));
img=im2double(img);
rgb=img;
for k=1:rows
    for j=1:cols
        % get h s i
        h=img(k,j,1);
        s=img(k,j,2);
        i=img(k,j,3);
        % set r g b
        r=0;
        g=0;
        b=0;
        if 0 <=h<=2/3*pi
            b=i*(1-s);
            r=i*(1+s*cos(h)/cos(pi/3-h));
            g=3*i-(b+r);
        elseif 2/3*pi<h<=4/3*pi
            r=i*(1-s);
            g=i*(1+s*cos(h-2/3*pi)/cos(pi-h));
            b=3*i-(r+g);
        elseif 4/3*pi<h<=5/3*pi
            g=i*(1-s);
            b=i*(1+s*cos(h-4/3*pi)/cos(5/3*pi-h));
            r=3*i-(g+b);
        end
        rgb(k,j,1)=r;
        rgb(k,j,2)=g;
        rgb(k,j,3)=b;
    end
end
rgb=im2uint8(rgb);

2.1.3 绘制灰度直方图

function frequency=drawgrayscalehistogram(img,filename)
% e.g.
%   img=imread('color.jpg');
%   filename='color';
%   frequency=drawgrayscalehistogram(img,filename);
img=im2uint8(img);
rows=length(img(:,1));
cols=length(img(1,:))/length(img(1,1,:));
frequency=zeros(1,256);
for i=1:rows
    for j=1:cols
        % Matlab number starts from 1
        temp_grayscale=abs(img(i,j,3))+1;
        frequency(temp_grayscale)=frequency(temp_grayscale)+1;
    end
end
ind=0:255;
h=figure;
bar(ind,frequency);
title([filename,' grayscale histogram']);
xlabel('grayscale');
ylabel('frequency');
saveas(h, [filename,'_grayscale_histogram'],'png');
close(h);

2.1.4 直方图均衡化

function histogramequalization(img,filename)
% e.g.
%   color=imread('color.jpg');
%   histogramequalization(color,'color');
img=im2uint8(img);
rows=length(img(:,1));
cols=length(img(1,:))/length(img(1,1,:));
iscolorful=false;
if img(1,1,1)~=img(1,1,2)
    iscolorful=true;
    img=rgb2hsi(img);
end
% (1)$n_j$
n_j=drawgrayscalehistogram(img,filename);
% (2)$P_f\left(f_j\right)=\frac{n_j}{n}$
p_f=n_j/(rows*cols);
% (3)$C_\left(f\right)$
c=p_f;
for i=2:length(p_f)
    c(i)=c(i) + c(i-1);
end
% (4)$ \left \lfloor 255C\left(f\right)+0.5 \right \rfloor $
g=round(255*c);
% (5)$n_i$
n_i=zeros(1,length(n_j));
for i=1:length(n_j)
    temp=abs(g(i)+1);
    n_i(temp)=n_i(temp)+n_j(i);
end
% (6)$P_g\left(g_i\right)=\frac{n_i}{n}$
p_g=n_i/(rows*cols);
% (7)$f_i\rightarrow g_i$
h=figure;
ind=0:255;
plot(ind, g, '-');
title(['Mapping relation: $f_i {\rightarrow} g_i$'],'Interpreter','LaTeX');
xlabel('$f_i$','Interpreter','LaTeX');
ylabel('$g_i$','Interpreter','LaTeX');
saveas(h,[filename, '_mapping_relation'],'png');
close();
% start histogram equalization
for i=1:rows
    for j=1:cols
        for k=1:3
            img(i,j,k)=g(abs(img(i,j,k)+1));
        end
    end
end
drawgrayscalehistogram(img,[filename,'_equalization']);
if iscolorful
    img=hsi2rgb(img);
end
imwrite(img,[filename,'_equalization.jpg']);
csvwrite([filename,'_equalization.csv'],[ind;n_j;c;g;n_i;p_g]);

2.1.5 直方图规定化

function histogrammatching(img,filename,refimg,refimgfilename)
% e.g.
%   img=imread('color.jpg');
%   filename='color';
%   refimg=imread('grey_equalization.jpg');
%   refimgfilename='grey_equalization';
%   histogrammatching(img,filename,refimg,refimgfilename);
img=im2uint8(img);
rows=length(img(:,1));
cols=length(img(1,:))/length(img(1,1,:));
iscolorful=false;
if img(1,1,1)~=img(1,1,2)
    iscolorful=true;
    img=rgb2hsi(img);
end
% (1)$n_j$
n_j=drawgrayscalehistogram(img,filename);
% (2)$P_f\left(f_j\right)=\frac{n_j}{n}$
p_f=n_j/(rows*cols);
% (3)$C_\left(f\right)$
c=p_f;
for i=2:length(p_f)
    c(i)=c(i) + c(i-1);
end
% (4)$ \left \lfloor 255C\left(f\right)+0.5 \right \rfloor $
g=round(255*c);
 
% rules
z_i=drawgrayscalehistogram(refimg,refimgfilename);
[rowss,colss,ss]=size(refimg);
p_z=z_i/(rowss*colss);
% (5)$ C\left(z\right)=\sum_{i=0}^{k} P_z\left(Z_i\right) $
c_z=p_z;
for i=2:length(p_z)
    c_z(i)=c_z(i)+c_z(i-1);
end
% (6)$ \left \lfloor 255C\left(f\right)+0.5 \right \rfloor $
y_n=round(255*c_z); 
% (7)$f_i \rightarrow Z_i$
mapping=zeros(1,256);
for i=1:length(g)
    gs=abs(g(i));
    % find a grayscale in y_n which is closest to gs
    [temp,mapping(i)]=findclosest(gs, y_n);
end
 
h=figure;
ind=0:255;
plot(ind, mapping, '-');
title(['Mapping relation (refers to ', refimgfilename, '): $f_i {\rightarrow} g_i$'],'Interpreter','LaTeX');
xlabel('$f_i$','Interpreter','LaTeX');
ylabel('$g_i$','Interpreter','LaTeX');
saveas(h,[filename,'_refer_',refimgfilename,'_mapping_relation'],'png');
close();
 
% (8) P_z\left(Z_i\right)
z_i_=zeros(1,length(z_i));
for i=1:length(z_i)
    temp=abs(mapping(i)+1);
    z_i_(temp)=z_i_(temp)+z_i(i);
end
p_z_z_i=z_i_/(rows*cols);
 
% start histogram matching
for i=1:rows
    for j=1:cols
        for k=1:3
            img(i,j,k)=mapping(abs(img(i,j,k)+1));
        end
    end
end
if iscolorful
    img=hsi2rgb(img);
end
imwrite(img,[filename,'_matching.jpg']);
drawgrayscalehistogram(img,[filename,'_matching']);
csvwrite([filename,'_matching.csv'],[ind;n_j;c;g;z_i;p_z;c_z;y_n;p_z_z_i]);

2.1.6 调用脚本

function main()
addpath('./functions');
format long;
color=imread('color.jpg');
histogramequalization(color,'color');
grey=imread('grey.png');
histogramequalization(grey,'grey');
 
grey_equalization=imread('grey_equalization.jpg');
histogrammatching(color,'color',grey,'grey');
histogrammatching(color,'color',grey_equalization,'grey_equalization');

2.2 Python实现

2.2.0 引入（import）

import csv
import matplotlib.pyplot as plt
import matplotlib.image as mping
import numpy as np
from numpy import cos, arccos, sqrt, power, pi

2.2.1 RGB转HSI

# FUNCTION  RGB转HSI
# INPUT     RGB图像数据
# OUTPUT    uint8格式HSI图像数据
def rgb2hsi(rgb):
    # 如果没有归一化处理，则需要进行归一化处理（传入的是[0,255]范围值）
    if rgb.dtype.type == np.uint8:
        rgb = rgb.astype('float64')/255.0
    for i in range(rgb.shape[0]):
        for j in range(rgb.shape[1]):
            r, g, b = rgb[i, j, 0], rgb[i, j, 1], rgb[i, j, 2]
            # 计算h
            num = 0.5 * ((r-g)+(r-b))
            den = sqrt(power(r-g, 2)+(r-b)*(g-b))
            theta = arccos(num/den) if den != 0 else 0
            rgb[i, j, 0] = theta if b <= g else (2*pi-theta)
            # 计算s
            rgb[i, j, 1] = (1 - 3 * min([r, g, b]) / (r+g+b)) if r+g+b != 0 else 0
            # 计算i
            rgb[i, j, 2] = 1 / 3 * (r+g+b)
    return (rgb * 255).astype('uint8')
    ```
### 2.2.2 HSI转RGB
```python
# FUNCTION  HSI转RGB
# INPUT     HSI图像数据
# OUTPUT    uint8格式RGB图像数据
def hsi2rgb(hsi):
    if hsi.dtype.type == np.uint8:
        hsi = (hsi).astype('float64') / 255.0
    for k in range(hsi.shape[0]):
        for j in range(hsi.shape[1]):
            h, s, i = hsi[k, j, 0], hsi[k, j, 1], hsi[k, j, 2]
            r, g, b = 0, 0, 0
            if 0 <= h < 2/3*pi:
                b = i * (1 - s)
                r = i * (1 + s * cos(h) / cos(pi/3-h))
                g = 3 * i - (b + r)
            elif 2/3*pi <= h < 4/3*pi:
                r = i * (1 - s)
                g = i * (1 + s * cos(h-2/3*pi) / cos(pi - h))
                b = 3 * i - (r + g)
            elif 4/3*pi <= h <= 5/3*pi:
                g = i * (1 - s)
                b = i * (1 + s * cos(h - 4/3*pi) / cos(5/3*pi - h))
                r = 3 * i - (g + b)
            hsi[k, j, 0], hsi[k, j, 1], hsi[k, j, 2] = r, g, b
    return (hsi * 255).astype('uint8')

2.2.3 绘制灰度直方图

# FUNCTION  绘制灰度直方图
# INPUT     图像数据、灰度直方图保存文件名
# OUTPUT    (灰度级，对应灰度级的频数)
def draw_grayscale_histogram(img, filename=''):
    # 给定的两幅图像jpg是uint8，png是float32
    # 如果传入的是png图像，需要转化为[0,255]
    if img.dtype.type != np.uint8:
        img = (img*255).astype(np.uint8)
    # 由于彩色图像的直方图均衡化是在亮度通道上进行，故暂认为灰度直方图也绘制在亮度通道上进行统计绘制
    arr = np.array([0]*256)
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            # H S I分别是0 1 2位置
            # 灰色PNG图像RGB三通道都一样，彩色图像转为HSI后取I通道
            arr[img[i, j, 2]] += 1
    ind = np.arange(256)
    plt.bar(ind, arr)
    filename = "_".join([filename, "grayscale", "histogram"])
    plt.title(filename)
    plt.savefig(filename+".png", dpi=72)
    plt.close()
    return ind, arr

2.2.4 直方图均衡化

# FUNCTION  直方图均衡化
# INPUT     图像数据
# OUTPUT    均衡化后的图像数据
def histogram_equalization(img, filename=''):
    is_jpg = False
    is_png = False
    # 如果输入图像是.jpg格式，需要转HSI
    if img.dtype.type == np.uint8:
        is_jpg = True
        img = rgb2hsi(img)
    # 如果输入图像是.png格式，为方便处理需要转[0,255]
    elif img.dtype.type != np.uint8:
        is_png = True
        img = (img*255).astype(np.uint8)
    # （1）统计$n_j$
    ind, n_j = draw_grayscale_histogram(img, filename)
    # （2）计算$P_f\left(f_j\right)=\frac{n_j}{n}$
    p_f = n_j / (img.shape[0]*img.shape[1])
    # （3）计算$C_\left(f\right)$
    # c = [p_f[i-1]+p_f[i] for i in range(1, len(p_f))]
    c = p_f
    for i in range(1, len(p_f)):
        c[i] += c[i-1]
    # （4）求$ \left \lfloor 255C\left(f\right)+0.5 \right \rfloor $
    g = [int(255*ele+0.5) for ele in c]
    # （5）计算$n_i$
    n_i = np.array([0]*len(n_j))
    for i in range(len(n_j)):
        n_i[g[i]] += n_j[i]
    # （6）计算$P_g\left(g_i\right)=\frac{n_i}{n}$
    p_g = n_i / (img.shape[0]*img.shape[1])
    # （7）映射关系$f_i\rightarrow g_i$
    plt.plot(ind, g, '-')
    plt.title('''Mapping relation: $f_i \\rightarrow g_i$''')
    plt.xlabel('''$f_i$''')
    plt.ylabel('''$g_i$''')
    plt.savefig('_'.join([filename, "mapping", "relation"]))
    plt.close()
    # 绘制直方图均衡化后的灰度直方图
    # 如果输入图像是.png格式，均衡化之前需要转化到[0,255]
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            for k in range(3):
                img[i, j, k] = g[img[i, j, k]]
    # 如果输入图像是.png格式，均衡化转换之后需要转化回[0,1]
    if is_png:
        img = (img/255).astype(np.float32)
    # 如果输入图像是.jpg格式，均衡化转换之后需要转换回RGB
    elif is_jpg:
        img = hsi2rgb(img)
    # 如果输入是.png，则输出也是".png"；如果输入是".jpg"，则输出也是".jpg"
    plt.imsave('_'.join([filename, "equalization"])+(".png" if is_png else ".jpg"), img)
    # 绘制直方图均衡化后的
    draw_grayscale_histogram(img, '_'.join([filename, "equalization"]))
    with open('_'.join([filename, "equalization"])+".csv", 'w+') as f:
        f_csv = csv.writer(f)
        f_csv.writerow(ind)
        f_csv.writerow(n_j)
        f_csv.writerow(c)
        f_csv.writerow(g)
        f_csv.writerow(n_i)
        f_csv.writerow(p_g)

2.2.5 调用脚本

# FUNCTION  调用测试代码
if __name__ == '__main__':
    png = mping.imread('./grey.png').copy()
    histogram_equalization(png, 'grey')

    jpg = mping.imread('./color.jpg').copy()
    histogram_equalization(jpg, 'color')

    # 测试RGB转HSI和HSI转RGB
    # jpg = mping.imread('./color.jpg').copy()
    # hsi = rgb2hsi(jpg)
    # plt.imsave('./color_hsi.jpg', hsi)
    # jpg_recover = hsi2rgb(hsi)
    # plt.imsave('./color_recover.jpg', jpg_recover)

3 实验结果

3.1 grey.png直方图均衡化处理结果

左、右分别为直方图均衡化前后。下一、下二分别为其灰度直方图。


下为直方图均衡化时所用的灰度级变换曲线图。

数值数据见实验文件生成的.csv文件。

3.2 color.jpg直方图均衡化处理结果

左、右分别为直方图均衡化前后。下一、下二分别为其（亮度通道）灰度直方图。

下为直方图均衡化时（对亮度通道）所用的灰度级变换曲线图。

数值数据见实验文件生成的.csv文件。

3.2 color.jpg直方图规定化处理结果

上图左、右分别为直方图均衡化前后。下图左为规定直方图，右为直方图规定化结果。

下图从上至下分别为原图、相对grey直方图规定化、相对grey直方图均衡化结果直方图规定化的灰度直方图。

从上至下分别为相对grey直方图规定化、相对grey直方图均衡化结果直方图规定化的灰度级变换曲线图。