close

Вход

Забыли?

вход по аккаунту

?

Veselov

код для вставкиСкачать
ВВЕДЕНИЕ
В ЦИФРОВУЮ ОБРАБОТКУ ИЗОБРАЖЕНИЙ.
МЕТОДЫ ФИЛЬТРАЦИИ И СЖАТИЯ ИЗОБРАЖЕНИЙ
ISBN: 9785808810792
9 785808 810792
Учебное издание
Гильмутдинов Марат Равилевич,
Веселов Антон Игоревич,
Ястребов Виктор Анатольевич,
Тюрликов Андрей Михайлович
МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ
РОССИЙСКОЙ ФЕДЕРАЦИИ
Федеральное государственное автономное образовательное
учреждение высшего образования
Санкт-Петербургский государственный университет аэрокосмического
приборостроения
ВВЕДЕНИЕ В ЦИФРОВУЮ ОБРАБОТКУ ИЗОБРАЖЕНИЙ.
МЕТОДЫ ФИЛЬТРАЦИИ И СЖАТИЯ ИЗОБРАЖЕНИЙ
Учебное пособие
ВВЕДЕНИЕ В ЦИФРОВУЮ ОБРАБОТКУ
ИЗОБРАЖЕНИЙ.
Методы фильтрации и сжатия изображений
Редактор В.П.Зуева
Учебное пособие
Сдано в набор 06.05.15. Подписано к печати 24.12.15. Формат бумаги
60 × 84 1/16. Бумага офсетная. Усл. печ. л. 4,4. Уч.-изд. л. 4,8. Тираж 100 экз.
Заказ № 581
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Санкт-Петербург
2015
УДК 004.627
ББК 30
В24
Рецензенты:
доктор технических наук, профессор кафедры вычислительной техники ИТМО
А. Ю. Тропченко;
1.3.
Обработка данных в базовом режиме . . .
1.4.
Дискретное косинусное преобразование . .
1.5.
Квантование спектральных коэффициентов
1.6.
Кодирование без потерь . . . . . . . . . . .
2.
Практические задания . . . . . . . . . . . . . .
3.
Задания повышенной сложности . . . . . . . .
4.
Вопросы для самопроверки . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
55
57
59
61
65
68
70
кандидат технических наук, доцент кафедры проблемно-ориентированных вычислительных
Библиографический список
комплексов ГУАП Е. А. Бакин
72
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
Гильмутдинов М. Р., Веселов А. И., Ястребов В. А., Тюрликов А. М.
В24
Введение в цифровую обработку изображений. Методы фильтрации и сжатия изобра-
жений: учебное пособие/ М. Р. Гильмутдинов, А. И. Веселов, В. А. Ястребов, А. М. Тюрликов;
СПб.:ГУАП, 2015. – 76 с.
ISBN 978-5-8088-1079-2
Рассматриваются методы повышения качества изображений, а также основные алгоритмы, используемые в современных стандартах сжатия JPEG и JPEG-LS. Рассмотрение теоретических вопросов сопровождается описанием практических работ, в ходе которых обучающиеся
получают практические навыки применения методов цифровой обработки изображений.
Учебное пособие может быть использовано студентами, обучающимися по направлениям “Информационная безопасность” и “Инфокоммуникационные технологии и системы связи”,
а также для самостоятельной работы аспирантов.
УДК 004.627
ББК 30
ISBN 978-5-8088-1079-2
© ГУАП, 2015
© М. Р. Гильмутдинов,
А. И. Веселов,
В. А. Ястребов,
А. М. Тюрликов, 2015
75
СОДЕРЖАНИЕ
ПРЕДИСЛОВИЕ
Предисловие
3
I. Анализ методов улучшения качества изображений
1.
Теоретические основы . . . . . . . . . . . . . . . . . . . . . .
1.1.
Основные положения . . . . . . . . . . . . . . . . . . . .
1.2.
Модели шума на изображении . . . . . . . . . . . . . . .
1.3.
Методы уменьшения шума на изображении . . . . . . . .
1.4.
Методы выделения контуров и повышения резкости изображения . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5.
Градационные преобразования для изображений . . . . .
2.
Практические задания . . . . . . . . . . . . . . . . . . . . . .
2.1.
Реализация методов шумоподавления . . . . . . . . . . .
2.2.
Реализация методов выделения контуров . . . . . . . . .
2.3.
Реализация градационных преобразований . . . . . . . .
3.
Вопросы для самопроверки . . . . . . . . . . . . . . . . . . .
II. Методы пиксельной обработки при сжатии изображений
1.
Теоретические основы . . . . . . . . . . . . . . . . . . . . . .
1.1.
Схема дифференциальной импульсной кодовой модуляции для сжатия изображений . . . . . . . . . . . . . . . .
1.2.
Методы предсказания данных . . . . . . . . . . . . . . . .
1.3.
Метод нелинейного предсказания из стандарта JPEG-LS .
1.4.
Квантование данных . . . . . . . . . . . . . . . . . . . . .
1.5.
Критерии оценки искажений . . . . . . . . . . . . . . . .
1.6.
Кодирование (сжатие) ошибок предсказания на основе кодов Голомба . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7.
Оценка эффективности сжатия . . . . . . . . . . . . . . .
2.
Практические задания . . . . . . . . . . . . . . . . . . . . . .
3.
Вопросы для самопроверки . . . . . . . . . . . . . . . . . . .
.
.
.
.
4
4
4
5
6
.
.
.
.
.
.
.
12
19
26
27
29
30
31
33
. 33
.
.
.
.
.
33
34
41
42
46
.
.
.
.
48
50
50
51
III. Методы блоковой обработки при сжатии с потерями на примере стандарта JPEG
1.
Теоретические основы . . . . . . . . . . . . . . . . . . . . . . .
1.1.
Характеристика стандарта JPEG . . . . . . . . . . . . . . .
1.2.
Требования к входным данным . . . . . . . . . . . . . . . .
74
Цифровая обработка изображений играет ключевую роль во многих современных задачах, связанных с хранением, отображением и пониманием
содержимого визуальных данных, представленных в цифровой форме. При
этом во многих прикладных областях в качестве целевого критерия используются взаимоисключающие понятия, такие как наиболее компактное представление и максимальное качество визуального восприятия. В связи с этим
понимание основных принципов компрессии и анализа особенностей изображений является основополагающим для будущих специалистов в данной
области.
Материал, изложенный в настоящем пособии, является логическим продолжением курса ранее опубликованного пособия «Цифровая обработка
изображений. Статистический анализ и квантование визуальных данных».
В пособии рассматриваются базовые задачи цифровой обработки изображений: фильтрация и сжатие. При рассмотрении методов фильтрации основное внимание уделено подходам, связанным с анализом особенностей изображений и повышением качества их визуального восприятия. Методы сжатия рассматриваются на примере подходов из популярных стандартов JPEG
и JPEG-LS, в которых используются два различных принципа компрессии:
алгоритмы сжатия с потерями, использующие блоковые преобразования, и
алгоритмы сжатия без потерь, основанные на пиксельной обработке.
53
53
53
54
3
I. АНАЛИЗ МЕТОДОВ УЛУЧШЕНИЯ
КАЧЕСТВА ИЗОБРАЖЕНИЙ
17. Кудряшов Б.Д. Теория информации. СПб.: Питер, 2009.
18. Гильмутдинов М.Р., Тюрликов А.М., Линский Е.М. Цифровая обработка изображений. Статистический анализ и квантование визуальных данных. СПб.: ГУАП, 2014.
19. Digital compression and coding of continuous-tone still images. 1991.
Основной задачей данного раздела учебного пособия является получение практических навыков использования алгоритмов улучшения качества
изображений, применяемых в пространственной области, на примерах методов шумоподавления, выделения контуров и градационных преобразований
(для упрощения анализа рассматриваются только однокомпонентные изображения).
20. ITU-T Recommendation H.263. Video Coding for Low Bitrate
Communication. 1996.
1. Теоретические основы
1.1.
Основные положения
В данном разделе будут рассмотрены общие принципы методов, улучшающих визуальное восприятие изображений. Эти методы также позволяют повысить эффективность применения алгоритмов видеоаналитики, таких
как распознавание образов, идентификация и сопровождение объектов.
Рассматриваемые методы можно разделить на несколько классов. В пособии предлагается рассматривать следующие методы:
– методы, использующие линейную и нелинейную фильтрацию, позволяющие уменьшать влияние шума, возникающего в процессе формирования
цифрового изображения;
– методы на базе дифференцирования, позволяющие выделять контуры,
специфические регионы на изображении, в которых происходит резкий
перепад яркости (такие регионы, как правило, отображают границы объектов фиксируемой на изображении сцены);
– методы на базе градационных преобразований, применяемые для повышения контраста изображения.
Будут затронуты также вопросы реализации рассматриваемых методов
с использованием арифметики с фиксированной точкой.
Изображение в цифровом формате рассматривается как двумерная функция Iy,x с дискретными аргументами, которая принимает дискретные значения в интервале от 0 до 2L − 1. L называют глубиной (bit depth). Поскольку принимаемые аргументами значения ограничены, удобнее рассматривать
4
73
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера,
2005.
изображение как двумерный массив или матрицу I[H×W ] , у которой H строк
и W столбцов. Отдельные позиции изображения называют пикселями. Значения, которые могут принимать пиксели, называют интенсивностями или
градациями цвета/яркости.
В качестве изображения рассматривается яркостная компонента с глубиной 8 бит. Для компонент RGB 24-битного файла BM P значение пикселя
яркостной компоненты на позиции (y, x) формируется независимо от других
пикселей в соответствии с формулой из стандарта JPEG/JFIF [1]:
3. Быков Р.Е. Цифровое преобразование изображений. М.: Горячая линияТелеком, 2003.
Iy,x = 0, 299Ry,x + 0, 587Gy,x + 0, 114By,x ,
1. ITU-T Recommendation T.81. Information technology – Digital
compression and coding of continuous-tone still images: JPEG File
Interchange Format (JFIF) . 2011.
4. Сойфер В.А. Методы компьютерной обработки изображений. М.: Физматлит, 2003.
5. Шапиро Л., Стокман Дж. Компьютерное зрение. М.: Бином, 2006.
6. ITU-T G.728 Coding of speech at 16 kbit/s using low-delay code excited
linear prediction. 1992.
7. Bocharova E. Compression for Multimedia. NY.: Cambrige University Press,
2010.
8. Красильников Н.Н. Цифровая обработка изображений. М.: Вузовская
книга, 2001.
9. CCITT Recommendation T.81. Information technology. Digital compression
and coding of continuous-tone still images. Requirements and guidelines.
1992.
10. Levinson N. The Wiener RMS error criterion in filter design and prediction //
J. Math. Phys. 1947. Т. 25. С. 261–278.
11. Durbin J. The fitting of time series models // Rev. Inst. Int. Stat. 1960. Т. 28.
С. 233–243.
12. Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. М.: Мир,
1989.
13. Сэломон Д. Сжатие данных, изображений и звука. М.: Техносфера,
2004.
14. Golomb S. Run-length encodings (Corresp.) // Information Theory, IEEE
Transactions on. 1966. Jul. Т. 12, № 3. С. 399–401.
(1.1)
где Iy,x , Ry,x , Gy,x , By,x — значение интенсивности пикселя на позиции
(y, x) для яркостной, красной, зеленой и голубой компонент соответственно.
В результате такого преобразования значения пикселей яркостной компоненты с глубиной 8 бит лежат в диапазоне [0, 255], однако при выполнении
некоторых других преобразований, результат может выйти за границы этого
интервала. В этом случае результат подвергается дополнительной операции
клиппирования, которая описывается формулой

min

, если Iy,x < I min

I
Clip(Iy,x , I min , I max ) = I max , если Iy,x > I max
,


I ,
min
max
если I
≤ Iy,x ≤ I
y,x
(1.2)
где I min и I max — минимальное и максимальное возможное значение пикселя рассматриваемой компоненты (0 и 255 соответственно для 8-битной компоненты).
1.2.
Модели шума на изображении
Рассмотрим модели аддитивного и импульсного шума. Модель формирования аддитивного шума для 8-битных значений интенсивностей пикселей изображения можно представить следующим образом:
15. ISO/IEC JTC 1/SC 29. Information technology – Lossless and near-lossless
compression of continuous-tone still images: Extended. 2003.
′
Iy,x
= Clip(Iy,x + Ny,x , 0, 255),
16. Joint Video Team of ITU-T and ISO/IEC JTC 1. Advanced video coding for
generic audiovisual services. Recommendation ITU-T H.264. 2009.
′
где Iy,x
— значение интенсивности компоненты I ′ , полученной в результате обработки исходной компоненты I; Iy,x — значение интенсивности, а
72
5
(1.3)
Ny,x — значение шума на позиции пикселя с координатами (y, x). Операция
Clip(·) соответствует формуле (1.2).
Модель формирования импульсного шума для 8-битных значений интенсивностей описывается следующей формулой:


с вероятностью pa

0,
′
Iy,x
=
255,
с вероятностью pb


