Предварительная обработка изображений для улучшения качества работы алгоритмов сжатия с потерями, основанных на всплесках

Одномерный случай

Покажем сначала, как произвольную функцию f = f (x), x ∈ [c, d], обладающую односторонними производными вплоть до k-го порядка в точке d, продолжить на (d, e] гладко до порядка k и с условием f (e) = 0 с помощью многочлена H n (x) подходящей степени, то есть так, чтобы выполнялись следующие условия:

H n (x) ∈ C k ((d, e)), H (i) n (d) = f (i) (d), H (i)

n (e) = 0, i = 0, 1, ..., k.

(3) Лемма 1 Явный вид многочлена H n (x), удовлетворяющего условиям (3), можно получить с помощью следующей формулы:

H 2k+1 (x) = (x − e) k+1 k i=0 f (i) (d) (x − d) i i! k−i j=0 (−1) j k + j k (d − e) −(j+k+1) (x − d) j . (4)

Доказательство Известно (как частный случай эрмитовой интерполяции), что задача нахождения многочлена H n (x) из условий (3) разрешима и, при этом имеет единственное решение (см., например [5], гл. 2, §2.7), если n = 2k + 1. H n (x) -многочлен степени n = 2k + 1, удовлетворяющий условиям (3) и называемый интерполяционным многочленом Эрмита, имеет вид:

H 2k+1 (x) = ω 2k+2 (x) (x − d) k+1 k i=0 f (i) (d) (x − d) i i! (x − d) k+1 ω 2k+2 (x) [k−i] [d] ,(5)

где

ω 2k+2 (x) = (x − d) k+1 (x − e) k+1 , а выражение (x−d) k+1 ω 2k+2 (x) [k−i] [d]

-есть совокупность первых членов разложения в ряд Тейлора по степеням (x−d) функции (x−d) k+1 ω 2k+2 (x) в точке d до (k −i)-го члена включительно. Путем несложных выкладок (5) преобразуется к виду (4).

Таким образом, с помощью леммы 1 мы можем построить многочлен

H n = H n (x), x ∈ (d, e], удо- влетворяющий всем условиям из (3). В качестве численного эксперимента экстраполируем с сохранением гладкости порядка k = 3 функцию α(x) = 3x 2 • cos x + 2x • sin x + 5 с отрезка [0, 4] на всю ось в виде функции с носителем [−2, 6] и проведем сначала всплеск-разложение так экстраполированной функции до уровня 5 с помощью всплесков семейства Добеши с 15-ю нулевыми моментами ∞ −∞ x n ψ(x)dx = 0, n = 0, 1, ...,] × [b 1 , b 2 ] и [a 1 , a 2 ] × (b 2 , B 2 ] так, чтобы выполнялись условия (2). На [A 1 , a 1 ) × [b 1 , b 2 ] и [a 1 , a 2 ] × [B 1 , b 1 ) продолжение будет строится аналогично. Используем формулу (4) из леммы 1 и построим g на (a 2 , A 2 ] × [b 1 , b 2 ] с помощью функции H a2 (x, y): H a2 (x, y) = (x − A 2 ) k+1 k i=0 f (i)

x (a 2 , y)

(x − a 2 ) i i! k−i j=0 (−1) j k + j k (a 2 − A 2 ) −(j+k+1) (x − a 2 ) j .(6)

Используя аналог формулы (4

) из леммы 1 и построим g на [a 1 , a 2 ] × (b 2 , B 2 ] с помощью функции H b2 (x, y): H b2 (x, y) = (y − B 2 ) k+1 k i=0 f (i) y (x, b 2 ) (y − b 2 ) i i! k−i j=0 (−1) j k + j k (b 2 − B 2 ) −(j+k+1) (y − b 2 ) j .(7)

Очевидно, что так построенная функция удовлетворяет условиям (2). Покажем теперь, как построить g на (a

2 , A 2 ] × (b 2 , B 2 ] в виде многочлена P (x, y) = 2k+1 i=1 2k+1 j=1 c i,j (x − A 2 ) i (y − B 2 ) j , (x, y) ∈ [a 2 , A 2 ] × [b 2 , B 2 ],

подчиненного условиям:

