数字图像处理MATLAB学习笔记（七）

数字图像处理MATLAB学习笔记（七）

Image Compression

1. Basic Background

图像压缩的两个结构块：

• 编码器：包含映射器、量化器、符号编码器
• 解码器：符号解码器、反量化器
  • 关于解码器输出的图像$\hat{f}(x,y)$，我们定义其与输入图像$f(x,y)$的误差为$e(x,y)$。则对任意位置的误差值为$e(x,y)=\hat{f}(x,y)-f(x,y)$。两幅图像间的总误差为$\sum^{{M-1}_{x=0}\sum}{N-1}_{y=0}[\hat{f}(x,y)-f(x,y)]\$
  • 我们设$f(x,y)$与$\hat{f}(x,y)$之间的***root-mean-square|rms均方根误差***为$e_{rms}$，也就是MxN数组的均方误差的平均值的平方根：$e_{rms}=\Big[\frac{1}{MN}\sum^{{M-1}_{x=0}\sum}{N-1}_{y=0}\big[\hat{f}(x,y)-f(x,y)\big]^2\Big]{1/2}\$

我们可以通过数字量化，令$n_1$和$n_2$分别表示原始图像和编码后图像所携带的信息单元的数量，那么压缩比为：$C_R=\frac{n_1}{n_2}\$。如果比率为10:1的压缩比，那么也就是说original image中的10个信息单元被压缩为1个信息单元。

我们也可以利用M-function计算压缩比：function imratio

function cr = imratio(f1, f2)
%IMRATIO Computes the ratio of the bytes in two images/variables.
%   CR = IMRATIO(F1, F2) returns the ratio of the number of bytes in
%   variables/ files F1 and F2. If F1 and F2 are an original and
%   compressed image, respectively, CR is the compression ratio.

error(nargchk(2, 2, nargin));   % check input arguments
cr = bytes(f1) / bytes(f2);

% ------------------------------------------------------------------------
function b = bytes(f)
% Return the number of bytes in input f. If f is a string, assume that it
% is an image filename; if not, it is an image variable.

if ischar(f)
    info = dir(f);
    b = info.bytes;
elseif isstruct(f)
    % MATLAB's whos function reports an extra 124 bytes of memory per
    % structure field because of the way MATLAB stores structures in
    % memory. Don't count this extra memory; instead, add up the memory
    % associated with each field.
    b = 0;
    fields = fieldnames(f);
    for k = 1:length(fields)
        elements = f.(fields{k});
        for m = 1:length(elements)
            b = b + bytes(elements(m));
        end
    end
else
    info = whos('f');
    b = info.bytes;
end

也可以自定义函数function compare来计算$e_{rms}$并声称对应的直方图：

function rmse = compare(f1, f2, scale)
%COMPARE Computes and displays the error between two matrices.
%   RMSE = COMPARE(F1, F2, SCALE) returns the root-mean-square error
%   between inputs F1 and F2, displays a histogram of the difference,
%   and displays a scaled difference image. When SCALE is omitted, a
%   scale factor of 1 is used.

% Check the input arguments and set defaults.
error(nargchk(2, 3, nargin));
if nargin < 3
    scale = 1;
end

% Compute the root-mean-square error
e = double(f1) - double(f2);
[m, n] = size(e);
rmse = sqrt(sum(e(:).^2) / (m*n));

% Output error image & histogram if an error (i.e. rmse ~= 0)
if rmse
    % Form error histogram.
    emax = max(abs(e(:)));
    [h, x] = hist(e(:), emax);
    if length(h) >= 1
        figure, bar(x, h, 'k');
        % Scale the error image symmertically and display.
        emax = emax / scale;
        e = mat2gray(e, [-emax, emax]);
        figure, imshow(e);
    end
end

2. Coding Redundancy 编码冗余

令具有概率$p_r(r_k)$的离散随机变量为$r_k$，其中$k=1,2,\dots ,L$，描述图像的灰度级，且$p_r(r_k)=\frac{n_k}{n}\$，$n_k$即为图像中出现第k级灰度的次数，$n$是图像的总像素数。

如果用于表示每个$r_k$值的比特数是$l(r_k)$，那么表示每个pixel的平均比特数为$L_{agv}=\sum_{k=1}^Ll(r_k)p_r(r_k)\$。那么一张MxN大小的图像的总比特数为$MNL_{agv}$

假定不丢失信息时，可以充分描绘一幅图像的最小数据量。我们设这个随机事件为$E$，发生该事件的概率为$P(E)$，信息单位为$I(E)=\log{\frac{1}{P(E)}}=-\log P(E)\$。

在离散的可能事件的集合 { a 1 , a 2 , … , a J } \{a_1,a_2,\dots,a_J\} { a1​,a2​,…,aJ​}中给定随机事件源，其对应的概率为 { P ( a 1 ) , P ( a 2 ) , … , P ( a J ) } \{P(a_1),P(a_2),\dots,P(a_J)\} { P(a1​),P(a2​),…,P(aJ​)}，每个输出信源的平均信息被称为信息源的熵entropy，即$H=-\sum_{j=1}^JP(a_j)\log P(a_j)\$。

如果一幅图像被看作发出自身灰度级信息源的样本，就可以用被观测图像的灰度级直方图来对信息源的符号概率建模，此时信息源的熵的一阶估计为$\tilde{H}=-\sum_{k=1}^Lp_r(r_k)\log{p_r(r_k)}\$

对应的M-Function计算代码如下：function ntrop

function h = ntrop(x, n)
%NTROP Computes a first-order estimate of the entropy of a matrix.
%   H = NTROP (X, N) returns the entropy of matrix X with N
%   symbols. N = 256 if omitted but it must be larger than the
%   number of unique values in X for accurate results. The estimate
%   assumes a statistically independent source characterized by the
%   relative frequency of occurrence of the e lements in X.
%   The estimate is a lower bound on the average number of bits per
%   unique value (or symbol) when coding without coding redundancy.

% Check the arguments
error(nargchk(1, 2, nargin));

if nargin < 2
     n = 256;
end

% Make input double
x = double(x);
% Compute N-bin historgam
xh = hist(x(:), n);
% Compute probalities
xh = xh / sum(xh(:));

% Make make to eliminate 0's since log2(0) = -inf
i = find(xh);
% Compute entropy
h = -sum(xh(i) .* log2(xh(i)));

Example Computing entropy：

>> f = [119 123 168 119; 123 119 168 168];
>> f = [f; 119 119 107 119; 107 107 119 119]
f =
   119   123   168   119
   123   119   168   168
   119   119   107   119
   107   107   119   119
>> p = hist(f(:), 8);
>> p = p / (sum(p))
p =
    0.1875    0.5000    0.1250         0         0         0         0    0.1875
>> h = ntrop(f)
h =
    1.7806

2.1 Huffman Code

离散数学----数据结构----霍夫曼编码

Huffman Tree算法流程：

1. 根据给定的n个权值 { w 1 , w 2 , … , w n } \{w_1,w_2,\dots,w_n\} { w