В статье разработан метод экстраполяции функций 2k-гладких в некоторой окрестности границы прямоугольника на R2 n с сохранением гладкости порядка k и со свойством компактности носителя. Также, в данной работе с помощью численных экспериментов показано, что применение алгоритмов сжатия с потерями, основанных на всплесках, для таким способом экстраполированных изображений дает улучшение основных показателей сжатия по сравнению с существующими методами.
а затем уже для новой функции, определённой на [2a b; 2b a], проделать тоже самое еще раз и так далее. Также, техника ¾сгибания¿ применяется для набора всплеск функций f j;kgj;k2Z и набора масштабирующих функций f'j;kgj;k2Z. Однако, может оказаться, что данная процедура не приведет нас к ортонормированным базисам всплесков, поэтому ¾сгибания¿ используют в основном для биортогональных всплесков.
Также существуют другие способы, такие как использование ¾left edge¿ , ¾right edge¿ , ¾interior¿ всплесков (предложены Мейером в [4]) или схожего с методом Мейера способа, предложенного Коэном, Добеши и Виалом в [2] (¾interior¿ и ¾edge¿ всплески). Оба варианта отличаются громоздкостью требуемых построений.
У всех вышеперечисленных методов имеются недостатки, например, при продолжении исходной функции одним из способов, описанных во втором абзаце данного раздела, серьезным недостатком является потеря точности аппроксимации из-за в общем случае нарушения гладкости исходной функции (см. результаты численного эксперимента в Таб. 1. Целью данной работы является устранение недостатков приведенных выше методов за счет экстраполяции исходной функции 2k-гладкой в некоторой окрестности границы прямоугольника на R2 n с сохранением гладкости порядка k и со свойством компактности носителя. Также, в данной работе проведено сравнение результатов применения основанных на всплесках алгоритмов сжатия изображений с потерями к экстраполированным изображениям с результатами применения тех же алгоритмов к исходным изображениям. 2
Экстраполяция 8> g 2 Ck(([A1; A2] [B1; B2]) n ); < gx(i)(aj; y) = fx(i)(aj; y); gx(i)(Aj; y) = 0; y 2 [b1; b2]; i = 0; :::; k; j = 1; 2; >: gy(i)(x; bj) = fy(i)(x; bj); gy(i)(x; Bj) = 0; x 2 [a1; a2]; i = 0; :::; k; j = 1; 2: 2.1
Одномерный случай Покажем сначала, как произвольную функцию f = f (x); x 2 [c; d], обладающую односторонними производными вплоть до k-го порядка в точке d, продолжить на (d; e] гладко до порядка k и с условием f (e) = 0 с помощью многочлена Hn(x) подходящей степени, то есть так, чтобы выполнялись следующие условия: Hn(x) 2 Ck((d; e));
Hn(i)(d) = f (i)(d); Hn(i)(e) = 0; i = 0; 1; :::; k:
Лемма 1 Явный вид многочлена Hn(x), удовлетворяющего условиям (
" k H2k+1(x) = (x e)k+1 X f (i)(d) i=0 (x d)i " k i ##
X(
j=0
Доказательство Известно (как частный случай эрмитовой интерполяции), что задача нахождения
многочлена Hn(x) из условий (
e)k+1, а выражение n (!x2k+d)2k(+x)1 o[[dk] i] – есть совокупность первых членов
Таким образом, с помощью леммы 1 мы можем построить многочлен Hn = Hn(x); x 2 (d; e],
удовлетворяющий всем условиям из (
В качестве среды для программной реализации экстраполяции исходной функции и выполнения всплескразложения был выбран пакет прикладных программ M atlab версии 8.4.0.150421 (R2014b) с расширениями W avelet T oolbox версии 4.14 и W aveLab версии 0.850. Всплеск-разложение исходной функции f и её последующая аппроксимация всплесками Добеши реализованы с помощью соответственно функций wavedec и wrcoef из библиотеки W avelet T oolbox . Разложение с использованием периодизированных и ¾interior¿, ¾edge¿ всплесков исходной функции f и её последующая аппроксимация реализованы с помощью соответственно функций M akeON F ilter, F W T _P O, M akeOBF ilter, F W T _CDJ V и U pDyadLo, CDJ V DyadU p из библиотеки W aveLab.
Как видно из данных численного эксперимента, представленных в Таб. 1, с точки зрения погрешности приближения при аппроксимации с помощью всплесков метод предварительной гладкой экстраполяции исходной функции с помощью многочленов Эрмита дает результаты лучше, чем другие сравниваемые методы. Эффект потери точности аппроксимации из-за нарушения гладкости продолжения исходной функции можно наблюдать при других уровнях разложения и для других семейств всплесков. 2.2
Двумерный случай
Покажем, как построить g = g(x; y) из (
##
A2) (j+k+1)(x a2)j :
" k
B2)k+1 X fy(i)(x; b2)
(y b2)i " k i
i! X(
##
B2) (j+k+1)(y b2)j :
тогда
Из (
l;i+k+1(y iP=0 ecb2 k l;i+k+1(x iP=0 eca2 k
l;i+k+1(y
iP=0 ecb2
Перепишем соответствующие условия из (
A2)i = Ppk=0 h Pik=0 (fx(i))(yl)(a2;b2) (
i! k
B2)i = Ppk=0 h Pik=0 (fy(i))(xl)(a2;b2) (
A2) (p i+k+1)i(x a2)p; B2) (p i+k+1)i(y b2)p; 2k+1 ec0a;2i = P ci;j(b2 j=1 2k+1 ecb02;j = P ci;j(a2 i=1
B2)j; ecla;2i = 2kP+1 (j j!l)! ci;j(b2
j=l A2)i; eclb;2j = 2ikP=+l1 (i i!l)! ci;j(a2
B2)j l; i = 1; :::; 2k + 1; A2)i l; j = 1; :::; 2k + 1;
l = 1; :::; k; Py(l)(x; b2) = 2ikP=+11 ecla;2i(x Px(l)(a2; y) = 2jkP=+11 eclb;2j(y
A2)i; B2)j; Обозначим . . . . . . n BB 0
B n+1 B 0
1 0 ... . . .
n ::: 2 1 ::: (n 1) n ::: 1 (n + 1) ::: 2 В полученной матрице вычтем из (n + 1)-й строки n-ю строку, умноженную на n!, n BB 0
B n+1 B 0
n ::: 2 1 ::: (n 1) 0 n ::: 1 (k 1) ::: (k n+1) k n C
1 ::: (n 1) C (k 1) ::: (k n) k n C
C
A n . . . 1 . . . n 1 1 0 ... . . . ...
n 1 0 n+1 . . . . .
. n+1 . . . . . .
n ::: 2 1 ::: (n 1) n+1 . . . 1 . . .
. . . . . .
. . . . . .
k . . . . .
. . . . . . .
1 где c4;4; c4;5; c5;4; c5;5 – свободные переменные. 2 ec1a;23 + 1 2c1;3 ea2 1
c1;3 2 ec0a;23 + 4 c3;3; c2;3 = eb2
2 1 2c0;3 ea2 2 2 2 c3;3; c1;3 = 2ecb02;3 c0;3 eb2 2
2 c3;3; c1;3 + 2 eb2 c3;3; (15) (16) Доказательство утверждения 2 аналогично доказательству утверждения 1. Теорема Система (14) совместна и ее решение имеет вид: при k = 0: c1;1 = eb2 ;
c0;1
0 2k+1 . . . . . .
j
Aj 1( ) j 11 Aj 2( ) . .
. (j 1) ::: (j k) Aj k 1( ) 1 ::: k
Bj( ) . .
.
Bj(k)( ) k+1
Ak( ) k Ak 1( ) 1 . . .
E 0 . . . 0
k+2 Ak+1( ) k+11 Ak( ) . .
. (k+1) ::: 2 A( ) 1 ::: k Bk+2( ) . .
.
Bk(k+)2( ) . . . . . . . . . . . .
2k+1
A2k( ) 1 21k A2k 1( ) C
... CCC 2k ::: (k+1) Ak( ) CCC 1 ::: k B2k+1( ) . . .
C C C C
A Несложно видеть, что ранг данной матрицы равен (k + 1)(2k + 1) + (k + 1)(2k + 1 k 1), т.е. (k + 1)(3k + 1). Для того, чтобы система (14) была совместна необходимо и достаточно, чтобы rank(Mk) = rank(Mk; eck). Найдем rank(Mk; ck). Из леммы 2 следует, что преобразования над матрицей Mk действуют на столбец e свободных членов из (14) следующим образом: (ecb02; :::; ecib2 1; :::; ecbk2)T
c0 i 1
( eb2 ; :::; X(
1 (i 1 j)!
k
j 1ci 1 j; :::; X(
1 (k j)! j 1ecbk2 j)T ; ecla2 = (0; 0; :::; 0; :::; 0; ecla;2k+1; cl;k+2; :::; cl;2k+1)T ea2 ea2 Докажем, что ec0a;2k+1 k + ::: + c0;2k+1 2k ea2 c0;k+1 k+1 + ::: + ecb02;2k+1 2k+1 eb2 = 0: (17) 0 1 B
B B
B 2 B
B . B .. B
B
B 2k+1 cl;i[ (i 1) ::: (i m+1) i m] m BBBB i=Pk+1 ea2 1 ::: (m 1) .. B . B
B
B k+1 BB
B k+2 BB . B . B . @ 2k+1 2k+1 cl;i[ (i 1) ::: (i k) i k 1] i=Pk+1 ea2 1 ::: k ecla;2k+1 k + ::: + cl;2k+1 2k
ea2 ecla;2k+1[ k1 k 1] + ::: + ecla;22k+1[ 21k 2k 1] j2=kPk++11 ecb12;1j (jj!l)! j l + j2=kPk++11 ecb02;j (2jj!l)! j l 2k+1 c0;j j=k+1 eb2 (jj!l)! j l
P
m 1 1
P (
C j l] CC
C; C C C C C C C C C C C C
A k ec0a;2k+1 k+1 + ::: + c0;2k+1 2k+1 = k+1 X Sa02;p(A2 a2)p 0 p ea2 0 k + ::: + 2k+1 X Sa02;p(A2 a2)p k p =
k p=k p=0 Рис. 1: Пример экстраполяции функции f (x; y) = 2 exp (x4) cos y + exp (y3) sin x + 13; (x; y) 2 [0; 1] [0; 1]
Применение экстраполяции для сжатия изображений с потерями
В качестве среды для программной реализации экстраполяции и применения алгоритмов сжатия
изображений на основе аппарата теории всплесков был выбран пакет прикладных программ Matlab с
расширением Wavelet Toolbox указанных во втором разделе версий. Пусть у нас есть цифровое растровое
изображение, которое можно рассматривать как функцию двух переменных X = X(i; j); (i; j) 2 Z2,
определенную в точках некоторого прямоугольника. Мы будем рассматривать только монохромные (с
глубиной цвета 8 бит на пиксель) изображения размера 2l 2l (l 2 N) пикселей, определенные в
целых точках [0; 2l 1] [0; 2l 1]. Экстраполируем исходное изображение X размера 2l 2l пикселей,
определенное на [0; 2l 1] [0; 2l 1], до изображения размера 2l+1 2l+1 пикселей, определенное на
[ 2l 1; 2l + 2l 1 1] [ 2l 1; 2l + 2l 1 1], с требуемой гладкостью k с помощью формул (
Рис. 2: Пример экстраполяции для изображения Lena.png После предварительной обработки исходного изображения будем применять к X и XE четыре основанных на всплесках алгоритма сжатия:
EZW (Embedded Zerotree W avelet, построил Шапиро в работе [6]); SP IHT (Set P artitioning in Hierarchical T rees, построили Саид и Перлман [7]); ST W (Spatial-orientation T ree W avelet, построили Саид и Перлман [8]);
W DR (W avelet Dif f erence Reduction, построили Тянь и Уэллс [9], [10], [11]). Отметим, что алгоритмы SP IHT и ST W являются измененными версиями алгоритма EZW . Как правило, в качестве характеристик, по которым сравниваются алгоритмы сжатия, выбирают показатели CR и P SN R:
CR (Compression Ratio) = SS0c (измеряется в процентах), где S0 – объём исходных данных, а Sc – объём сжатых данных; PSNR (P eak Signal-to-N oise Ratio) = 10 log10 M25S5E2 (измеряется в децибелах), где MSE (M ean Squared m n Error) = m1n iP=1 j=1
P jX(i; j) Xcomp(i; j)j2, X = X(i; j) – исходное изображение, Xcomp = Xcomp(i; j) – восстановленное после сжатия изображение, m
n – размер изображений X и Xcomp. Однако, так как само XE изначально нигде не хранится и содержит незначимые пиксели, нужные лишь в процессе сжатия и отбрасываемые после восстановления изображения, то показатель CR не подходит для сравнения алгоритмов сжатия изображений с предварительной обработкой и без неё. Поэтому вместо CR мы будем использовать показатель F S (F ile Size) – размер файла сжатого изображения с предварительной обработкой (F SE) и без предварительной обработки (F S). MSE для алгоритмов с предварительной обработкой изображений будем вычислять следующим образом: MSEE := 2l1+1 P2l P2l jX(i; j) XEcomp(i; j)j2, i=1 j=1 где XEcomp – изображение, полученное восстановлением после сжатия изображения XE и взятия только той его части, которая задана на [0; 2l 1] [0; 2l 1] (то есть изображения 2l 2l пикселей). Для алгоритмов 2l 2l без предварительной обработки MSE := 2l1+1 P P jX(i; j) Xcomp(i; j)j2. Тогда PSNRE = 10 log10 M2S55E2E , i=1 j=1 PSNR = 10 log10 M25S5E2 .
Для того, чтобы сравнить перечисленные алгоритмы сжатия в двух вариациях: без предварительной обработки и с предварительной обработкой, будем подбирать такие значения параметров в самих алгоритмах, чтобы при FSE FS вариант с предварительной обработкой давал большее PSNR (то есть лучшее качество изображения). Все четыре алгоритма сжатия реализованы с помощью функции wcompress из библиотеки Wavelet Toolbox. Основным параметром функции wcompress является максимальное число итераций соответствующего алгоритма сжатия – maxloop. Его увеличение ведет к улучшению качества изображения (к большему PSNR), но к увеличению размера файла изображения (к большему FS). Еще один важный параметр – level, который определяет уровень, до которого происходит всплеск-разложение.
Все численные эксперименты проведены для монохромного изображения Lena:png размера 512 512 пикселей с размером файла 167034 байт, а также для его экстраполированной (k = 1) версии размера 1024 1024. В приведенных алгоритмах сжатия на этапе всплеск-разложения были использованы все базисы всплесков доступные в расширении Wavelet Toolbox: дискретные всплески Мейера (dmey), всплески Добеши (db1 – db45), койфлеты (coif 1 – coif 5), симмлеты (sym1 – sym45), биортогональные всплески (bior1:1, bior1:3, bior1:5, bior2:2, bior2:4, bior2:6, bior2:8, bior3:1, bior3:3, bior3:5, bior3:7, bior3:9, bior4:4, bior5:5, bior6:8), обратные биортогональные всплески (rbio1:1, rbio1:3, rbio1:5, rbio2:2, rbio2:4, rbio2:6, rbio2:8, rbio3:1, rbio3:3, rbio3:5, rbio3:7, rbio3:9, rbio4:4, rbio5:5, rbio6:8). Далее будут приведены результаты только для тех базисов всплесков, которые дали улучшение в выбранных показателях сжатия в варианте с предварительной обработкой изображения.
В Таб. 2 приведены результаты численных экспериментов для алгоритма сжатия EZW (левая часть таблицы – вариант без предварительной обработки, правая часть – с предварительной обработкой). Для дискретных всплесков Мейера (dmey) всплеск-разложение проводилось до уровня level = 5, число итераций алгоритма сжатия без предварительной обработки maxloop = 12, с предварительной обработкой maxloopE = 11. Для остальных базисов всплесков (db3, db4, db10, db12, coif5, sym3, sym4, sym8, sym10, sym12, sym13, sym23, sym25, sym27, bior4.4, bior5.5, rbio4.4) level = 4, maxloop = 11, maxloopE = 10.
В Таб. 3 приведены результаты численных экспериментов для алгоритма сжатия SPIHT (левая часть таблицы – вариант без предварительной обработки, правая часть – с предварительной обработкой). Для дискретных всплесков Мейера (dmey) level = 5, maxloop = 13, maxloopE = 12. Для остальных базисов всплесков (db3, db4, db10, db12, coif5, sym3, sym4, sym8, sym10, sym12, sym13, sym23, sym25, bior4.4, bior5.5, rbio4.4) level = 4, maxloop = 12, maxloopE = 11.
В Таб. 4 приведены результаты численных экспериментов для алгоритма сжатия WDR (левая часть таблицы – вариант без предварительной обработки, правая часть – с предварительной обработкой). Для дискретных всплесков Мейера (dmey) level = 5, maxloop = 11, maxloopE = 10. Для остальных базисов всплесков (coif5, bior4.4, bior5.5) level = 4, maxloop = 10, maxloopE = 9.
Для алгоритма сжатия STW ни один из базисов всплесков при всевозможных параметрах не дал положительных результатов в варианте сжатия с предварительной обработкой изображения. Из Таб. 2, 3, 4 видно, что наибольшее значение P SN RE P SN R для алгоритмов EZW, SPIHT, WDR получено с использованием биортогональных всплесков bior5.5, а наибольшее значение F SE F S для тех же алгоритмов получено с использованием дискретных всплесков Мейера dmey. Таблица 2: Результаты численных экспериментов для алгоритма сжатия EZW
Показатели сжатия F SE (байт) P SN RE (дБ) 63700 38.8236 61629 38.6860 61027 38.6562 60774 38.2070
Работа выполнена за счет гранта Российского научного фонда (проект 14-11-00702). Список литературы
Society for Industrial and Applied Mathematics (SIAM), [11] J. Tian, R.O. Wells, Jr. Image data processing in the compressed wavelet domain. 3rd International Conference on Signal Processing Proc., Yuan, B. and Tang, X., Eds., Beijing, China, 978–981, 1996.
Image preprocessing for improving performance of lossy compression algorithms based on wavelets
Dmitriy A. Yamkovoy Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia) Keywords: smooth extrapolation, image compression, lossy compression.
In this paper we developed extrapolation method of 2k-smooth functions in a neighborhood of the rectangle to the R2 n with k-smooth-preserving and finite support properties. Also, advantages of new construction applying to lossy image compression using wavelet theory in comparison with other methods are discussed according to basic characteristics of image compression.