P (i) y (x, B 2 ) = 0, P (i) y (x, b 2 ) = (H a2 ) (i) y (x, b 2 ), x ∈ (a 2 , A 2 ], P (i) x (A 2 , y) = 0, P (i) x (a 2 , y) = (H b2 ) (i) x (a 2 , y), y ∈ (b 2 , B 2 ], i = 0, ..., k,(8)где c i,j -неизвестные коэффициенты. На [A 1 , a 1 )×[B 1 , b 1 ) и [A 1 , a 1 )×(b 2 , B 2 ], (a 2 , A 2 ]×[B 1 , b 1 ) продолжение строится аналогично. Покажем, как найти все c i,j . Обозначим c 0,i a2 = 2k+1 j=1 c i,j (b 2 − B 2 ) j , c l,i a2 = 2k+1 j=l j! (j−l)! c i,j (b 2 − B 2 ) j−l , i = 1, ..., 2k + 1, c 0,j b2 = 2k+1 i=1 c i,j (a 2 − A 2 ) i , c l,j b2 = 2k+1 i=l i! (i−l)! c i,j (a 2 − A 2 ) i−l , j = 1, ..., 2k + 1, l = 1, ..., k,(9)

тогда

P (l) y (x, b 2 ) = 2k+1 i=1 c l,i a2 (x − A 2 ) i , P(l)

x (a 2 , y)

= 2k+1 j=1 c l,j b2 (y − B 2 ) j , l = 0, ..., k.(10)

Из ( 6)-( 10) можно видеть, что c l,i a2 = 0, i = 1, ..., k, c l,j b2 = 0, j = 1, ..., k, l = 0, ..., k.

Перепишем соответствующие условия из (8) с учетом ( 11)

k i=0 c l,i+k+1 a2 (x − A 2 ) i = k i=0 (f (i) x ) (l) y (a2,b2) i! k−i j=0 (−1) j k+j k (a 2 − A 2 ) −(j+k+1) (x − a 2 ) i+j , k i=0 c l,i+k+1 b2 (y − B 2 ) i = k i=0 (f (i) y ) (l) x (a2,b2) i! k−i j=0 (−1) j k+j k (b 2 − B 2 ) −(j+k+1) (y − b 2 ) i+j , l = 0, ..., k. В правых частях двух последних уравнений сделаем замену j = p − i, поменяем порядок суммирования и получим k i=0 c l,i+k+1 a2 (x − A 2 ) i = k p=0 k i=0 (f (i) x ) (l) y (a2,b2) i! (−1) p−i k+p−i k (a 2 − A 2 ) −(p−i+k+1) (x − a 2 ) p , k i=0 c l,i+k+1 b2 (y − B 2 ) i = k p=0 k i=0 (f (i) y ) (l) x (a2,b2) i! (−1) p−i k+p−i k (b 2 − B 2 ) −(p−i+k+1) (y − b 2 ) p , l = 0, ..., k. (12) Обозначим S l,p a2 = k i=0 (f (i) x ) (l) y (a2,b2) i! (−1) p−i k+p−i k (a 2 − A 2 ) −(p−i+k+1) , S l,p b2 = k i=0 (f (i) y ) (l) x (a2,b2) i! (−1) p−i k+p−i k (b 2 − B 2 ) −(p−i+k+1) , p = 0, ..., k, l = 0, ..., k. Тогда дифференцируя соответствующие выражения из (12) в точках A 2 и B 2 соответственно получим c l,i+k+1 a2 = k p=i S l,p a2 (A 2 − a 2 ) p−i p i , c l,i+k+1 b2 = k p=i S l,p b2 (B 2 − b 2 ) p−i p i , i = 0, ..., k, l = 0, ..., k. (13) Обозначим α = a 2 − A 2 , β = b 2 − B 2 . Пусть c i

a2 , c i b2 , i = 0, ..., k -вектор-столбцы, координаты которых вычислены по формулам (11) и (13), тогда мы можем переписать (9) в виде системы линейных уравнений

M k • c = c k (14)

относительно столбца неизвестных c, где

M k =           A(α) A 2 (α) A 3 (α) • • • A 2k+1 (α) . . . . . . . . . . . . . . . A (k) (α) (A 2 (α)) (k) (A 3 (α)) (k) • • • (A 2k+1 (α)) (k) B 1 (β) B 2 (β) B 3 (β) • • • B 2k+1 (β) . . . . . . . . . . . . . . . B (k) 1 (β) B (k) 2 (β) B (k) 3 (β) • • • B (k) 2k+1 (β)           -матрица размерности 2(k + 1)(2k + 1) × (2k + 1) 2 ; A(α) = diag{α, α ..., α} -диагональная матрица раз- мерности (2k + 1) × (2k + 1); B i (β) -матрица размерности (2k + 1) × (2k + 1): (B i (β))(j, ) = (β β 2 • • • β 2k+1 ), j = i, (0 0 • • • 0), j = i, c = (c 1,1 , ..., c 1,2k+1 , ..., c 2k+1,1 , ..., c 2k+1,2k+1 ) T , c k = (( c 0 b2 ) T , ..., ( c k b2 ) T , ( c 0 a2 ) T , ..., ( c k a2 ) T ) T .

