JPEG2000图像压缩标准研究.doc

通信工程一班 孙大江 邵军立 潘亚辉JPEG2000图像压缩标准研究 作者 孙大江20092420128 (算法) 邵军立20092420126 (程序) 潘亚辉20092420124 (仿真) 组长 邵军立摘要 JPEG 2000 正式名称为:ISO 15444，是由 JPEG 组织负责制定。该标准是由联合摄影专家组于 1997 年开始征集提案，把它作为 JPEG 标准的一个更新换代标准。它的目标是进一步改进目前压缩算法的性能，以适应低带宽、高噪声的环境，以及医疗图像、电子图书馆、传真、Internet 网上服务和保安等方面的应用。本文详细、全面的深入研究了 JPEG2000 图像压缩标准；系统研究了 JPEG2000压缩标准构成的理论基础；按模块分析了的 JPEG2000 整个的编码、解码系统；最后选择图像进行了matlab仿真,效果良好.关键词：图像压缩，JPEG2000，小波变换Abstract The Official name of JPEG2000 is ISO 15444.it is developed by the organization of JPEG. The Joint Photographic Expert Group began collecting sponsors from 1997.They want to make it as a replacement standards of the JPEG standard and its goal is to further improve the performance of the current compression algorithm so that it can adapt to the application of low-bandwidth, high noise environment, and medical images, electronic library, fax, Internet and other online services and security aspects. This article studyed the JPEG2000 coding standard comprehensively; studyed the based theory of the JPEG2000 compression standards constitute; analysis the entire JPEG2000 encoding and decoding system by the modules.finally,the simulation of the pictures got perfect result. Key Words: Image compression,JPEG2000, wavelet transform1. 引言 近年来,随着人们对图像压缩处理的要求越来越高,旧的图像处理标准(如JPEG)已经难以适应最新的形势,在这种背景下JPEG2000作为一个新的标准,以其对静止图像更优秀的压缩性能,和处理图像过程中更大的灵活性和伸缩性,在各种具体情况中得到了越来越多的应用.本文主要介绍了JPEG2000图像压缩标准的主要压缩算法及其编程实现,并给出了具体实例加以说明,二.JPEG2000标准概述JPEG2000 与传统 JPEG 最大的不同，在于它放弃了 JPEG 所采用的以离散余弦变换(Discrete Cosine Transform) 为主的区块编码方式，而采用以小波转换(Wavelet Transform) 为主的多解析编码方式。离散子波变换算法是现代谱分析工具，在包括压缩在内的图像处理与图像分析领域正得到越来越广泛的应用。此外 JPEG2000 还将彩色静态画面采用的 JPEG 编码方式与 2 值图像采用的 JBIG 编码方式统一起来，成为对应各种图像的通用编码方式。简单原理如下图所示JPEG200 标准提供了一套新的特征，这些特征对于一些新产品（如数码相机）和应用（如互联网）是非常重要的。它把 JPEG 的四种模式（顺序模式、渐进模式、无损模式和分层模式）集成在一个标准之中。在编码端以最大的压缩质量（包括无失真压缩）和最大的图像分辨率压缩图像，在解码端可以从码流中以任意的图像质量和分辨率解压图像，最大可达到编码时的图像质量和分辨率。JPEG2000应用的领域包括互联网、彩色传真、打印、扫描、数字摄像、遥感、移动通信、医疗图像和电子商务等等。标准提供了一套全新的特征，它的最主要的特征包括:1.高压缩率：由于在离散子波变换算法中，图像可以转换成一系列可更加有效存储像素模块的“子波”，因此，JPEG2000 格式的图片压缩比可在现在的 JPEG基础上再提高 1030，而且压缩后的图像显得更加细腻平滑，这一特征在互联网和遥感等图像传输领域有着广泛的应用。2.无损压缩和有损压缩：JPEG2000 提供无损和有损两种压缩方式，无损压缩在许多领域是必须的，例如医学图像中有时有损压缩是不能忍受的，再如图像档案中为了保存重要的信息较高的图像质量是必然的要求。同时 JPEG2000 提供的是嵌入式码流，允许从有损到无损的渐进解压。3.渐进传输：现在网络上的 JPEG 图像下载时是按“块”传输的，因此只能一行一行地显示，而采用 JPEG2000 格式的图像支持渐进传输(ProgressiveTransmission)。所谓的渐进传输就是先传输图像轮廓数据，然后再逐步传输其他数据来不断提高图像质量。互联网、打印机和图像文档是这一特性的主要应用场合。4.感兴趣区域压缩：可以指定图片上感兴趣区域（Region of Interest），然后在压缩时对这些区域指定压缩质量，或在恢复时指定某些区域的解压缩要求。这是因为子波在空间和频率域上具有局域性，要完全恢复图像中的某个局部，并不需要所有编码都被精确保留，只要对应它的一部分编码没有误差就可以了。5.码流的随机访问和处理：这一特征允许用户在图像中随机地定义感兴趣区域，使得这一区域的的图像质量高于其它图像区域.码流的随机处理允许用户进行旋转、移动、滤波和特征提取等操作。6.容错性：在码流中提供容错性有时是必要的，例如在无线等传输误码很高的通信信道中传输图像时，没有容错性是让人不能接受的。7.开放的框架结构：为了在不同的图像类型和应用领域优化编码系统，提供一个开放的框架结构在是必须的，在这种开放的结构中编码器只实现核心的工具算法和码流的解析，如果需要解码器可以要求数据源发送未知的工具算法。8.基于内容的描述：图像文档、图像索引和搜索在图像处理中是一个重要的领域，MPEG-7 就是支持用户对其感兴趣的各种“资料”进行快速、有效的检索的一个国际标准。基于内容的描述在 JPEG2000 中是压缩系统的特性之一。这些优秀的特征，使得 JPEG2000 成为新世纪的静止图像压缩标准。3. JPEG2000的编码过程及算法1. 编码过程JPEG2000 编码器的结构框图如图上部分所示，首先对源图像数据进行离散小波变换，然后对变换后的小波系数进行量化，接着对量化后的数据熵编码，最后形成输出码流。解码器是编码器的逆过程，如图下部分所示，首先对码流进行熵解码，然后解量化和小波反变换，最后生成重建图像数据。 JPEG2000 的处理对象不是整幅图像，而是把图像分成若干图像片（image tiles），对每一个图像片进行独立的编解码操作。术语“图像片”（tiling）是指原始图像被分成互不重叠的矩形块，对每一个图像片进行独立的编解码处理。在对每个图像片进行小波变换之前，通过减去一个相同的数量值对所有的图像片进行水平移位，如下图所示。 编码器的最后使用了算术编码器，在 JPEG2000 中使用的是 MQ 编码器，MQ编码器在本质上与 JPEG 中的 QM 编码器很相似。整个 JPEG2000 的编码过程可概括如下：1.把原图像分解成各个成分（亮度信号和色度信号）。2.把图像和它的各个成分分解成矩形图像片。图像片是原始图像和重建图像的基本处理单元。3.对每个图像片实施小波变换。4.对分解后的小波系数进行量化并组成矩形的编码块（code-block）。5.对在编码块中的系数“位平面”熵编码。6.为使码流具有容错性，在码流中添加相应的标识符（Maker）。7.可选的文件格式用来描述图像和它的各个成分的意义。2.JPEG2000的算法(1).小波变换和离散小波变换定义设函数 满足条件 (1)其中()是(t )的傅立叶变换,则称(t )为基本小波或母小波.上式*(1)称为“允许条件”，可以看出，能用作基本小波的函数 (t )必须满足否则C会在=0处趋于无穷大。与允许条件相应的时域要求是： (2)由式(2)可知基本小波必须是振荡的且平均值为 0。原则上讲，满足允许条件的函数(t )便可用作基本小波，不过往往要求更高些，即 (t )还要满足“正规条件”，以保证()在频域上表现出较好的局域性能。为达到此目的，要求(t )的前 n 阶原点矩为零，且 n 的值越大越好，即： (3)式(3)即“正规条件”。之所以要求 n 越大越好，是因为要尽可能地消除 f (t )的多项式中 (p n)项在小波变换中的作用，以便突出信号的高阶起伏和高阶导数中可能存在的奇点，即让小波变换能充分反映信号的高阶(细节)变化。引入尺度因子( 0)和平移因子 b,将基本小波进行伸缩和平移，得到函数族：称为分析小波。那么对于函数f (t ) L ( R)，其小波变换定义如下:由于 a,b 是连续变量，故成为连续小波变换(continue wavelet transform,CWT)。但是,在实际应用中，尤其是在计算机上实现时，连续小波必须加以离散化。在连续小波中，若对a取离散值，,j 为倍频程次数，进一步取 b 为 a 的倍数，为方便起见，设a 1，则所对应的离散小波函数可写为：而离散小波变换系数则可表示为：其重构公式为：其中，C 是一个与信号无关的常数。在实际应用中，常设a = 2、b = 1，由此得到的小波变换称为二进小波(2).小波提升算法提升算法(Lifting Scheme)的概念首先是由 Wim Sweldens 提出来的，是一种更快更有效的小波变换方法，它不依赖傅立叶变换，可以在空间域直接计算小波系数提升算法主要有两方面的应用。首先，它可以用来实现已经存在的小波。Daubechies和 Sweldens 已经证明，任何离散小波变换或者具有有限长滤波器的双通道子带滤波都可以分解成为一系列的提升步骤，所有能够用 Mallat 算法实现的小波换，都可以用提升算法来实现。因此，从这个层面来讲，提升算法只是一种新的实现小波变换的方法信号在小波变换后的性质只取决于所使用的小波，而与提升算法本身无关。其次，提升算法能够构造新的小波。虽然，某个具体的设计可能会使提升算法等同于某个第一代小波变换，但从本质而言，提升算法这一层次的运用属于第二代小波变换的范畴。由于提升算法不依赖于傅立叶变换，不需要伸缩和平移就可以构造小波，因此这样的小波称为第二代小波。与第一代小波相比，提升算法具有如下优点：第一，可以更快速地实现小波变换。第二，允许完全在空域上计算小波变换。换句话说，就是原始信号的数据可以由小波变换后的数据直接代替而不需要额外的存储器。第三，传统情况下，不能从小波反变换的形式上直接看出它是小波正变换的逆过程，只有通过傅立叶变换才能证明它的完全重建性质。而提升算法则很直观，只要进行与正变换相反的操作就可立即得到小波的反变换形式。(3).EBCOT 算法的研究及实现EBCOT 算法是 JPEG2000 标准中的核心算法之一，它能实现对图像的有效压缩，同时产生码流具有分辨率可伸缩性，信噪比可伸缩性，随机访问和处理等非常好的特性。EBCOT 算法与早期的嵌入式算法有相似的地方，都是采用小波变换进行子带采样，然后对小波系数进行量化和编码。与 EZW 和 SPIHT 算法不同的是，EBCOT 算法并没有使用零树结构而是使用编码块，同时提出了一种称为“压缩后率失真优化”的算法，基于这种思想产生的码流具有完全的嵌入特性。编码思想EBCOT 算法的基本思想是将小波变换以后的子带划分为大小固定的码块，对码块系数量化，按照二进制位分层的方法，从高有效位平面开始，依次对每个位平面上的所有小波系数位进行三个通道扫描建模(位平面编码)，生成上下文和 0、1 符号对，然后对这些上下文和符号对进行上下文算术编码，形成码块码流，完成第一阶段编码块编码(tier-1)；最后根据一定参数指标如码率、失真度，按率失真最优原则在每个独立码块码流中截取合适的位流组装成最终的图像压缩码流，完成第二阶段码流组装过程(tier-2)。算法框图如图所示。 第一阶段块编码(tier-1)位平面编码又称上下文建模(Context Formation)，其作用是对码块进行逐个位平面的扫描，对每个Bit位形成一个上下文标志(Context Label)，然后把这个上下文标志交给后面的算术编码模块进行处理编码框图如下图所示。JPEG2000 所用的 MQ 算术编码器是一种特殊的二进制算术编码器，属于自适应二进制算术编码器。所谓自适应算术编码是指编码系统用来划分区间的当前符号概率估计，是可以根据已经传输和编码的信息串进行调整的。一个自适应二进制算术编码需要使用统计模型，以便用来选择编码区间划分时所用的条件概率估计。MQ 算术编码器概率估计依赖于编码的某些“特征”，故又称为基于上下文的二进制算术编码。其主流程图如下:第二阶段码流组装(tier-2)第一阶段块编码得到的仅仅是各个独立码块的码流，为了使解码得到的图像具有不同的特性，必须对这些码流进行有效的组织。从上面的介绍可以知道以码块为单位的码流是按照块的不同失真度组织的，随着块码流的增加，失真度减少。为了使得全图像在一定码长下的失真度最小，就要从每块中裁剪部分码流组织在一起，这个过程称为打包过程，也就是第二阶段编码。这一过程实现了一定保真度和分辨率的码率可伸缩性和渐进性。 在JPEG2000 编解码系统中，EBCOT 算法是其重要的组成部分。而 EBCOT 算法当中的第一阶段块编码又是整个算法的核心，它占用了大量的编码时间，具体如下表所示。由表内容可以看出，无论是无损压缩还是有损压缩，EBCOT 算法都占整个编码器耗时的 70%以上，而其中的位平面编码时间，更是占到整个编码耗时的 50%以上。所以，自从 EBCOT 算法提出后，由于第一阶段块编码的运算量比较大、编码速度较慢，针对这种情况的优化研究很有必要。目前很多学者提出了相应的改进方法，比较有代表性的有样点省略法和群列省略法。这里不加以详细论述.块编码算法改进在进行位平面编码时虽然每个系数仅在一次通道扫描中编码输出，但是三次通道扫描都需要扫描该系数，而且每次扫描都有可能需要收集上下文信息。例如编码一个3232比特，深度为 N 的码块，由高到低需要编码的位平面 N 个，则需要扫描的次数高达3232( N3 2)次，与此同时还需要收集大量的上下文。所以位平面编码的编码时间占整个编码器编码时间 50%以上。而上下文形成的特点使三个编码通道中的系数的数量在不同位平面间变化很大。如下图所示，此图为barbara图像(256256)在位平面编码时三个通道系数编码数量变化示意图。图中通道1表示重要性通道，通道2表示幅度精练通道，通道3表示清除通道。由图可以看出，在最高位平面MSB，所有系数都只在清除通道进行编码。重要性传播通道中编码的系数数量先是增加，然后由于在重要性传播通道中的系数已变为重要的，所以在重要性传播通道中编码的系数数量又逐渐减少。在低位平面(0、1、2)大部分系数在幅度精练通道中编码，只有一少部分在重要通道中编码，没有系数在清除通道中编码。在整个扫描编码过程中，三次扫描要形成三次上下文并判断其所属编码通道，这样就会使编码时间大幅度增加。根据对图中的数据分析研究，目前存在两种针对于位平面编码的改进方法。即低位平面清除通道编码省略法和高位平面一次扫描法.这里也不再加以详细讲解.只是给出二者编码框图低位平面清除通道编码省略法编码框图高位平面一次扫描法编码框图4. 程序1.以下程序保存为.m文件到matlab安装文件中的work文件夹中，供主程序运行时调用。 imratiofunction cr = imratio(f1, f2) error(nargchk(2, 2, nargin); % Check input argumentscr = bytes(f1) / bytes(f2); % Compute the ratiofunction b = bytes(f) if ischar(f) info = dir(f); b = info.bytes;elseif isstruct(f) b = 0; fields = fieldnames(f); for k = 1:length(fields) b = b + bytes(f.(fieldsk); endelse info = whos(f); b = info.bytes;end wavefast function c, s = wavefast(x, n, varargin)error(nargchk(3, 4, nargin); if nargin = 3 if ischar(varargin1) lp, hp = wavefilter(varargin1, d); else error(Missing wavelet name.); endelse lp = varargin1; hp = varargin2; end fl = length(lp); sx = size(x); if (ndims(x) = 2) | (min(sx) 2) | isreal(x) | isnumeric(x) error(X must be a real, numeric matrix.); end if (ndims(lp) = 2) | isreal(lp) | isnumeric(lp) . | (ndims(hp) = 2) | isreal(hp) | isnumeric(hp) . | (fl = length(hp) | rem(fl, 2) = 0 error(LP and HP must be even and equal length real, . numeric filter vectors.); end if isreal(n) | isnumeric(n) | (n log2(max(sx) error(N must be a real scalar between 1 and . log2(max(size(X).); endc = ; s = sx; app = double(x);for i = 1:n app, keep = symextend(app, fl); rows = symconv(app, hp, row, fl, keep); coefs = symconv(rows, hp, col, fl, keep); c = coefs(:) c; s = size(coefs); s; coefs = symconv(rows, lp, col, fl, keep); c = coefs(:) c; rows = symconv(app, lp, row, fl, keep); coefs = symconv(rows, hp, col, fl, keep); c = coefs(:) c; app = symconv(rows, lp, col, fl, keep);endc = app(:) c; s = size(app); s;function y, keep = symextend(x, fl)keep = floor(fl + size(x) - 1) / 2);y = padarray(x, (fl - 1) (fl - 1), symmetric, both);function y = symconv(x, h, type, fl, keep)if strcmp(type, row) y = conv2(x, h); y = y(:, 1:2:end); y = y(:, fl / 2 + 1:fl / 2 + keep(2);else y = conv2(x, h); y = y(1:2:end, :); y = y(fl / 2 + 1:fl / 2 + keep(1), :);Endhuff2matfunction x = huff2mat(y)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)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 1) | isreal(p) | isnumeric(p) error(P must be a real numeric vector.); endglobal CODECODE = cell(length(p), 1); % Init the global cell arrayif 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);s = cell(length(p), 1);for i = 1:length(p) si = 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 s2 = s1, s2; % and merge & prune its nodes s(1) = ; % to match the probabilitiesendfunction makecode(sc, codeword)global CODE if isa(sc, cell) % For cell array nodes, makecode(sc1, codeword 0); % add a 0 if the 1st element makecode(sc2, codeword 1); % or a 1 if the 2ndelse % For leaf (numeric) nodes, CODEsc = char(0 + codeword); % create a char code stringendwave*function nc = wavepaste(type, c, s, n, x)error(nargchk(5, 5, nargin)nc = wavework(paste, type, c, s, n, x);function varargout = waveback(c, s, varargin)error(nargchk(3, 5, nargin);error(nargchk(1, 2, nargout);if (ndims(c) = 2) | (size(c, 1) = 1) error(C must be a row vector.); end if (ndims(s) = 2) | isreal(s) | isnumeric(s) | (size(s,2) = 2) error(S must be a real, numeric two-column array.); end elements = prod(s, 2);if (length(c) = elements(end) error(C S must be a standard wavelet . decomposition structure.); endnmax = size(s, 1) - 2; wname = varargin1; filterchk = 0; nchk = 0; switch nargincase 3 if ischar(wname) lp, hp = wavefilter(wname, r); n = nmax; else error(Undefined filter.); end if nargout = 1 error(Wron