Предварительная обработка изображений
для улучшения качества работы алгоритмов сжатия
с потерями, основанных на всплесках
Д.А. Ямковой
dmitriiyamkovoi@bk.ru
ИММ УрО РАН (Екатеринбург)
Аннотация
В статье разработан метод экстраполяции функций 2k-гладких в
некоторой окрестности границы прямоугольника Ω на R2 \ Ω с со-
хранением гладкости порядка k и со свойством компактности носи-
теля. Также, в данной работе с помощью численных экспериментов
показано, что применение алгоритмов сжатия с потерями, основан-
ных на всплесках, для таким способом экстраполированных изобра-
жений дает улучшение основных показателей сжатия по сравнению
с существующими методами.
1 Введение
Рассмотрим задачу аппроксимации функций, заданных на отрезке (в данном разделе будем в основном
говорить только об одномерном случае, так как в двумерном случае возникают те же проблемы) с помо-
щью аппарата теории всплесков. Такая задача может возникнуть, например, при сжатии и дальнейшем
восстановлении некоторых достаточно гладких данных. При решении поставленной задачи, однако, воз-
никают некоторые трудности – аппарат классических базисных масштабирующих функций пространств
кратномасштабного анализа (КМА) (см., например, [1], гл. 5), используемый для аппроксимации, пло-
хо приспособлен для приближения функций, определенных на отрезке. Существуют различные способы,
описанные в [2], адаптации базиса L2 (R) к случаю L2 ([a, b]).
Один из способов решения данной проблемы заключается в том, чтобы продолжить исходную функцию
нулем на всю вещественную ось и использовать стандартный аппарат КМА. Также, можно сначала про-
должать исходную функцию на более широкий промежуток симметрично, периодически или непрерывно,
а затем уже положить функцию нулем вне новой области определения.
Довольно часто применяют периодизированные всплески (см., например, [1], гл. 9, §9.3). Суть дан-
ного метода заключается в том, чтобы периодизировать каждую функцию из набора всплеск функций
{ψj,k }j,k∈Z и набора масштабирующих функций {ϕj,k }j,k∈Z , получив таким образом ортонормированные
базисы подпространств Vjпер пространства j∈Z Vjпер = L2 ([a, b]).
S
Еще один способ решить возникающую проблему – метод так называемых «folding» всплесков, опреде-
ленных в [3]. Основная идея состоит в том, чтобы исходную функцию f симметрично «согнуть» относи-
Copyright c by the paper’s authors. Copying permitted for private and academic purposes.
In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the International Youth School-conference «SoProMat-2017», Yekaterinburg,
Russia, 06-Feb-2017, published at http://ceur-ws.org
183
тельно правого и левого края промежутка [a, b]:
f (−x + 2b), x ∈ [b, 2b − a],
f (x) =
f (−x + 2a), x ∈ [2a − b, a],
а затем уже для новой функции, определённой на [2a − b, 2b − a], проделать тоже самое еще раз и так
далее. Также, техника «сгибания» применяется для набора всплеск функций {ψj,k }j,k∈Z и набора мас-
штабирующих функций {ϕj,k }j,k∈Z . Однако, может оказаться, что данная процедура не приведет нас к
ортонормированным базисам всплесков, поэтому «сгибания» используют в основном для биортогональ-
ных всплесков.
Также существуют другие способы, такие как использование «left edge» , «right edge» , «interior» всплес-
ков (предложены Мейером в [4]) или схожего с методом Мейера способа, предложенного Коэном, Добеши
и Виалом в [2] («interior» и «edge» всплески). Оба варианта отличаются громоздкостью требуемых постро-
ений.
У всех вышеперечисленных методов имеются недостатки, например, при продолжении исходной функ-
ции одним из способов, описанных во втором абзаце данного раздела, серьезным недостатком является
потеря точности аппроксимации из-за в общем случае нарушения гладкости исходной функции (см. ре-
зультаты численного эксперимента в Таб. 1. Целью данной работы является устранение недостатков при-
веденных выше методов за счет экстраполяции исходной функции 2k-гладкой в некоторой окрестности
границы прямоугольника Ω на R2 \ Ω с сохранением гладкости порядка k и со свойством компактности
носителя. Также, в данной работе проведено сравнение результатов применения основанных на всплес-
ках алгоритмов сжатия изображений с потерями к экстраполированным изображениям с результатами
применения тех же алгоритмов к исходным изображениям.
2 Экстраполяция
Пусть f : Ω → R – 2k-гладкая функция в некоторой окрестности границы Ω (то есть существуют и
k
непрерывны все частные производные ∂xk∂1 ∂y f
k2 функции f = f (x, y), где k ∈ {0} ∪ N, k1 , k2 ∈ [0, k] ∩ Z,
k1 + k2 = k), Ω = [a1 , a2 ] × [b1 , b2 ]. Пусть A1 < a1 < a2 < A2 , B1 < b1 < b2 < B2 . Требуется построить
функцию
f (x, y), (x, y) ∈ Ω
F (x, y) := g(x, y), (x, y) ∈ ([A1 , A2 ] × [B1 , B2 ]) \ Ω (1)
0, (x, y) 6∈ [A1 , A2 ] × [B1 , B2 ]
такую, чтобы выполнялись условия
k
g ∈ C (([A1 , A2 ] × [B1 , B2 ]) \ Ω),
(i) (i) (i)
gx (aj , y) = fx (aj , y), gx (Aj , y) = 0, y ∈ [b1 , b2 ], i = 0, ..., k, j = 1, 2, (2)
g (i) (x, b ) = f (i) (x, b ), g (i) (x, B ) = 0, x ∈ [a , a ], i = 0, ..., k, j = 1, 2.
y j y j y j 1 2
2.1 Одномерный случай
Покажем сначала, как произвольную функцию f = f (x), x ∈ [c, d], обладающую односторонними произ-
водными вплоть до k-го порядка в точке d, продолжить на (d, e] гладко до порядка k и с условием f (e) = 0
с помощью многочлена Hn (x) подходящей степени, то есть так, чтобы выполнялись следующие условия:
Hn (x) ∈ C k ((d, e)),
(i) (i) (3)
Hn (d) = f (i) (d), Hn (e) = 0, i = 0, 1, ..., k.
Лемма 1 Явный вид многочлена Hn (x), удовлетворяющего условиям (3), можно получить с помощью
следующей формулы:
" k " k−i ##
i X
X (x − d) k + j
H2k+1 (x) = (x − e)k+1 f (i) (d) (−1)j (d − e)−(j+k+1) (x − d)j . (4)
i=0
i! j=0
k
184
Доказательство Известно (как частный случай эрмитовой интерполяции), что задача нахождения мно-
гочлена Hn (x) из условий (3) разрешима и, при этом имеет единственное решение (см., например [5], гл. 2,
§2.7), если n = 2k + 1. Hn (x) – многочлен степени n = 2k + 1, удовлетворяющий условиям (3) и называемый
интерполяционным многочленом Эрмита, имеет вид:
" k [k−i] #
ω2k+2 (x) X (i) (x − d)i (x − d)k+1
H2k+1 (x) = f (d) , (5)
(x − d)k+1 i=0 i! ω2k+2 (x) [d]
n o[k−i]
(x−d)k+1
где ω2k+2 (x) = (x − d)k+1 (x − e)k+1 , а выражение ω2k+2 (x) – есть совокупность первых членов
[d]
k+1
разложения в ряд Тейлора по степеням (x−d) функции (x−d)
ω2k+2 (x) в точке d до (k −i)-го члена включительно.
Путем несложных выкладок (5) преобразуется к виду (4).
Таким образом, с помощью леммы 1 мы можем построить многочлен Hn = Hn (x), x ∈ (d, e], удо-
влетворяющий всем условиям из (3). В качестве численного эксперимента экстраполируем с сохранением
гладкости порядка k = 3 функцию α(x) = 3x2 · cos x + 2x · sin x + 5 с отрезка [0, 4] на всю ось в виде функции
с носителем [−2, 6] и проведем сначала всплеск-разложение так экстраполированной
R∞ функции до уровня 5
с помощью всплесков семейства Добеши с 15-ю нулевыми моментами −∞ xn ψ(x)dx = 0, n = 0, 1, ..., 14 (это
соответствует длине носителей масштабирующей функции ϕ и всплеска ψ равной 29). Далее, применим
разложение до того же уровня по тем же всплескам к функции α, продолженной описанными во втором
абзаце введения способами. Применим также к исходной функции описанную процедуру, но основанную
на периодизированных и «interior», «edge» всплесках, упомянутых в третьем и пятом абзацах введения.
Будем сравнивать время выполнения (в секундах) всплеск-разложения до выбранного уровня и погреш-
ности приближения на отрезке [0, 4] в нормах C и L2 при аппроксимации с помощью всплесков Добеши с
использованием указанных выше способов продолжения исходной функции, а также периодизированных
и «interior», «edge» всплесков.
В качестве среды для программной реализации экстраполяции исходной функции и выполнения всплеск-
разложения был выбран пакет прикладных программ 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 _CDJV и U pDyadLo, CDJV DyadU p
из библиотеки W aveLab.
Как видно из данных численного эксперимента, представленных в Таб. 1, с точки зрения погрешности
приближения при аппроксимации с помощью всплесков метод предварительной гладкой экстраполяции ис-
ходной функции с помощью многочленов Эрмита дает результаты лучше, чем другие сравниваемые мето-
ды. Эффект потери точности аппроксимации из-за нарушения гладкости продолжения исходной функции
можно наблюдать при других уровнях разложения и для других семейств всплесков.
Таблица 1: Результаты описанного численного эксперимента в одномерном случае
Числовые характеристики
Продолжение Время(с) Погрешность в норме C Погрешность в норме L2
1. Нулем 0.0144 0.1680 · 102 0.4645
2. Периодическое 0.0128 0.1780 · 102 0.4671
3. Симметричное 0.0148 0.1613 · 10−1 0.3119 · 10−3
4. Непрерывное 0.0142 0.1595 · 10−1 0.2864 · 10−3
−4
5. Эрмитово (k = 3) 0.0150 0.1558 · 10 0.4121 · 10−6
Числовые характеристики
Всплески Время(с) Погрешность в норме C Погрешность в норме L2
1. Периодизированные 0.0048 0.1899 · 102 0.7407
2. «interior», «edge» 0.0825 0.2378 · 10−2 0.2421 · 10−3
185
2.2 Двумерный случай
Покажем, как построить g = g(x, y) из (1) на (a2 , A2 ] × [b1 , b2 ] и [a1 , a2 ] × (b2 , B2 ] так, чтобы выполнялись
условия (2). На [A1 , a1 ) × [b1 , b2 ] и [a1 , a2 ] × [B1 , b1 ) продолжение будет строится аналогично. Используем
формулу (4) из леммы 1 и построим g на (a2 , A2 ] × [b1 , b2 ] с помощью функции Ha2 (x, y):
" k " k−i ##
(x − a2 )i X
j k+j
X
k+1 (i) −(j+k+1) j
Ha2 (x, y) = (x − A2 ) fx (a2 , y) (−1) (a2 − A2 ) (x − a2 ) . (6)
i=0
i! j=0
k
Используя аналог формулы (4) из леммы 1 и построим g на [a1 , a2 ] × (b2 , B2 ] с помощью функции Hb2 (x, y):
" k " k−i ##
(y − b2 )i X
j k+j
X
k+1 (i) −(j+k+1) j
Hb2 (x, y) = (y − B2 ) fy (x, b2 ) (−1) (b2 − B2 ) (y − b2 ) . (7)
i=0
i! j=0
k
Очевидно, что так построенная функция удовлетворяет условиям (2).
Покажем теперь, как построить g на (a2 , A2 ] × (b2 , B2 ] в виде многочлена
2k+1
X 2k+1
X
P (x, y) = ci,j (x − A2 )i (y − B2 )j , (x, y) ∈ [a2 , A2 ] × [b2 , B2 ],
i=1 j=1
подчиненного условиям:
(
(i) (i) (i)
Py (x, B2 ) = 0, Py (x, b2 ) = (Ha2 )y (x, b2 ), x ∈ (a2 , A2 ],
(i) (i) (i) i = 0, ..., k, (8)
Px (A2 , y) = 0, Px (a2 , y) = (Hb2 )x (a2 , y), y ∈ (b2 , B2 ],
где ci,j – неизвестные коэффициенты. На [A1 , a1 )×[B1 , b1 ) и [A1 , a1 )×(b2 , B2 ], (a2 , A2 ]×[B1 , b1 ) продолжение
строится аналогично. Покажем, как найти все ci,j . Обозначим
2k+1 2k+1
c0,i cl,i
ci,j (b2 − B2 )j , e j! j−l
P P
ea2 = a2 = (j−l)! ci,j (b2 − B2 ) , i = 1, ..., 2k + 1,
j=1 j=l
2k+1 2k+1
l = 1, ..., k, (9)
c0,j cl,j
ci,j (a2 − A2 )i , e i! i−l
P P
eb2 = b2 = (i−l)! ci,j (a2 − A2 ) , j = 1, ..., 2k + 1,
i=1 i=l
тогда
2k+1
(l)
cl,i i
P
Py (x, b2 ) = ea2 (x − A2 ) ,
i=1
2k+1 l = 0, ..., k. (10)
(l)
cl,j j
P
Px (a2 , y) = eb2 (y − B2 ) ,
j=1
Из (6)–(10) можно видеть, что
cl,i
a2 = 0, i = 1, ..., k, l = 0, ..., k. (11)
e
cl,j
eb2 = 0, j = 1, ..., k,
Перепишем соответствующие условия из (8) с учетом (11)
k Pk (fx(i) )y
(l)
(a2 ,b2 )
hP i
l,i+k+1 k−i j k+j
−(j+k+1)
(x − A2 )i = (x − a2 )i+j
P
ca2
e i=0 i! j=0 (−1) k (a2 − A2 ) ,
i=0
k Pk (fy(i) )x
(l)
(a2 ,b2 )
hP i
k−i
cl,i+k+1 j k+j −(j+k+1)
(y − B2 )i = (y − b2 )i+j
P
eb2 i=0 i! j=0 (−1) k (b2 − B2 ) ,
i=0
l = 0, ..., k.
В правых частях двух последних уравнений сделаем замену j = p − i, поменяем порядок суммирования и
получим
k Pk h Pk (fx(i) )(l) i
y (a2 ,b2 )
cl,i+k+1 p−i k+p−i −(p−i+k+1)
(x − A2 )i = p=0 (x − a2 )p ,
P
ea2 i=0 i! (−1) k (a 2 − A 2 )
i=0
k Pk h Pk (fy(i) )x(l) (a2 ,b2 ) i (12)
cl,i+k+1 p−i k+p−i −(p−i+k+1)
i
(y − b2 )p ,
P
e b2 (y − B 2 ) = p=0 i=0 i! (−1) k (b 2 − B 2 )
i=0
l = 0, ..., k.
186
Обозначим
l,p (fx(i) )(l)
Pk y (a2 ,b2 )
(−1)p−i k+p−i (a2 − A2 )−(p−i+k+1) ,
Sa2 = i=0 i! k
(i) (l)
l,p k (f ) (a2 ,b2 )
(−1)p−i k+p−i (b2 − B2 )−(p−i+k+1) ,
P
Sb2 = i=0 y xi! k
p = 0, ..., k, l = 0, ..., k.
Тогда дифференцируя соответствующие выражения из (12) в точках A2 и B2 соответственно получим
l,i+k+1 Pk l,p
(A2 − a2 )p−i pi ,
ca2
e = p=i Sa2
l,i+k+1 Pk l,p
(B2 − b2 )p−i pi ,
cb2
e = p=i Sb2 (13)
i = 0, ..., k, l = 0, ..., k.
cia2 , e
Обозначим α = a2 − A2 , β = b2 − B2 . Пусть e cib2 , i = 0, ..., k – вектор-столбцы, координаты которых
вычислены по формулам (11) и (13), тогда мы можем переписать (9) в виде системы линейных уравнений
Mk · c = e
ck (14)
относительно столбца неизвестных c, где
A(α) A2 (α) A3 (α) ··· A2k+1 (α)
.. .. .. .. ..
.
. . . .
A(k) (α) (A2 (α))(k) (A3 (α))(k) (A2k+1 (α))(k)
···
Mk =
B1 (β) B2 (β) B3 (β) ··· B2k+1 (β)
.. .. .. ..
..
. . . . .
(k) (k) (k) (k)
B1 (β) B2 (β) B3 (β) ··· B2k+1 (β)
– матрица размерности 2(k + 1)(2k + 1) × (2k + 1)2 ; A(α) = diag{α, α ..., α} – диагональная матрица раз-
мерности (2k + 1) × (2k + 1); Bi (β) – матрица размерности (2k + 1) × (2k + 1):
(β β 2 · · · β 2k+1 ), j = i,
(Bi (β))(j, ) =
(0 0 · · · 0), j 6= i,
c = (c1,1 , ..., c1,2k+1 , ..., c2k+1,1 , ..., c2k+1,2k+1 )T , e c0b2 )T , ..., (e
ck = ((e ckb2 )T , (e
c0a2 )T , ..., (e
cka2 )T )T . Решив систему (14),
мы определим все неизвестные коэффициенты ci,j многочлена P (x, y).
Лемма 2 Справедливы следующие утверждения.
1. Матрицу m × k
1 2 3 ··· j ··· k
α2 α3 αj αk
α ··· ···
1
2α 3α2 ··· jαj−1 ··· kαk−1
0
2·1 3 · 2α ··· j(j − 1)αj−2 ··· k(k − 1)αk−2
.. .. .. .. .. .. ..
. . . . . . .
0 0 0 ··· j · ... · (j − m + 2)αj−m+1 ··· k · ... · (k − m + 2)αk−m+1
эквивалентными преобразованиями можно привести к виду
1 2 3 ··· j ··· k
α2 αj−1 αk−1
1 α ··· ···
j−1 j−2 k−1 k−2
0
1 2α ··· 1 α ··· 1 α
(j−1)(j−2) j−3 (k−1)(k−2) k−3
0
0 1 ··· 1·2 α ··· 1·2 α
. .. .. .. .. .. ..
. . .
. . . . .
(j−1)·...·(j−m+1) j−m (k−1)·...·(k−m+1) k−m
0 0 0 ··· 1·...·(m−1) α ··· 1·...·(m−1) α
2. Вектор-столбец (c1 , ..., ci , ..., cm )T преобразованиями, осуществляющими пункт 1, приводится к виду
i−1 m−1
( cα1 , ..., 1
(−1)j (i−1−j)! α−j−1 ci−j , ..., 1
(−1)j (m−1−j)! α−j−1 cm−j )T .
P P
j=0 j=0
187
Доказательство Докажем утверждение 1 с помощью индукции по n – количество первых строк исходной
матрицы (n < m).
База индукции при n = 1 очевидна. Предположим, что исходную матрицу можно эквивалентными преоб-
разованиями привести к искомому виду
1 ··· n n+1 ··· k
. .. .. ..
.. . . . ..
. . . .
(k−1)·...·(k−n+1) k−n
n·...·2
n 0 ··· 1 1·...·(n−1) α ··· 1·...·(n−1) α
.
k−n
n+1 0 · · · n · ... · 1 (n + 1) · ... · 2α ··· k · ... · (k − n + 1)α
.. . . .. .. .. ..
. . . . . .
В полученной матрице вычтем из (n + 1)-й строки n-ю строку, умноженную на n!,
1 ··· n n+1 ··· k
. .. .. ..
.. ..
.. . . . . .
(k−1)·...·(k−n+1) k−n
n·...·2
n 0
··· 1 1·...·(n−1) α ··· 1·...·(n−1) α
∼
k−n
n+1 0
··· 0 n · ... · 1α ··· (k − 1) · ... · (k − n)α
.. .. .. .. .. ..
. . . . . .
1 ··· n n+1 ··· k
.. .. .. ..
.. ..
. . . . . .
n·...·2 (k−1)·...·(k−n+1) k−n
n 0
··· 1 1·...·(n−1) α ··· 1·...·(n−1) α
.
(k−1)·...·(k−n) k−n−1
n+1 0 ··· 0 1 ··· 1·...·n α
.. .. .. ..
.. ..
. . . . . .
Доказательство утверждения 2 аналогично доказательству утверждения 1.
Теорема Система (14) совместна и ее решение имеет вид:
c0,1
при k = 0: c1,1 = bα2 ;
e
при k = 1:
c1,3 c0,3 c0,3
2e
a2
c3,2 = 1·β − βa22 − 2βc3,3 , c3,1 = βa2 − ec1,3 2
e e
a2 + β c3,3 ,
1,2
c0,2 c1,3 c0,3
0,3
4αec c 2αec1,3 2αec0,3
c2,1 = − βa2 + 2αe c1,3 2 b2 b2 a2
a2 − 2αβ c3,3 , c2,2 = 1·α − α2 − 1·β + β 2
a2 b2
− αb22 − 2αc3,3 ,
e e e e
+ 4αβc3,3 , c2,3 = 1·α
c0,2 c0,3
2 0,3
c = 2α eca2 − α2 e 2e α2 e
c1,3 α2 e
c0,3 2e
c1,3 + α2 β 2 c , c = b2 − e c1,2 + a2 a2
− 2α2 βc , c = b2 − e c1,3 + α2 c ,
1,1 β a2 3,3 1,2 − α b2 1·β β2 3,3 1,3 α b2 3,3
(15)
где c3,3 – свободная переменная;
при k = 2:
2,5 1,5
c0,5
ca ca a2 3·2 4·3
c5,3 = 2·1·β − 1·β 2 + β 3 − 2·1 c5,4 − 2·1 c5,5 ,
e e e
2 2
1,5 0,5
ca c
c5,2 = 1·β − βa22 − 2βc5,3 − 3β 2 c5,4 − 4β 3 c5,5 ,
e e
2
0,5
c
c5,1 = aβ2 − βc5,2 − β 2 c5,3 − β 3 c5,4 − β 4 c5,5 ,
e
2,4 1,4 0,4
c4,3 = eca2 − eca22 + eca32 − 3·2 c4,4 − 4·3 c4,5 ,
2·1·β 1·β β 2·1 2·1
1,4 0,4
ca ca
c4,2 = 1·β − β 2 − 2βc4,3 − 3β c4,4 − 4β 3 c4,5 ,
e
2
e
2 2 (16)
c0,4
c4,1 = aβ2 − βc4,2 − β 2 c4,3 − β 3 c4,4 − β 4 c4,5 ,
e
c2,i c1,i c0,i
3·2 4·3 2
b2 b2 b2
− − −
e e e
c3,i = 2·1 αc4,i 2·1 α c5,i + 2·1·α 1·α 2 + α3 , i = 1, ..., 5,
1,i 0,i
cb c
c2,i = −2αc3,i − 3α2 c4,i − 4α3 c5,i + 1·α − αb22 , i = 1, ..., 5,
e
2
e
c0,i
c1,i = −αc2,i − α2 c3,i − α3 c4,i − α4 c5,i + bα2 , i = 1, ..., 5,
e
где c4,4 , c4,5 , c5,4 , c5,5 – свободные переменные.
188
Доказательство С помощью леммы 2 эквивалентными преобразованиями приведем матрицу Mk из (14)
к виду
1 ··· j ··· 2k+1
Aj−1 (α) A2k (α)
1 E ··· ···
j−1 j−2 2k 2k−1
2
0 ··· 1 A (α) ··· 1 A (α)
.. .. .. .. .. ..
.
. . . . .
(j−1)·...·(j−k) j−k−1 2k·...·(k+1) k
0 ··· A (α) ··· A (α) ∼
k+1 1·...·k 1·...·k
k+2 B1 (β) ··· Bj (β) ··· B2k+1 (β)
.. .
.. .. .. .. ..
. . . . .
(k) (k) (k)
2k+1 B1 (β) · · · Bj (β) ··· B2k+1 (β)
··· ···
1 k+1 k+2 2k+1
1 E ··· Ak (α) Ak+1 (α) ··· A2k (α)
k k−1 k+1 k 2k 2k−1
2 0
··· 1A (α) 1 A (α) ··· 1 A (α)
.. .. .. .. .. .. ..
. . . . . . .
(k+1)·...·2 2k·...·(k+1) k
k+1 0 ··· E 1·...·k A(α) ··· A (α) .
1·...·k
k+2 0 ··· 0 Bk+2 (β) ··· B2k+1 (β)
..
. .. .. .. .. ..
. .. . . . . .
(k) (k)
2k+1 0 ··· 0 Bk+2 (β) ··· B2k+1 (β)
Несложно видеть, что ранг данной матрицы равен (k + 1)(2k + 1) + (k + 1)(2k + 1 − k − 1), т.е. (k + 1)(3k + 1).
Для того, чтобы система (14) была совместна необходимо и достаточно, чтобы rank(Mk ) = rank(Mk , e ck ).
Найдем rank(Mk , e
ck ). Из леммы 2 следует, что преобразования над матрицей Mk действуют на столбец
свободных членов из (14) следующим образом:
i−1 k
c0b2 X 1 −j−1 i−1−j
X 1
ci−1
c0b2 , ..., e k T j
(−1)j ck−j
α−j−1 e T
e
(e b2 , ..., c b2 ) ∼ ( , ..., (−1) α c b2 , ..., b2 ) ,
− − −
e e
α j=0
(i 1 j)! j=0
(k j)!
cla2 = (0, 0, ..., 0, ..., 0, e
e cl,k+1
a2 cl,k+2
,ea2 cl,2k+1
, ..., ea2 )T ∼
2k+1
c0,j j! j−l
P
eb (j−l)! β
2
j=k+1
1
cl,k+1
e a2 αk + ... + e cl,2k+1
a2 α2k − α
2k+1 2k+1
c1,j j! j−l 0,j j! j−l
P P
eb2 (j−l)! β c
eb2 (j−l)! β
k k−1 l,2k+1 2k 2k−1 j=k+1 j=k+1
cl,k+1 −
2 ea2 [ 1 α ] + ... + c
ea2 [ 1 α ] 1·α + α 2
.. ..
.
.
2k+1
P l,i (i−1)·...·(i−m+1) i−m m−1
P i 1 −i−1
2k+1
P m−1−i,j j! j−l
m c
ea2 [ 1·...·(m−1) α ] − (−1) (m−1−i)! α [ c
e b2 (j−l)! β ]
,
i=k+1 i=0 j=k+1
.. ..
. .
2k+1
P l,i (i−1)·...·(i−k) i−k−1 k 2k+1
P i 1 −i−1
P k−i,j j! j−l
ca2 [ α ] − (−1) α [ c β ]
k+1 e 1·...·k (k−i)! e b (j−l)!
2
i=k+1 i=0 j=k+1
l,k+2
k+2 ca2
e
.. ..
. .
2k+1 cl,2k+1
ea2
l = 0, ..., k.
Докажем, что
c0,k+1
b2 c0,2k+1
β k+1 + ... + eb2 β 2k+1
c0,k+1 αk + ... + e
c0,2k+1 α2k −
e
ea2 a2 = 0. (17)
α
189
Используя (13) получаем
k k
X 0,p p X 0,p p
c0,k+1
ea2 αk+1 + ... + e
c0,2k+1
a2 α2k+1 = αk+1 Sa2 (A2 − a2 )p−0 + ... + α2k+1 Sa2 (A2 − a2 )p−k =
p=0
0 k
p=k
0,0 k+1 0! 0,1 1! 1! 0,k k!
Sa2 α + Sa2 [(A2 − a2 )αk+1 + αk+2 ] + ... + Sa2 [(A2 − a2 )k αk+1 +
(0 − 0)!0! (1 − 0)!0! (1 − 1)!1! (k − 0)!0!
k! 0,0 k+1 0,1 k+2 0,k 2k+1 0,0 k+1
... + (A2 − a2 )0 α2k+1 ] = Sa2 α + Sa2 α (−1 + 1)1 + ... + Sa2 α (−1 + 1)k = Sa2 α .
(k − k)!k!
Аналогично доказывается, что e c0,k+1
b2 c0,2k+1
β k+1 + ... + eb2
0,0 k+1
β 2k+1 = Sb2 β . Таким образом, (17) верно тогда
0,0 k+1 0,0 k+1
и только тогда, когда Sa2 α = Sb2 β , но
k (i) k (i)
X fx (a2 , b2 ) −i k−i −(−i+k+1) k+1
X fy (a2 , b2 ) −i k − i
(−1) (a2 −A2 ) α = (−1) (b2 −B2 )−(−i+k+1) β k+1 ⇔
i=0
i! k i=0
i! k
f (a2 , b2 )(a2 − A2 )−(k+1) αk+1 = f (a2 , b2 )(b2 − B2 )−(k+1) β k+1 .
Аналогично доказывается, что
2k+1 m−1 2k+1
X (i − 1) · ... · (i − m + 1) i−m X 1 X m−1−i,j j!
cl,i
a2 [ α ]− (−1)i α−i−1 [ cb2 β j−l ] = 0
1 · ... · (m − 1) − − −
e e
i=0
(m 1 i)! (j l)!
i=k+1 j=k+1
для любого m = 1, ..., k + 1 и любого l = 0, ..., k. Таким образом, rank(Mk ) = rank(Mk , e
ck ).
c0,1
При k = 0 сразу же получаем c1,1 = bα2 . Выпишем теперь решение (14) при k = 1 в явном виде. Можно
e
видеть, что rank(M1 ) равен 8; в качестве свободной переменной можем взять c3,3 . Выразим оставшиеся
переменные через свободную и получим (15). Также, выпишем решение (14) при k = 2 в явном виде.
Можно видеть, что rank(M2 ) равен 21; в качестве свободных переменных можем взять c4,4 , c4,5 , c5,4 , c5,5 .
Выразим оставшиеся переменные через свободные и получим (16).
Таким образом, с помощью формул (6), (7) и теоремы 1 мы можем построить требуемую функцию g из
(1), удовлетворяющую всем условиям из (2). Пример экстраполяции с [0, 1]×[0, 1] на [−1, 2]×[−1, 2] (границы
экстраполяции выбраны произвольным образом) для функции f (x, y) = 2·exp (x4 )·cos y +exp (y 3 )·sin x+13
при k = 3 изображен на Рис. 1.
Рис. 1: Пример экстраполяции функции f (x, y) = 2 · exp (x4 ) · cos y + exp (y 3 ) · sin x + 13, (x, y) ∈ [0, 1] × [0, 1]
190
3 Применение экстраполяции для сжатия изображений с потерями
В качестве среды для программной реализации экстраполяции и применения алгоритмов сжатия изоб-
ражений на основе аппарата теории всплесков был выбран пакет прикладных программ Matlab с рас-
ширением Wavelet Toolbox указанных во втором разделе версий. Пусть у нас есть цифровое растровое
изображение, которое можно рассматривать как функцию двух переменных X = X(i, j), (i, j) ∈ Z2 ,
определенную в точках некоторого прямоугольника. Мы будем рассматривать только монохромные (с
глубиной цвета 8 бит на пиксель) изображения размера 2l × 2l (l ∈ 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 с помощью формул (6), (7) и теоремы
1 и получим новое изображение XE = XE (i, j), (i, j) ∈ Z2 . В терминах второго раздела данной статьи
a1 = b1 = 0, a2 = b2 = 2l − 1, A1 = B1 = −2l−1 , A2 = B2 = 2l + 2l−1 − 1. На Рис. 2 показан пример
экстраполяции при k = 1, где из цифрового растрового монохромного изображения X размера 512 × 512
пикселей (область, ограниченная красной рамкой) получено новое изображение XE размера 1024 × 1024
пикселей.
Рис. 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 – объём
сжатых данных;
2
255
• PSNR (P eak Signal-to-N oise Ratio) = 10·log10 M SE (измеряется в децибелах), где MSE (M ean Squared
m n
1
|X(i, j) − X comp (i, j)|2 , X = X(i, j) – исходное изображение, X comp = X comp (i, j) –
P P
Error) = mn
i=1 j=1
восстановленное после сжатия изображение, m × n – размер изображений X и X comp .
191
Однако, так как само XE изначально нигде не хранится и содержит незначимые пиксели, нужные лишь
в процессе сжатия и отбрасываемые после восстановления изображения, то показатель CR не подходит для
сравнения алгоритмов сжатия изображений с предварительной обработкой и без неё. Поэтому вместо CR
мы будем использовать показатель F S (F ile Size) – размер файла сжатого изображения с предваритель-
ной обработкой (F SE ) и без предварительной обработки (F S). MSE для алгоритмов с предварительной
2l P
2l
1 comp
(i, j)|2 ,
P
обработкой изображений будем вычислять следующим образом: MSEE := 2l+1 |X(i, j)−XE
i=1 j=1
comp
где XE – изображение, полученное восстановлением после сжатия изображения XE и взятия только той
его части, которая задана на [0, 2l − 1] × [0, 2l − 1] (то есть изображения 2l × 2l пикселей). Для алгоритмов
2l P
2l 2
1
|X(i, j) − X comp (i, j)|2 . Тогда PSNRE = 10 · log10 M255
P
без предварительной обработки MSE := 2l+1 SEE ,
i=1 j=1
2
255
PSNR = 10 · log10 M SE .
Для того, чтобы сравнить перечисленные алгоритмы сжатия в двух вариациях: без предварительной
обработки и с предварительной обработкой, будем подбирать такие значения параметров в самих алгорит-
мах, чтобы при 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.
192
Таблица 2: Результаты численных экспериментов для алгоритма сжатия EZW
Показатели сжатия Показатели сжатия
Всплески F S (байт) P SN R (дБ) Всплески F SE (байт) P SN RE (дБ)
db3 157041 40.8522 db3 121002 42.7873
db4 156024 36.8092 db4 120931 37.4881
db10 156293 40.6419 db10 122333 42.5396
db12 156942 40.7429 db12 122830 42.7604
coif5 153755 38.3765 coif5 118125 42.6963
sym3 157043 40.8522 sym3 121004 42.7873
sym4 156004 40.7542 sym4 120142 42.6876
sym8 154422 40.6459 sym8 118491 42.5209
sym10 154312 40.7670 sym10 118207 41.4837
sym12 154322 40.7153 sym12 118269 42.5855
sym13 154213 40.7887 sym13 118600 42.6128
sym23 154423 39.9929 sym23 119873 42.5573
sym25 154899 39.4138 sym25 119071 41.4616
sym27 155281 40.8059 sym27 120551 42.6504
dmey 154371 35.6957 dmey 115462 42.4252
bior4.4 151223 37.9625 bior4.4 114871 42.5363
bior5.5 154840 33.8310 bior5.5 118983 42.1763
rbio4.4 164147 40.7066 rbio4.4 129227 42.3120
Таблица 3: Результаты численных экспериментов для алгоритма сжатия SPIHT
Показатели сжатия Показатели сжатия
Всплески F S (байт) P SN R (дБ) Всплески F SE (байт) P SN RE (дБ)
db3 105561 40.6507 db3 91923 41.9744
db4 104675 36.7293 db4 91884 37.2345
db10 104408 40.4485 db10 92269 41.7706
db12 104771 40.5452 db12 92673 41.9442
coif5 103330 38.2626 coif5 90068 41.9338
sym3 105563 40.6507 sym3 91924 41.9744
sym4 104680 40.5580 sym4 91076 41.9092
sym8 103448 40.4553 sym8 89888 41.7851
sym10 103348 40.5700 sym10 89747 40.8933
sym12 103394 40.5214 sym12 89760 41.8377
sym13 103260 40.5918 sym13 89646 41.8677
sym23 103474 39.8281 sym23 90492 41.8076
sym25 103468 39.2689 sym25 89887 40.8666
dmey 101152 35.6334 dmey 81015 41.6875
bior4.4 102195 37.8512 bior4.4 88231 41.7458
bior5.5 105452 33.7903 bior5.5 92536 41.3213
rbio4.4 108951 40.5387 rbio4.4 96576 41.6475
Таблица 4: Результаты численных экспериментов для алгоритма сжатия WDR
Показатели сжатия Показатели сжатия
Всплески F S (байт) P SN R (дБ) Всплески F SE (байт) P SN RE (дБ)
coif5 101088 37.5026 coif5 63700 38.8236
dmey 102063 35.2019 dmey 61629 38.6860
bior4.4 98239 37.1329 bior4.4 61027 38.6562
bior5.5 101805 33.4843 bior5.5 60774 38.2070
193
Благодарности
Работа выполнена за счет гранта Российского научного фонда (проект 14-11-00702).
Список литературы
[1] I. Daubechies. Ten lectures on wavelets. Society for Industrial and Applied Mathematics (SIAM),
Philadelphia, PA, 1992.
[2] A. Cohen, I. Daubechies, P. Vial. Wavelets on the interval and fast wavelet transforms. Appl. Comp. Harm.
Anal., 1(1):54–81, 1992.
[3] A. Cohen, I. Daubechies, J.-C. Feauveau. Biorthogonal bases of compactly supported wavelets.
Communications on Pure and Applied Mathematics, 45(5):485–560, 1992.
[4] Y. Meyer. Wavelets on the interval. Revista Matemática Iberoamericana, 7(2):115–133, 1991.
[5] A.A. Privalov. Theory of approximation of functions. Izd-vo Saratovskogo gos-go un-ta, Saratov, 1990 (in
Russian). = А.А. Привалов. Теория приближения функций. Изд-во Саратовского гос-го ун-та, Саратов,
1990.
[6] J.M. Shapiro. Embedded image coding using zerotrees of wavelet coefficients. IEEE Trans. on Signal
Processing, 41:3445–3462, 1993.
[7] A. Said, W.A. Pearlman. A new, fast, and efficient image codec based on set partitioning in hierarchical
trees. IEEE Trans. on Circuits and Systems for Video Technology, 6(3):243–250, 1996.
[8] A. Said, W.A. Pearlman. Image compression using the spatial-orientation tree. IEEE Int. Symp. on Circuits
and Systems, Chicago, IL, 279–282, 1993.
[9] J. Tian, R.O. Wells, Jr. A lossy image codec based on index coding. IEEE Data Compression Conference,
DCC’96, 456, 1996.
[10] J. Tian, R.O. Wells, Jr. Embedded image coding using wavelet-difference reduction. Wavelet Image and
Video Compression, Topiwala, P., Ed., Kluwer Academic, Norwell, MA, 289–301, 1998.
[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.
194
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 \ Ω 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.
195