Решив систему (14), мы определим все неизвестные коэффициенты c i,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 1 α α 2 • • • α j−1 • • • α k−1 0 1 2α • • • j−1 1 α j−2 • • • k−1 1 α k−2 0 0 1 • • • (j−1)(j−2) 1•2 α j−3 • • • (k−1)(k−2) 1•2 α k−3 . . . . . . . . . . . . . . . . . . . . . 0 0 0 • • • (j−1)•...•(j−m+1) 1•...•(m−1) α j−m • • • (k−1)•...•(k−m+1) 1•...•(m−1) α k−m         2. Вектор-столбец (c 1 , ..., c i , ..., c m ) T преобразованиями, осуществляющими пункт 1, приводится к виду ( c1 α , ..., i−1 j=0 (−1) j 1 (i−1−j)! α −j−1 c i−j , ..., m−1 j=0 (−1) j 1 (m−1−j)! α −j−1 c m−j ) T .

Доказательство Докажем утверждение 1 с помощью индукции по n -количество первых строк исходной матрицы (n < m).

База индукции при n = 1 очевидна. Предположим, что исходную матрицу можно эквивалентными преобразованиями привести к искомому виду

      1 ••• n n+1 ••• k . . . . . . . . . . . . . . . . . . n 0 • • • 1 n•...•2 1•...•(n−1) α • • • (k−1)•...•(k−n+1) 1•...•(n−1) α k−n n+1 0 • • • n • ... • 1 (n + 1) • ... • 2α • • • k • ... • (k − n + 1)α k−n . . . . . . . . . . . . . . . . . .       . В полученной матрице вычтем из (n + 1)-й строки n-ю строку, умноженную на n!,       1 ••• n n+1 ••• k . . . . . . . . . . . . . . . . . . n 0 • • • 1 n•...•2 1•...•(n−1) α • • • (k−1)•...•(k−n+1) 1•...•(n−1) α k−n n+1 0 • • • 0 n • ... • 1α • • • (k − 1) • ... • (k − n)α k−n . . . . . . . . . . . . . . . . . .       ∼        1 ••• n n+1 ••• k . . . . . . . . . . . . . . . . . . n 0 • • • 1 n•...•2 1•...•(n−1) α • • • (k−1)•...•(k−n+1) 1•...•(n−1) α k−n n+1 0 • • • 0 1 • • • (k−1)•...•(k−n) 1•...•n α k−n−1 . . . . . . . . . . . . . . . . .        .

Доказательство утверждения 2 аналогично доказательству утверждения 1.

Теорема Система (14) совместна и ее решение имеет вид:

