CC BY
99
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
Сер. 10. 2010. Вып. 2
УДК 004.627:004.932.2 В. В. Окунев
ОБ ОДНОМ МЕТОДЕ ОПТИМИЗАЦИИ
ФРАКТАЛЬНОГО АЛГОРИТМА СЖАТИЯ ИЗОБРАЖЕНИЙ
1. Введение. С тех пор как Барнсли в 1991 г. представил основанную на IFS (системах итерируемых функций) схему кодирования цифровых графических изображений как схему с потенциально высокой степенью компрессии [1], фрактальное сжатие значительно привлекло внимание исследователей в области обработки изображений. Чуть позже Жакен предложил алгоритм, базирующийся на PIFS (системах итерируемых в определенных областях функций) [2]. Несмотря на то, что соотношение качества изображения и степени сжатия для данного алгоритма в целом было признано ниже, чем для других схем сжатия (например, JPEG) [3], новизна и простота подхода вызвали интерес в научном сообществе.
2. Классический алгоритм. Рассмотрим классический подход к фрактальному сжатию изображений. Его существо состоит в поиске подобных участков (блоков) изображения и вычислении параметров преобразований, позволяющих переводить одни блоки в другие внутри найденных пар (отсюда происходит название «фрактальный алгоритм» - самоподобие является ключевым свойством, характеризующим фракталы). В сжатом файле вместо информации о пикселях хранится информация о найденных отображениях. Для декодирования эти преобразования итеративно применяются к любому начальному изображению до тех пор, пока изображение, полученное на текущем шаге, не будет отличаться от предыдущего (с учетом квантования значений пикселей).
Ранговыми блоками назовем участки изображения, которые покрывают его полностью, при этом взаимно не перекрываясь; доменными блоками (или доменами) -некоторым образом выбранные участки изображения, которые могут перекрываться. Классический алгоритм включает в себя разбиение изображения f на ранговые блоки Ri, покрытие изображения f доменными блоками Dj и нахождение для каждого рангового блока соответствующего домена одновременно с настройкой параметров преобразования (контрастности и яркости) для наилучшего соответствия. Таким образом, основной вычислительный шаг алгоритма - это сравнение доменного и рангового блоков. Для каждого рангового блока сравниваются варианты преобразования всех доменов к этому ранговому блоку, такие как аффинные преобразования (параллельный перенос, сжатие и один из восьми вариантов ориентации). Поиск соответствия между доменными и ранговыми блоками ведется следующим образом: к выбранному домену применяется один из восьми базовых поворотов/отражений, затем доменный блок сжимается, чтобы соответствовать размеру рангового блока, и в конечном
Окунев Вадим Вячеславович — аспирант кафедры компьютерных технологий и систем факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Научный руководитель: доктор физико-математических наук, проф. Е. И. Веремей. Количество опубликованных работ: 3. Научное направление: алгоритмы анализа и сжатия цифровых графических изображений. E-mail: vadik-okunev@yandex.ru.
© В. В. Окунев, 2010
итоге методом наименьших квадратов вычисляются оптимальные значения яркости и контрастности. Для расчета оптимальной контрастности (с) и яркости (Ь) необходимо определить соответствующие переменные, которые бы минимизировали выражение У~] У~] \cdij + Ь — ]2. Здесь и - значения пикселей рангового и доменного (сжа-
г j
того до размеров рангового) блоков соответственно. Эти пиксели находятся в прямоугольных массивах с М строками и N столбцами. Решение представляется в следующем виде [4]:
а а -
с=0’ Ь = 'Г-13ё’
где а = ЕЕ(<1,3-3){ггз -г); /? = ЕЕ “ ^)2; Л=Жм ЕЕ<% г=шЕЕп3-
г j г j г j г j
Фрактальное сжатие требует большого объема вычислений, с точки зрения как поиска подобных блоков (для каждого рангового блока необходим перебор достаточно большого количества доменных блоков), так и их сравнения. Первые реализации фрактальных методов компрессии были очень медленными: на самых мощных рабочих станциях вычисления обычно занимали много часов, а иногда и дней. По сравнению с существующими методами фрактальный алгоритм компрессии до настоящего времени достаточно медленный. Это препятствует использованию фрактальных методов сжатия на практике. Попытки увеличить скорость кодирования делались, в частности, в следующих направлениях:
• аппаратное распараллеливание вычислений;
• классификация блоков;
• выделение характеристических особенностей блоков.
Целью данной статьи является разработка метода оптимизации фрактального алгоритма, обеспечивающего существенное сокращение необходимого времени вычислений при приемлемых потерях в качестве декодированного изображения.
3. Оптимизация классического алгоритма. Для компенсации временных затрат при кодировании изображения рассмотрим классификационный метод, основанный на применении DCT-преобразований и вейвлет-коэффициентов Хаара. В целях классификации всего множества блоков изображения для каждого из них вычисляется характеристический вектор, состоящий из двух компонент. Первую компоненту вектора условно можно назвать вейвлет-индексом, а вторую - БОТ-индексом. Перед тем, как описать схему вычисления компонент характеристического вектора, введем следующие обозначения. Матрицу пикселей исходного изображения обозначим через ], а ее элементы - через fij, i, j = 1, п. Под изменением яркости будем понимать изменение значения всех пикселей блока на константу: f'ij = fij + а; под изменением контрастности - умножение значений всех пикселей на константу: f'j = afij, 0 < а < 1.
Пусть Г - матрица прямого DCT-преобразования матрицы f. Ее элементы обозначим через i,j = 1, п. В [5] приведены равенства, которые задают правила изменения матрицы Г при следующих преобразованиях, применяемых к матрице f:
Г^ = ( — 1)i+^Fij, f'ij = ,
= ( — 1)j+1Fij, ,
Г^ = Щг, fij = fji,
где f' - матрица преобразованного в соответствии с указанными равенствами исходного блока; Р' - матрица БСТ-преобразования блока , г,^ = 1, п. Воспользовавшись
этими правилами, можно существенно уменьшить количество операций, необходимых для вычисления матрицы Е для блоков, полученных в результате применения преобразований симметрии (таких как отражение или поворот) к исходным блокам.
Приведем алгоритм вычисления DCT-индекса с использованием полученных DCT-матриц (рис. 1):
1) рассматривается «верхний левый угол» матрицы Е, соответствующей данному блоку / (подматрица Р = {-^Л^цз);
2) поочередно оцениваются элементы этой подматрицы от 1-го до г-го (где г - выбранная размерность индекса, г Е [1; 8]) в порядке, показанном на рис. 1; если текущий элемент меньше нуля, то соответствующему разряду индекса в двоичном представлении, начиная слева, присваивается значение 0, если больше или равен нулю - 1;
3) полученное число преобразуется к десятичному виду.
1 5
2 4 6
3 7 00
Рис. 1. Порядок выбора элементов для вычисления DCT-индекса
Пример 1. Рассмотрим блок размером 8 х 8 пикселей (рис. 2) и вычислим для него DCT-индекс.
Рис. 2. Блок изображения размером 8 х 8 пикселей Для этого приведем матрицу значений самого блока
Г :
100 76 77 90 133 136 9 9 5 1 1
92 62 88 118 135 133 111 117
83 47 74 117 110 114 110 126
81 48 69 96 99 110 115 132
/ = 86 56 84 129 136 124 135 157 ,
88 71 108 165 158 126 152 179
100 93 118 174 148 128 161 186
130 123 118 173 145 126 174 180
зе DCT-преобразования
930, 50 -170, 88 -36, 71 1 4 8 0 О 36, 02 0, 87 -3, 55
-128, 05 22, 68 -11, 02 62, 32 -23, 45 -17, 76 30, 89 7, 48
64, 40 33, 23 -0, 78 10, 28 -14, 68 -5, 30 - -16, 53 -0, 37
19,81 - -14, 99 -19, 55 5, 76 2, 74 -11, 06 0, 66 11, 33
-17, 50 11,41 24, 68 20, 80 -15, 00 -5, 33 -5, 08 3, 27
-22, 25 -0, 62 10,13 1 6 7, - 7 2 2, 1 0, 16 1, 44 -1, 78
0, 65 1,40 0, 97 0, 55 3, 10 2, 48 -5, 72 0, 69
1,04 -2, 72 -1, 57 0, 55 2, 83 1, 61 -0, 34 -2, 60
и «верхний левый угол» матрицы Г
Г =
930,50 -170,88 -36,71
-128, 05 22, 68 -11,02
64,40 33,23 -0,78
Пусть размерность DCT-индекса равна 5. Приведем табл. 1, поясняющую вычисление значений бит DCT-индекса. В результате двоичная запись индекса будет выглядеть как 001102, что при переводе в десятичную систему счисления даст число 6 - оно и будет искомым.
Таблица 1. Формирование значения ЮСТ-индекса
Номер Значение Знак Значение бита
элемента элемента элемента БСТ-индекса
1 -170,88 - 0
2 -128,05 - 0
3 64,40 + 1
4 22,68 + 1
5 -36,71 - 0
Для вычисления вейвлет-индекса используются вейвлет-коэффициенты Хаара е2, ез, е4 и е5. Их можно найти следующим образом. Пусть имеется изображение А размером п х п пикселей, п = 2т, т € N. Коэффициент е\ будет равен среднему значению всех пикселей:
1 П П
61 =
П i=1 3 = 1
Для того чтобы найти оставшиеся коэффициенты, изображение делится на 4 равные части так, как показано на рис. 3. Определив среднее значение каждого из получившихся блоков и вычитая его из е1 , найдем искомые величины коэффициентов:
1 т т 1 т п
= е1 _ ез = е1_^Е Е оУ’
г=1 ^=1 г=1 ^=т+1
1 пт 1 п п
е4 = е1 ~ ^2 X Еа^ е5 = е1 _ ^2 X X
г=т+1 ^=1 г=т+1 ^=т+1
После вычисления они выстраиваются в порядке возрастания значения. Возможны всего 4! = 24 варианта перестановки четырех индексов коэффициентов. Им можно поставить в соответствие числа от 1 до 24, которые и будут искомыми значениями вейвлет-индекса.
ез
е4 е5
Рис. 3. Вейвлет-коэффициенты Хаара
Пример 2. Вычислим вейвлет-индекс для блока изображения, приведенного на рис. 2. Для этого выделим четыре подматрицы из матрицы / значений пикселей
блока: /2 = {/у}г^=1“4) /з = {/и'}*=тд, ^=5|8> /4 = {/у}*=5|8, ^=М> ^"5 = {/у}*,^=5|8-Чтобы получить коэффициент е1, рассчитаем среднее значение матрицы /. Для нахождения коэффициентов е2-е5 определим средние значения матриц /2-/5 и вычтем их из е1 соответственно. В результате имеем: е1 = 115.3125; е2 = 115.3125 — 82.375 = 32.9375; е3 = 115.3125 — 118.4375 = —3.125; е4 = 115.3125 — 109.5 = 5.8125; е5 = 115.3125 — 150.9375 = —35.625. Расположив данные коэффициенты в порядке возрастания, получим последовательность коэффициентов е5 ез е4 в2 и соответствующую ей последовательность индексов коэффициентов 5 3 4 2. Выберем способ нумерации последовательностей индексов коэффициентов так, как показано в табл. 2, тогда выявленная последовательность будет соответствовать числу 22 - это и будет искомый вейвлет-индекс.
Таким образом, характеристический вектор для блока изображения, приведенного на рис. 2, будет равен ^ «_
Такой метод позволяет разбить все множество блоков на классы, количество которых в пределе равно 24 • 2Г, где г - количество бит DCT-индекса, используемых для его формирования, и впоследствии вести поиск подобных блоков только внутри соответствующих классов, что существенно сокращает время перебора блоков -основного вычислительного шага во фрактальном сжатии. Однако при применении алгоритма может возникать следующая проблема: для отдельно взятого рангового блока может не найтись ни одного домена с полностью совпадающим значением характеристического вектора. В целях преодоления таких ситуаций перед началом
Последовательность индексов вейвлет-коэффициентов Соответствующее значение вейвлет-индекса
2 3 4 5 1
2 3 5 4 2
2 4 3 5 3
2 5 4 3 6
5 4 3 2 24
поиска подобных блоков задается критерий близости характеристических векторов по DCT-индексу (вейвлет-индекс считается более грубым критерием подобия, и если для конкретного рангового блока не нашлось ни одного доменного блока с совпадающим значением этого индекса, то поиск ведется заново по всему множеству доменов). Степень близости в смысле DCT-индекса определяется так: к DCT-компонентам двух характеристических векторов применяется логическая операция «исключающее ИЛИ», и чем больше результат операции, тем менее близкими считаются векторы.
Итоговый алгоритм выглядит следующим образом:
1. Строится разбиение изображения на ранговые и доменные блоки.
2. Для каждого блока вычисляется характеристический вектор.
3. Для первого рангового блока ведется поиск подобного ему домена из множества таких доменов, для которых вейвлет-индекс совпадает с вейвлет-индексом рангового.
4. При нахождении такого блока сравниваются DCT-индексы. В случае их совпадения проводится сопоставление (сравнение блоков по схеме, приведенной в п. 2, а также запись координат блоков и полученных значений в файл), в случае несовпадения - дополнительный поиск доменов с совпадающим значением вейвлет-индекса, затем среди них выбирается наиболее близкий домен в смысле близости DCT-индексов и с найденным блоком проводится сопоставление.
5. Если не нашлось доменов с совпадающей величиной вейвлет-индекса, поиск ведется заново по всему множеству доменов, пока не будет установлен домен, доставляющий минимальное значение норме разности блоков (с учетом сжатия домена и настроек яркости и контрастности).
6. Для второго и последующих ранговых блоков повторяются шаги 3-5, пока все ранговые блоки не окажутся «покрытыми».
4. Заключение. Сравнительные результаты компрессии приведены в табл. 3, в которой FrCl - классический алгоритм фрактального сжатия, FrOpt - оптимизированный алгоритм с применением классификации с помощью характеристических векторов, JPEG - алгоритм JPEG.
Таблица 3. Сравнительные результаты компрессии
а б
Алгоритм Время компрессии, с PSNR, дБ Алгоритм Степень сжатия, % PSNR, дБ
FrCl 476 31.0309 JPEG 8.2147 32.6294
FrOpt 213 30.1055 FrOpt 4.3929 30.1055
ІІ9
На рис. 4 приведено оригинальное изображение, на рис. 5 - восстановленные изображения, соответствующие табл. 3, а, на рис. 6 - соответствующие табл. 3, б.
Рис. 6. Изображения, восстановленные с помощью алгоритмов JPEG (а) и FrOpt (б)
Оптимизация классического алгоритма с использованием классификационного метода, основанного на применении DCT-преобразований, была ранее реализована в среде MATLAB [6]. Последующая модификация алгоритма путем внедрения характеристических векторов, а также реализация итогового алгоритма на языке C++ дали возможность существенно улучшить результаты по сравнению с предыдущими показателями. Полученные результаты позволяют сделать вывод о том, что использование предложенной оптимизированной схемы компрессии обеспечивает преимущество по сравнению с классическими методами, состоящее в существенном снижении временных затрат при допустимых потерях в качестве.
Литература
1. Fisher Y., Jacobs E. W., Boss R. D. Fractal image compression using iterated transform // NOSC: Technical Report 1408. 1991.
2. Jacquin A. Image coding based on a fractal theory of iterated contractive image transformations // IEEE Transactions on Image Processing. 1992. N 1. P. 18-30.
3. Ghosh S. K., Mukherjee J., Das P. P. Fractal image compression: a randomized approach // Pattern Recognition Letters. 2004. Vol. 25. P. 1013-1024.
4. Welstead S. Fractal and wavelet image compression techniques. Bellingham (Wash.): SPIE Optical Engineering Press, 1999. XV, 232 p.
5. Ватолин Д. С. Использование ДКП для ускорения фрактального сжатия изображений // Программирование. 1999. № 3. С. 51-57.
6. Окунев В. В. Оптимизация фрактального алгоритма сжатия изображений // Процессы управления и устойчивость : Труды 38-й Междунар. науч. конференции аспирантов и студентов / под ред. А. В. Платонова, Н. В. Смирнова. СПб.: Изд-во С.-Петерб. ун-та, 2007. С. 423-428.
Статья рекомендована к печати проф. Е. И. Веремеем.
Статья принята к печати 24 декабря 2009 г.
