Реализация алгоритма сжатия изображений на основе чебышевских преобразований (алгоритм GDCT) Текст научной статьи по специальности «Компьютерные и информационные науки»

РЕАЛИЗАЦИЯ АЛГОРИТМА СЖАТИЯ ИЗОБРАЖЕНИИ НА ОСНОВЕ ЧЕБЫШЕВСКИХ ПРЕОБРАЗОВАНИИ (АЛГОРИТМ GDCT)

Ю.С. Радченко (rad@sendmail.ru)

Воронежский госуниверситет

Предложен метод сжатия изображения, основанный на разложении по полиномам Чебышева и применении квадратурных формул Гаусса - Чебышева (алгоритм GDCT). Показано, что алгоритм GDCT является обобщением дискретного косинусного преобразования, используемого в алгоритмах JPEG и MPEG. Алгоритм имеет четыре ступени сжатия: прореживание исходного изображения, усечение спектра, обнаружение изменяющихся фрагментов, определение вектора движения фрагментов. На примерах обработки искусственных текстур и натурных изображений продемонстрированы возможности алгоритма GDCT .

Системы связи нового поколения ориентированы на передачу информационных потоков различного вида, то есть являются мультимедийными телекоммуникационными системами [1,2]. В связи с этим передача телевизионных и компьютерных изображений, видеоконференции и видеосотовая связь предполагают применение процедур сжатия изображений. В последнее время выполнен ряд работ, в которых предлагается применять разложения по базису классических ортогональных полиномов Эрмита, Якоби (Лежандра, Чебышева) [3,4]. Эти преобразования обладают весьма компактными спектрами, которые чувствительны к сдвигу анализируемого фрагмента. В данной работе приводятся результаты исследования одного из вариантов полиномиальных преобразований - разложение и синтез по базису полиномов Чебышева I и II рода. Вычисление спектральных коэффициентов производится с помощью квадратурной формулы Гаусса - Чебышева, обладающей наивысшей степенью алгебраической точности при заданном числе отсчетов. Показано, что при этом получаются преобразования, которые являются обобщением дискретного косинусного и синусного преобразований. 1. Разложение сигналов по базису ортогональных полиномов.

Пусть в подобласти (x,y}eQ наблюдается поле s(x,y,T), представляющее собой фрагмент u(x,y)IQ(x,y) пространственного сигнала. Здесь IQ (x, y) - индикаторная функция подобласти, T=(tx, Ty) параметры сдвига фрагмента в данном кадре. Базисные функции 9mk(x,y), используемые для дискретного представления сигнала s(x,y,T), определяются аналитическими свойствами самого сигнала и геометрической формой подобласти Q. Однако, реализация процедуры сжатия, особенно для полиномиальных базисов, существенно упрощается, если имеет место факторизация функций 9mk(x,y)= 9m(x) 9k(y).Здесь 9m(x), 9k(y) - одномерные функции, основанные на ортогональных полиномах. Тогда для полезного сигнала s(x,y,T) имеет место пара преобразований

s(x,y,т) = XXCmk( T)9m(x)9k(yX Cmk(T) = Ц s(x,y, Т)ф m(x^k(y)dxdy. (1)

m k Q

Для разложения сигналов и изображений можно применять следующие классические ортогональные многочлены: а) Эрмита; б) Лежандра; в) полиномы Чебышева I и II рода. Если обозначить ax,ay - характерные размеры подобласти Q, z1=x/ax, z2=y/ay, то (1) можно переписать с использованием ортогональных полиномов в виде

s(x,y) = X CmkPm(x/ax)Pk(y/ay)

m,k

Сшк = -1 Ц ^у^М^ )Pш(z1 )P(z2 )Pk(z2 )dz1dz2 =

О* (2)

= ^ш^-1 JP(z1 )Pш(z1 )dz11 ^у^М^Рк^^.

Здесь dш - норма ортогонального с весом р^) полинома рш^). Поскольку для разложимых базисных функций фшк(х,у)= фш(х) фк(у) вычисление спектральных коэффициентов производится последовательным интегрированием по координатам (х,у), то для упрощения анализа рассмотрим сначала одномерные преобразования, а затем обобщим их на двумерный случай.

2. Алгоритм преобразований Чебышева.

Пусть имеется процесс Дг), квадратично интегрируемый с весом р^). Тогда можно записать квадратурную формулу гауссовского типа наивысшей алгебраической степени

N

точности порядка 2№1 как 1 f(z)р(z)dz = ^А,ПД^П) . Здесь zn - нули полиномов рN

п=1