при k = 0: c 1,1 = c 0,1 b 2 α ; при k = 1:          c 3,2 = c 1,3 a 2 1•β − c 0,3 a 2 β 2 − 2βc 3,3 , c 3,1 = 2 c 0,3 a 2 β − c 1,3 a2 + β 2 c 3,3 , c 2,1 = − 4α c 0,3 a 2 β + 2α c 1,3 a2 − 2αβ 2 c 3,3 , c 2,2 = c 1,2 b 2 1•α − c 0,2 b 2 α 2 − 2α c 1,3 a 2 1•β + 2α c 0,3 a 2 β 2 + 4αβc 3,3 , c 2,3 = c 1,3 b 2 1•α − c 0,3 b 2 α 2 − 2αc 3,3 , c 1,1 = 2α 2 c 0,3 a 2 β − α 2 c 1,3 a2 + α 2 β 2 c 3,3 , c 1,2 = 2 c 0,2 b 2 α − c 1,2 b2 + α 2 c 1,3 a 2 1•β − α 2 c 0,3 a 2 β 2 − 2α 2 βc 3,3 , c 1,3 = 2 c 0,3 b 2 α − c 1,3 b2 + α 2 c 3,3 , (15) где c 3,3 -свободная переменная; при k = 2:                                            c 5,3 = Доказательство С помощью леммы 2 эквивалентными преобразованиями приведем матрицу M k из (14) к виду              1 ••• j ••• 2k+1 1 E • • • A j−1 (α) • • • A 2k (α) 2 0 • • • j−1 1 A j−2 (α) • • • 2k 1 A 2k−1 (α) . . . . . . . . . . . . . . . . . . k+1 0 • • • (j−1)•...•(j−k) 1•...•k A j−k−1 (α) • • • 2k•...•(k+1) 1•...•k A k (α) k+2 B 1 (β) • • • B j (β) • • • B 2k+1 (β) . . . . . . . . . . . . . . . . . . 2k+1 B (k) 1 (β) • • • B (k) j (β) • • • B (k) 2k+1 (β)              ∼              1 ••• k+1 k+2 ••• 2k+1 1 E • • • A k (α) A k+1 (α) • • • A 2k (α) 2 0 • • • k 1 A k−1 (α) k+1 1 A k (α) • • • 2k 1 A 2k−1 (α) . . . . . . . . . . . . . . . . . . . . . k+1 0 • • • E (k+1)•...•2 1•...•k A(α) • • • 2k•...•(k+1) 1•...•k A k (α) k+2 0 • • • 0 B k+2 (β) • • • B 2k+1 (β) . . . . . . . . . . . . . . . . . . . . . 2k+1 0 • • • 0 B (k) k+2 (β) • • • B (k) 2k+1 (β)              . Несложно видеть,( c 0 b2 , ..., c i−1 b2 , ..., c k b2 ) T ∼ ( c 0 b2 α , ..., i−1 j=0 (−1) j 1 (i − 1 − j)! α −j−1 c i−1−j b2 , ..., k j=0 (−1) j 1 (k − j)! α −j−1 c k−j b2 ) T ,

c l a2 = (0, 0, ..., 0, ..., 0, c l,k+1 a2 , c l,k+2 a2 , ..., c l,2k+1 a2

) T ∼                            1 c l,k+1 a2 α k + ... + c l,2k+1 a2 α 2k − 2k+1 j=k+1 c 0,j b 2 j! (j−l)! β j−l α 2 c l,k+1 a2 [ k 1 α k−1 ] + ... + c l,2k+1 a2 [ 2k 1 α 2k−1 ] − 2k+1 j=k+1 c 1,j b 2 j! (j−l)! β j−l 1•α + 2k+1 j=k+1 c 0,j b 2 j! (j−l)! β j−l α 2 . . . . . . m 2k+1 i=k+1 c l,i a2 [ (i−1)•...•(i−m+1) 1•...•(m−1) α i−m ] − m−1 i=0 (−1) i 1 (m−1−i)! α −i−1 [ 2k+1 j=k+1 c m−1−i,j b2 j! (j−l)! β j−l ] . . . . . . k+1 2k+1 i=k+1 c l,i a2 [ (i−1)•...•(i−k) 1•...•k α i−k−1 ] − k i=0 (−1) i 1 (k−i)! α −i−1 [ 2k+1 j=k+1 c k−i,j b2 j! (j−l)! β j−l ] k+2 c l,k+2 a2 . . . . . . 2k+1 c l,2k+1 a2                            , l = 0, ..., k. Докажем, что c 0,k+1 a2 α k + ... + c 0,2k+1 a2 α 2k − c 0,k+1 b2 β k+1 + ... + c 0,2k+1b2β 2k+1 α = 0. (17)

Используя (13) получаем

c 0,k+1 a2 α k+1 + ... + c 0,2k+1 a2 α 2k+1 = α k+1 k p=0 S 0,p a2 (A 2 − a 2 ) p−0 p 0 + ... + α 2k+1 k p=k S 0,p a2 (A 2 − a 2 ) p−k p k = S 0,0 a2 α k+1 0! (0 − 0)!0! + S 0,1 a2 [(A 2 − a 2 )α k+1 1! (1 − 0)!0! + α k+2 1! (1 − 1)!1! ] + ... + S 0,k a2 [(A 2 − a 2 ) k α k+1 k! (k − 0)!0! + ... + (A 2 − a 2 ) 0 α 2k+1 k! (k − k)!k! ] = S 0,0 a2 α k+1 + S 0,1 a2 α k+2 (−1 + 1) 1 + ... + S 0,k a2 α 2k+1 (−1 + 1) k = S 0,0 a2 α k+1 .

Аналогично доказывается, что c 0,k+1 b2 β k+1 + ... + c 0,2k+1 b2 β 2k+1 = S 0,0 b2 β k+1 . Таким образом, (17) верно тогда и только тогда, когда S 0,0 a2 α k+1 = S 0,0 b2 β k+1 , но