I
,
(1.4)
с вероятностью 1 − pa − pb
y,x ,
7. Можно ли использовать кодирование длин серий для сжатия разностей
коэффициентов DC?
1.3.1 Идея линейной фильтрации изображения
В процессе выполнения фильтрации входного (зашумленного) изобра′
жения I[H×W ] формируется новое изображение I[H×W
] с теми же размера′
ми. Каждый пиксель Iy,x формируется в результате применения некоторого
оператора к пикселю Iy,x и подмножеству соседних с ним пикселей, образующих так называемую (апертуру) фильтра. Применяемый оператор может
быть как линейным, так и нелинейным. При линейной фильтрации используется следующее правило для расчета значений отфильтрованных пикселей:
R
1 ∑
Z
R
∑
wk,m Iy+k,x+m ,
4. Для чего делается перегруппировка AC коэффициентов в соответствии с
зигзагообразным обходом?
6. Раскройте специфику формирования значений шагов квантования в зависимости от квантуемой полосы.
Методы уменьшения шума на изображении
′
Iy,x
=
3. С какой целью в дискретном косинусном преобразовании используются
нормирующие коэффициенты?
5. С какой целью выполняется квантование спектральных коэффициентов?
где pa , pb — параметры модели импульсного шума.
1.3.
2. Сколько будет отрицательных и сколько положительных элементов в базисном векторе ДКП длины N при частоте f ? Значения N и f получить у
преподавателя.
8. Пусть сжимается 8-битное изображение. Чему равна битовая категория
квантованного коэффициента AC, равного 8071?
(1.5)
k=−R m=−R
где R — радиус фильтра, определяющий апертуру, wk,m — весовые коэффициенты фильтра, Z — коэффициент нормировки, рассчитывающийся как
Z=
R
∑
R
∑
wk,m .
(1.6)
k=−R m=−R
Следует отметить, что использование термина ”радиус” является общепринятым, хотя в реальности апертура имеет форму квадрата со стороной 2R+1.
На краях массива I, где апертура выходит за границы изображения, недостающее значение, как правило, заменяется интенсивностью ближайшего по
Евклидову расстоянию пикселя.
Значение и положение весовых коэффициентов определяется так называемым ядром фильтра. Ядро фильтра удобно представлять в виде маски
6
71
i = 0, ..., 7, j = 0, ..., 7
Знаки спектральных коэффициентов кодируются отдельно равномерным кодом (1 бит на знак). При оценке размера битового потока
обратите внимание, что в битовом потоке нет необходимости хранить знаки для нулевых значений.
3 Оценить распределение ошибок квантования ACi,j в каждой позиции (i, j).
– определить три матрицы квантования, каждая из которых во
всех позициях содержит одно и то же значение q1 , q2 и q3 , которые определяют «мягкое», «умеренное» и «грубое» квантование соответственно;
– выполнить процедуру квантования, пользуясь формулами (3.5) и (3.6), в процессе обработки изображений из
тестового набора;
– построить гистограммы ошибок квантования ACi,j , набрав
статистику по ошибкам в отдельной позиции (i, j);
– сделать выводы о распределении ошибок квантования в каждой позиции и предложить альтернативу кодированию по битовым категориям при сжатии коэффициентов переменного
тока.
4. Кодирование без потерь
— матрицы весовых коэффициентов. Пример такого представления можно
увидеть на рис. 1.6.
Операции с фиксированной точкой в вычислительных системах, как правило, выполняются быстрее, чем операции с плавающей точкой. Поэтому
весовые коэффициенты wk,m в выражении (1.5) на практике представляют
в формате с фиксированной точкой. Способы подобной аппроксимации весовых коэффициентов будут рассмотрены далее. Следует отметить, что при
проведении операций с фиксированной точкой в случае применения округления или клиппирования процедура фильтрации может стать нелинейной.
В данной работе рассматриваются два линейных фильтра [2, 3, 4, 5]:
– фильтр на основе скользящего среднего (усредняющий фильтр);
– фильтр Гаусса.
1.3.2 Метод скользящего среднего
В методе скользящего среднего значения весовых коэффициентов постоянны и не зависят от расстояния до центрального пикселя. Все весовые коэффициенты равны единице, в связи с этим выход фильтра рассчитывается
по следующей формуле:
′
Ii,j
=
R
∑
1
(2R + 1)2
R
∑
Ii+k,j+m .
(1.7)
k=−R m=−R
4.1 Исследование различных вариантов зигзагообразного обхода.
q
Необходимо провести статистический анализ значений ACi,j
и
предложить альтернативный зигзагобразному обходу (рис. 3.4)
вариант перегруппировки квантованных коэффициентов, который
бы обеспечил лучшее сжатие.
4.2 Исследование предсказания коэффициентов постоянного тока.
Предложить вариант предсказания значений DC q , который бы давал лучшее сжатие по сравнению с разностным кодированием в
стандарте JPEG.
4. Вопросы для самопроверки
1. Какой динамический диапазон у спектральных коэффициентов, рассчитанных для блока с размерами 4 × 4, 8 × 8 , 16 × 16 , 32 × 32 пикселя, если
входные данные занимают 8 или 16 бит?
70
Графическое представление импульсного отклика на единичный импульс приведено на рис. 1.1б. Как видно из рисунка, отклик имеет форму
квадрата, состоящего из одинаковых интенсивностей, которые по своему
значению в (2R + 1)2 раз меньше интенсивности пикселя входного изображения, представляющего собой одиночный импульс — один белый пиксель
на черном фоне. Для удобства отображения значения отклика умножены на
шкалирующий коэффициент.
Как уже было сказано, фильтр скользящего среднего вносит специфическое размытие в области резких перепадов интенсивностей на изображении.
Это хорошо заметно на рис. 1.2,б, где представлен результат фильтрации белого круга (рис. 1.2,а). Края круга после фильтрации становятся более гладкими, а сам круг увеличивается в размерах.
7
а)
б)
Рис. 1.1. Результаты применения метода скользящего среднего при использовании маски
7х7 (радиус 3). Шкалирующий коэффициент повышения значений интенсивностей отфильтрованного изображения равен 25: а — исходное изображение; б — отфильтрованное изображение
1.3.3 Фильтр Гаусса
Поскольку на фотореалистичном изображении соседние пиксели сильно коррелируют, значения весовых коэффициентов, представленных в выражении 1.5, как правило, рассчитываются на базе некоторой функции от
расстояния до центральной позиции Ii,j , значения которой убывают по мере удаления позиции коэффициента от центральной. Однако такой подход
имеет негативный эффект с точки зрения визуального восприятия. В окрестности контуров (резких перепадов) применение линейных фильтров приводит к ”размытию”контура. Поэтому данный класс фильтров также называют
фильтрами размытия (в англоязычной литературе Blur).
Функции расстояния для расчета значений весовых коэффициентов формируется на базе функции Гаусса от двух переменных (отсюда и название
фильтра):
(
)
−(k 2 + m2 )
wk,m = exp
,
(1.8)
2σ2
где σ — параметр фильтра, определяющий скорость убывания коэффициентов wk,m по мере удаления от позиции центрального пикселя.
Импульсный отклик (рис. 1.3,б), в отличие от скользящего среднего, да-
8
2.2 То же, что и в предыдущем задании, но для двумерного преобразования ДКП 4 × 4.
2.3 То же, что и в предыдущем задании, но для 12-битных входных значений.
2.4 Оценить число бит, необходимых для представления спектральных
коэффициентов на выходе ДКП, которое бы позволило организовать
преобразование без потерь данных.
2.5 Модифицировать реализации кодера и декодера для использования
ДКП 16×16.
2.6 Модифицировать реализации кодера и декодера для использования
ДКП 4×4.
3. Квантование
1 Оценить влияние позиции спектрального коэффициента в блоке 8×
8 на искажения, вносимые квантованием. Для яркостных компонент
(luma)
набора тестовых изображений построить графики PSNR(q(i,j) )
(luma)
при фиксированных значениях (i, j). Значения q(i,j) должны изменяться в интервале [1, 50]. Все элементы матрицы квантования,
отличные от (i, j), приравнять к единице. Рассмотреть следующие
пары значений (i, j): (0, 0), (1, 0), (1, 1), (4, 0), (0, 4), (4, 4), (7, 7).
Сделать вывод о влиянии позиции квантуемого спектрального коэффициента на вносимые искажения.
2 Реализация упрощенной процедуры квантования в соответствии со
стандартом ITU-T H.263 [20]. В стандарте H.263 вместо матрицы
квантования в каждой полосе используется, кроме полосы коэффициентов постоянного тока, один и тот же шаг квантования, рассчитываемый по значению QP . Процедуры квантования и деквантования необходимо реализуются по следующим формулам:
(
)
|DC| + 2
q
|DC | = round
,
4
)
(
|ACi,j |
q
, i = 0, ..., 7, j = 0, ..., 7 ,
|ACi,j
| = round
2QP
dq
|ACi,j
|



0,
=
DC dq = 4DC q ,
q
если ACi,j
= 0;
q
QP (2|ACi,j
|


QP (2|AC q |
i,j
+ 1),
если QP — нечетный;
+ 1) − 1, если QP — четный,
69
3.
Задания повышенной сложности
а)
б)
Большинство индивидуальных заданий предполагают модификацию одного или нескольких методов стандарта JPEG. В качестве результата индивидуального задания необходимо привести и проанализировать кривые
«скорость-искажение» для оригинального и модифицированного методов.
Все индивидуальные задания необходимо выполнить для всех изображений
из тестового набора.
1. Предварительная обработка изображений
1.1 Проанализировать влияние процедуры прореживания (децимации)
цветоразностных компонент на характеристики сжатия и вносимые
искажения. Для анализа необходимо построить и сравнить графики «скорость-искажение» для компонент R, G, B, полученных для
входных данных в форматах 4:4:4 и 4:2:0. Процедуру децимации
осуществлять путем исключения строк и столбцов с нечетными номерами (нумерация начинается с нуля). Восстановление исходного
размера путем копирования данных из ближайшей позиции, находящейся слева или сверху.
1.2 То же, что и в предыдущем задании, но прореживание выполнять
путем замены четырех смежных пикселей на их среднее значение.
1.3 То же, что и в предыдущем задании, но вместо прореживания использовать обработку цветоразностных компонент линейным фильтром на основе скользящего среднего. Для вычисления нового значения компоненты на позиции (i, j) использовать следующую формулу:
I ′ (i, j) =
R
∑
1
(2R + 1)2
R
∑
I(i + m, j + n) ,
m=−R n=−R
где R — радиус фильтра (использовать значения 1 и 2); I(i + m, j +
n) — значение исходной компоненты на позиции (i + m, j + n).
2. Дискретное косинусное преобразование
2.1 Определить минимальное и максимальное значения для спектральных коэффициентов, которые возможно получить на выходе одномерного ДКП длины 4 при условии, что входные данные являются
целыми числами в 8-битном беззнаковом представлении.
68
Рис. 1.2. Результаты применения метода скользящего среднего при использовании маски
7х7 (радиус 3): а — исходное изображение; б — отфильтрованное изображение
ет форму более близкую к точке. Аналогично уменьшаются искажения при
фильтрации изображения круга (рис. 1.4). Еще один пример влияния усредняющего фильтра и фильтра Гаусса на исходное изображение приведен на
рис. 1.5.
1.3.4 Переход к арифметике с фиксированной точкой
Изначально большинство рассмотренных алгоритмов разрабатывалось
для арифметики с фиксированной точкой. В процессе перехода от представления в формате с плавающей точкой для формирования целочисленной маски фильтра можно выделить следующие основные этапы:
– выбор радиуса фильтра R;
– определение коэффициентов фильтра wk,m и коэффициента нормировки
Z, как суммы всех весовых коэффициентов;
– вычисление коэффициентов Wk,m в формате с фиксированной точкой:
⌊
⌋
1
Wk,m = 2b wk,m ,
(1.9)
Z
где b определяет число бит после запятой, которое должно использоваться
при целочисленных вычислениях;
9
а)
б)
Описанный алгоритм применяется для всех пикселей исходного изображения I[H×W ] . На краях изображения, отсутствующие в апертуре значения, как
и в случае линейной фильтрации, заполняются некоторым образом, например, значением ближайшего по Евклидову расстоянию существующего пикселя.
графики PSNR(Q(luma) ) для яркостной и PSNR(Q(chroma) ) для обеих цветоразностных компонент. Обе матрицы квантования построить по формуле (3.7) с одинаковыми значениями параметра R. По
построенным графикам сделать выводы об эффективности квантования яркостной и цветоразностных компонент при одних и тех же
матрицах квантования.
3. Сжатие без потерь
3.1 Реализовать процедуру кодирования квантованных коэффициентов
постоянного тока DC, включающую в себя:
3.1.1. Вычисление разностей ∆DC в соответствии с формулой (3.8).
3.1.2. Разбиение
значений
∆DC
на
пару
(
)
BC(∆DC ) , M agnitude(∆DC ) .
3.2 Оценить эффективность использования разностного кодирования
для коэффициентов постоянного тока.
3.2.1. Построить гистограммы частот f (DC q ) и f (∆DC ), полученных для блоков 8 × 8 яркостной составляющей тестового изображения.
3.2.2. Вычислить оценки энтропии, воспользовавшись данными из
множеств {DC q } и {∆DC }.
3.3 Реализовать процедуру кодирования длинами серий (Run-Length
Encoding — RLE), которая для каждого блока 8 × 8 из квантованных коэффициентов переменного тока AC включает в себя:
3.3.1. Перегруппировку AC q коэффициентов в соответствии с зигзагообразной последовательностью.
3.3.2. Замену сформированной последовательности на последовательность пар (Run, Level).
3.3.3. Формирование для каждой пары (Run, Level) на основе значения Level новой пары — (BC(Level), M agnitude(Level)).
3.3.4. Формирование троек ((Run, BC(Level)) , M agnitude(Level)).
3.4 Определить соотношение размеров в сжатом битовом потоке для
следующих данных:
– BC(∆DC ) ,
– M agnitude(∆DC ) ,
– (Run, BC(Level)),
– M agnitude(Level).
Соотношение определить в процентах к общей длине оценки битового потока.
10
67
Рис. 1.3. Результаты применения фильтра Гаусса при использовании маски 7х7 (радиус 3).
Шкалирующий коэффициент повышения значений интенсивностей отфильтрованного изображения равен 25: а — исходное изображение; б — отфильтрованное изображение
– обновление нормирующего коэффициента Z, как суммы всех весовых коэффициентов Wk,m . Нормирующий коэффициент должен быть равен 2b и
деление на него реализуется через сдвиг на b бит вправо.
Пример целочисленной маски фильтра Гаусса приведен на рис. 1.6.
1.3.5 Фильтры, основанные на порядковых статистиках
Фильтры, основанные на порядковых статистиках относятся к одному из
частных случаев нелинейной фильтрации. Рассмотрим медианный фильтр,
наиболее распространенный в обработке изображений алгоритм данного
класса. Алгоритм медианной фильтрации для пикселя Ii,j можно представить в виде следующих шагов:
– формирование одномерного массива A из апертуры пикселя Ii,j ;
– взятие медианы одномерного массива A.
– выполнить все пункты из списка общих заданий, представленных далее;
а)
б)
– выполнить один или несколько пунктов из списка индивидуальных заданий, предварительно согласовав их с преподавателем;
– в качестве меры искажений использовать PSNR для компонент Y , Cb, Cr
при построении графиков «скорость-искажение», если иное не уточняется в задании (PSNR необходимо вычислять между входными значениями
компонент и восстановленными после обратного ДКП);
– в качестве оценки размера битового потока вместо кодирования по Хаффману использовать оценку одномерной энтропии, выполняемую по формуле
∑
H=−
p̂(x) · log2 p̂(x) ,
(3.9)
x
а также формулой
Ix = ⌈− log2 p̂(x)⌉
(3.10)
Рис. 1.4. Результаты применения фильтра Гаусса при использовании маски 7х7 (радиус 3):
а — исходное изображение, б — отфильтрованное изображение
для верхней оценки числа бит кодового слова кода Хаффмана.
1. Дискретное косинусное преобразование (ДКП)
1.1 Реализовать процедуру прямого и обратного ДКП для блоков N ×N .
Выход обоих преобразований сделать целочисленным с помощью
операции округления.
1.2 Оценить искажения, вносимые ДКП.
1.2.1. Пиксели яркостной компоненты тестового изображения разбить на непересекающиеся блоки 8 × 8.
1.2.2. Над каждым блоком выполнить преобразование ДКП и получить целочисленные значения спектральных коэффициентов.
1.2.3. Над полученными спектральными коэффициентами выполнить обратное ДКП с последующим округлением выходных
данных.
1.2.4. По исходной и восстановленной после обратного ДКП яркостным компонентам вычислить значение PSNR.
2. Квантование спектральных коэффициентов
2.1 Реализовать процедуру квантования и деквантования в соответствии с формулами (3.5) и 3.6.
2.2 Оценить влияние искажений, вносимых квантованием, на компоненты изображений. Для набора тестовых изображений построить
Рассмотрим принципы нахождения медианы на примере произвольного
одномерного массива A = a1 , ..., aN , где ai ∈ R. Можно выделить два типовых способа вычисления медианы median(A). При использовании первого
способа следует выполнить следующие шаги:
– произвести сортировку элементов одномерного массива A (неважно по
возрастанию или по убыванию);
′
– сохранить в Ii,j
элемент, находящийся в середине отсортированного массива A. Для апертуры радиуса
(⌊R необходимо
⌋
) взять из отсортированного
2
массива элемент с индексом (2R+1)
+
1
.
2
Рассмотрим второй способ нахождения медианы. Пусть a[N /2] — N /2
порядковая статистика массива A. Из свойств медианы известно, что
N
∑
|ai − a[N /2] | ≤
i=1
|ai − x|, ∀x ∈ ℜ.
(1.10)
i=1
Таким образом, срединное значение median(A) для рассматриваемого массива A можно получить путем решения следующей оптимизационной задачи:
median(A) = a[N /2] = min
x∈R
66
N
∑
N
∑
∥x − ai ∥1 = min
x∈R
i=1
11
N
∑
i=1
|x − ai |,
(1.11)
а)
б)
в)
Рис. 1.5.
Результаты применения усредняющего фильтра, а также фильтра Гаусса при
использовании маски 7х7 (радиус 3): а — исходное изображение; б — изображение после применения усредняющего фильтра; в — изображение после применения фильтра Гаусса
Рис. 1.6. Пример целочисленного ядра фильтра Гаусса (нормирующий коэффициент равен
16)
где ∥ · ∥1 — обозначает l1 норму.
Следует отметить, что первый способ является предпочтительным, в
случае когда размер исходного массива невелик. Применение второго способа более эффективно при поиске медианы на фильтрах больших размеров.
1.4. Методы выделения контуров и повышения резкости изображения
1.4.1 Понятие первой и второй производной функции целочисленного
аргумента
Принцип действия типовых фильтров рассматриваемого класса алгоритмов обработки изображений основан на использовании свойств первой и второй производных. Рассмотрим понятие производной функции от дискретного аргумента. Производную подобных функций принято определять в терминах разностей, которые могут быть аппроксимированы различными способами. К результату аппроксимации предъявляются следующие требования.
Первая производная должна быть:
12
Значение ∆DC представляется в форме битовой категории и амплитуды.
Значение битовой категории кодируется по Хаффману. Код Хаффмана для
битовой категории ∆DC должен сохранятся в заголовочной части файла.
2. Перед кодированием коэффициенты AC q переупорядочиваются в соответствии с зигзагообразной последовательностью, показанной на рис. 3.4.
При условии, что энергия блока сконцентрирована в левом верхнем углу
блока, новая последовательность будет обладать следующим свойством
— значимые (ненулевые) коэффициенты будут преобладать в ее начале, а
ближе к концу последовательности в основном будут находиться нулевые
значения.
3. Сформированная на предыдущем шаге последовательность коэффициентов AC q кодируется длинами серий. Результатом будет являться новая
последовательность, состоящая из пар (Run/Level). Здесь Run определяет число нулевых AC q коэффициентов в серии. Каждая серия завершается Level — ненулевым значением AC q коэффициента. Обработка завершается, когда закодирован последний ненулевой AC q коэффициент. При
этом в поток вставляется дополнительная специальная пара (0, 0), которая
носит название End Of Block (EOB). Если длина серии, которая не является
последней в блоке, превышает 15 нулей, то она разбивается на несколько
подсерий, каждая из которых содержит не более 16 нулей. Подпоследовательность из шестнадцати нулей кодируется символом (15, 0). Такая пара
носит специальное название Zero Run Length (ZRL).
4. Значение Level, если оно не относится к паре ZRL или EOB, кодируется по
битовым категориям, формируя пару (BC(Level), M agnitude(Level)).
Дальнейшему кодированию подвергается по Хаффману пара
(Run, BC(Level)), за которой в битовый поток записывается значение амплитуды M agnitude(Level), состоящее из BC(Level) бит.
2. Практические задания
В данном разделе входной информацией для кодера являются компоненты цветного изображения в формате Y CbCr 4:4:4 (без децимации) по 8 бит
на пиксель в каждой компоненте. В качестве тестовых изображений необходимо использовать lena.bmp, peppers.bmp и baboon.bmp. В ходе выполнения
задания необходимо:
– реализовать упрощенную модель кодера и декодера по стандарту JPEG в
соответствии со списком общих заданий;
65
стандарта JPEG при кодировании квантованных спектральных коэффициентов используется второй подход, поскольку наиболее часто повторяющимся
элементом в кодируемой последовательности является нулевое значение.
1.6.4 Основные этапы кодирования
Кодирование каждой компоненты осуществляется независимо. Обход
блоков 8 × 8 каждой компоненты изображения осуществляется в сканирующем порядке (рис. 3.2). Для каждого блока выполняются следующие действия.
– равной нулю на плоских участках(областях постоянной яркости);
– ненулевой в начале и в конце ступеньки или склона яркости (участках изменения яркости);
– ненулевой на склонах яркости.
Аналогично, вторая производная должна быть:
– равной нулю на плоских участках;
– ненулевой в начале и в конце ступеньки или склона яркости;
– равной нулю на склонах постоянной крутизны.
Будем полагать, что приращение аргумента рассматриваемой дискретной функции равно единице. Тогда производная одномерной функции dI(x)
dx
совпадает с приращением функции dI(x) и определяется как разность значений соседних элементов [2]:
dIx = Ix+1 − Ix .
(1.12)
Аналогично, вторая производная определяется как разность соседних
значений первой производной:
Рис. 3.4. Зигзагообразный обход при обработке AC коэффициентов блока 8 × 8
1. Кодирование коэффициента постоянного тока DC q .
2. Перегруппировка 63 коэффициентов переменного тока AC q и формирование одномерного массива в соответствии с зигзагообразной последовательностью.
3. Применение кодирования длин серий для последовательности из 63 AC q
коэффициентов.
4. Кодирование пар (Run, Level).
Рассмотрим подробно каждый из перечисленных шагов.
1. Для кодирования DC q используется разностный метод. Дальнейшей обработке подвергается разность DC q коэффициента текущего и предыдущего обрабатываемого блоков:
d2 Ix = Ix+1 − Ix − (Ix − Ix−1 ) = Ix+1 − 2Ix + Ix−1 .
(1.13)
Один из видов аппроксимации первой производной — разность элементов, расположенных на расстоянии одного пикселя от позиции, для которой
ищется производная:
dIx = Ix+1 − Ix−1 .
(1.14)
Подобный подход, в частности, используется при формировании фильтров
Превитт и Собеля [2, 5]. Фильтр Собеля будет рассмотрен далее.
а)
б)
Рис. 1.7. Пример масок для расчета аппроксимации первой производной функции целочисленного аргумента: а — маска для шага в один пиксель; б — маска для шага в два пикселя
q
∆DC = DCiq − DCi−1
,
(3.8)
где i — номер текущего обрабатываемого блока. Для первого блока значение ∆DC вычисляется как разность DC0q и среднего по всем возможным
значениям DC q .
64
Легко показать, что процедуру вычисления производных по формулам (1.12) и (1.14) можно реализовать в виде фильтров, маски которых приведены на рис. 1.7.
13
1.4.2 Оператор Лапласа
В задаче увеличения резкости изображения требуется подчеркнуть контуры изображения. Для ее выполнения удобнее пользоваться второй производной. Рассмотрим применение двумерной второй производной в задаче
улучшения изображений. Подход сводится к построению второй производной для дискретного случая и к последующему построению соответствующего ядра фильтра. Простейшим оператором, основанным на второй производной, является оператор Лапласа (Лапласиан). В случае обработки функции двух переменных Ix,y он определяется следующим образом:
∇2 Iy,x =
∂ 2 Iy,x
∂ 2 Iy,x
+
.
2
∂x
∂y 2
(1.15)
Поскольку теперь имеются две переменные, вторые частные производные рассчитывают по следующим формулам:
∂ 2 Iy,x
= Iy,x+1 − 2Iy,x + Iy,x−1 ;
∂x2
(1.16)
∂ 2 Iy,x
= Iy+1,x − 2Iy,x + Iy−1,x .
∂y 2
(1.17)
Путем объединения (1.16), (1.17) и (1.15) можно определить двумерный
Лапласиан в дискретной форме:
(
)
∇2 Iy,x = Iy,x+1 + Iy,x−1 + Iy+1,x + Iy−1,x − 4Iy,x .
Амплитудой фактически является само кодируемое число. Такое представление может показаться избыточным, однако оно необходимо для упрощения последующего формирования и декодирования битового потока. С
помощью битовой категории легко определить минимальное количество значимых бит, которое необходимо для представления произвольного числа в
двоичном формате. Типовой алгоритм формирования битового потока для
кодирования по битовым категориям формулируется следующим образом:
1. Сначала в битовый поток укладывается код битовой категории BC. В качестве кода для битовой категории в базовом режиме стандарта JPEG используется код Хаффмана. Сам код Хаффмана хранится в заголовочной
части файла.
2. Вслед за кодом битовой категории в поток укладывается битовое представление амплитуды M G с помощью равномерного кода. Для представления амплитуды затрачивается ровно BC бит. Если значение M G отрицательное, то для представления амплитуды используется обратный код
(биты представления амплитуды инвертируются).
Декодер легко может восстановить исходное значение, прочитав из битового потока сначала код Хаффмана для битовой категории BC . Далее продекодировав код Хаффмана, декодер извлекает дополнительно BC бит, которые интерпретирует как амплитуду M G. Если старший бит M G имеет нулевое значение, то равномерный код должен быть проинвертирован, чтобы
получить окончательное значение коэффициента.
(1.18)
Маски фильтров, реализующих правило (1.18), представлены на рис. 1.8.
Нормировка для таких фильтров не требуется.
1.6.3 Кодирование длин серий
Результат применения фильтра на базе оператора Лапласа будет отличен
от нуля только в окрестностях контуров (существенных перепадов интенсивностей) на изображении. Применение Лапласиана для усиления перепадов
Кодирование длин серий (Run Length Encoding — RLE) применяется в случае, когда символы на выходе кодируемого источника часто повторяются. Тогда исходная последовательность символов источника s1 , s2 , s3 , ..., sN заменяется на последовательность пар вида
(r(1) , s(1) ), (r(2) , s(2) ), ..., (r(M ) , s(M ) ), где M ≤ N , r(i) — число, указывающее количество подряд идущих символов, а затем код самого повторяющегося символа s(i) . Количество повторов также называют серией.
В случае, когда на выходе источника часто повторяется только один символ S из алфавита, исходная последовательность заменяется на последовательность пар (r(i) , l(i) ), где r(i) — длина i-й серии символа S, а l(i) — код
символа, который прервал серию из символов S. Заметьте, что длины серий символа S в этой реализации могут быть нулевыми. В базовом режиме
14
63
Рис. 1.8. Маски фильтра Лапласа
и кодирование длин серий. Далее представлено описание последовательности шагов кодирования.
1.6.2 Кодирование битовых категорий
Кодирование битовых категорий предполагает представление целого
числа x в виде пары чисел: битовой категории BC (bit category) и амплитуды M G (magnitude). Битовая категория числа x вычисляется по следующей
формуле:
BC(x) = ⌈log2 (|x| + 1)⌉ .
Определить к какой битовой категории относится обрабатываемое значение
можно и с помощью табл. 3.1, строки которой легко построить с помощью
правила:
{BC == i} : [−2i + 1, −2i−1 ], [2i−1 , 2i − 1] ,
исходного изображения I сводится к следующему действию:
{
I + ∇2 I, если w0,0 > 0
′
I =
,
I − ∇2 I, иначе
где w0,0 — значение центрального элемента маски Лапласиана.
Рис. 1.9. Примеры масок фильтров подчеркивания контуров на базе оператора Лапласа
а)
где i — номер строки.
(1.19)
б)
Таблица 3.1.
Битовые категории целых чисел
Категория
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Амплитуда
—
-1, 1
-3, -2, 2, 3
7,…, -4, 4,…7
-15,…, -8, 8,…15
-31,…-16, 16,…31
-63,…-32, 32,…63
-127,…-64, 64,…127
-255,…-128, 128,…255
-511,…-256, 256,…511
-1023,…-512, 512,…1023
-2047,…-1024, 1024,…2047
-4095,…-2048, 2048,…4095
-8191,…-4096, 4096,…8191
-16383,…-8192, 8192,…16383
-32767,…-16384, 16384,…32767
-32768
62
Рис. 1.10. Результат применения оператора Лапласа на изображении ’lena’: а — исходное
изображение; б — изображение после фильтрации
Применение метода обработки, соответствующего формуле (1.19), можно представить в виде суммы двух линейных операторов, применяемых к
изображению. Первый из них — тождественный оператор, размер ядра которого совпадает с ядром оператора Лапласа. Второе же слагаемое ∇2 I можно
рассматривать как исходное изображение I, обработанное с помощью фильтра на базе оператора Лапласа. Благодаря свойству линейности данную процедуру можно представить в виде одного фильтра, ядро которого является
15
а)
б)
изображения. Некоторые источники [13] для упрощения процедуры построения матриц квантования приводят следующую формулу:
Y
qi,j
(R) = 1 + (i + j) · R,
Рис. 1.11. Результат применения оператора Лапласа на изображении ’peppers’: а — исходное изображение, б — изображение после фильтрации
суммой рассматриваемых ядер. Примеры масок таких фильтров показаны
на рис. 1.9. Примеры обработки изображений с помощью оператора Лапласа представлены на рис. 1.10 и 1.11.
¯
I s = I − I,
(1.20)
(3.7)
где R — целочисленный параметр, управляющий качеством обработки. Очевидно, чем больше значение R, тем больше формируемые шаги квантования и хуже качество восстановленного изображения, но выше степень сжатия. Стандарт CCITT T.81 [9]) содержит два примера матриц квантования,
один — для яркостной составляющей и второй — для цветоразностных:
16 11 10 16
24
40
51
61 12 12 14 19
26
58
60
55 14 13 16 24
40
57
69
56 51
87
80
62 14 17 22 29
Q(luma) = 18 22 37 56
68 109 103
77 24 36 55 64
81 104 113
92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103
99 1.4.3 Фильтрация с подъемом высоких частот
Процедура, заключающаяся в вычитании из изображения его расфокусированной копии, называется нерезким маскированием [2]. Такой подход
применялся в фотографии и фотопечати до появления цифровых камер. Эту
процедуру можно описать следующим образом:
i = 0, ..., 7, j = 0, ..., 7 ,
Q(chroma)
=
17
18
24
47
99
99
99
99
18
21
26
66
99
99
99
99
24
26
56
99
99
99
99
99
47
66
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
99
1.6. Кодирование без потерь
s
где I — изображение с повышенной резкостью, полученное в результате
нерезкого маскирования, I¯ — расфокусированная копия исходного изображения I.
Некоторым обобщением нерезкого маскирования можно считать фильтрацию с подъемом высоких частот [2]. Изображение I ′ , полученное в результате применения фильтра с подъемом высоких частот, можно описать
следующим образом:
¯
I ′ = αI − I,
(1.21)
16
1.6.1 Общие сведения
Кодирование без потерь является завершающим этапом работы кодера
JPEG, на котором происходит формирование битового потока. Этот этап
в англоязычной литературе (в частности в стандартах JPEG [19] и CCITT
T.81 [9]) также называют энтропийное кодирование (entropy coding). Прежде
чем приступить к описанию этапов кодирования без потерь, рассмотрим два
из применяемых вспомогательных методов: кодирование битовых категорий
61
для двух цветоразностных компонент. Qluma , как правило, имеет индекс 0,
а Qchroma — индекс 1.
Идея выбора значений шагов квантования опирается на неодинаковое
психовизуальное восприятие искажений в разных полосах спектра, которое
можно описать следующей последовательностью утверждений.
1. Для фотореалистичных изображений энергия спектра блока, как правило,
концентрируется в левом верхнем углу матрицы спектральных коэффициентов, где определены полосы низких частот спектра ДКП. Таким образом среди коэффициентов расположенных вблизи левого верхнего угла
будут преобладать большие значения. В области высоких частот (ближе
к правому нижнему углу блока), как правило, значения спектральных коэффициентов близки к нулю.
2. Грубое (с большим шагом) квантование будет приводить к обнулению
большинства высокочастотных элементов спектра блока.
где α ≥ 1. Преобразуем представленное выражение:
¯
I ′ = (α − 1)I + I − I.
Используя уравнение (1.20), можно получить итоговую формулу для вычисления результата фильтрации изображения с подъемом высоких частот:
I ′ = (α − 1)I + I s .
(c)
где шаг квантования qi,j является элементом соответствующей матрицы Q(c)
в зависимости от текущей обрабатываемой компоненты. В дальнейшем коq
дироваться будут именно номера квантов Yi,j
, а аппроксимирующие значеdq
ния каждого кванта Yi,j будут вычисляться при декодировании как произведение номера кванта и шага квантования:
(c)
dq
q
Yi,j
= Yi,j
· qi,j , i = 0, ..., 7, j = 0, ..., 7,
0 -1 0
-1 +4 -1
0 -1 0
Рис. 1.12. Маска фильтра изображения с подъемом высоких частот
Фильтрация с подъемом высоких частот также может быть реализована
за один проход при использовании маски, приведенной на рис. 1.12. Параметр α, как правило, варьируется в интервале [1, 1.5]. Следует отметить, что
при α = 1 фильтрация с подъемом высоких частот сводится к задаче повышения резкости с помощью Лапласиана. С ростом параметра α эффект
увеличения резкости падает.
Результаты фильтрации изображений с подъемом высоких частот представлены на рис. 1.13 и 1.14.
(3.6)
Эту процедуру также называют деквантованием и обычно обозначают Q−1 ,
как операцию, обратную процедуре квантования Q.
Построение эффективных матриц квантования, обеспечивающих высокую степень сжатия при минимуме визуальных искажений, является нетривиальной задачей. Результат также зависит от специфики обрабатываемого
60
(1.23)
В случае использования Лапласиана в качестве способа получения резкого изображения, уравнение (1.23) можно представить в следующем виде:
{
(α − 1)I + ∇2 I, если w0,0 > 0
′
I =
,
(1.24)
(α − 1)I − ∇2 I, иначе
3. В соответствии с известными наблюдениями, искажения в низкочастотных полосах визуально более заметны, чем искажения той же интенсивности в высокочастотной части спектра. Поэтому для полос, соответствующих низким частотам, шаги квантования выбираются небольшими.
Формально процедуру квантования спектральных коэффициентов Yi,j
можно определить следующим образом:
)
(
Yi,j
q
Yi,j = round
, i = 0, ..., 7, j = 0, ..., 7 ,
(3.5)
(c)
qi,j
(1.22)
1.4.4 Оператор Собеля
Оператор Собеля является ключевым во многих алгоритмах анализа контуров изображения. В основе оператора Собеля лежит расчет двух производных: по вертикали и по горизонтали. Соответствующие этим операциям
маски фильтров приведены на рис. 1.15.
17
а)
б)
1.4.3 Особенности дискретного косинусного преобразования в
стандарте JPEG
Рис. 1.13. Результат фильтрации с подъемом высоких частот (при параметре α = 1.5)
изображения ’lena’: а — исходное изображение; б — изображение после фильтрации
В простейшем случае для определения наличия контура на позиции
(y, x) значение длины рассчитанного градиента |∇Iy,x | сравнивается с некоторым порогом thr. При превышении значения порога на позиции рассматриваемого пикселя детектируется контур. Рассмотрим алгоритм расчета значения длины градиента |∇Iy,x | и его направления θy,x для пикселя Iy,x .
1. Рассчитать значения откликов Ghy,x и Gvy,x , используя фильтры с максами,
показанными на рис. 1.15 .
2. Рассчитать величину силы (длины) контура |∇Iy,x | по формуле
√
|∇Iy,x | = (Ghy,x )2 + (Gvy,x )2 .
(1.25)
3. Рассчитать направление градиента θy,x по формуле
θy,x = arctg
Gvy,x
.
Ghy,x
(1.26)
Результат сравнения значений градиентов с порогом наглядно представить в виде бинарной карты, пример которой изображен на рис. 1.16 и 1.17.
Здесь белым отмечены позиции, отнесенные к контурам на изображении.
18
В стандарте JPEG размерность двумерного преобразования N ×N равна
8. Результатом ДКП является матрица Y размером 8×8. Ее элементы принято
называть спектральными коэффициентами. Коэффициент на позиции (0, 0)
принято называть коэффициентом постоянного тока. Его обозначают DC
от английского Direct Current, поскольку его значение получено в результате использования двух функций косинусов нулевой частоты: по горизонтали
и по вертикали. Остальные спектральные коэффициенты называют коэффициентами переменного тока и обозначают AC от английского Alternating
Current.
Каждая позиция (k, l) в матрице Y соответствует своей полосе (band)
спектра. Другими словами, значения спектральных коэффициентов всех блоков в компоненте, находящихся на одной и той же позиции (k, l), образуют
полосу (k, l) спектра.
1.5. Квантование спектральных коэффициентов
В стандарте JPEG используется равномерное скалярное квантование.
Подробнее с методами квантования можно ознакомиться в учебном пособии
[18]. В данном разделе излагаются некоторые аспекты, касающиеся выбора
шагов квантования и технической реализации этой процедуры в стандарте
JPEG.
Стандарт JPEG предполагает использование индивидуального шага
квантования для каждой полосы (k, l). Шаги квантования для всех полос
объединяются в матрицу квантования Q, которая также имеет размерность
8 × 8. В битовом потоке может содержаться до четырех различных таблиц
квантования, каждая из которых задается в специальной секции, именуемой
DQT. Каждая секция DQT начинается со специального двухбайтового кода
’FFDB’. Каждая таблица квантования имеет свой уникальный индекс (селектор) Tq , числовое значение от 0 до 3. Это значение используется для ссылки
на таблицу в заголовках кадра, компонент или сканов. Стандарт не предусматривает таблиц квантования по умолчанию, т. е. в заголовках битового
потока должна быть определена хотя бы одна таблица квантования. Поскольку статистики спектров цветоразностных компонент в целом похожи и сильно отличаются от спектра яркостной компоненты, на практике используют
две таблицы квантования: Qluma для яркостной составляющей и Qchroma
59
в пространстве (времени) t-го отсчета, для которого вычисляется значение
косинуса:
(2t + 1)π
θt =
.
2N
√
Нормирующий коэффициент Cf зависит от частоты f , для которой он вычисляется:
{
1
N , если f = 0
Cf =
.
(3.1)
2
N , иначе
а)
б)
Таким образом, значения матрицы преобразования T вычисляются по
формуле
(
)
√
(2t + 1)π
tf,t = Cf · cos
f , f = 0, ..., N − 1, t = 0, ..., N − 1. (3.2)
2N
Обратное преобразование в соответствии со свойством ортонормированности выполняeтся путем умножения на T T , транспонированную матрицу прямого ДКП.
Прямое и обратное ДКП для двумерного случая, в соответствии со свойствами ортонормированности и сепарабельности, в матричной форме выглядит следующим образом:
Рис. 1.14. Результат фильтрации с подъемом высоких частот (при параметре α = 1.5)
изображения ’peppers’: а — исходное изображение; б — изображение после фильтрации
а)
б)
Y = (T · X) · T T ;
(
)
X = T T · Y · T,
где X и Y являются матрицами размерности N × N . Первое произведение
следует интерпретировать как выполнение преобразования по столбцам, в
то время как второе — преобразование по строкам. Окончательные формулы прямого и обратного преобразований, используемые в стандарте JPEG
выглядят следующим образом:
(
)
(
)
N
−1 N
−1
∑
∑
√ √
(2i + 1)π
(2j + 1)π
yk,l = Ck Cl ·
xi,j cos
k cos
l , (3.3)
2N
2N
i=0 j=0
k = 0, ..., N − 1; l = 0, ..., N − 1;
xi,j =
N
−1 N
−1 √
∑
∑
k=0 l=0
Ck
)
(
)
(
√
(2j + 1)π
(2i + 1)π
k cos
l , (3.4)
Cl · yk,l cos
2N
2N
Рис. 1.15. Маски фильтров, используемые в операторе Собеля для аппроксимации частных
производных: а — по горизонтали; б — по вертикали
1.5. Градационные преобразования для изображений
1.5.1 Понятие градационного преобразования
Назовем множество пикселей, формирующих изображение, пространственной областью. Пространственные методы повышения качества изображений оперируют непосредственно значениями этих пикселей и могут
описываться следующим выражением:
I ′ = T [I],
i = 0, ..., N − 1; j = 0, ..., N − 1.
58
19
(1.27)
а)
б)
преобразования спектральные коэффициенты подвергаются процедуре равномерного скалярного квантования, после чего кодируются в соответствии с
процедурой, описанной в подразд. 1.6.
1.4. Дискретное косинусное преобразование
1.4.1 Свойства дискретного косинусного преобразования
1. Ортонормированность. Если T — ортогональная матрица преобразования, то должно выполнятся следующее равенство T T T = T T T = E, где
E — обозначает единичную матрицу.
Рис. 1.16.
Результат применения оператора Собеля для обнаружения границ на изображении ’lena’ при thr = 127: а — исходное изображение; б — изображение после фильтрации
где I — входные данные, I ′ — обработанные (улучшенные) данные, T [·] —
оператор над I, определенный в некоторой окрестности точки с координатами (y, x). Один из основных подходов к определению окрестности (апертуры) вокруг точки (y, x) заключается в использовании квадратной или прямоугольной области на изображении, центрированной в точке (y, x). Простейшую форму оператора T можно получить в случае использования вырожденной окрестности, когда преобразование использует только интенсивность
пикселя на текущей, обрабатываемой позиции. В этом случае T становится
оператором градационного преобразования, которое определяется функцией
отображения:
′
Iy,x
= T (Iy,x ).
(1.28)
В простейшем случае функция отображения не зависит от местоположения обрабатываемого пикселя и неизменна для всего изображения. Градационное преобразование можно представить в виде “черного” ящика с отсутствующей памятью, который для конкретного входного значения формирует соответствующее выходное на основе заранее определенной функции
(рис 1.18).
Значения функции отображения принято хранить в виде так называемой
таблицы отображения (Look-up table), представленной в виде одномерного
массива. Размер такого массива определяется числом градаций, которые ис-
20
2. Локализация энергии. Преобразование локализует энергию сигнала в низкочастотных коэффициентах результата преобразования. Для двумерного
случая — это верхний левый угол,
3. Декорреляция коэффициентов. Коэффициенты ДКП меньше коррелируют между собой, чем исходные значения пикселей в пространственной
области,
4. Сепарабельность преобразования. Двумерное преобразование сводится
к последовательному выполнению одномерного преобразования по строкам и по столбцам (или в обратном порядке). Данное свойство позволяет
снизить вычислительную сложность двумерного преобразования.
Также следует отметить, что в случае выполнения преобразования над вещественными данными, результат останется вещественным.
1.4.2 Построение матрицы дискретного косинусного преобразования
Процедуру выполнения ДКП нагляднее описывать в форме матричного умножения. Одномерный случай ДКП размерности N представляется как
умножение вектора-столбца x̄ длины N , описывающего сигнал в пространственной (временной) области, на матрицу преобразования T размерности
N × N:
ȳ = T · x̄.
Матрицу преобразования также принято называть ядром (kernel). Строки
√
матрицы T состоят из векторов, образованных значениями косинусов Cf ·
cos (θt · f ). Здесь f соответствует номеру строки и является значение частоты косинуса; Cf является нормирующим коэффициентом, а θt — положение
57
1. Без перемежения (non-interleaved). В этом режиме компоненты обрабатываются последовательно, одна за другой: сначала все блоки первой компоненты, потом второй и т.д. Для режима non-interleaved MCU — это блок
8×8 пикселей. MCU обрабатываются в порядке, показанном на рис. 3.2.
а)
б)
2. С перемежением (interleaved). В этом режиме MCU объединяет несколько
блоков, относящихся к разным компонентам. Структуры MCU для форматов Y CbCr 4:2:2 и 4:2:0 показаны на рис. 3.3,а и 3.3,б соответственно.
Рис. 1.17. Результат применения оператора Собеля для обнаружения границ на изображении ’peppers’ при thr = 127: а — исходное изображение; б — изображение после фильтрации
Рис. 3.2.
Обход блоков 8 × 8 при энтропийном кодировании. Последовательный режим
(non-interleave)
а)
б)
пользуются для представления интенсивностей пикселей. Например, для 8битного представления таблица отображения будет состоять из 256 элементов.
В данном подразделе будут рассмотрены следующие простые методы
градационных преобразований:
– с функцией преобразования на базе опорных точек;
– гамма преобразование;
– выравнивание гистограмм.
1.5.2 Метод на базе опорных точек
Рис. 3.3. Перемежающийся (interleaved) режим обработки данных: а — формат 4:2:2; б —
Одним из ключевых алгоритмов, применяемых в блоковых кодеках, является алгоритм спектрального преобразования. В стандарте JPEG используется разновидность дискретного косинусного преобразования (ДКП), которая получила название DCT-IV. Это преобразование будет рассмотрено в
следующем подразделе. Дискретное косинусное преобразование в стандарте JPEG также, как и блок, имеет размер 8 × 8. Полученные в результате
Один из простейших способов улучшения качества изображений на основе анализа и видоизменения гистограммы основывается на использовании
метода на базе опорных точек. Следует отметить, что для его применения
требуется участие человека, поскольку сами опорные точки задаются вручную. Процедуры, основанные на выравнивании гистограмм, а также гаммапреобразовании выполняются в автоматическом режиме.
Пример функции преобразования с помощью двух опорных точек для 8битного изображения приведен на рис. 1.19. Промежуточные значения функции преобразования получаются путем линейной интерполяции с использо-
56
21
формат 4:2:0;
а)
Входное
значение
Градационное
преобразование
чем для яркостной компоненты Y . Информация о размерах компонент и их
соотношении сохраняется в заголовке в виде следующих переменных:
Выходное
значение
X = W, Y = H;
H0 = 2, V0 = 2;
б)
H1 = 1, V1 = 1;
250
H2 = 1, V2 = 1.
200
Здесь W — ширина, а H — высота исходного изображения (в пикселях).
Нулевой индекс соответствует компоненте Y , первый — Cb, а второй — Cr.
150
100
1.3.
Обработка данных в базовом режиме
50
0
0
50
100
150
200
Принципиальная схема сжатия в базовом режиме стандарта JPEG показана на рис. 3.1.
250
Рис. 1.18. Градационное преобразование: а — “черный ящик“; б — пример функции преоб-
Cs1,
, CsN
Разбиение на
блоки 8*8
разования
Кодирование без потерь
ДКП
(DCT)
ванием двух фиксированных точек (0, 0), (255, 255) и промежуточных точек,
задаваемых вручную.
Матрицы
квантования
Cs 1,
1.5.3 Гамма-преобразование
Квантование
(Q)
, Cs N
Зигзагообр.
обход
Кодирование
длин серий
(RLE)
Кодирование по
Хаффману
Таблицы кода
Хаффмана
Битовый
поток
Формирование
компоненты
Декодирование без потерь
Обратное
ДКП
Степенные преобразования выполняются по следующей формуле:
′
= c(Iy,x )γ ,
Iy,x
(1.29)
(IDCT)
Деквантование
(Q-1)
Формирование
блока 8x8
Декодирование по
Хаффману
Рис. 3.1. Сжатие в базовом режиме стандарта JPEG
где c и γ являются положительными величинами. На практике процедура получила широкое распространение благодаря гамма-коррекции, выполняемой
в системах отображения с помощью электронно-лучевых трубок, поскольку
гамма-коррекция имеет аналогичную передаточную функцию. При этом перед выполнением расчета по формуле (1.29) входное значение Iy,x следует
′
привести в диапазон [0, 1]. Выходное значение Iy,x
возвращают в диапазон
[0, 255] путем умножения на 255.
Очевидно, что кривые, полученные со значениями γ > 1 и γ < 1,
имеют различное поведение. Кривые степенных зависимостей при малых γ
Каждая компонента разбивается на непересекающиеся блоки размером
8×8 пикселей, которые в стандарте также называют элементом данных (data
unit). Если какой-либо из размеров W или H не кратен 8, то изображение
расширяется до необходимых размеров в соответствующем направлении. В
стандарте вводится понятие минимального элемента данных для кодирования (minimal coded unit — MCU). Его конфигурация зависит от того, какой из
двух режимов обработки выбран.
22
55
1. Возможность независимой обработки компонент.
2. Независимая обработка блоков (за исключением коэффициентов постоянного тока) внутри одной компоненты.
3. Сжатие изображений с потерями.
1.2. Требования к входным данным
Стандарт предполагает использование до 4 входных компонент с 8битными (для базового режима) или 12-битными (для расширений и для прогрессивного режима) значениями. Режим кодирования без потерь позволяет
использовать входные данные в интервале от 2 до 16 бит. В иерархическом
режиме входные данные могут иметь знак.
Стандарт позволяет использовать децимированные данные. Для определения размеров каждой компоненты в заголовке сжатого файла должны
присутствовать следующие параметры:
Рис. 1.19.
Пример функции градационного преобразования на базе двух опорных точек
(65, 40) и (180, 210)
– Y — максимальная высота изображения среди всех компонент,
– X — максимальная ширина изображения среди всех компонент,
– Vi — коэффициент децимации по высоте i-й компоненты,
– Hi — коэффициент децимации по ширине i-й компоненты.
Размеры компонент вычисляются по следующим формулам:
⌈
⌉
Hi
xi = X ×
;
Hmax
⌈
⌉
Vi
yi = Y ×
,
Vmax
где Vmax и Hmax — максимальные значения среди Vi и Hi .
Манипулировать исходными размерами компонент можно с помощью
алгоритмов пространственной интерполяции. Также перед выполнением
процедуры обработки к исходному набору компонент может быть применено декоррелирующее преобразование. Примером предварительной подготовки исходных данных может быть преобразование из формата RGB в
YCbCr с последующим прореживанием цветоразностных компонент. В данном разделе будем предполагать, что алгоритм сжатия в базовом режиме
стандарта JPEG обрабатывает данные в формате Y CbCr 4:2:0. Прореживание цветоразностных компонент выполняется таким образом, чтобы размеры ширины и высоты цветоразностных Cb и Cr были в два раза меньше,
Рис. 1.20. Графики степенных преобразований для γ = 0.1, 0.6, 1.5 и 6.0 соответственно
(параметр c = 1)
54
23
отображают узкий диапазон малых входных значений в широкий диапазон
выходных значений, однако для больших входных значений верно обратное
утверждение. Следует также отметить, что выбор γ = 1 превращает выражение (1.29) в тождественное преобразование.
1.5.4 Метод выравнивания гистограмм
Видоизменение (выравнивание) гистограммы может быть использовано в задаче улучшения качества изображений. Данный подход может увеличить контрастность обрабатываемого изображения. Использование гистограмм для обработки данных в режиме реального времени является широко
используемым, так как их вычисление является достаточно простым как для
программной, так и для аппаратной реализации.
Гистограммой интенсивности цифрового изображения с уровнями яркости в диапазоне [0, 2L − 1] называется дискретная функция h(rk ) = nk ,
где rk — k-я градация яркости, а nk — число пикселей на изображении, имеющих яркость rk . Как правило, производится нормирование гистограммы
путем деления каждого из ее значений на общее число пикселей в изображении и формируется оценка вероятности появления на изображении пикселя
с интенсивностью rk :
p(rk ) = nk /N,
(1.30)
где N — общее число пикселей на изображении. Сумма всех значений нормированной гистограммы равна единице. Для выполнения операции выравнивания гистограмм следует воспользоваться формулой на базе эмпирической функции распределения:
III. МЕТОДЫ БЛОКОВОЙ
ОБРАБОТКИ ПРИ СЖАТИИ С
ПОТЕРЯМИ НА ПРИМЕРЕ
СТАНДАРТА JPEG
Основная задача данного раздела учебного пособия — изучение алгоритмов, используемых в базовом (baseline) режиме стандарта JPEG, анализ
статистических свойств, используемых при сжатии коэффициентов дискретного косинусного преобразования (ДКП), а также получение практических
навыков разработки методов блоковой обработки при сжатии изображений
с потерями.
1. Теоретические основы
1.1.
Характеристика стандарта JPEG
где 2 − 1 — максимально допустимое число градаций интенсивности на
изображении (255 для 8-битного изображения); sk — значение интенсивности в выходном изображении, в которое будет преобразовываться rk .
На основе формулы (1.31) строится таблица поиска (Lookup-table). В
этой таблице каждому входному числу (интенсивности пикселя) поставлено
в соответствие выходное. В дальнейшем сформированная таблица используется в качестве градационного преобразования. В результате применения метода выравнивания гистограмм полученная функция распределения должна
приобрести характер близкий к функции распределения равномерного закона.
В методе выравнивания гистограмм следует выделить три основных этапа:
– аппроксимация распределения вероятностей градаций интенсивности
пикселей изображения с помощью гистограммы;
Стандарт JPEG получил название по аббревиатуре группы Joint
Photographic Expert Group, которая образовывала совместный комитет от двух международных организаций по стандартизации — CCITT
(Consultative Committee for International Telephony and Telegraphy) и
ISO/IEC (International Organization for Standardization / International
Electrotechnical Commission). CCITT в 1992 г. была преобразована в ITU-T(
International Telecommunications Union - Telecommunication Standardization
Sector). Стандарт именуется этими организациями соответственно ITU-T
Recommendation T.81 (с 1992 г.) и ISO/IEC 10918-1 (с 1994г.).
Стандарт JPEG предполагает наличие четырех основных режимов кодирования:
1. Базового (baseline или sequential).
2. Прогрессивного (progressive).
3. Иерархического (hierarchical).
4. Кодирования без потерь (lossless).
Ключевыми особенностями обработки многокомпонентных (цветных)
изображений в базовом режиме JPEG являются:
24
53
sk = (2L − 1)
k
∑
j=0
2L − 1 ∑
nj , k = 0, 1, 2, ..., 2L − 1
N j=0
k
pr (rj ) =
L
(1.31)
b) ∃i : x′i < xi ;
c) ∀i :
x′i
x′i
– переход к эмпирической функции распределения на основе анализа построенной гистограммы;
– использование эмпирической функции распределения в качестве функции
преобразования.
К недостатку данного метода можно отнести нечувствительность к локальным особенностям изображения, что в ряде случаев может негативно
сказаться на качестве обработанного изображения.
< xi ;
d) ∀i : > xi ;
n
n
∑
∑
e)
x′i >
xi ;
f)
g)
i=1
n
∑
i=1
n
∑
|x′i |
i=1
n
∑
>
(x′i )2 >
i=1
|xi |;
i=1
n
∑
(xi )2 .
а)
i=1
б)
4. Является ли целевая функция в задаче подбора оптимальных параметров
кода Голомба унимодальной?
5. Какой контур детектируется предсказанием MED для изложенных далее ситуаций (ответ сформулировать на качественном уровне в терминах ”вертикальный/горизонтальный перепад интенсивности с высокой/низкой на низкую/высокую в окрестности предсказываемого пикселя”):
a) C ⩽ min(A, B) и min(A, B) = A;
b) C ⩽ min(A, B) и min(A, B) = B;
c) C ⩾ max(A, B) и max(A, B) = A;
d) C ⩾ max(A, B) и max(A, B) = B.
6. Рассматривается схема сжатия, основанная на дифференциальной импульсной кодовой модуляции с MED предсказанием и равномерным скалярным квантованием с шагом Q = 1. Допустим, сжимаются 8-битные
данные. Будем считать, что данные на входе можно рассматривать как сообщения некоторого источника без памяти. Сравните максимальное значение энтропии квантованных ошибок предсказания и входных данных.
Объясните, за счет чего достигается сжатие в подобной схеме.
Рис. 1.21.
Результат применения метода выравнивания гистограмм: а — исходное (затемненное) изображение; б — изображение после выравнивания
Результат применения метода выравнивания гистограмм к затемненному
изображению показан на рис. 1.21 и 1.22.
1.5.5 Построение бинарного изображения с помощью градационного
преобразования
Рассмотрим процесс построения бинарного изображения, основанный
на использовании градационного преобразования, функция преобразования
которого показана на рис. 1.23.
В зависимости от выбора порога T выходной пиксель бинарного значения соответствует белому (интенсивность соответствующего пикселя исходного изображения превышает порог) или черному цвету (иначе).
Примеры соответствующих бинарных изображений показаны на
рис. 1.24 и 1.25.
52
25
а)
б)
1.2 Методы скалярного квантования
1.2.1. равномерное квантование типа mid-rise.
1.2.2. равномерное квантование типа mid-tread.
1.2.3. квантование с расширенной нулевой зоной (как пример неравномерного квантования).
1.3 Кодирование и декодирование на основе кода Голомба (параметры
кода должны задаваться как параметр).
2. Проанализировать зависимость ошибки предсказания от порядка авторегрессионной модели n.
Рис. 1.22.
Результаты применения метода выравнивания гистограмм: а — исходная гистограмма; б — гистограмма после выравнивания
3. Сравнить результаты, полученные для схемы с предсказанием на основе двумерной авторегрессионной модели, с результатами, полученными с
использованием предсказания MED.
4. Сравнить эффективность реализованных способов квантования, зафиксировав способ предсказания. Требуемый способ предсказания определяет
преподаватель.
5. Меняя шаг квантования ∆, построить зависимость PSNR между исходным и декодированным изображениями от размера битового потока на
выходе кодера.
6. Подобрать оптимальные параметры кода Голомба для каждого шага квантования ∆ в схеме с зафиксированными способами предсказания и квантования. Требуемые способы определяет преподаватель.
7. Оценить эффективность применения кодов Голомба, сравнив усредненный по количеству кодируемых символов размер полученного битового
потока с оценкой энтропии, вычисленной по формуле (2.27).
3. Вопросы для самопроверки
Рис. 1.23.
Пример функции преобразования для построения бинарного изображения при
1. Как изменится ошибка предсказания, если в DPCM кодеке убрать обратную связь?
T = 127
2. Какова сложность алгоритма Левинсона-Дарбина?
2. Практические задания
В процессе выполнения заданий необходимо работать с яркостной компонентой Y изображения, которую можно выделить, воспользовавшись формулой (1.1). Для визуальной оценки выходное изображение следует формировать в градациях серого цвета, распространяя значения яркостной состав-
26
3. Пусть для одних и тех же входных данных получены два набора ошибок
предсказания: x = {x1 , x2 , ..., xn } — ошибки на выходе MED предсказания; x′ = {x′1 , x′2 , ..., x′n } — ошибки на выходе авторегрессионной модели. Укажите, какие из следующих утверждений верны, объясните на качественном уровне почему:
a) ∃i : x′i > xi ;
51
{
e′ =
а)
2e, если e ≥ 0
2|e| − 1, иначе
.
б)
(2.26)
Данное правило выполняет отображение отрицательных чисел на нечетные, а положительных чисел — на четные положительные числа.
1.7.
Оценка эффективности сжатия
Эффективность сжатия обычно оценивается с применением формулы
энтропии, как среднего числа бит, затрачиваемых на представление одного
символа кодируемого источника. Наиболее простым способом анализа эффективности сжатия является оценка энтропии при поэлементном независимом сжатии. В этом случае кодируемые значения интерпретируются как сообщения некоторого постоянного источника X (стационарного источника,
на выходе которого сообщения независимы). Алфавит источника X формируется по всем возможным кодируемым значениям . Величина p̂(x) является
оценкой вероятности сообщения x, которая вычисляется:
p̂(x) =
nx
,
n
где nx — количество пикселей в данной компоненте, которые имеют значение x, а n — общее число пикселей в компоненте. Тогда оценка энтропии
при поэлементном независимом сжатии вычисляется по следующей формуле [17]:
∑
Ĥ(X) = −
p̂(x) log2 p̂(x).
(2.27)
x
2. Практические задания
1. Написать программу, реализующую следующие методы схемы дифференциальной импульсной кодовой модуляции:
1.1 Предсказание
1.1.1. Двухпроходное предсказание на основе двумерной авторегрессионной модели порядка n. Порядок модели и конфигурацию соседних пикселей, используемых для предсказания,
определяет преподаватель.
1.1.2. Схема MED из стандарта JPEG-LS
50
Рис. 1.24.
Построения бинарного изображения с помощью градационного преобразования на примере изображения ’lena’ при T = 127: а — исходное изображение; б — выходное
бинарное изображение
ляющей пикселя на соответствующие позиции компонент R, G и B.
2.1. Реализация методов шумоподавления
1. Реализовать модель аддитивного шума в соответствии с формулой 1.3, используя Гауссовскую случайную величину с распределением N (0, σ2 ).
2. Реализовать модель импульсного шума в соответсвии с формулой (1.4).
3. Построить графики PSNR(σ) для модели аддитивного шума и
PSNR(pa , pb ) для модели импульсного шума. Значение P SN R вычислять
по оригинальной и зашумленным компонентам Y . Преподаватель задает
диапазон для параметров σ, pa и pb .
4. Обработка изображений с Гауссовским шумом.
4.1 Реализовать метод скользящего среднего.
4.2 Подобрать размер окна R для метода скользящего среднего, который бы позволил максимизировать значение PSNR, вычисляемое по
оригинальному изображению и зашумленному изображению после
фильтрации.
4.3 Реализовать метод Гауссовской фильтрации.
27
а)
б)
Приведем алгоритм расчета кода Голомба с параметром m для целого
неотрицательного числа n.
Кодирование:
1. Расчет параметра c:
c = ⌈log2 m⌉.
(2.23)
2. Расчет целой части q и остатка r:
⌊n⌋
q=
, r = n − qm.
m
(2.24)
3. Кодирование целой части q унарным кодом: унарный код состоит из q единиц и завершающего нуля.
Рис. 1.25. Построения бинарного изображения с помощью градационного преобразования
на примере изображения ’peppers’ при T = 127: а — исходное изображение; б — выходное
бинарное изображение
4.4 Подобрать параметр σ для метода Гауссовской фильтрации, который бы позволил максимизировать значение PSNR, вычисляемое по
оригинальному изображению и зашумленному изображению после
фильтрации.
4.5 Построить графики PSNR(σ) для нескольких фиксированных значений R, определяемых преподавателем.
4.6 Сформировать выходные изображения для наилучших с точки зрения PSNR вариантов использования фильтра Гаусса и скользящего
среднего. Провести визуальное сравнение. Указать искажения, специфические для рассмотренных методов.
4.7 Реализовать метод медианной фильтрации с размерами окон (2R +
1) × (2R + 1), где R определяет радиус фильтра.
4.8 Подобрать размер окна R для метода медианной фильтрации, который бы позволил максимизировать значение PSNR, вычисляемое по
оригинальному изображению и зашумленному изображению после
фильтрации.
4.9 Проанализировать значения PSNR, рассчитанные для зашумленных
изображений после фильтрации.
4.10 Произвести анализ полученных значений PSNR после применения
различных фильтров.
28
4. Кодирование остатка r равномерным кодом. При этом если r < 2c − m,
то равномерный код занимает c − 1 бит, иначе код занимает c бит.
5. Кодовое слово для числа n получается конкатенацией унарного кода целой
части q и равномерного кода для остатка r.
Следует отметить, что код r строится таким образом, чтобы максимальный остаток имел код, состоящий из всех единиц. Приведем простое правило, позволяющее строить такие коды. Обозначим наибольшее значение, кодируемое c − 1 битом, как b. Код для первого c-битного значения получается
конкатенацией двоичного представления числа b + 1 с нулевым битом.
Декодирование:
1. Сначала декодируется значение q как число единиц до первого нуля.
2. Декодирование n
{
n=
n = m ∗ q + R, если R < 2c − m
n = m ∗ q + R′ − (2c − m), иначе
,
(2.25)
где R и R′ — значения, хранящиеся в c − 1 и c битах после завершающего
нуля унарного кода соответственно.
Следует отметить, что в том случае, когда параметр кода Голомба m является степенью двойки, коды для всех остатков занимают одинаковое число
бит c. Процедуры кодирования и декодирования в таком случае упрощаются.
Таже отметим, что так как ошибки предсказания могут быть отрицательными, перед кодированием кодом Голомба их надо отобразить на множество
неотрицательных чисел с использованием следующего правила:
49
ношению сигнал/шум, (SNR — signal-to-noise ratio):
H ∑
W
∑
2
(Iy,x )
ES
y=1 x=1
SNR = 10 lg
= 10 lg H W
,
∑ ∑
EN
2
(Iy,x − Iˆy,x )
(2.21)
y=1 x=1
где ES определяет энергию исходного сигнала; EN — энергию шума. Энергия дискретного по времени сигнала определяется как сумма квадратов всех
его значений.
И для аудиоданных, и для изображений шум интерпретируется как несоответствие исходного и обработанного (восстановленного) сигнала, а энергия шума — как сумма квадратов разностей отсчетов исходного и обработанного сигналов, аналогично формуле (2.18). Отличие PSNR от SNR заключается в том, что все значения исходного сигнала в числителе дроби заменяются на максимальное значение интенсивности пикселя в данной компоненте. Однако для упрощения вычислений максимальное значение сигнала, как
правило, заменяют на максимально возможное значение, которое может принимать интенсивность пикселя. Таким образом, в используемой разрядной
сетке L бит значение PSNR вычисляется по следующей формуле:
PSNR = 10 lg
W H(2L − 1)2
H ∑
W
∑
.
(2.22)
2
(Iy,x − Iˆy,x )
y=1 x=1
Более подробно о критериях оценки искажений можно узнать в работах [8, 13].
1.6. Кодирование (сжатие) ошибок предсказания на основе кодов
Голомба
Кодами Голомба называют семейство кодов для сжатия информации без
потерь, предложенное в 1960-х г. Соломоном Голомбом [14]. Данные коды
можно рассматривать как эффективные коды Хаффмана, оптимальные для
источников, описываемых геометрическим распределением. Также одним из
существенных преимуществ кодов Голомба является то, что при кодировании и декодировании нет необходимости хранить таблицу кодовых слов.
Различные модификации кодов Голомба используются в современных
стандартах сжатия изображений и видеоданных, в частности, в стандартах
JPEG-LS [15] и H.264 [16].
48
5. Обработка изображений с импульсным шумом.
5.1 Наложить на исходное изображение импульсный шум в соответствии с моделью (1.4). Подготовить несколько вариантов зашумленного изображения, для которых доля искаженных пикселей составит 5%, 10%, 25%, 50%.
5.2 Вычислить значения PSNR по исходному и зашумленным изображениям.
5.3 Реализовать метод медианной фильтрации с размерами окон (2R +
1) × (2R + 1), где R определяет радиус фильтра.
5.4 Построить зависимость PSNR(R) для всех вариантов зашумленных
изображений. Провести визуальное сравнение. Указать искажения,
специфические для рассмотренных методов. Выбрать наилучший
параметр фильтра для каждого варианта зашумленного изображения.
2.2. Реализация методов выделения контуров
1. Реализация оператора Лапласа
1.1 Применить оператор Лапласа I ′ = L(I) к изображению, воспользовавшись формулой (1.5). В качестве весов ωi,j использовать элементы матрицы из рис. 1.8.
1.2 Сформировать изображение по отклику оператора Лапласа I ′ . Значение пикселя на позиции (y, x) вычислить по следующей формуле:
′
Clip(Iy,x
+ 128, 0, 255)
(1.32)
1.3 Синтезировать изображение с усилением высоких частот, используя
разность исходного изображения I и отклика оператора Лапласа I ′ :
′′
′
Iy,x
= Clip(Iy,x + Iy,x
, 0, 255),
(1.33)
применив маску (рис. 1.12).
1.4 Синтезировать изображения, применив маску, приведенной на
рис. 1.12 для разных значений параметра α. Параметр α следует варьировать в интервале [1, 1.5] с шагом 0.1. Сделать вывод о влиянии
параметра α на выходное изображение.
1.5 Определить значение параметра α, при котором синтезированное
изображение совпадет с I ′′ . Обосновать определенное значение α.
29
1.6 Рассчитать среднее значение яркости для каждого полученного
изображения. Сделать выводы о причинах изменения этого значения.
1.7 Построить и проанализировать гистограммы для исходного изображения, а также изображений, полученных в результате выполнения
пункта 14.
2. Реализация оператора Собеля
2.1 Применить оператор Собеля к исходному изображению I.
2.2 Рассчитать значения |∇Iy,x | и θy,x , воспользовавшись формулами
(1.25) и (1.26).
2.3 Подобрать порог thr для формирования наиболее точной бинарной
карты градиентов на изображении по следующему правилу: если
значение |∇Iy,x | для рассматриваемого пикселя превысило порог
thr, то установить значение его интенсивности в 0, иначе — 255.
2.4 Сформировать файлы для рассмотренных порогов thr и выбрать
наилучший. Сделать выводы о влиянии порога thr на бинарную
карту градиентов.
2.5 Построить 4-цветную карту направлений градиентов на полученных изображениях. Каждому из четырех квадрантов в соответствие
ставится определенный цвет: красный, зеленый, синий или белый.
Пиксель карты направлений окрашивается в цвет, соответствующий квадранту, в котором располагается вектор, характеризующий
направление контура в рассматриваемом пикселе.
2.3.
Реализация градационных преобразований
а Iˆy,x — интенсивность пикселя восстановленного изображения на той же
позиции. Размеры изображения будем полагать равными W по ширине и H
по высоте. Тогда шумом на позиции (y, x) можно считать разность между
исходным и восстановленным пикселем:
Ny,x = Iy,x − Iˆy,x .
Наиболее простыми с вычислительной точки зрения критериями оценки искажений при обработке изображений с потерями являются сумма абсолютных значений разностей между исходными и восстановленными после обработки пикселями (Sum of Absolute Differences — SAD) и сумма квадратов разностей между исходными и восстановленными пикселями (Sum of
Squared Differences — SSD):
SAD =
2. Применить к ним градационное преобразование на базе двух опорных точек для повышения качества исходных изображений.
|Iy,x − Iˆy,x |;
(2.17)
2
(Iy,x − Iˆy,x ) .
(2.18)
y=1 x=1
SSD =
H ∑
W
∑
y=1 x=1
Вместо слова разность (Difference) иногда используют слово ошибка (Error).
Данные критерии часто используются для поиска и сравнения отдельных фрагментов изображения. Для сравнения изображений разной размерности или фрагментов изображений используют усредненные значения SAD
(Mean of Absolute Errors — MAE) и SSD (Mean of Squared Errors — MSE):
H W
1 ∑∑
|Iy,x − Iˆy,x |,
W H y=1 x=1
(2.19)
H W
2
1 ∑∑
(Iy,x − Iˆy,x ) ,
W H y=1 x=1
(2.20)
MAE =
2.3.1 Градационные преобразования на базе опорных точек
1. Синтезировать засвеченное, затемненное, а также сбалансированное
изображения.
H ∑
W
∑
MSE =
3. Сформировать гистограммы для исходных и полученных изображений.
Сделать выводы о влиянии алгоритма на исходные изображения.
где W является шириной, а H — высотой фрагмента или изображения.
Более распространенным критерием оценки искажений, вносимых в
изображение, является пиковое отношение сигнал/шум (PSNR — peak signalto-noise ratio). Эта мера аналогична используемому при обработке звука от-
30
47
а)
б)
2.3.2 Гамма преобразования
′
1. Получить графики уравнения Iy,x
= c(Iy,x )γ для различных значений
параметра γ для каждого из трех входных изображений. Значение c следует положить равным 1, а параметр γ следует перебирать в интервале
[0.04, 25].
2. Сделать выводы о влиянии параметра γ на исходные изображения.
3. Сформировать гистограммы для исходных и полученных изображений.
Сделать выводы о влиянии алгоритмов на исходные изображения.
2.3.3 Методы выравнивания гистограмм
Рис. 2.5. Примеры равномерного и неравномерного квантования скалярного квантования:
а — равномерное скалярное квантование; б — квантование с расширенной нулевой зоной (квант
нулевой зоны больше, чем остальные)
1. Синтезировать засвеченное, затемненное, а также сбалансированное
изображения.
2. Применить к ним алгоритм выравнивания гистограмм в соответствии с
формулой (1.31).
{
sign(x) =
1, если x ≥ 0,
−1, иначе
.
Скалярное квантование вещественных чисел с шагом ∆ и расширенной
нулевой зоной можно определить с помощью формулы:
⌊
⌋
|x| + f (∆)
qext_dz =
∗ sign(x),
(2.15)
∆
где функция f (∆) контроля ширины нулевой зоны определяется как f (∆) =
a
a
1
a
b ∆. Для отношения b должно выполняться условие b < 2 . Вычисление аппроксимирующего значения для обоих случаев выполняется по следующей
формуле:
y = qx ∆.
(2.16)
1.5. Критерии оценки искажений
3. Сформировать гистограммы для исходных и полученных изображений.
Сделать выводы о влиянии алгоритма на исходные изображения.
2.3.4 Методы построения карты контуров на основе градационных
преобразований
1. Синтезировать градационную функцию преобразования для построения
бинарного изображения. Порог T следует перебирать в интервале [16, 240]
с фиксированным шагом.
2. Подобрать порог T для формирования наиболее наглядного бинарного
изображения по следующему правилу: если значение s для рассматриваемого пикселя превысило порог T , то установить значение его интенсивности в 0, иначе — 255.
3. Сформировать файлы для рассмотренных порогов T и выбрать наилучший. Сделать выводы о влиянии порога T на полученные бинарные изображения.
Для оценки искажений, вносимых при обработке с потерями, часто используют понятие шум, который как бы был внесен в процессе применения
цепочки «преобразование — обратное преобразование». Для определения
количественных характеристик шума обозначим Iy,x — интенсивность пикселя исходного изображения на позиции (y, x) (в строке y и в столбце x ),
1. Какой размер имеет ядро двумерного фильтра при использовании апертуры радиуса R?
46
31
3. Вопросы для самопроверки
2. При выполнении фильтрации с подъемом высоких частот средняя яркость
изображения увеличивается. Как это можно объяснить?
3. Опишите ключевые особенности функции передачи для градационного
преобразования, которое увеличивает яркость темных областей изображения.
4. Для алгоритма построения бинарного изображения из подразд. 1.5.5 предложите процедуру автоматического выбора порога T .
зывают mid-tread. При этом ноль, как правило, является аппроксимирующим значением. По этой причине такой тип равномерного квантования
также называют квантованием с нулевой зоной (dead-zone). Примером такого вида скалярного квантования является арифметическое округление.
а)
б)
Рис. 2.4. Виды скалярного квантования: а — mid-rise – ноль на границе кванта; б — midtread – ноль внутри кванта
Как уже было сказано ранее, в случае равномерного скалярного квантования (рис. 2.5,а) интервалы всех квантов равны между собой и равны значению ∆, который является параметром алгоритма. При неравномерном скалярном квантовании хотя бы один квант должен отличаться по длине от всех
остальных. Примером такого вида квантования является квантование с расширенной нулевой зоной (рис. 2.5,б), для которого длина кванта, в который
попадает ноль, больше, чем размер всех остальных квантов, равных между
собой.
Равномерное скалярное квантование вещественных чисел с шагом ∆
можно описать следующей формулой:
⌊
⌋
|x| + 12 ∆
quni =
∗ sign(x),
(2.14)
∆
где с помощью операции ⌊·⌋ обозначается деление с отбрасыванием дробной части (целочисленное деление), функция sign определена следующим
образом:
32
45
II. МЕТОДЫ ПИКСЕЛЬНОЙ
ОБРАБОТКИ ПРИ СЖАТИИ
ИЗОБРАЖЕНИЙ
a)
xMIN
t0
y1
y2
t1
y3
t2
0
y4
t3
xMAX
t4
X
б)
Изучение алгоритмов, используемых при пиксельной обработке изображений на примерах разностных схем и стандарта сжатия JPEG-LS, а также
получение практических навыков разработки методов пиксельной обработки при сжатии изображений. Для упрощения анализа в данном разделе рассматриваются только однокомпонентные изображения.
y4
xMAX
xMIN
t0
t1
t2
y3
t3
t4
X
1. Теоретические основы
y2
1.1. Схема дифференциальной импульсной кодовой модуляции для
сжатия изображений
y1
Сжатие с использованием пиксельной обработки, как правило, базируется на алгоритме дифференциальной импульсной кодовой модуляции
(DPCM, Difference Pulse Coding Modulation). Типовая схема этого алгоритма показана на рис. 2.1.
Рис. 2.3. Пример скалярного квантования N = 4: а — пример разбиения числовой оси; б —
ступенчатая функция
Здесь на оси аргумента отложены элементы исходного множества X, а
по оси ординат — значения аппроксимирующего множества Y . Если множество X является неограниченным, то t0 = −∞, а tN = ∞. Для ограниченного множества X границы крайних квантов и определяются предельными
элементами множества X.
В зависимости от положения нулевого значения различают два вида скалярного квантования:
Сжатие
x
Восстановление
e  x  xˆ
i, j
i, j
+
i, j
e'
i, j
x'
i, j
Кодер
Q
Декодер
Q-1
i, j
+
__
Q-1
xˆ
i, j
e'
i, j
xˆ
1. Нулевое значение расположено на границе одного из квантов (рис. 2.4,а).
В англоязычной литературе этот вид квантования называют mid-rise. При
этом ноль не является аппроксимирующим значением.
+
Предсказание
i, j
Предсказание
x'
i, j
2. Нулевое значение расположено внутри одного из квантов, часто в его середине (рис. 2.4,а). В англоязычной литературе этот вид квантования на-
Рис. 2.1. Принципиальная схема кодирования на базе алгоритма импульсной кодовой моду-
44
33
ляции
Один из ключевых блоков этой схемы — предсказание, целью которого является декорреляция исходных данных, которые для изображений обладают двумерной пространственной избыточностью. Результат предсказания
x̂i,j для текущего пикселя вычитается из значения этого пикселя xi,j , и полученная таким образом ошибка предсказания ei,j подвергается дальнейшей
обработке.
Если при сжатии допустимы потери информации, то значения ошибок
предсказания подвергаются квантованию. В данном разделе будут рассмотрены методы скалярного квантования.
Окончательное формирование битового потока осуществляется путем
кодирования ошибок предсказания (для случая сжатия без потерь) или номеров квантов (для случая сжатия с потерями) с помощью статистических
методов кодирования источников. Кодирование осуществляется в блоке кодер, а декодирование при восстановлении в блоке декодер (рис. 2.1). Примерами методов кодирования источников являются коды Голомба, Хаффмана
или алгоритм арифметического кодирования. В данном разделе будут рассмотрены методы сжатия на базе кодов Голомба.
В случае применения в схеме DPCM квантования предсказание должно осуществляться с использованием восстановленных данных. Для этого в
схеме предусмотрена обратная связь, в которой выполняется восстановление
данных, в процессе которого осуществляется деквантование (определение
аппроксимирующего значения кванта) квантованной ошибки предсказания
и сложение аппроксимированной ошибки предсказания с предсказанием,
выполненным для данного пикселя. Аппроксимированное значение ошибки квантования на рис. 2.1 отмечено как e′i,j , а результат восстановления
пикселя — x′i,j . В случае отсутствия процедур квантования-деквантования
x′i,j = xi,j и обратная связь на стороне кодера не нужна, так как кодер может
оперировать входными данными.
Методы предсказания, анализируемые в данном разделе, будут рассмотрены в подразд. 1.2. Используемые методы квантования рассмотрены в подразд. 1.4. Сжатие на основе кодов Голомба представлено в подразд. 1.6.
1.2. Методы предсказания данных
3. Векторное.
Квантование называют скалярным, если элементами множеств X и Y являются числовые (скалярные) величины. Если элементы обоих множеств являются векторами, то квантование называют векторным. В этом случае при
отображении по формуле (2.13) длины всех векторов x и y должны совпадать. В данном подразделе рассматривается только скалярное квантование.
1.4.2 Скалярное квантование
Поскольку элементами множества X и Y для скалярного квантования
являются числа, для простоты будем далее полагать, что множество X есть
подмножество вещественных чисел. Тогда кванты будут представлять собой
интервалы вещественной оси, а определение процедуры квантования можно интерпретировать как задание способа разбиения числовой оси на эти
интервалы и выбора в каждом из них единственного аппроксимирующего
значения. Таким образом, для того чтобы определить процедуру скалярного
квантования, необходимо задать следующие величины:
1. Число квантов N . Кванты при этом нумеруются от 1 до N . Для упрощения
кодирования номера кванта с помощью равномерного кода, число квантов
выбирают как N = 2R , где R определяет число двоичных разрядов равномерного кода.
2. Границы квантов ti , i = 0, ..., N ; j-й квант определяет полуоткрытый интервал [tj−1 ; tj ). В общем случае неважно, находится закрытая граница
’[’ справа или слева. В данной работе будем полагать, что такая граница
всегда будет располагаться слева. Исключение могут составлять первый
и последний кванты. Если множество X не ограничено слева, то первый
квант (−∞, t1 ) будет открытым. Если множество X ограничено справа
значением tN , то последний квант [tN −1 , tN ] будет закрытым.
3. Аппроксимирующие значения yi ∈ [t(i−1) ; ti ), i = 1, ..., N .
Сама процедура скалярного квантования сводится к замене вещественного
числа x на аппроксимирующее значение кванта, в который оно попадает:
x ∈ [t(i−1) , ti ) 7−→ Q(x) = yi .
Для простоты изложения в этом подразделе будет рассмотрен одномерный случай, который может быть легко обобщен на двумерный. Для последо-
Заметьте, что yi также являются вещественными числами. Процедуру отображения при скалярном квантовании часто изображают графически в виде
ступенчатой функции в декартовой системе координат, как это показано на
рис. 2.3,textitб.
34
43
1.2.1 Постановка задачи предсказания
Обработка в
сканирующем
порядке
вательности значений x0 , ..., xt−M , xt−M +1 , ..., xt−1 требуется определить
способ, который бы позволил найти значение оценки x
bt значения следующего элемента последовательности как функции от M предыдущих значений:
C
B
D
A
X
Y
H
Рис. 2.2. Алгоритм нелинейного предсказания MED из стандарта JPEG-LS
и снизу. Таким образом, если C ⩾ max(A, B) или C ⩽ min(A, B), уместно
предположить, что в окрестности текущего обрабатываемого пикселя проходит контур.
1.4.
Квантование данных
1.4.1 Общие сведения
x
bt = F (xt−M , xt−M +1 , ..., xt−1 ) .
(2.1)
Такой способ будем называть предсказанием. Алгоритмы предсказания можно разделить на два основных класса.
1. Линейные методы используют для предсказания линейную комбинацию
значений ранее обработанных пикселей. Примером линейного предсказания, который будет рассмотрен в данном разделе, является метод статистического анализа данных на основе авторегрессионных моделей.
2. Нелинейные методы используют для предсказания нелинейные операции
или правила. В качестве примера нелинейного предсказания в разделе будет рассмотрен метод Median Edge Detector (MED), включенный в стандарт JPEG-LS.
1.2.2 Оценка эффективности предсказания
В широком смысле под квантованием понимают правило отображения
элементов x из некоторого исходного множества X в элементы другого множества Y :
y = Q(x),
(2.13)
причем процедура отображения такова, что несколько элементов из X (для
непрерывного множества — это подмножество значений) могут отображаться в один и тот же элемент из Y . Подмножество из X, все элементы которого
отображаются в один и тот же элемент из Y , образуют квант. Множество
Y принято называть аппроксимирующим множеством. С точки зрения теории множеств, множество X мощнее Y . Как правило, множество Y является
конечным или счетным.
Конечный результат процедуры квантования зависит от конкретной прикладной задачи, в которой она применяется. Это может быть аппроксимирующий элемент y или номер кванта, в который попадает x. Результат представляется в виде индекса или кода, удобного для дальнейшего хранения или
передачи в цифровой форме.
Методы квантования можно разделить на следующие основные классы.
1. Скалярное равномерное квантование.
2. Скалярное неравномерное квантование.
42
В качестве меры точности способа предсказания будем использовать
сумму квадратов ошибок предсказания, формируемых алгоритмом по всей
последовательности:
∑
∑
2
E=
e2t =
(xt − x
bt ) .
(2.2)
t
t
1.2.3 Поиск оптимальных коэффициентов авторегрессионной модели в
задаче линейного предсказания данных
Примером линейных методов предсказания является авторегрессионная
модель порядка M , для которой предсказание очередного значения x
bt вычисляется как линейная комбинация значений интенсивностей соседних элементов, которые были обработаны ранее:
x
bt =
M
∑
ai xt−i .
(2.3)
i=1
Наибольшую популярность метод линейного предсказания получил в
обработке звука [6], [7]. Применительно к изображениям частный случай
35
для M = 1 рассматривается в работе [8]. Линейное предсказание на базе
авторегрессионной модели первого порядка также используется в стандарте
JPEG [9] для сжатия коэффициентов постоянного потока (этот метод будет
рассмотрен в следующем разделе).
Для конечной выборки кодируемых элементов существует процедура,
которая позволяет определить оптимальные с точки зрения ошибки предсказания значения коэффициентов линейной комбинации. Процедура базируется на составлении и решении системы уравнений Юла-Уокера (иногда
их еще называют уравнениями Винера-Хопфа). Если последовательность xt
является значениями детерминированного процесса, то для определения оптимальных значений ai необходимо минимизировать значение суммы квадратов ошибок (рис. 2.2):
E=
∑
t
e2t
=
∑
(
xt −
t
M
∑
)2
ai xt−i
i=1
−→ min .
ai
(2.4)
Для этого используют метод наименьших квадратов, который сводит
поиск оптимальных значений ai к решению системы из M однородных уравнений. Каждое уравнение системы представляет собой частную производную по искомому параметру aj приравненную к нулю для нахождения точки
экстремума (в нашем случае это должен быть минимум):
...
(i+1)
aj
=
(i)
aj
(i)
+ λai−j+1 ,
...
(i+1)
ai
(i)
(i)
= ai + λa1 ,
(i+1)
ai+1 = λ,
i+1
∑
(i)
aj Ri−j+1 + λS (i) = 0.
j=0
Таким образом, из последнего тождества находится значение λ, с помо(i+1)
щью которого рассчитываются значения aj
и S (i+1) .
1.3.
Метод нелинейного предсказания из стандарта JPEG-LS
Решение полученной системы зависит от диапазона значений параметра
t. В случае, когда −∞ < t < ∞, применяется автокорреляционный метод.
Когда 0 < t < T − 1, применяется ковариационный метод. Рассмотрим
автокорреляционный метод.
Способ учета корреляции между соседями, реализованный в стандарте JPEG-LS, учитывает наличие контуров на изображении. Процедура получила название Median Edge Detector (MED). Входными данными алгоритма MED являются восстановленные значения интенсивностей соседних пикселей, которые уже были обработаны кодером и, соответственно, будут доступны декодеру, который выполняет аналогичный алгоритм. Расположение
соседних пикселей, чьи восстановленные значения будут равны A, B и C,
показаны на рис. 2.2. Описание алгоритма MED можно представить в виде
следующего псевдокода.
IF C ⩾ max(A, B)
THEN X̂ = min(A, B)
ELSE
IF C ⩽ min(A, B)
THEN X̂ = max(A, B)
ELSE X̂ = A + B − C
Как уже было сказано ранее, алгоритм MED направлен на анализ наличия контура на изображении в окрестности предсказываемого пикселя X.
Под контуром понимается перепад между значениями интенсивностей смежных пикселей, т.е один из пикселей имеет интенсивность, значительно превышающую интенсивность соседа. Перепад может быть вертикальным, если существенно различаются интенсивности соседа слева и справа, или горизонтальным, если существенно различаются интенсивности соседей сверху
36
41
∂E
= 0; j = 1, ..., M.
∂aj
Учитывая, что
[ (
)2 ]
M
∑
∑
∂
xt −
ai xt−i
M
∑
∑
∑
t
i=1
∂E
=
=
ai
xt−j xt−i −
xt xt−j ,
∂aj
∂aj
t
t
i=1
получаем систему уравнений следующего вида:
M
∑
i=1
ai
∑
t
xt−j xt−i =
∑
xt xt−j ;
j = 1, ..., M.
(2.5)
t
Следующим шагом алгоритма Левинсона-Дарбина является формирование
дополнительной системы на базе уравнений (2.10). В новой системе элементы обоих векторов-столбцов записаны в обратном порядке. Матрица R(i+1)
при этом остается неизменной.
 i+1


 
∑ (i)
0
a
R
R0
R1
...
Ri
Ri+1
i−j+1 
j

 R
 (i)  j=0

R0
... Ri−1
Ri 
 1
 ai  

0

  ...  

R1
... Ri−2 Ri−1    
 R2

...

  (i)  = 
 . (2.11)
 ...

...
...
...
...  a2  

  

0

 Ri
 
Ri−1 ... R0
R1  a(i)
1


0
Ri+1
Ri
... R1
R0
−1
(i)
S
Введем определение автокорреляционной функции дискретного сигнала
xt :
∞
∑
Ri ≜
∞
∑
xt xt−i =
t=−∞
xt xt+i ,
t=−∞
для которой справедливо свойство симметрии:
Ri = R−i .
Уравнения (2.5) можно переписать в виде:
M
∑
ai R|j−i| = Rj ;
j = 1, ..., M.
i=1
Поскольку матрица R(i+1) является симметричной и Теплицевой, обе системы (2.10) и (2.11) будут совпадать с точностью до положения коэффициентов.
Далее домножим обе части системы (2.11) на скаляр λ и сложим системы (2.10) и (2.11) друг с другом. В результате получим систему следующего
вида:



−1
R0
R1
...
Ri
Ri+1
(i) 
 (i)
 R
R0
... Ri−1
Ri 
 1
 a1 + λai 



...
R1
... Ri−2 Ri−1  

 R2
=

  (i)

 ...
...
...
...
...   ai−1 λa(i)
2 

  (i)

 Ri
Ri−1 ... R0
R1  ai + λa(i)
1
Ri+1
Ri
... R1
R0
λ


i+1
∑
(i)
(i)
.
(2.12)
S + λ j=0 aj Ri−j+1 




0




...


=



0




0


 i+1

∑ (i)
aj Ri−j+1 + λS (i)
j=0
Для того чтобы система (2.12) соответствовала первоначальной, определенной в (2.7), необходимо выполнение следующих тождеств, обеспечива(i+1)
ющих определение коэффициентов {aj
} на итерации (i + 1):
(i+1)
a1
(i)
(i)
= a1 + λai ,
40
В матричной форме:

R0
R1
R2
 R1
R
R1
0

 ...
...
...
RM −1 RM −2 RM −3
... RM −2
... RM −3
...
...
...
R1
  

RM −1
a1
R1
  

RM −2 
  a2  =  R2  . (2.6)




...
...
... 
R0
aM
RM
Решение системы линейных уравнений произвольного вида имеет сложность M 3 . Однако сложность решения можно свести к M 2 , если воспользоваться тем фактом, что матрица автокорреляционных коэффициентов в системе (2.6) является симметричной и теплицевой. Для решения такой системы уравнений рассмотрим итеративный алгоритм, предложенный Левинсоном и Дарбином [10, 11] (метод также описан в [12]), который описывает
итеративную процедуру определения неизвестных параметров aj .
1.2.4 Алгоритм Левинсона-Дарбина
Метод предполагает нахождение вектора неизвестных значений вектора
a длины M , последовательно определяя неизвестные для систем меньшего
размера от 1 до M − 1. Перед тем как приступить к обоснованию алгоритма
преобразуем, систему (2.6) за счет переноса в левую часть вектора свободных
членов:

   
R1
R0
R1
. . . RM −2 RM −1
−1
0
 R2




a
0
R
R
.
.
.
R
R
1
0
M −3
M −2   1 

=   . (2.7)
 ...
...
...
...
...
. . .   . . .  . . .
RM
RM −1
RM −2
...
R1
37
R0
aM
0
Теперь столбец свободных членов стал нулевым, однако матрица автокорреляционных коэффициентов (2.7) перестала быть симметричной. Для того,
чтобы вернуть матрице (2.7) это свойство, добавим к ней сверху еще одну
строку. Это также приведет к добавлению ненулевого элемента в первую позицию вектора-столбца свободных членов. Обозначим этот элемент как SM :

  

R0
R1
R2
... RM −1 RM
−1
SM
R
  

R0
R1
... RM −2 RM −1 
 1
  a1   0 

  

R1
R0
... RM −3 RM −2   a2  =  0  . (2.8)
 R2

  

 ...
...
...
...
...
...   ...   ... 
RM RM −1 RM −2 ...
R1
R0
aM
0
Обозначим вектор-столбец коэффициентов ai через a = (a1 , ..., aM )T , а
матрицу автокорреляционных коэффициентов Ri через R(M ) , где M определяет число элементов в матрице.
Алгоритм Левинсона-Дарбина представляет собой итеративную процедуру решения системы уравнений (2.8). Опишем первую и (i+1)-ю итерации
алгоритма.
Решение системы уравнений на итерации 1. На первой итерации рас(1)
сматривается система(2.8), содержащая неизвестные коэффициент a1 и
(1)
значение S . Верхний индекс указывает на номер итерации, на которой
осуществляется вычисление неизвестного. В матричном виде система (2.8)
выглядит следующим образом:
) ( )
(
)(
−1
R0 R1
S1
=
.
(1)
R1 R0
0
a1

R0
R
 1

R2

 ...
Ri
R1
R0
R1
...
Ri−1
R2
R1
R0
...
Ri−2
R3
R2
R1
...
Ri−3
... Ri−1
... Ri−2
... Ri−3
...
...
... R1
  

1
Ri
S (i)
 (i)  

Ri−1 
 a 1   0 
  (i)  

Ri−2  a2  =  0  . (2.9)
  

...   ...   ... 
(i)
R0
0
ai
Для перехода к следующей итерации (i + 1) необходимо выполнить следующие преобразования:
1. Добавить справа к матрице R(i) дополнительный столбец коэффициентов
автокорреляции (Ri+1 , Ri , Ri−1 , ...R1 )T , состоящий из (i + 1) элемента.
2. Чтобы соблюсти размерность матричного умножения необходимо добавить к вектору-столбцу коэффициентов a(i) один нулевой элемент (верхний индекс i показывает, что данный вектор коэффициентов был получен на итерации i). На данном этапе решение модифицированной системы
совпадает с решением системы (2.9).
3. Добавить к модифицированной на первом шаге матрице R(i) строку
(Ri+1 , Ri , Ri−1 , ...R0 ), содержащую (i + 2) элемента.
4. Для компенсации этой вставки и получения новой корректной системы
уравнений необходимо добавить к вектору-столбцу свободных членов
элемент, который будет равен результату произведения добавленной строки и модифицированного вектора коэффициентов a(i) :
i+1
∑
(i)
aj Ri−j+1 ,
j=0
(1)
Отсюда можно найти неизвестные a1 и S (1) :
(1)
(i)
(i)
где a0 = −1, а a(i+1) = 0.
R1
;
R0
На этом первая итерация считается законченной.
Решение системы уравнений на итерации i+1. При выполнении данной
итерации предыдущие шаги от 1 до i выполнены и мы уже имеем решение
для следующей системы:
В результате выполненных преобразований получим новую систему вида:



 
S (i)
−1
R0
R1
...
Ri
Ri+1


0
 R
 (i)  

R0
... Ri−1
Ri 
 1
 a 1  

0

  (i)  

R1
... Ri−2 Ri−1  a2  
 R2

...

  = 
 . (2.10)
 ...

...
...
...
...   ...  

  (i)  

0

 Ri
Ri−1 ... R0
R 1  a i  
 i+1

∑ (i)
a
R
Ri+1
Ri
... R1
R0
i−j+1
0
j
38
39
a1 =
(1)
S (1) = −R0 + R1 a1 = −R0 + R1
R1
R2 − R02
= 1
.
R0
R0
j=0
Документ
Категория
Без категории
Просмотров
0
Размер файла
3 724 Кб
Теги
veselov
1/--страниц
Пожаловаться на содержимое документа