ортогональных с весом р^); Ап - числа Кристоффеля. Узлы и веса {Ап} однозначно определяются видом ортогонального с весом р^) полинома рш^). Для полиномов Чебышева получается формула Меллера (Гаусса-Чебышева)

1 1 s(z)Tm(z)J п ^К ^ , л

Г / 2 dz = -— Xs(Zn)Tш(Zn). (3)

dш -1 л/1 - Z2 dшNn=0

Здесь zn=cos(п(2n+1)/2N) -нули полинома Чебышева I рода Т^); Ап=л/К - то есть все весовые коэффициенты одинаковы, норма полинома Чебышева dm=п/2, если и dm=п,

если ш=0. Учитывая, что Tm(zn)=cos(m•arccos(zn))=cos(пm(n+0.5)/N), приходим к выражению для прямого и обратного преобразований

Г^-1 п + 05 Г1N-1

Сш = л — X s(zn)c0s(пш С0 =лЬ X^^

V N п=0 N Шп=0 (4)

Г2 м Г2 м

SM(z) = ёшд N XСшТш(z)= ёшд N XCшc0s(ш• arccos(z))

V ^ ш=0 V ^ ш = 0

0.5 , ш = 0

Здесь = < .В формуле (4) использована симметричная нормировка

1 , ш > 0

базисных функций, удобная для практической реализации преобразований в матричном виде. Выражение (4) весьма похоже на так называемое четное дискретное косинусное преобразование (ДКП) [2]. Но имеет три важных отличия: 1)Точки отсчета zn = C0s(п(п + 0.5) /^ сигнала s(z) берутся неравномерно. 2)Синтез сигнала SM(z)

выполняется в произвольной точке ze[-1,1], а не в дискретном наборе точек отсчета, как в ДКП. 3)Точность формулы (4) при преобразовании достаточно гладких функций существенно выше, чем у ДКП. Поэтому число отсчетов можно взять значительно меньше.

Если при восстановлении взять неравномерную сетку отсчетов по закону Zj=cos(п(i+0.5)/L), ]=0,.. Х-1, то обратное преобразование принимает вид

ад^ = ёшд 2 X Сш ^^^^г05). (5)

^^ = 0 L

Преобразование (5) выглядит аналогично ДКП, однако спектральные коэффициенты

отличаются от спектральных коэффициентов обычного ДКП. В том случае, если нас

интересует равномерная сетка отсчетов Zj=2j/(L-1) -1 +8, где 8 - некоторый сдвиг, ]=0,.. Х-1, то

1~2 м

= Втл — XСт С0^т • arccOS(2j/(L - 1) - 1 + 8). (6)

^т-0

Алгоритм преобразования (4), (6) назовем обобщенным дискретным косинусным преобразованием (ОБСТ) . В формулах (5), (6) L число точек в блоке восстановленного изображения может быть произвольным. В таком случае возникает эффект масштабирования восстановленного изображения по сравнению с размерами исходного блока с изображением. Если использовать разложение по полиномам Чебышева II рода, то получаем

2 1

(7)

Ст - N1 -z2s(z)Um(z)dz = п-1

2 V / ч ■ / (к + 1)(т+1Х . , к +1,

N +1 к=

к=0 м

N +1

N +1

= X CmUm(z).

т=0

Здесь учтено, что zk=cos(п(k+1)/(N+1)) и

ит^к)

= sin((M +1) arccos(zk)) = sin(п(k + 1)(т +1) /(N +1))

1 - zí

sin(п(k + 1)/(N +1))

(8)

Алгоритм преобразований (7) сходен с дискретным синусным преобразованием

Для использования квадратурных формул типа Гаусса (3) при переходе от интегральных преобразований к дискретным требуется определить число узлов N обеспечивающее вычисление необходимого числа спектральных коэффициентов Ст с заданной точностью. Как уже указывалось выше, формулы типа Гаусса являются точными для всех полиномов степени 2№1, а порядок ошибки формулы (3) Однако

эта оценка чисто качественная. Кроме того, ряд моделей полезных сигналов имеет ряд особых точек, в которых функция не дифференцируема. Поэтому необходимы расчеты для типовых моделей сигналов, с помощью которых можно оценить необходимое число отсчетов N. Вопросы точности аппроксимации сигналов и выбора числа узлов квадратурной формулы исследовались на примере следующих тестовых моделей: непрерывных дифференцируемых и недифференцируемых на множестве точек, финитных, разрывных (импульсных).

Например, если в качестве модельного сигнала, взять s(z) — V1 — Z2 , то точные значения спектральных коэффициентов имеют значения

СТ --С0 =

м

^^ - 2

