数学形态学及图像压缩

                            数学形态学及图像压缩

1.基本概念：

  （1）数学形态学一般是使用二值图像，进行边界提取，骨架提取，孔洞填充，角点提取，图像重建。

   基本的算法: 膨胀腐蚀，开操作，闭操作，击中击不中变换。

 （2）形态学结构元素选择：

    结构元素的选择非常重要, 对同一类型的形态学操作, 选择不同形状的结构元素会得到不同的结果。结构元素的大小也会影响最终效果。如果结构元素太小, 对噪声的滤除作用不明显;若结构元素太大 ,滤除噪声的同时又会对目标本身造成较大影响, 故结构元素对形态学操作来说是关键。

    结构元素形状的影响：在进行形态学操作时,结构元素形状的选择是一个重要因素。对同一图像,当结构元素形状不同时,形态学操作结果也不一样 。

   

    从图3可以看出, 对原始图像 ,膨胀操作整体上使之变得“饱满”。 图(b)中的结构元素A1是边长为4的正方形结构素 ,由于原始图像为矩形,因此膨胀后的图像仍为矩形; 图(c)中的结构元素B1是半径为2的菱形,膨胀后的图像为一个 近似的八边形。对于同一幅原始图像,采用不同的结构元素进行膨胀, 得到的结果有明显差异, 因此, 结构元素的形状 影响数学形态学操作的结果。
   结构元素大小的影响：结构元素形状相同时, 结构元素大小影响形态学操作的结果。 选择两个正方形的结构元素 A2和B2, A2的半径为5, B2的半径为9,分别对图2中的图(a)进行腐蚀操作, 结果如图4所示:

   

    从图4可以看出, 对原始图像, 腐蚀操作使之整体上变“瘦”。图(b)中, 结构元素A2半径为5, 左侧的圆和右的矩 形变小,中间的横线明显变细; 随着结构元素半径的增大, 当B2半径为9时, 如图(c)所示,不仅左侧的圆和右侧的矩形变得更小,中间的横线甚至完全消失。 当结构元素形状相同而半径不同时, 腐蚀的效果有较大差异,因此,结构元素的大小影响数学形态学操作的结果。

 （3）JPEG的压缩原理：

   JPEG是联合图象专家组(Joint Picture Expert Group)的英文缩写，是国际标准化组织(ISO)和CCITT联合制定的静态图象的压缩编码标准。和相同图象质量的其它常用文件格式(如GIF，TIFF，PCX)相比，JPEG是目前静态图象中压缩比最高的。

   正是由于JPEG的高压缩比，使得它广泛地应用于多媒体和网络程序中，例如HTML语法中选用的图象格式之一就是JPEG(另一种是GIF)。这是显然的，因为网络的带宽非常宝贵，选用一种高压缩比的文件格式是十分必要的。

   JPEG有几种模式，其中最常用的是基于DCT变换的顺序型模式，又称为基线系统(Baseline)。

   JPEG编码器流程：

   解码器流程：

   
   本处引用自：JPEG编码 http://web.gdut.edu.cn/~dj/jxsb/cai/mp82.htm

  （4）Huffman编码是一种编码方式，是一种用于无损数据压缩的熵编码（权编码）算法。

  Huffman编码具体步骤：

     ① 概率统计（如对一幅图像，或m幅同种类型图像作灰度信号统计），得到n个不同概率的信息符号。　      ② 将n个信源信息符号的n个概率，按概率大小排序。　      ③ 将n个概率中，最后两个小概率相加，这时概率个数减为n-1个。　      ④ 将n-1个概率，按大小重新排序。　      ⑤ 重复③，将新排序后的最后两个小概率再相加，相加和与其余概率再排序。　 　　 ⑥ 如此反复重复n-2次，得到只剩两个概率序列。　 　　 ⑦ 以二进制码元（0.1）赋值，构成哈夫曼码字。编码结束。　 

    哈夫曼码字长度和信息符号出现概率大小次序正好相反，即大概信息符号分配码字长度短，小概率信息符号分配码字长度长。 

