灰度图像的Huffman编码 MATLAB代码

clear all;
close all;
clc;

% src = imread('kids.tif');
src = imread('lenna_RGB.tif');
src_d = rgb2gray(src);

[m, n] = size(src);
src_size = m * n;
gray_level = 256; % 灰度级;

% prob数组保存图像中各灰度出现的概率
prob = [];
for src_value = 0 : (gray_level - 1)
    index = find(src_d == src_value);
    i = src_value + 1;
    prob(i) = length(index) / src_size;
end

% 画出直方图
% stem(0: gray_level-1, prob);
% xlabel('灰度值');
% ylabel('概率');
% title('灰度直方图');

% huffman编码
p = prob;
n = length(p);
q = p;
m = zeros(n-1, n);
for i = 1 : n - 1 
    [q, l] = sort(q);
    m(i,:) = [l(1 : n - i + 1), zeros(1, i - 1)];
    q = [q(1) + q(2), q(3 : n), 1];
end

bre = zeros(n-1, n); 
bre(n-1, 1) = 0 + j; % 虚部表示当前的二进制数的位数，以下类似
bre(n-1, 2) = 1 + j; % 给bre最后一行的前两列赋值
% 通过当前行给上一行的bre赋值，前两列通过m等于为1的值的虚实部赋值，后面等于m大于1得值
% m值为1的数都是上一行两个最小的概率相加得到的
for time = 1 : n-2   % time = 1 : 7
    loc_1 = find(real(m(n-time, :)) == 1); % 找m的行实部为1的返回指数,从最后一行往上找
    prebit = bre(n-time,loc_1);   % prebit等于bre当前行，实部为1的列的值)
    % 给bre往上一行一二列赋值，第一列为prebit实部*2，prebit虚部加1
    % 第二列为prebit实部*2加1，prebit虚部加1
    % 实部乘二代表多加一位，最低位一个赋0，一个赋1；虚部为位数，加一
    bre(n-time-1,1) = (real(prebit)*2 + 0) + j *(imag(prebit)+1); 
    bre(n-time-1,2) = (real(prebit)*2 + 1) + j *(imag(prebit)+1); 
    % 找m的当前行实部大于1的返回指数
    loc_not1 = find(real(m(n-time, :)) > 1);
    % 给bre的上一行，第三列及之后列赋值
    % 等于当前行的实部大于1的值
    bre(n-time-1,3:3+time-1) = bre(n-time,loc_not1);
end

[m1,index] = sort(m(1,:));  % 对m进行排序得到m1(按照灰度值排序)，m是对概率进行排序
code = bre(1,index);  % code等于bre按照index排序的第一行
code_data = real(code);
code_bits = imag(code);
disp(['gray level',' ', 'huffman code']);
for i = 1 : length(code)
    disp([num2str(i-1), ' ', num2str(dec2bin(code_data(i)))]);
    disp([num2str(i-1), ' ', num2str(dec2bin(code_data(i),code_bits(i)))]);
end
code_binary = dec2bin(code_data);

%逐点编码
out = [];
for row = 1:m
    for vol = 1:n
        now_gray = src_d(row,vol);
        now_code = code_binary(now_gray+1,:);
        now_bits = code_bits(now_gray+1);
        now_code = now_code(end-now_bits+1:end);
        out = [out, now_code];
    end
end

%计算压缩比
real_bitnum = length(out);
bitnum_no_huffman = src_size * nextpow2(gray_level);
comp_ratio = bitnum_no_huffman / real_bitnum;
Lavg = real_bitnum / src_size;
% Hshannon = (-1) * prob * (log2(prob))';
disp(['Lavg = ', num2str(Lavg)]);
disp(['normal bit num = ', num2str(nextpow2(gray_level))]);
disp(['comp_ratio = ', num2str(comp_ratio)]);
% disp(['Hshannon = ', num2str(Hshannon)]);