2

п—1л/1 - Z2 п

2 1 = 1

СТт

К л

-1

1 - z2

7Г

т+1

Тт +1 <Ж>"

т -1

Тт-1Ф

(9)

-1

1

1

1

- 4

--- при четном ш > 0,

п(ш2 -1)

0 при нечетном ш.

При вычислении по формуле (4) соответствующие коэффициенты имеют значения:

1 1 . , к + 0.5,

С0 =— X sin(п-

к=0

N

),

2 ^ • , к + 0.5, /пш(к + 0.5),

Сш = ~ X ^^Т")c0s(-^->

N г=0 N N

(10)

Т

В таблице ¡приведены значения Сш и Сш при различном числе квадратурных точек N (N=4,8,12,16).

Как видно из таблицы достаточно точно вычисляются коэффициенты с номером ш<3М4.

Таблица 1

ш 0 2 4 6 8 10 12

С 0.673 -0.424 -0.085 -0.036 -0.02 -0.013 -0.0089

Cm(N=4) 0.653 -0.383 -0 0.383 -1.307 0.383

Cm(N=8) 0.641 -0.416 -0.075 -0.023 0 0.023 0.075

^(N=12) 0.638 -0.421 -0.081 -0.032 -0.015 -0.0065 0

Cm(N=16) 0.638 -0.422 -0.083 -0.034 -0.018 -0.01 -0.0055

Были исследованы ошибки восстановления сигналов с использованием метрик

dp(Sм,S) = { |Sм(x) - S(x)|pdx}1 /p

в пространстве Lp измеримых и интегрируемых функций для значений p=1,2, ^ . Расчеты для различных классов сигналов показали достаточность 6-8 точек для использования чебышевских преобразований. 3. Блочные преобразования (ОБСТ).

Перейдем теперь к двумерному варианту преобразований сигнала s(x,y). Пусть в пределах блока из ШхШ точек, берутся отсчеты сигнала NxN штук по закону