k i=0 f (i) x (a 2 , b 2 ) i! (−1) −i k − i k (a 2 −A 2 ) −(−i+k+1) α k+1 = k i=0 f (i) y (a 2 , b 2 ) i! (−1) −i k − i k (b 2 −B 2 ) −(−i+k+1) β k+1 ⇔ f (a 2 , b 2 )(a 2 − A 2 ) −(k+1) α k+1 = f (a 2 , b 2 )(b 2 − B 2 ) −(k+1) β k+1 .

Аналогично доказывается, что 3 Применение экстраполяции для сжатия изображений с потерями В качестве среды для программной реализации экстраполяции и применения алгоритмов сжатия изображений на основе аппарата теории всплесков был выбран пакет прикладных программ Matlab с расширением Wavelet Toolbox указанных во втором разделе версий. Пусть у нас есть цифровое растровое изображение, которое можно рассматривать как функцию двух переменных X = X(i, j), (i, j) ∈ Z 2 , определенную в точках некоторого прямоугольника. Мы будем рассматривать только монохромные (с глубиной цвета 8 бит на пиксель) изображения размера 6), ( 7) и теоремы 1 и получим новое изображение X E = X E (i, j), (i, j) ∈ Z 2 . В терминах второго раздела данной статьи После предварительной обработки исходного изображения будем применять к X и X E четыре основанных на всплесках алгоритма сжатия:

2k+1 i=k+1 c l,i a2 [ (i − 1) • ... • (i − m + 1) 1 • ... • (m − 1) α i−m ] − m−1 i=0 (−1) i 1 (m − 1 − i)! α −i−1 [ 2k+1 j=k+1 c m−1−i,j b2 j! (j − l)! β j−l ] = 0 для любого m = 1, ..., k + 1 и любого l = 0, ..., k. Таким образом, rank(M k ) = rank(M k , c k ). При k = 0 сразу же получаем c 1,1 = c 0,2 l × 2 l (l ∈ N) пикселей, определенные в це- лых точках [0, 2 l − 1] × [0, 2 l − 1]. Экстраполируем исходное изображение X размера 2 l × 2 l пикселей, определенное на [0, 2 l − 1] × [0, 2 l − 1], до изображения размера 2 l+1 × 2 l+1 пикселей, определенное на [−2 l−1 , 2 l + 2 l−1 − 1] × [−2 l−1 , 2 l + 2 l−1 − 1], с требуемой гладкостью k с помощью формул (a 1 = b 1 = 0, a 2 = b 2 = 2 l − 1, A 1 = B 1 = −2 l−1 , A 2 = B 2 = 2 l + 2 l−1 − 1. На Рис.

• 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) = Sc S0 (измеряется в процентах), где S 0 -объём исходных данных, а S c -объём сжатых данных;

• PSNR (P eak Signal-to-N oise Ratio) = 10•log 10

255 2 M SE (измеряется в децибелах), где MSE (M ean Squared Error) = 1 mn m i=1 n j=1 |X(i, j) − X comp (i, j)| 2 , X = X(i, j) -исходное изображение, X comp = X comp (i, j) - восстановленное после сжатия изображение, m × n -размер изображений X и X comp .

Однако, так как само X E изначально нигде не хранится и содержит незначимые пиксели, нужные лишь в процессе сжатия и отбрасываемые после восстановления изображения, то показатель CR не подходит для сравнения алгоритмов сжатия изображений с предварительной обработкой и без неё. Поэтому вместо CR мы будем использовать показатель F S (F ile Size) -размер файла сжатого изображения с предварительной обработкой (F S E ) и без предварительной обработки (F S). MSE для алгоритмов с предварительной обработкой изображений будем вычислять следующим образом: MSE In this paper we developed extrapolation method of 2k-smooth functions in a neighborhood of the rectangle Ω to the R 2 \ Ω 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.

E := 1 2 l+1 2 l i=1 2 l j=1 |X(i, j)−X comp E (i, j)| 2 , где X comp E -изображение, полученное восстановлением после сжатия изображения X E и взятия только той его части, которая задана на [0, 2 l − 1] × [0, 2 l − 1] (то есть изображения 2 l × 2 l пикселей). Для алгоритмов без предварительной обработки MSE := 1 2 l+1 2 l i=1 2 l j=1 |X(i, j) − X comp (i, j)| 2 . Тогда PSNR E = 10 •