2.函数说明：

   （1）形态学用到matlab自带的函数：

     imdilate ：膨胀     imerode：腐蚀     strel：生成结构元素     imopen: 开操作     imclose: 闭操作     bwhitmiss: 击中击不中     bwmorph: 实现多种形态学操作

    （2）图像压缩编码用到的matlab函数：

   huffman 编解码：mat2huff，huff2mat   JPEG 编解码：im2jpeg，jpeg2im       JPEG2000编解码：im2jpeg2k，jpeg2k2im

 3.MATLAB源码：

  （1）compare.m

    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. %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004.%   $Revision: 1.3 $  $Date: 2003/04/18 05:07:33 $% Check 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 symmetrically and display      emax = emax / scale;      e = mat2gray(e, [-emax, emax]);      figure;   imshow(e);   endend

   (2) huff2mat.m

function x = huff2mat(y)%HUFF2MAT Decodes a Huffman encoded matrix.%   X = HUFF2MAT(Y) decodes a Huffman encoded structure Y with uint16%   fields: %      Y.min    Minimum value of X plus 32768%      Y.size   Size of X%      Y.hist   Histogram of X%      Y.code   Huffman code%%   The output X is of class double.%%   See also MAT2HUFF.%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004%   $Revision: 1.5 $  $Date: 2003/11/21 13:17:50 $if ~isstruct(y) | ~isfield(y, 'min') | ~isfield(y, 'size') | ...       ~isfield(y, 'hist') | ~isfield(y, 'code')   error('The input must be a structure as returned by MAT2HUFF.');endsz = double(y.size);   m = sz(1);   n = sz(2);xmin = double(y.min) - 32768;            % Get X minimummap = huffman(double(y.hist));           % Get Huffman code (cell)% Create a binary search table for the Huffman decoding process.% 'code' contains source symbol strings corresponding to 'link'% nodes, while 'link' contains the addresses (+) to node pairs for% node symbol strings plus '0' and '1' or addresses (-) to decoded% Huffman codewords in 'map'. Array 'left' is a list of nodes yet to% be processed for 'link' entries.code = cellstr(char('', '0', '1'));     % Set starting conditions aslink = [2; 0; 0];   left = [2 3];       % 3 nodes w/2 unprocessedfound = 0;   tofind = length(map);      % Tracking variableswhile length(left) & (found < tofind)   look = find(strcmp(map, code{left(1)}));    % Is string in map?   if look                            % Yes      link(left(1)) = -look;          % Point to Huffman map      left = left(2:end);             % Delete current node      found = found + 1;              % Increment codes found         else                               % No, add 2 nodes & pointers      len = length(code);             % Put pointers in node      link(left(1)) = len + 1;            link = [link; 0; 0];            % Add unprocessed nodes      code{end + 1} = strcat(code{left(1)}, '0');      code{end + 1} = strcat(code{left(1)}, '1');            left = left(2:end);             % Remove processed node      left = [left len + 1 len + 2];  % Add 2 unprocessed nodes   endendx = unravel(y.code', link, m * n);    % Decode using C 'unravel'x = x + xmin - 1;                     % X minimum offset adjustx = reshape(x, m, n);                 % Make vector an array

   (3) huffman.m

function CODE = huffman(p)%HUFFMAN Builds a variable-length Huffman code for a symbol source.%   CODE = HUFFMAN(P) returns a Huffman code as binary strings in%   cell array CODE for input symbol probability vector P. Each word%   in CODE corresponds to a symbol whose probability is at the%   corresponding index of P. %%   Based on huffman5 by Sean Danaher, University of Northumbria,%   Newcastle UK. Available at the MATLAB Central File Exchange:%   Category General DSP in Signal Processing and Communications. %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004%   $Revision: 1.5 $  $Date: 2003/10/26 18:37:16 $% Check the input arguments for reasonableness.error(nargchk(1, 1, nargin));if (ndims(p) ~= 2) | (min(size(p)) > 1) | ~isreal(p) | ~isnumeric(p)   error('P must be a real numeric vector.');     end   % Global variable surviving all recursions of function 'makecode'global CODECODE = cell(length(p), 1);  % Init the global cell array                            if length(p) > 1            % When more than one symbol ...   p = p / sum(p);          % Normalize the input probabilities   s = reduce(p);           % Do Huffman source symbol reductions   makecode(s, []);         % Recursively generate the codeelse     CODE = {'1'};            % Else, trivial one symbol case!end;   %-------------------------------------------------------------------%function s = reduce(p);% Create a Huffman source reduction tree in a MATLAB cell structure% by performing source symbol reductions until there are only two% reduced symbols remainings = cell(length(p), 1);% Generate a starting tree with symbol nodes 1, 2, 3, ... to % reference the symbol probabilities.for i = 1:length(p)   s{i} = i; endwhile numel(s) > 2   [p, i] = sort(p);    % Sort the symbol probabilities   p(2) = p(1) + p(2);  % Merge the 2 lowest probabilities   p(1) = [];           % and prune the lowest one      s = s(i);            % Reorder tree for new probabilities   s{2} = {s{1}, s{2}}; % and merge & prune its nodes   s(1) = [];           % to match the probabilitiesend%-------------------------------------------------------------------%function makecode(sc, codeword)% Scan the nodes of a Huffman source reduction tree recursively to% generate the indicated variable length code words.% Global variable surviving all recursive callsglobal CODE                           if isa(sc, 'cell')                   % For cell array nodes,   makecode(sc{1}, [codeword 0]);    % add a 0 if the 1st element   makecode(sc{2}, [codeword 1]);    % or a 1 if the 2ndelse                                 % For leaf (numeric) nodes,   CODE{sc} = char('0' + codeword);  % create a char code stringend

   (4) im2jpeg2k.m

function y=im2jpeg2k(x,n,q)% IM2JPEG2K Compresses an image wsing a JPEG 2000 approximation.%  Y=IM2JPEG2K(X,N,Q) compresses image X using an N-scale JPEG%  2K wavelet transform,implicit or explicit coefficient%  quantization, and Huffman sybol coding augmented by zero%  run-length coding. IF quantization vector Q contains two%  elementers, they are assumed to be implicit quantization%  parameters; else, it is assumed to contain explicit subband step%  sizes. Y is an encoding structure containing Huffman-encoded%  data and additional parameters needed by JPEG2K2IM for decoding%%  See also JPEG2K2IM.global RUNSerror(nargchk(3, 3, nargin));      %Check input argumentsif ndims(x) ~= 2 | ~isreal(x) | ~isnumeric(x) | ~isa(x, 'uint8')    error('The input must be a UINT8 image.');endif length(q) ~= 2 & length(q) ~= 3 * n + 1    error('The quantization step size vector is bad.');end% Level shift the input and compute its wavelet transformx = double(x) - 128;[c, s] = wavefast(x, n,'jpeg9.7');% Quantize the wavelet coefficients.q = stepsize(n, q);sgn = sign(c);    sgn(find(sgn == 0)) = 1;    c = abs(c);for k = 1:n    qi = 3 * k - 2;    c = wavepaste('h', c, s, k, wavecopy('h', c, s, k) / q(qi));    c = wavepaste('v', c, s, k, wavecopy('v', c, s, k) / q(qi + 1));    c = wavepaste('d', c, s, k, wavecopy('d', c, s, k) / q(qi + 2));endc = wavepaste('a', c, s, k, wavecopy('a', c, s, k) / q(qi + 3));c = floor(c);     c = c .* sgn;% Run-length code zero runs of more than 10. Begin by creating% a special code for 0 runs ('zrc') and end-of code ('eoc') and% making a run-length table.zrc = min(c(:)) - 1;   eoc = zrc - 1;   RUNS = [65535];%  Find the run transition points: 'plus' contains the index of the%  start of a zero run; the correspongding 'minus' is its end + 1.z = c ==0;                z = z - [0 z(1:end - 1)];puls = find(z == 1);       minus = find(z == -1);% Remove any terminating zero run from 'c'.if length(plus) ~= length(minus)    c(plus(end) :end) = [];  c = [c eoc];end% Remove all other zero runs (based on 'puls' end 'minus') from 'c'.for i = length(minus) :-1:1    run = minus(i) - plus(i);    if run > 10        ovrflo = floor(run / 65535);      run = run - ovrflo * 65535;        c = [c(1:plus(i) - 1) repmat([zrc 1], 1, ovrflo) zrc ...                runcode(run) c(minus(i):end)];    endend% Huffman encode and add misc. information fordecoding.y.runs     = uint16(RUNS);y.s        = uint16(s(:));y.zrc      = uint16(-zrc);y.q        = uint16(100 * q);y.n        = uint16(n);y.huffman  = mat2huff(c);%------------------------------------------------------------------------%function y = runcode(x)% Find a zero run in the run-length table. If lot found, create a% new entry in the table.Return the index of the run.global RUNSy = find(RUNS == x);if length(y) ~= 1    RUNS = [RUNS; x];    y = length(RUNS);end%-------------------------------------------------------------------------%function q = stepsize(n, p)% Create a subband quantization array of stepsizes ordered by% decomposition (first to last) and subband (herizontal, vertical'% diagonal, and for final decomposition the approximation subband).if length(p) == 2            % Implicit Quantization    q = [];    qn = 2 ^ (8 - p(2) + n) * (1+ p(1) / 2 ^ 11);    for k = l:n        qk = 2 ^ - k * qn;        q = [q (2 * qk) (2 * qk) (4 * qk)];    end    q = [q qk];else                         % Explicit Quantization          q = p;endq = round(q * 100) / 100;    % Round to 1/100th placeif any(100 * q > 65535)    error('The quantizing steps are not UNIT16 representable.');endif any(q == 0)    error('A quantizing step of 0 is not allowed.');end 

   (5) imratio.m

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. %   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004%   $Revision: 1.4 $  $Date: 2003/11/21 13:11:15 $error(nargchk(2, 2, nargin));         % Check input argumentscr = bytes(f1) / bytes(f2);           % Compute the ratio%-------------------------------------------------------------------%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)      b = b + bytes(f.(fields{k}));   endelse   info = whos('f');     b = info.bytes;end

   (6) jpeg2k2im.m

function x = jpeg2k2im(y)%JPEG2K2IM Decodes an IM2JPEG2K compressed image.% x = JPEG2K2IM(Y) decodes compressed image Y,reconstructing an % approximation of the original image X. Y is an encoding% strucure returned by IM2JPEG2K%% See also IM2JPEG2K.error(nargchk(1, 1, nargin));       % Check input arguments% Get decoding parameters: scale, quantization vector, run-length% table size,zero run code, end-of-data code, wavelet bookkeeping% array, and run-length table.n = double(y.n);q = double(y.q) / 100;runs = double(y,runs);rlen = length(runs);zrc = -double(y.zrc);eoc = zrc - 1;s = double(y.s);s = reshape(s, n + 2, 2);% Compute the size of the wavelet transform.c1 = prod(s(1, :));for i = 2:n + 1    c1 = c1 +3 * prod(s(i, :));end% perform Huffman decoding followed by zero run decoding.r = huff2mat(y.huffman);c = [];  zi = find(r == zrc);   i = 1;for j = 1:length(zi)    c = [c r(i:zi(j) - 1) zeros(1, runs(r(zi(j) + 1)))];    i = zi(j) + 2;end zi = find(r == eoc);               % Undo terminating zero runif length(zi) == 1    c = [c r(i:zi - 1)];    c = [c zeros(1, c1 - length(c))];else    c = [c r(i:end)];end % Denotmalize the coefficients.c = c + (c > 0) - (c < 0);for k = 1:n    qi = 3 * k - 2;    c = wavepaste('h', c, s, k, wavecopy('h', c, s, k) * q(qi));    c = wavepaste('v', c, s, k, wavecopy('v', c, s, k) * q(qi + 1));    c = wavepaste('d', c, s, k, wavecopy('d', c, s, k) * q(qi + 2));endc = wavepaste('a', c, s, k, wavecopy('a', c, s, k) * q(qi + 3));% Compute the inverse wavelect transform and level shift.x = waveback(c, s, 'jpeg9.7', n);x = uint8(x + 128);

   (7) mat2huff.m

function y = mat2huff(x)% MAT2HUFF Huffman encode a matrix.%   Y = MAT2HUFF(X) Huffman encodes matrix X using symbol%   probabilities in unit-width histogram bins between X's minimum%   and maximum values. The encode data is returned as a structure%   Y:%       Y.code      The Huffman-encoded values of X, stored in%                   a uint16 vector. The other fields of Y contain%                   additional decoding informaition, including:%       Y.min       The minimum value of X plus 32768%       Y.size      The size of X%       Y.hist      The histogram of X%%   If X is logical, uint8, uint16, uint32, int8, int16, or double,%   with integer values, it can be input directly to MAT2HUFF. The%   minimum value of X must be representable as an int16.%%   If X is double with non-integer values---for example, an image%   with values between 0 and 1---first scale X to an appropriate%   integer range before the call. For example, use Y =%   MAT2HUFF(255*X) for 256 gray level encoding.%%   NOTE: the number of Huffman code words is round(max(X(:))) -%   round(min(X(:))) + 1. You may need to scale input X to generate%   codes of reasonable length. The maximum row or colum dimension%   of X is 65535.%%   See also HUFF2MATif ndims(x) ~= 2 | ~isreal(x) | (~isnumeric(x) & ~islogical(x))    error('X must be a 2-D real nemeric or logical matrix.');end% Store the size of input x.y.size = uint32(size(x));% Find the range of x values and store its minimum value biased% by +32768 as a UINT16.x = round(double(x));xmin = min(x(:));xmax = max(x(:));pmin = double(int16(xmin));pmin = uint16(pmin + 32768);    y.min = pmin;% Compute the input histogram between xmin and xmax with unit% width bins, scale to UINT16, and store.x = x(:)';h = histc(x , xmin:xmax);if max(h) > 65535    h = 65535 * h / max(h);endh = uint16(h);  y.hist = h;% Code the input matrix and store the result.map = huffman(double(h));           % Make Huffman code maphx = map(x(:) - xmin + 1);          % Map imagehx = char(hx)';                     % Convert to char arrayhx = hx(:)';hx(hx == ' ') = [];                 % Remove blanksysize = ceil(length(hx) / 16);      % Compute encoded sizehx16 = repmat('0' , 1 , ysize*16);  % Pre-allocate modulo-16 vectorhx16(1 : length(hx)) = hx;          % Make hx modulo-16 in lengthhx16 = reshape(hx16 , 16 , ysize);  % Reshape to 16-character wordshx16 = hx16' - '0';                 % Convert binary string to decimaltwos = pow2(15:-1:0);y.code = uint16(sum(hx16 .* twos(ones(ysize , 1) , :) , 2))';

   (8) Compressiontest.m

% f=imread('lena.bmp');% c=mat2huff(f);% cr1=imratio(f,c)% save lena c;% % load lena.mat;% g=huff2mat(c);% f=imread('lena.bmp');% rmse=compare(f,g);%  f=imread('lena.bmp');%  c1=im2jpeg(f,2);%  f1=jpeg2im(c1);%  imshow(f1);%  imratio(f,c1) %  f=load('lena.mat');   f=imread('worldcup.jpg'); c1=im2jpeg2k(f,7,[8 8.5]); f1=jpeg2k2im(c1); imshow(f1); imratio(f,c1)

4.附录：

 测试图片：