хп = Я0ШБ(0.5 • (N1 -1) • (1 + ^(п(П + 0.5) /

ук = Я0иКО(0.5 • (N1 -1) • (1 + ^((к + 0.5)/N))) (11)

Отсчеты сигнала образуют матрицу 8 = ^Ц = ^(х^ ук)||. Затем эта матрица преобразуется в

матрицу спектральных коэффициентов С размером МхМ. При обратном преобразовании может использоваться прямоугольная матрица размером LxM. То есть восстановленный блок имеет размеры LxL. Прямое и обратное преобразование Чебышева (ОБСТ) в матричном виде определяются операциями

С = С

ш,1

■ ф 8фТ

8 = ¥ ТС¥

(12)

Матрицы прямого и обратного преобразования ОБСТ имеют вид

N2

ф =||ф ш(п)|

11КМ

cos(пш

0.5

(П + 0.5)

N

¥

cos(m • arccos(Zj)

LM

cos(ш • arccos(

L -1

1 + 5))

(13)

В частном случае можно полагать 5=0. Дискретизация аргументов {x,y} по правилу (11) приводит к появлению погрешности квадратурной формулы и неточному вычислению спектральных коэффициентов. Для уточнения расчетов можно произвести интерполяцию значения поля s(x,y) в точках расположения нулей полиномов Чебышева по ближайшим отсчетам.. Например, для этого можно использовать четырехточечную интерполяционную формулу Бесселя

s(xn,yk)=0.25(s00+s01+s10+sn)+0.5(u-0.5)(s10-s00+s11-s01)+

+0.5(v-0.5)(s01-s00+s11-s10)+(u-0.5)(v-0.5)(sn-s10-s01+s00), (14)

где u=xn-x00, v=yk-y00, Sij (i,j =0.. 1 ) - значения сигнала в ближайших четырех пикселях, окружающих точку xn,yk. Применение интерполяции улучшает передачу мелкомасштабной текстуры изображений при реализации процедуры сжатия - восстановления сигналов.

4. Результаты исследования алгоритма GDCT.

Итак, в блоке изображений размером N1xN1 элементов выбирается матрица NxN отсчетов (N<N1). Так происходит сжатие на ступени I процедуры. Затем с помощью преобразования GDCT вычисляется матрица MxM из наиболее значимых спектральных коэффициентов, которые в дальнейшем прореживаются или при помощи матрицы квантования, или путем регулируемой низкочастотной фильтрации. Так происходит сжатие изображения на II ступени процедуры. При восстановлении изображения формируется блок размером LxL элементов. При восстановлении размер блока LxL может быть меньше, равен или больше размеров исходного блока N1xN1.

Для экспериментов был выбран набор искусственных изображений и натурных фотографий. Цветные изображения были представлены в 24 битовом формате BMP. С этими изображениями осуществлялись следующие операции. 1) Перекодировка 8-битовых R, G, B -компонент цвета в Y,U,V по правилу: Y=int(0.299R+ 0.587G +0.114), U=int(R-Y=0.701R-0.587G-0.114B), V= B-Y=int(-0.299R-0.587G+0.886B). 2) Субдискретизация матриц Y,U,V по одному из законов: 4:4:4, 4:2:2, 4:2:0. 3) Спектральное преобразование матриц Y,U,V на основе GDCT. В зависимости от типа изображения на первой ступени сжатия проводилось прореживание в пропорции N1/N=16/6, 16/8, 8/6. Затем проводилось преобразование спектральной матрицы C для реализации второй ступени сжатия изображения. 4) Для упрощения сопоставления исходного и восстановленного изображений проводилось восстановление блока размером N1xN1, то есть было выбрано N1=L. Затем производилась визуальная оценка качества изображения. На фигурах 1 -4 показаны результаты преобразований объекта типа «портрет» с помощью алгоритмов GDCT и JPEG. Степень сжатия за счет усечения спектров в алгоритмах JPEG и GDCT взята 5.6 .Но в алгоритме GDCT проведено дополнительно сжатие первой ступени в 1.8 (фигура 2) и 4 раза (фигура 3).

Поскольку алгоритм сжатия и восстановления на основе чебышевских преобразований обладают свойством масштабирования восстановленных изображений по сравнению с оригиналами, было проведено исследование такой возможности GDCT. В зависимости от степени гладкости исходного изображения коэффициент пространственного увеличения (L/N1) может доходить до 20. Исследование преобразований типовых натурных изображений с резкими границами и мелкомасштабной текстурой показало, что при коэффициенте увеличения (L/N1)>4 начинает проявляться блочная структура изображений. Выводы

Исследования показали, что изображения с плавными градациями яркости и цвета могут быть подвергнуты сильному сжатию уже на I ступени процедуры. Насыщенные мелкими деталями текстуры и натурные изображения испытывают искажения в GDCT, и рекомендуемая степень сжатия I ступени порядка 30-50%. Скорость убывания спектров в GDCT такого же порядка, как и в классическом дискретном косинусном преобразовании. Поэтому сжатие на II ступени такого же порядка, как и в алгоритме JPEG. Так как положение узлов полиномов Чебышева образуют дробные числа, то были исследованы

две реализации GDCT - с округлением положения узлов до ближайшего пикселя, и с межпиксельной линейной интерполяцией. Эксперименты показали, что интерполяция улучшает передачу тонких деталей и наклонных краев изображений. Выполнено сравнение в статическом режиме алгоритмов сжатия GDCT и классического JPEG для различных типов изображений. Показано, что в зависимости от типа изображений алгоритм GDCT в статическом режиме позволяет увеличить степень сжатия до 3 раз по сравнению с алгоритмом JPEG.

Автор благодарен М.Ю. Радченко и к.т.н. А.Ю. Савинкову за помощь при разработке программного кодека алгоритма GDCT на языке С++.

Литература

1. А.В. Дворкович, В.П. Дворкович, Б.Н. Мохин и др. Единые принципы сжатия цветных динамических изображений различного разрешения. Цифровая обработка сигналов, №1, 1999, с. 27-35.

2. Цифровая обработка телевизионных и компьютерных изображений./ Под ред. Ю.Б. Зубарева и В.П. Дворковича. М.: МЦНТИ, 1997, - 217 с.

3. Ю.С. Радченко, М.Ю. Радченко Оптимальные быстрые алгоритмы представления изображений в базисе ортогональных полиномов. Труды I международной конференции «Цифровая обработка сигналов и ее применения» DSPA'98, 1998, Москва, т. III, с. 163-166.

4. Ю.С. Радченко, А.Ю. Кожин, М.Ю. Радченко Обнаружение и оценка параметра сдвига сжатых с помощью ортогональных полиномов сигналов. Радиотехника, №6, 1999, с. 17-19.

5. Ю.С. Радченко, А.Ю. Савинков Исследование алгоритмов сжатия изображений на основе полиномиальных преобразований (алгоритм GDKT). Труды III конференции «Цифровая обработка сигналов и ее применения» DSPA'2000, 2000, Москва, т. II, с. 89-91

Объект-портет

Фигура 1. Исходный Фигура 2. GDCT 8^6

Фигура 3. GDCT 12^6

Фигура 4. JPEG -5.6