Fil 17 декабря 2013 в 14:09

Изобретаем JPEG

• Алгоритмы ,
• Обработка изображений

• Tutorial

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

Знание алгоритма JPEG очень полезно не только для сжатия изображений. В нем используется теория из цифровой обработки сигналов, математического анализа, линейной алгебры, теории информации, в частности, преобразование Фурье, кодирование без потерь и др. Поэтому полученные знания могут пригодиться где угодно.

Если есть желание, то предлагаю пройти те же этапы самостоятельно параллельно со статьей. Проверить, насколько приведенные рассуждения подходят для разных изображений, попытаться внести свои модификации в алгоритм. Это очень интересно. В качестве инструмента могу порекомендовать замечательную связку Python + NumPy + Matplotlib + PIL(Pillow). Почти вся моя работа (в т. ч. графики и анимация), была произведена с помощью них.

Внимание, трафик! Много иллюстраций, графиков и анимаций (~ 10Мб). По иронии судьбы, в статье про JPEG всего 2 изображения с этим форматом из полусотни.

Каков бы ни был алгоритм сжатия информации, его принцип всегда будет один - нахождение и описание закономерностей. Чем больше закономерностей, тем больше избыточности, тем меньше информации. Архиваторы и кодеры обычно «заточены» под конкретный тип информации, и знают где можно их найти. В некоторых случаях закономерность видна сразу, например картина голубого неба. Каждый ряд его цифрового представления можно довольно точно описать прямой.

Будем тренироваться на кошках енотах. В качестве примера взято серое изображение, приведенное выше. Оно хорошо совмещает как однородные области, так и контрастные. А если мы научимся сжимать серое, то и с цветным не будет проблем.

Векторное представление

Для начала проверим насколько зависимы два соседних пикселя. Логично предположить, что скорее всего, они будут очень похожи. Проверим это для всех пар изображения. Отметим их на координатной плоскости точками так, что значение точки по оси X - значение первого пикселя, по оси Y - второго. Для нашего изображения размером 256 на 256 получим 256*256/2 точек:

Предсказуемо, что большинство точек находится на или рядом с прямой y=x (а их там еще больше, чем видно на рисунке, так как они многократно накладываются друг на друга, и, к тому же, они полупрозрачные). А раз так, то было бы проще работать, повернув их на 45°. Для этого нужно выразить их в другой системе координат.

Базисные вектора новой системы, очевидно, такие: . Вынуждены делить на корень из двойки, чтобы получить ортонормированную систему (длины базисных векторов равны единичке). Здесь показано, что некоторая точка p = (x, y) в новой системе будет представлена как точка (a 0 , a 1). Зная новые коэффициенты, мы легко можем получить старые обратным поворотом. Очевидно, первая (новая) координата является средним, а вторая - разностью x и y (но деленные на корень из 2). Представьте, что вам предложено оставить только одно из значений: либо a 0 , либо a 1 (то есть другое приравнять нулю). Лучше выбрать a 0 , так как значение a 1 и так, скорее всего, будет около нуля. Вот, что получится, если мы восстановим изображение только по a 0:

Увеличение в 4 раза:

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

Это один и тот же график, но с разных точек зрения. Красные линии - оси, которые напрашивались сами собой. Им соответствуют вектора: . Напоминаю, что приходится делить на некоторые константы, чтобы длины векторов стали равны единице. Таким образом, разложив по такому базису, мы получим 3 значения a 0 , a 1 , a 2 , причем a 0 важнее a 1 , а a 1 важнее a 2 . Если мы выбросим a 2 , то график «сплющится» в направлении вектора e 2 . Этот и так довольно не толстый трехмерный лист станет плоским. Он потеряет не так много, зато мы избавимся от трети значений. Сравним изображения, восстановленные по тройкам: (a 0 , 0, 0), (a 1 , a 2 , 0) и (a 0 , a 1 , a 2). В последнем варианте мы ничего не выбросили, поэтому получим оригинал.

Увеличение в 4 раза:

Второй рисунок уже хорош. Резкие участки немного сгладились, но в целом картинка сохранилась очень неплохо. А теперь, точно так же поделим на четверки и визуально определим базис в четырехмерном пространстве… А, ну да. Но можно догадаться, каким будет один из векторов базиса, это: (1,1,1,1)/2. Поэтому можно посмотреть проекцию четырехмерного пространства на пространство, перпендикулярное вектору (1,1,1,1), чтобы выявить другие. Но это не лучший путь.
Наша цель - научиться преобразовывать (x 0 , x 1 , ..., x n-1) к (a 0 , a 1 , ..., a n-1) так, что каждое значение a i все менее важно, чем предыдущие. Если мы сможем так делать, то, возможно, последние значения последовательности вообще можно будет выбросить. Вышеприведенные опыты намекают, что можно. Но без математического аппарата не обойтись.
Итак, нужно преобразовать точки к новому базису. Но сначала необходимо найти подходящий базис. Вернемся к первому эксперименту разбиения на пары. Будем считать обобщенно. Мы определили базисные векторы:

Выразили через них вектор p :

или в координатах:

Чтобы найти a 0 и a 1 нужно спроецировать p на e 0 и e 1 соответственно. А для этого нужно найти скалярное произведение

аналогично:

В координатах:

Часто бывает удобнее проводить преобразование в матричной форме.

Тогда A = EX и X = E T A. Это красивая и удобная форма. Матрица E называется матрицей преобразования и является ортогональной, с ней мы еще встретимся.

Переход от векторов к функциям.

С векторами малых размерностей работать удобно. Однако мы хотим находить закономерности в бОльших блоках, поэтому вместо N-мерных векторов удобнее оперировать последовательностями, которыми представлено изображение. Такие последовательности я буду называть дискретными функциями, так как следующие рассуждения применимы и к непрерывным функциям.
Возвращаясь к нашему примеру, представим такую функцию f(i), которая определена всего в двух точках: f(0)=x и f(1)=y. Аналогично зададим базисные функции e 0 (i) и e 1 (i) на основе базисов e 0 и e 1 . Получим:

Это очень важный вывод. Теперь во фразе «разложение вектора по ортонормированным векторам» мы можем заменить слово «вектор» на «функция» и получить вполне корректное выражение «разложение функции по ортонормированным функциям». Не беда, что мы получили такую куцую функцию, так как такие же рассуждения работают и для N-мерного вектора, который можно представить как дискретную функцию с N значениями. А работа с функциями нагляднее, чем с N-мерными векторами. Можно и наоборот, представить такую функцию как вектор. Более того, обычную непрерывную функцию можно представить бесконечномерным вектором, правда уже не в евклидовом, а гильбертовом пространстве. Но мы туда не пойдем, нас будут интересовать только дискретные функции.
А наша задача нахождения базиса превращается в задачу нахождения подходящей системы ортонормированных функций. В следующих рассуждениях предполагается, что мы уже как-то определили набор базисных функций, по которым и будем раскладывать.
Допустим, у нас есть некоторая функция (представленная, например, значениями), которую мы хотим представить в виде суммы других. Можно представлять этот процесс в векторном виде. Для разложения функции нужно «спроецировать» ее на базисные функции по очереди. В векторном смысле вычисление проекции дает минимальное сближение исходного вектора к другому по расстоянию. Помня о том, что расстояние вычисляется с помощью теоремы Пифагора, то аналогичное представление в виде функций дает наилучшее среднеквадратичное приближение функции к другой. Таким образом, каждый коэффициент (k) определяет «близость» функции. Более формально, k*e(x) - лучшее среднеквадратичное приближение к f(x) среди l*e(x).
В следующем примере показан процесс приближения функции только по двум точкам. Справа - векторное представление.

Применительно к нашему эксперименту разбивания на пары, можно сказать, что эти две точки (0 и 1 по абсцисс) - пара соседних пикселей (x, y).
То же самое, но с анимацией:

Если мы возьмем 3 точки, то нужно рассматривать трехмерные вектора, однако приближение будет точнее. А для дискретной функции с N значениями нужно рассматривать N-мерные вектора.
Имея набор полученных коэффициентов, можно легко получить исходную функцию, просуммировав базисные функции, взятые с соответствующими коэффициентами. Анализ этих коэффициентов может дать много полезной информации (в зависимости от базиса). Частным случаем этих соображений является принцип разложения в ряд Фурье. Ведь наши рассуждения применимы к любому базису, а при разложении в ряд Фурье берется вполне конкретный.

Дискретные преобразования Фурье (ДПФ)

В предыдущей части мы пришли к выводу, что неплохо было бы разлагать функцию на составные. В начале 19 века Фурье тоже задумался над этим. Правда картинки с енотом у него не было, поэтому ему пришлось исследовать распределение тепла по металлическому кольцу. Тогда он выяснил, что очень удобно выражать температуру (и ее изменение) в каждой точке кольца как сумму синусоид с разными периодами. «Фурье установил (рекомендую к прочтению , интересно), что вторая гармоника затухает в 4 раза быстрее, чем первая, а гармоники более высоких порядков затухают с ещё большей скоростью».
В общем, вскоре оказалось, что периодичные функции замечательно раскладываются на сумму синусоид. А так как в природе существует много объектов и процессов, описывающимися периодичными функциями, то появился мощный инструмент их анализа.
Пожалуй, один из самых наглядных периодических процессов - это звук.

• 1-й график - чистый тон частотой 2500 герц.
• 2-й - белый шум. Т. е. шум c равномерно распределенными частотами по всему диапазону.
• 3-й - сумма первых двух.

Если бы мне дали значения последней функции на тот момент, когда я не знал про ряды Фурье, и попросили проанализировать их, то я бы точно растерялся и не смог бы сказать ничего путного. Ну, да, какая-то функция, но как понять, что там есть что-то упорядоченное? Но если бы я догадался прослушать последнюю функцию, то ухо уловило бы чистый тон среди шума. Хотя и не очень хорошо, так как я специально при генерации подобрал такие параметры, чтобы на суммарном графике сигнал визуально растворился в шуме. Как я понял, до сих пор точно не уставлено, как слуховой аппарат делает это. Однако, недавно стало ясно, что он не раскладывает звук на синусоиды. Возможно, когда-нибудь мы поймем как это происходит, и появятся более совершенные алгоритмы. Ну, а мы пока по старинке.
Почему бы не попробовать взять синусоиды в качестве базиса? На самом деле мы фактически уже сделали это. Вспомним наше разложение на 3 базисных вектора и представим их на графике:

Да-да, знаю, это выглядит как подгонка, но с тремя векторами трудно ожидать большего. Зато теперь понятно, как получить, например, 8 базисных векторов:


Не очень сложная проверка показывает, что эти вектора попарно перпендикулярны, т. е. ортогональны. Это значит, их можно использовать как базис. Преобразование по такому базису широко известно, и называется дискретным косинусным преобразованием (DCT). Думаю, из приведенных графиков понятно как получается формула DCT-преобразования:

Это все та же формула: A = EX с подставленным базисом. Базисные вектора указанного DCT (они же векторы-строки матрицы E) ортогональны, но не ортонормированы. Это следует помнить при обратном преобразовании (не буду останавливаться на этом, но, кому интересно - у inverse DCT появляется слагаемое 0.5*a 0 , так как нулевой вектор базиса больше остальных).
На следующем примере показан процесс приближения промежуточных сумм к исходным значениям. На каждой итерации очередной базис умножается на очередной коэффициент и прибавляется к промежуточной сумме (то есть так же, как и в ранних опытах над енотом - треть значений, две трети).

Но, все-таки, несмотря на некоторые доводы о целесообразности выбора такого базиса, реальных аргументов пока нет. Действительно, в отличие от звука, целесообразность разложения изображения на периодические функции гораздо менее очевидна. Впрочем, изображение действительно может быть слишком непредсказуемым даже на небольшом участке. Поэтому, картинку делят на достаточно маленькие кусочки, но не совсем крохотные, чтобы разложение имело смысл. В JPEG изображение «нарезается» на квадраты 8x8. В пределах такого кусочка фотографии обычно очень однородны: фон плюс небольшие колебания. Такие области шикарно приближаются синусоидами.
Ну, допустим, этот факт более или менее понятен интуитивно. Но появляется нехорошее предчувствие насчет резких цветовых переходов, ведь медленно изменяющиеся функции нас не спасут. Приходится добавлять разные высокочастотные функции, которые справляются со своей работой, но побочно проявляются на однородном фоне. Возьмем изображение 256x256 с двумя контрастными областями:

Разложим каждую строку с помощью DCT, получив, таким образом, по 256 коэффициентов на строку.
Затем оставим только первые n коэффициентов, а остальные приравняем нулю, и, поэтому, изображение будет представлено в виде суммы только первых гармоник:

Число на картинке - количество оставленных коэффициентов. На первом изображении осталось только среднее значение. На второй уже добавилась одна низкочастотная синусоида, и т. д. Кстати, обратите внимание на границу - несмотря на все лучшее приближение, рядом с диагональю хорошо заметны 2 полоски, одна светлее, другая темнее. Часть последнего изображения увеличенного в 4 раза:

И вообще, если вдали от границы мы видим первоначальный равномерный фон, то при приближении к ней, амплитуда начинает расти, наконец достигает минимального значения, а затем резко становится максимальным. Это явление известно как эффект Гиббса.

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


При изучении DCT может сложиться ложное впечатление, что всегда вполне достаточно всего нескольких первых (низкочастотных) коэффициентов. Это верно для многих кусочков фотографий, тех, чьи значения не меняются резко. Однако, на границе контрастных участков значения будут резво «скакать» и даже последние коэффициенты будут велики. Поэтому, когда слышите о свойстве сохранения энергии DCT, делайте поправку на то, что оно относится ко многим видам встречаемых сигналов, но не ко всем. Для примера подумайте, как будет выглядеть дискретная функция, коэффициенты разложения которой равны нулю, кроме последнего. Подсказка: представьте разложение в векторном виде.
Несмотря на недостатки, выбранный базис является одним из лучших на реальных фотографиях. Чуть позже мы увидим небольшое сравнение с другими.

DCT vs все остальное

Когда я изучал вопрос ортогональных преобразований, то, честно говоря, меня не очень убеждали доводы о том, что все вокруг - это сумма гармонических колебаний, поэтому нужно и картинки раскладывать на синусоиды. А может быть лучше подойдут какие-нибудь ступенчатые функции? Поэтому искал результаты исследований об оптимальности DCT на реальных изображениях. То, что «Именно DCT чаще всего встречается в практических приложениях благодаря свойству «уплотнения энергии»» написано везде. Это свойство означает, что максимальное количество информации заключено в первых коэффициентах. А почему? Нетрудно провести исследование: вооружаемся кучей разных картинок, различными известными базисами и начинаем считать среднеквадратичное отклонение от реального изображения для разного количества коэффициентов. Нашел небольшое исследование в статье (использованные изображения ) по этой методике. В ней приведены графики зависимости сохраненной энергии от количества первых коэффициентов разложений по разным базисам. Если вы просмотрели графики, то убедились, что DCT стабильно занимает почетное… эмм… 3-место. Как же так? Что еще за KLT преобразование? Я восхвалял DCT, а тут…

KLT

Все преобразования, кроме KLT, являются преобразованиями с постоянным базисом. А в KLT (преобразование Карунена-Лоэва) вычисляется самый оптимальный базис для нескольких векторов. Он вычисляется таким образом, что первые коэффициенты дадут наименьшую среднеквадратичную погрешность суммарно для всех векторов. Похожую работу мы проводили ранее вручную, визуально определяя базис. Сначала кажется, что это здравая идея. Мы могли бы, например, разбивать изображение на небольшие секции и для каждой вычислять свой базис. Но мало того, что появляется забота хранения этого базиса, так еще и операция его вычисления достаточно трудоемкая. А DCT проигрывает лишь немного, и к тому же у DCT существуют алгоритмы быстрого преобразования.

DFT

DFT (Discrete Fourier Transform) - дискретное преобразование Фурье. Под этим названием иногда упоминается не только конкретная трансформация, но и весь класс дискретных трансформаций (DCT, DST...). Посмотрим на формулу DFT:

Как вы догадываетесь, это ортогональное преобразование с каким-то комплексным базисом. Так как подобная комплексная форма встречается чуть чаще, чем всегда, то имеет смысл изучить ее вывод.
Может сложится впечатление, что любой чистый гармонический сигнал (с целой частотой) при DCT разложении будет давать только один ненулевой коэффициент, соответствующий этой гармонике. Это не так, поскольку помимо частоты, важна и фаза этого сигнала. Например, разложение синуса по косинусам (подобным образом и в дискретном разложении) будет таким:

Вот вам и чистая гармоника. Она наплодила кучу других. На анимации показаны коэффициенты DCT синусоиды в разных фазах.


Если вам показалось, что столбики вращаются вокруг оси, то вам не показалось.
Значит теперь будем раскладывать функцию на сумму синусоид не просто разных частот, но еще и смещенных по какой-то фазе. Будет удобнее рассмотреть сдвиг фаз на примере косинуса:

Простое тригонометрическое тождество дает важный результат: сдвиг по фазе заменяется суммой синуса и косинуса, взятых с коэффициентами cos(b) и sin(b). Значит, можно раскладывать функции на сумму синусов и косинусов (без всяких фаз). Это распространенная тригонометрическая форма. Однако, гораздо чаще используется комплексная. Для ее получения нужно воспользоваться формулой Эйлера . Просто подставим производные формулы для синуса и косинуса, получим:


Теперь небольшая замена. Верхнее подчеркивание - сопряженное число.

Получим итоговое равенство:

c - комплексный коэффициент, действительная часть которого равна косинусному коэффициенту, а мнимая - синусному. А множество точек (cos(b), sin(b)) является окружностью. В такой записи каждая гармоника входит в разложение и с положительной и с отрицательной частотой. Поэтому в различных формулах Фурье-анализа обычно происходит суммирование или интегрирование от минус до плюс бесконечности. Производить вычисления часто бывает удобнее именно в такой комплексной форме.
Преобразование раскладывает сигнал на гармоники с частотами от одного до N колебаний на области сигнала. Но частота дискретизации составляет N на области сигнала. А по теореме Котельникова (aka теорема Найквиста - Шеннона) частота дискретизации должна по крайней мере в два раза превышать частоту сигнала. Если это не так, то получается эффект появления сигнала с ложной частотой:

Пунктирной линий показан неверно восстановленный сигнал. С таким явлением вы часто сталкивались в жизни. Например, забавное движение колес автомобиля на видео, или муаровый эффект.
Это приводит к тому, что вторая половина из N комплексных амплитуд как будто состоит из других частот. Эти ложные гармоники второй половины являются зеркальным отображением первой и не несут дополнительной информации. Таким образом, у нас остается N/2 косинусов и N/2 синусов (образующих ортогональный базис).
Ладно, базис есть. Его составляющие - гармоники с целым числом колебаний на области сигнала, а значит, крайние значения гармоник равны. Точнее почти равны, так как последнее значение берется не совсем с края. Более того - каждая гармоника почти зеркально симметрична относительно своего центра. Все эти явления особенно сильны на низких частотах, которые нам и важны при кодировании. Это плохо еще и тем, что на сжатом изображении будут заметны границы блоков. Проиллюстрирую DFT-базис с N=8. Первые 2 ряда - косинусные составляющие, последние - синусные:

Обратите внимание на появление дублей составляющих при повышении частоты.

Можете мысленно подумать, как мог бы быть разложен сигнал, значения которого плавно уменьшаются с максимального значения в начале до минимального в конце. Более-менее адекватное приближение смогли бы сделать лишь гармоники ближе к концу, что для нас не очень здорово. На рисунке слева приближение несимметричного сигнала. Справа - симметричного:


С первым дела крайне плохи.
Так может быть сделать как в DCT - уменьшить частоты в 2 или другое количество раз, чтобы количество некоторых колебаний было дробным и границы находились в разных фазах? Тогда составляющие будут неортогональны. И ничего тут не поделать.

DST

Что если вместо косинусов в DCT использовать синусы? Мы получим Discrete Sine Transform (DST). Но для нашей задачи все они неинтересны, так как и целые и половинки периодов синусов близки к нулю на границах. То есть мы получим примерно такое же неподходящее разложение, как и у DFT.

Возвращаясь к DCT

Как у него дела на границах? Хорошо. Есть противофазы и нет нулей.

Все остальное

Не-Фурье преобразования. Не буду описывать.
WHT - матрица состоит только из ступенчатых составляющих со значениями -1 и 1.
Haar - по совместительству ортогональное вейвлет-преобразование.
Они уступают DCT, но легче для вычислений.

Итак, вас посетила мысль придумать свое преобразование. Помните вот что:

1. Базис должен быть ортогонален.
2. С фиксированным базисом вы не сможете превзойти KLT по качеству сжатия. Между тем, на реальных фотографиях DCT почти не уступает.
3. На примере DFT и DST нужно помнить про границы.
4. И помнить, что у DCT есть еще хорошее преимущество - вблизи границ составляющих их производные равны нулю, а значит, переход между соседними блоками будет довольно плавным.
5. У преобразований Фурье существуют быстрые алгоритмы со сложностью O(N*logN), в отличие от вычисления в лоб: O(N 2).

Будет непросто, правда? Впрочем, для некоторых типов изображений можно подобрать лучший базис, чем у DCT.

Двумерные преобразования

Сейчас попробуем провести такой эксперимент. Возьмем, для примера, кусочек изображения.

Его 3D график:

Пройдемся DCT(N=32) по каждой строке:

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

Коэффициент (0,0) получился слишком большим, поэтому на графике он уменьшен в 4 раза.
Итак, что получилось?
Левый верхний угол - самые значащие коэффициенты разложения самых значащих коэффициентов.
Левый нижний угол - самые незначащие коэффициенты разложения самых значащих коэффициентов.
Правый верхний угол - самые значащие коэффициенты разложения самых незначащих коэффициентов.
Правый нижний угол - самые незначащие коэффициенты разложения самых незначащих коэффициентов.
Понятно, что значимость коэффициентов уменьшается, если двигаться по диагонали из левого верхнего угла в правый нижний. А какой важнее: (0, 7) или (7, 0)? Что они вообще означают?
Сначала по строкам: A 0 = (EX T) T = XE T (транспонировали, так как формула A=EX для столбцов), затем по столбцам: A=EA 0 = EXE T . Если аккуратно посчитать, то получится формула:

Таким образом, если вектор раскладывается на синусоиды, то матрица на функции вида cos(ax)*cos(by). Каждый блок 8x8 в JPEG представляется в виде суммы 64-х функций вида:


В Википедии и других источниках такие функции представлены в более удобной форме:

Поэтому коэффициенты (0, 7) или (7, 0) одинаково полезны.
Впрочем, фактически это обычное одномерное разложение на 64 64-мерных базиса. Все вышесказанное применимо не только к DCT, но и к любому ортогональному разложению. Действуя по аналогии, в общем случае получаем N-мерное ортогональное преобразование.
А вот уже 2-мерное преобразование енота (DCT 256x256). Опять же с обнуленными значениями. Числа - количество необнуленных коэффициентов из всех (оставлялись самые значимые значения, находящиеся в треугольной области в левом верхнем углу).

Запомните, что коэффициент (0, 0) называется DC, остальные 63 - AC.

Выбор размера блока

Товарищ спрашивает : почему в JPEG используется разбиение именно 8x8. Из заплюсованного ответа:

The DCT treats the block as if it were periodic and has to reconstruct the resulting jump at the boundaries. If you take 64x64 blocks, you"ll most likely have a huge jump at the boundaries, and you"ll need lots of high-frequency components to reconstruct that to a satisfactory precision

Мол, DCT работает хорошо только на периодических функциях, и если вы возьмете большой размер, то, скорее всего, получите гигантский скачок на границах блока и понадобится много высокочастотных компонентов для его покрытия. Это неверно! Такое объяснение очень похоже на DFT, но не на DCT, так как оно отлично покрывает такие скачки уже первыми составляющими.
На той же странице приводится ответ из MPEG FAQ, с основными аргументами против больших блоков:

• Мало прибыли при разбиении на большие блоки.
• Увеличение вычислительной сложности.
• Высокая вероятность большого количества резких границ в одном блоке, что вызовет эффект Гиббса.

Предлагаю самостоятельно исследовать это. Начнем с первого .


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

Большие блоки не показаны, так как визуально почти неотличимы от 32x32. Теперь посмотрим на абсолютную разность с исходным изображением (усиленную в 2 раза, иначе ничего толком не видно):

8x8 дает лучший результат, чем 4x4. Дальнейшее увеличение размера уже не дает хорошо заметного преимущества. Хотя я всерьез бы задумался над 16x16, вместо 8x8: увеличение сложности на 33% (о сложности в следующем абзаце), дает небольшое, но все-таки видимое улучшение при одинаковом количестве коэффициентов. Однако, выбор 8x8 выглядит достаточно обоснованным и, возможно, является золотой серединой. JPEG был опубликован в 1991. Думаю, что такое сжатие являлось очень сложным для процессоров того времени.

Второй аргумент. Нужно помнить, что при увеличении размера блока потребуется больше вычислений. Давайте оценим насколько. Сложность преобразования в лоб, как мы уже вполне умеем: O(N 2), так как каждый коэффициент состоит из N слагаемых. Но на практике используется эффективный алгоритм быстрого преобразования Фурье (БПФ, Fast Fourier Transform, FFT). Его описание выходит за рамки статьи. Его сложность: O(N*logN). Для двумерного разложения нужно воспользоваться им дважды по N раз. Таким образом, сложность 2D DCT - O(N 2 logN). Теперь сравним сложности вычисления изображения одним блоком и несколькими маленькими:

• Одним блоком (kN)x(kN): O((kN) 2 log(kN)) = O(k 2 N 2 log(kN))
• k*k блоками N*N: O(k 2 N 2 logN)

Это значит, что, например, вычисление при разбиении на 64x64 в два раза сложнее, чем на 8x8.

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


Хотя искажения блоков 16x16 простираются дальше, чем у 8x8, но надпись более плавная. Поэтому меня убедили только первые два аргумента. Но мне что-то больше нравится разделение на 16x16.

Квантование

На данном этапе мы имеем кучу матриц 8x8 с коэффициентами косинусного преобразования. Пришло время избавляться от малозначащих коэффициентов. Существует более элегантное решение, чем просто обнулять последние коэффициенты, как мы делали выше. Нас не устраивает этот способ, так как необнуленные значения хранятся с избыточной точностью, а среди тех, кому не повезло, могли оказаться достаточно важные. Выход - нужно использовать матрицу квантования. Потери происходят именно на это этапе. Каждый Фурье-коэффициент делится на соответствующее число в матрице квантования. Рассмотрим на примере. Возьмем первый блок от нашего енота и произведем квантование. В спецификации JPEG приводится стандартная матрица:


Стандартная матрица соответствует 50% качеству в FastStone и IrfanView. Такая таблица была выбрана с точки зрения баланса качества и степени сжатия. Думаю, что значение для DC-коэффициента больше соседних из-за того, что DCT ненормализовано и первое значение получается больше, чем следовало бы. Высокочастотные коэффициенты огрубляются сильнее из-за их меньшей важности. Думаю, сейчас такие матрицы используются редко, так как ухудшение качества хорошо заметно. Никто не запрещает использовать свою таблицу (со значениями от 1 до 255)
При декодировании происходит обратный процесс - квантованные коэффициенты почленно умножаются на значения матрицы квантования. Но так как мы округляли значения, то не сможем точно восстановить исходные коэффициенты Фурье. Чем больше число квантования, тем больше погрешность. Таким образом, восстановленный коэффициент является лишь ближайшим кратным.
Еще пример:

И на десерт, рассмотрим качество 5% (при кодировании в Fast Stone).


При восстановлении этого блока мы получим только усредненное значение плюс вертикальный градиент (из-за сохранившегося значения -1). Зато для него хранится всего два значения: 7 и -1. C другими блоками ситуация не лучше, вот восстановленная картинка:

Кстати, насчет 100% качества. Как вы догадываетесь, в этом случае матрица квантования состоит полностью из единиц, то есть квантования не происходит. Однако, из-за округления коэффициентов до целого, мы не можем в точности восстановить исходную картинку. Например, енот сохранил 96% пикселей точно, а 4% отличались на 1/256. Разумеется, такие «искажения» невозможно заметить визуально.
А можете посмотреть матрицы квантования различных фотоаппаратов.

Кодирование

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

Пример 0 (для разминки)
Представьте такую ситуацию, что ваш знакомый забыл у вас дома листочек со списком и теперь просит продиктовать его по телефону (других способов связи нет).
Список:

• d9rg3
• wfr43gt
• wfr43gt
• d9rg3
• d9rg3
• d9rg3
• wfr43gt
• d9rg3

Как бы вы облегчили свою задачу? Особого желания мучительно диктовать все эти слова у вас нет. Но их всего два и они повторяются. Поэтому вы просто как-нибудь диктуете первые два слова и договариваетесь, что далее «d9rg3» будете называть первым словом, а «wfr43gt» - вторым. Тогда достаточно будет продиктовать: 1, 2, 2, 1, 1, 1, 2, 1.

Подобные слова мы будем обозначать как A, B, C..., и называть их символами. Причем под символом может скрываться что угодно: буква алфавита, слово или бегемот в зоопарке. Главное, что одинаковым символам соответствуют одинаковые понятия, а разным - разные. Так как наша задача - эффективное кодирование (сжатие), то будем работать с битами, так как это наименьшие единицы представления информации. Поэтому, запишем список как ABBAAABA. Вместо «первое слово» и «второе слово» можно использовать биты 0 и 1. Тогда ABBAAABA закодируется как 01100010 (8 бит = 1 байт).

Пример 1
Закодировать ABC.
3-м разным символам (A, B, C) никак нельзя сопоставить 2 возможных значений бита (0 и 1). А раз так, то можно использовать по 2 бита на символ. Например:

• A: 00
• B: 01
• C: 10

Последовательность битов, сопоставленная символу, будем называть кодом. ABC будет кодироваться так: 000110.

Пример 2
Закодировать AAAAAABC.
Использовать по 2 бита на символ A кажется немного расточительным. Что, если попробовать так:

• C: 00

Закодированная последовательность: 000000100.
Очевидно, этот вариант не подходит, так как непонятно, как декодировать первые два бита этой последовательности: как AA или как C? Использовать какой-нибудь разделитель между кодами очень расточительно, будем думать как по-другому обойти это препятствие. Итак, неудача произошла из-за того, что код C начинается с кода A. Но мы полны решимости кодировать A одним битом, пусть даже B и С будут по два. Исходя из такого пожелания, A дадим код 0. Тогда коды B и C не могут начинаться на 0. Но могут на 1:

• B: 10
• C: 11

Последовательность закодируется так: 0000001011. Попробуйте мысленно декодировать ее. Вы сможете сделать это только одним способом.
Мы выработали два требования к кодированию:

1. Чем больше вес символа, тем короче должен быть его код. И наоборот.
2. Для однозначного декодирования код символа не может начинаться с кода любого другого символа.

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

Пример 3
Рассмотрим общий случай для 4-х символов с любыми весами.

• A: pa
• B: pb
• C: pc
• D: pd

Без потери общности, положим pa ≥ pb ≥ pc ≥ pd. Существуют всего два принципиально разных по длинам кодов варианта:


Какое из них предпочтительнее? Для этого нужно вычислить получаемые длины закодированных сообщений:
W1 = 2*pa + 2*pb + 2*pc + 2*pd
W2 = pa + 2*pb + 3*pc + 3*pd
Если W1 меньше W2 (W1-W2<0), то лучше использовать первый вариант:
W1-W2 = pa - (pc+pd) < 0 => pa < pc+pd.
Если C и D вместе встречаются чаще других, то их общая вершина получает самый короткий код из одного бита. В противном случае, один бит достается символу A. Значит, объединение символов ведет себя как самостоятельный символ и имеет вес равный сумме входящих символов.
Вообще, если p - вес символа представленный долей его вхождения (от 0 до 1), то лучшая длина кода s=-log 2 p.
Рассмотрим это на простом случае (его легко представить в виде дерева). Итак, нужно закодировать 2 s символов с равными весами (1/2 s). Из-за равенства весов длины кодов будут одинаковыми. Каждому символу потребуется s бит. Значит, если вес символа 1/2 s , то его длина s. Если вес заменить заменить на p, то получим длину кода s=-log 2 p . Значит, если один символ встречается в два раза реже другого, то длина его кода будет на бит длиннее. Впрочем такой вывод легко сделать, если вспомнить, что добавление одного бита позволяет в два раза увеличить количество возможных вариантов.
И еще одно наблюдение - два символа с наименьшими весами всегда имеют наибольшие, но равные длины кодов. Более того, их биты, кроме последнего, совпадают. Если бы это было неверно, то, по крайней мере, один код можно было бы укоротить на 1 бит, не нарушая префиксности. Значит, два символа с наименьшими весами в кодовом дереве имеют общего родителя уровнем выше. Вы можете видеть это на примере С и D выше.

Пример 4
Попробуем решить следующий пример, по выводам, полученным в предыдущем примере.

1. Все символы сортируются в порядке убывания весов.
2. Два последних символа объединяются в группу. Этой группе присваивается вес, равный сумме весов этих элементов. Эта группа участвует в алгоритме наравне с символами и другими группами.

Шаги повторяются, пока не останется только одна группа. В каждой группе одному символу (или подгруппе) присваивается бит 0, а другому бит 1.
Этот алгоритм называется кодированием Хаффмана.
На иллюстрации приведен пример с 5-ю символами (A: 8, B: 6, C: 5, D: 4, E: 3). Справа указан вес символа (или группы).

Кодируем коэффициенты

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

Обратите внимание - шкала логарифмическая! Сможете объяснить причину появления скопления значений превышающих 200? Это DC-коэффициенты. Так как они сильно отличаются от остальных, то неудивительно, что их кодируют отдельно. Вот только DC:

Обратите внимание, что форма графика напоминает форму графиков из самих ранних экспериментов деления на пары и тройки пикселей
Вообще, значения DC-коэффициентов могут меняться от 0 до 2047 (точнее от -1024 до 1023, так как в JPEG производится вычитание 128 из всех исходных значений, что соответствует вычитанию 1024 из DC) и распределяться довольно равномерно с небольшими пиками. Поэтому кодирование Хаффмана здесь не очень-то поможет. А еще представьте, каким большим будет дерево кодирования! И во время декодирования придется искать в нем значения. Это очень затратно. Думаем дальше.
DC-коэффициент - усредненное значение блока 8x8. Представим градиентный переход (пусть не идеальный), который часто встречается в фотографиях. Сами DC значения будут разными, но они будут представлять арифметическую прогрессию. Значит, их разность будет более-менее постоянна. Построим гистограмму разностей:

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

То есть положительные значения прямо кодируются их двоичным представлением, а отрицательные - так же, но с заменой ведущей 1 на 0. Осталось решить, как кодировать длины. Так как их 12 возможных значений, то можно использовать 4 бита для хранения длины. Но вот тут-то как раз лучше использовать кодирование Хаффмана.

Значений с длинами 4 и 6 больше всего, поэтому им достались самые короткие коды (00 и 01).

Может возникнуть вопрос: почему на примере у значения 9 код 1111110, а не 1111111? Ведь можно смело поднять «9» на уровень выше, рядом с «0»? Дело в том, что в JPEG нельзя использовать код, состоящий только из единиц - такой код зарезервирован.
Есть еще одна особенность. Коды, полученные описанным алгоритмом Хаффмана могут не совпасть по битам с кодами в JPEG, хотя их длины будут одинаковыми. Используя алгоритм Хаффмана, получают длины кодов, а сами коды генерируются (алгоритм прост - начинают с коротких кодов и добавляют их по очереди в дерево как можно левее, сохраняя свойство префиксности). Например, для дерева выше хранится список: 0,2,3,1,1,1,1,1. И, разумеется, хранится список значений: 4,6,3,5,7,2,8,1,0,9. При декодировании коды генерируются таким же способом.

Теперь порядок. Мы разобрались как хранятся DC:
[код Хаффмана для длины DC diff (в битах)]
где DC diff = DC текущее - DC предыдущее

Смотрим AC:

Так как график очень похож на график для разностей DC, то принцип тот же: [код Хаффмана для длины AC (в битах)]. Но не совсем! Так как на графике шкала логарифмическая, то не сразу заметно, что нулевых значений примерно в 10 раз больше, чем значения 2 - следующего по частоте. Это понятно - не все пережили квантование. Вернемся к матрице значений, полученной на этапе квантования (используя матрицу квантования FastStone, 90%).

Так как встречается много групп подряд идущих нулей, то появляется идея - записывать только количество нулей в группе. Такой алгоритм сжатия называется RLE (Run-length encoding, кодирование повторами). Осталось выяснить направление обхода «подряд идущих» - кто за кем? Выписать слева направо и сверху вниз - не очень эффективно, так как ненулевые коэффициенты концентрируются около левого верхнего угла, а чем ближе к правому нижнему - тем больше нулей.

Поэтому, в JPEG используется порядок, называемый «Zig-zag», он показан на левом рисунке. Такой способ хорошо выделяет группы нулей. На правом рисунке - альтернативный способ обхода, не относящийся к JPEG, зато с любопытным названием (пруф). Он может использоваться в MPEG при сжатии видео с чересстрочной разверткой. Выбор алгоритма обхода не влияет на качество изображения, но может увеличить количество кодируемых групп нулей, что в итоге может отразиться на размере файла.
Модифицируем нашу запись. Для каждого ненулевого AC - коэффициента:
[Количество нулей перед AC][код Хаффмана для длины AC (в битах)]
Думаю, что вы сразу скажете - количество нулей тоже отлично закодируется Хаффманом! Это очень близкий и неплохой ответ. Но можно немного оптимизировать. Представьте, что имеем некоторый коэффициент AC, перед которым было 7 нулей (разумеется, если выписывать в зигзагообразном порядке). Эти нули - дух значений, которые не выдержали квантования. Скорее всего, наш коэффициент тоже сильно потрепало и он стал маленьким, а, значит, его длина - короткой. Значит, количество нулей перед AC и длина AC - зависимые величины. Поэтому записывают так:
[код Хаффмана для (Количество нулей перед AC, длина AC (в битах)]
Алгоритм кодирования остается тем же: те пары (количество нулей перед AC, длина AC), которые встречаются часто, получат короткие коды и наоборот.

Строим гистограмму зависимости количества по этим парам и дерево Хаффмана.


Длинный «горный хребет» подтверждает наше предположение.

Особенности реализации в JPEG:
Такая пара занимает 1 байт: 4 бита на количество нулей и 4 бита на длину AC. 4 бита - это значения от 0 до 15. Для длины AC хватит с избытком, но ведь нулей может быть больше 15? Тогда используется больше пар. Например, для 20 нулей: (15, 0)(5, AC). То есть, 16-й ноль кодируется как ненулевой коэффициент. Так как ближе к концу блока всегда полно нулей, то после последнего ненулевого коэффициента используется пара (0,0). Если она встретится при декодировании, значит оставшиеся значения равны 0.

Выяснили, что каждый блок закодирован хранится в файле так:
[код Хаффмана для длины DC diff ]
[код Хаффмана для (количество нулей перед AC 1 , длина AC 1 ]
…
[код Хаффмана для (количество нулей перед AC n , длина AC n ]
Где AC i - ненулевые AC коэффициенты.

Цветное изображение

Способ представления цветного изображения зависит от выбранной цветовой модели. Простое решение - использовать RGB и кодировать каждый цветовой канал изображения по отдельности. Тогда кодирование не будет отличаться от кодирования серого изображения, только работы в 3 раза больше. Но сжатие изображения можно увеличить, если вспомнить, что глаз более чувствительнее к изменению яркости, чем цвета. Это значит, что цвет можно хранить с бОльшими потерями, чем яркость. У RGB нет отдельного канала яркости. Она зависит от суммы значений каждого канала. Поэтому, RGB-куб (это представление всех возможных значений) просто «ставят» на диагональ - чем выше, тем ярче. Но на этом не ограничиваются - куб немного поджимают с боков, и получается скорее параллелепипед, но это лишь для учета особенностей глаза. Например, он более восприимчив к зеленому, чем синему. Так появилась модель YCbCr.

(Изображение с Intel.com)
Y - компонента яркости, Cb и Cr являются синей и красной цветоразностными компонентами. Поэтому, если хотят сильнее сжать изображение, то RGB переводят в YCbCr, и каналы Cb и Cr прореживают. То есть разбивают на небольшие блоки, например 2x2, 4x2, 1x2, и усредняют все значения одного блока. Или, другими словами, уменьшают размер изображения для этого канала в 2 или 4 раза по вертикали и/или горизонтали.

Каждый блок 8x8 кодируется (DCT + Хаффман), и закодированные последовательности записываются в таком порядке:

Любопытно, что спецификация JPEG не ограничивает в выборе модели, то есть реализация кодировщика может как угодно разделить изображение по цветовым компонентам (каналам) и каждый будет сохранен по отдельности. Мне известно об использовании Grayscale (1 канал), YCbCr (3), RGB (3), YCbCrK (4), CMYK (4). Первые три поддерживаются почти всеми, а вот с последними 4-канальными бывают проблемы. FastStone, GIMP поддерживают их корректно, а штатные программы Windows, paint.net корректно извлекают всю информацию, но потом выбрасывают 4 черный канал, поэтому ( сказал, что не выбрасывают, читайте его комментарии) показывают более светлое изображение. Слева - классический YCbCr JPEG, справа CMYK JPEG:



Если они различаются по цветам, или видна только одна картинка, то, скорее всего, у вас IE (любой версии) (UPD. в комментариях говорят «или Safari»). Можете попробовать открыть статью в разных браузерах.

И еще кое-что

В двух словах о дополнительных возможностях.

Progressive mode

Разложим полученные таблицы коэффициентов DCT на сумму таблиц (примерно так (DC, -19, -22, 2, 1) = (DC, 0, 0, 0, 0) + (0, -20, -20, 0, 0) + (0, 1, -2, 2, 1)). Сначала закодируем все первые слагаемые (как мы уже научились: Хаффман и обход зигзагом), затем вторые и т. д. Такой трюк полезен при медленном интернете, так как сперва загружаются только DC коэффициенты, по которым строится грубая картинка c «пикселями» 8x8. Затем округленные AC коэффициенты, позволяющие уточнить рисунок. Затем грубые поправки к ним, затем более точные. Ну и так далее. Коэффициенты округляются, так как на ранних этапах загрузки точность не столь важна, зато округление положительно сказывается на длине кодов, так как для каждого этапа используется своя таблица Хаффмана.

Lossless mode

Сжатие без потерь. DCT нет. Используется предсказание 4-й точки по трем соседним. Ошибки предсказания кодируются Хаффманом. По-моему, используется чуть чаще, чем никогда.

Hierarhical mode

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

Арифметическое кодирование

Алгоритм Хаффмана создает оптимальные коды по весу символов, но это верно только для фиксированного соответствия символов с кодами. Арифметическое не имеет такой жесткой привязки, что позволяет использовать коды как бы с дробным числом бит. Утверждается, что оно уменьшает размер файла в среднем на 10% по сравнению с Хаффманом. Не распространено из-за проблем с патентом, поддерживается не всеми.

Я надеюсь, что теперь вам понятен алгоритм JPEG интуитивно. Спасибо за прочтение!

UPD
предложил указать использованное ПО. С удовольствием сообщаю, что все доступны и бесплатны:

• Python + NumPy + Matplotlib + PIL(Pillow) . Основной инструмент. Нашелся по выдаче «Matlab free alternative». Рекомендую! Даже если вам не знаком Python, то уже через пару часов научитесь производить расчеты и строить красивые графики.
• JpegSnoop . Показывает подробную информацию о jpeg-файле.
• yEd . Редактор графов.
• Inkscape . Делал в нем иллюстрации, такие как пример алгоритма Хаффмана. Прочитал несколько уроков, оказалось очень здорово.
• Daum Equation Editor . Искал визуальный редактор формул, так как с Latex-ом не очень дружу. Daum Equation - плагин к Хрому, мне показался очень удобен. Помимо мышкотыкания, можно редактировать Latex.
• FastStone . Думаю, его представлять не надо.
• PicPick . Бесплатная альтернатива SnagIt. Сидит в трее, скриншотит что скажут куда скажут. Плюс всякие плюшки, типа линейки, пипетки, угломера и пр.

Теги:

• jpeg
• dct
• dft
• фурье
• хаффман

Добавить метки

JPEG-сжатие цифрового изображения

Одним из наиболее полных и популярных стандартов сжатия изображений является стандарт JPEG.

Сам процесс сжатия состоит из трех последовательных шагов:

а) Вычисление дискретного косинусного преобразования (ДКП) для матриц 8*8-блоков, полученных после стандартного разбиения матрицы ЦИ;

б) квантование коэффициентов ДКП;

в) кодирование неравномерным кодом.

Сначала ЦИ разбивается на отдельные блоки размером 8*8 элементов, которые обрабатываются последовательно слева направо и сверху вниз. Обработка каждого блока начинается со сдвига по яркости значений всех его 64 элементов, что достигается вычитанием величины , где - максимальное число уровней яркости. Затем вычисляется двумерное ДКП элементов блока. Полученные значения коэффициентов квантуются в соответствии с формулой:

,

где - результат квантования значения коэффициента ДКП , а - соответствующий элемент матрицы коэффициентов квантования:

.

(Необходимо отметить, что перед тем, как квантованные коэффициенты ДКП могут быть подвергнуты обратному ДКП для восстановления блока изображения, они должны быть умножены на :

. (2.5)

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

Отквантованные значения коэффициентов переупорядочиваются зигзаг-преобразованием согласно:

где показана очередность, в которой выбираются коэффициенты. Результатом является одномерная последовательность квантованных коэффициентов.

Одномерный массив, полученный после зигзаг-преобразования, упорядочивается по возрастанию пространственной частоты, при этом, как правило, возникают длинные последовательности нулей, что эффективно используется процедурой JPEG-кодирования. Рекомендованная JPEG матрица квантования имеет следующий вид:

Пример . Последовательное кодирование и декодирование JPEG . Рассмотрим сжатие и восстановление следующего блока 8*8 элементов согласно стандарту последовательного кодирования JPEG:

Исходные пиксели могут иметь 256 или 2 8 уровней яркости, так что процесс кодирования начинается со сдвига диапазона значений – вычитания из значений пикселей величины 2 7 или 128. В результате получается массив:

который после прямого ДКП будет иметь вид:

Если для квантования полученных данных используется приведенная выше матрица квантования, то после квантования коэффициенты примут вид:

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

(-26 -31 -3 -2 -6 2 -4 1 -4 1 1 5 0 2 0 0 -1 2 0 0 0 0 0 -1 -1 КБ)

Кодовое слово КБ означает конец блока, указывает на то, что все оставшиеся коэффициенты в переупорядоченной последовательности равны 0. Для кодирования полученного массива используются стандартные коды Хаффмана, преобразующие массив в непрерывный поток битов.

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

Полностью восстановленный блок получается после выполнения обратного ДКП полученного массива:

и обратного сдвига диапазона значений на +2 7 =+128. В результате получаем:

Все отличия значений элементов исходного и восстановленного блоков возникают вследствие самой природы сжатия с потерями, являющегося сутью JPEG процедур сжатия и восстановления. В данном примере ошибки восстановления находятся в диапазоне от -14 до 11 и распределены следующим образом:

Характерные особенности сингулярных чисел блоков матрицы цифрового изображения при JPEG-сжатии. Пусть исходное ЦИ в градациях серого, хранящееся в некотором формате без потерь, например, в формате TIF, матрица которого имеет размеры , разбивается стандартным образом на блоки . Если для каждого блока ЦИ определить множество всех СНЧ (сингулярный спектр), то оказывается, что в среднем лишь 2.40% от общего числа блоков (ОЧБ) имеют нулевые СНЧ.

Данный факт не случаен. Ранг любой матрицы определяется количеством ее ненулевых СНЧ, а значит наличие нулей в сингулярном спектре будет говорить о том, что число ее линейно независимых строк (столбцов) меньше размера. Однако, для произвольного реального ЦИ, даже с учетом коррелированности значений яркости пикселей, вероятность того, что строки (столбцы) очередного блока окажутся линейно зависимыми, невелика.

Квантование коэффициентов DCT, которое происходит в процессе сохранения ЦИ в формате JPEG (с потерями), является необратимой процедурой и приводит к некоторым особенностям возмущений СНЧ блоков.

Пусть исходное ЦИ подверглось JPEG-сжатию. Проведем для него операцию частичного восстановления (ЧВ), которая включает в себя: 1) энтропийное декодирование; 2) умножение полученных коэффициентов на соответствующие элементы массива нормализации (матрицы квантования); 3) применение обратного DCT, но без последующего округления.

У полученной матрицы практически все блоки содержат нулевые СНЧ, причем таких значений в блоках будет достаточно много (табл.2.1). Такая ситуация закономерна. После квантования и округления коэффициентов DCT блоков многие из них, отвечающие высоким и средним частотам, обнулятся, оставаясь нулями после ЧВ, что, учитывая соответствие между коэффициентами дискретного преобразования Фурье и сингулярными тройками матрицы изображения, где - СНЧ и отвечающие ему левый и правый СНВ соответственно, приведет к обнулению наименьших (а возможно и средних по величине) СНЧ матриц блоков.

Табл.2.1. Результаты сингулярного разложения блоков частично восстановленных изображений

	ОЧБ		Кол-во блоков, у кот-х нулевых СНЧ больше 2-х, по отношению к ОЧБ (в %)
m=8	m=7	m=6	m=5	m=4	m=3	m=2	m=1	m=0
POUT
CAMERAMAN
TIRE
MOON											99.8
CELL

Заметим, что, чем меньше нулевых СНЧ в рассматриваемом блоке, тем больше линий контура он содержит. Действительно, наличие контуров в блоке говорит о значительной высокочастотной составляющей в сигнале, отвечающем этому блоку. Тогда коэффициенты DCT, соответствующие высоким и средним частотам, будут сравнительно большими и могут остаться ненулевыми после квантования и ЧВ, а значит внесут свой вклад не только в максимальные СНЧ.

Для наглядного представления справедливости вышесказанного рассмотрим изображение СELL.TIF (рис.2.5(а)). На рис.2.5(б) представлена матрица нулевых СНЧ блоков (МНСЧБ) размерности ЧВ-изображения, каждый элемент которой равен количеству нулевых СНЧ в соответствующем блоке. На рисунке выделены элементы, имеющие наименьшие значения, что позволяет наглядно увидеть соответствие между контурами исходного ЦИ и блоками, содержащими наименьшее количество нулевых СНЧ.

Пусть исходное изображение, подвергшееся JPEG-сжатию, восстанавливается полностью. Это означает, что после ЧВ все значения яркости пикселей округляются до целых и вводятся в диапазон . Это действие возмутит матрицу изображения, полученную после ЧВ, определенным образом изменится количество нулевых СНЧ в блоках (табл.2.2). Там, где после ЧВ не было элементов, значительно меньших 0 или больших 255, возмущение матрицы будет небольшим. В соответствии с соотношением

, (2.6)

имеющим место для произвольной матрицы, где - СНЧ исходной и возмущенной матриц соответственно, - матрица возмущений блока, - спектральная матричная норма, СНЧ являются нечувствительными к возмущающим воздействиям. Если некоторые из нулевых СНЧ блоков матрицы ЧВ-изображения станут ненулями после полного восстановления (ПВ), то их значения будут сравнимы с погрешностью округления, что не характерно для блоков исходного ЦИ.

Рис.2.5. Исходное изображение СELL.TIF (а); МНСЧБ после ЧВ (б); МНСЧБ после полного восстановления (в)

Наиболее заметным различие между совокупным исходным изображением и полностью восстановленным после JPEG-сжатия будет при сравнении их МНСЧБ. Типичная картина представлена на рис.2.5(в), при этом МНСЧБ для CELL.TIF имела только нулевые значения.

Таблица 2.2. Результаты сингулярного разложения блоков полностью восстановленных изображений

Изображение в формате без потерь (TIF)	ОЧБ	Количества блоков, имеющих нулевых СНЧ	Кол-во блоков, у которых нулевых СЧ больше двух, по отношению к ОЧБ (%)
m=8	m=7	m=6	m=5	m=4	m=3	m=2	m=1	m=0
POUT
CAMERAMAN
TIRE
MOON
CELL

Вопросы

1. Что означает сжатие даннях? Что такое избыточность даннях?
2. Основные виды избыточности данных.
3. Как реализуется сжатие посредством квантования?
4. Что такое малоранговая оппроксимация изображения? Как реализуется сжатие посредством использования малоранговых аппроксимаций изображения?
5. Что такое сингулярное разложение матрицы?
6. Что такое спектральное разложение матрицы?
7. Соответствие между параметрами цифрового изображения в пространственной и частотной областях.
8. Основные шаги JPEG-сжатие цифрового изображения. Матрицы квантования.
9. Характерные особенности сингулярных чисел блоков матрицы цифрового изображения при JPEG-сжатии.
10. Частичное и полное восстановление цифрового изображения после сжатия.

Литература

1. Кобозева А.А. Анализ информационной безопасности / А.А.Кобозева, В.А.Хорошко. – К.: Изд. ГУИКТ, 2009. – 251 с.
2. Деммель Дж. Вычислительная линейная алгебра / Дж.Деммель; пер.с англ. Х.Д.Икрамова. - М.: Мир, 2001. - 430 с.
3. Бахвалов Н.С. Численные методы / Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. - М.: БИНОМ. Лаборатория знаний, 2006. - 636 с.
4. Гонсалес Р. Цифровая обработка изображений / Р.Гонсалес, Р.Вудс; пер. с англ. под ред. П.А.Чочиа. - М.: Техносфера, 2005. - 1072 с.
5. Каханер Д. Численные методы и программное обеспечение / Д.Каханер, К.Моулер, С.Нэш; пер. с англ. Х.Д.Икрамова. - М.: Мир, 2001. - 575 с.
6. Гантмахер Ф.Р. Теория матриц / Ф.Р.Гантмахер. - М.: Наука, 1988. - 552 с.
7. Воеводин В.В. Вычислительные основы линейной алгебры / В.В.Воеводин. - М.: Наука. Гл.ред.физ.-мат.лит., 1977. - 304 с.

Область применения

Алгоритм JPEG в наибольшей степени пригоден для сжатия фотографий и картин, содержащих реалистичные сцены с плавными переходами яркости и цвета. Наибольшее распространение JPEG получил в цифровой фотографии и для хранения и передачи изображений с использованием сети Интернет .

С другой стороны, JPEG малопригоден для сжатия чертежей, текстовой и знаковой графики, где резкий контраст между соседними пикселами приводит к появлению заметных артефактов . Такие изображения целесообразно сохранять в форматах без потерь, таких как TIFF , GIF или PNG .

JPEG (как и другие методы искажающего сжатия) не подходит для сжатия изображений при многоступенчатой обработке, так как искажения в изображения будут вноситься каждый раз при сохранении промежуточных результатов обработки.

JPEG не должен использоваться и в тех случаях, когда недопустимы даже минимальные потери, например, при сжатии астрономических или медицинских изображений. В таких случаях может быть рекомендован предусмотренный стандартом JPEG режим сжатия Lossless JPEG (который, однако, не поддерживается большинством популярных кодеков) или стандарт сжатия JPEG-LS .

Сжатие

При сжатии изображение преобразуется из цветового пространства RGB в YCbCr (YUV). Следует отметить, что стандарт JPEG (ISO/IEC 10918-1) никак не регламентирует выбор именно YCbCr, допуская и другие виды преобразования (например, с числом компонентов , отличным от трёх), и сжатие без преобразования (непосредственно в RGB), однако спецификация JFIF (JPEG File Interchange Format, предложенная в 1991 году специалистами компании C-Cube Microsystems, и ставшая в настоящее время стандартом де-факто) предполагает использование преобразования RGB->YCbCr.

После преобразования RGB->YCbCr для каналов изображения Cb и Cr, отвечающих за цвет, может выполняться «прореживание» (subsampling ), которое заключается в том, что каждому блоку из 4 пикселов (2х2) яркостного канала Y ставятся в соответствие усреднённые значения Cb и Cr (схема прореживания «4:2:0» ). При этом для каждого блока 2х2 вместо 12 значений (4 Y, 4 Cb и 4 Cr) используется всего 6 (4 Y и по одному усреднённому Cb и Cr). Если к качеству восстановленного после сжатия изображения предъявляются повышенные требования, прореживание может выполняться лишь в каком-то одном направлении - по вертикали (схема «4:4:0») или по горизонтали («4:2:2»), или не выполняться вовсе («4:4:4»).

Стандарт допускает также прореживание с усреднением Cb и Cr не для блока 2х2, а для четырёх расположенных последовательно (по вертикали или по горизонтали) пикселов, то есть для блоков 1х4, 4х1 (схема «4:1:1»), а также 2х4 и 4х2 (схема «4:1:0»). Допускается также использование различных типов прореживания для Cb и Cr, но на практике такие схемы применяются исключительно редко.

Далее яркостный компонент Y и отвечающие за цвет компоненты Cb и Cr разбиваются на блоки 8х8 пикселов. Каждый такой блок подвергается дискретному косинусному преобразованию (ДКП) . Полученные коэффициенты ДКП квантуются (для Y, Cb и Cr в общем случае используются разные матрицы квантования) и пакуются с использованием кодирования серий и кодов Хаффмана . Стандарт JPEG допускает также использование значительно более эффективного арифметического кодирования , однако из-за патентных ограничений (патент на описанный в стандарте JPEG арифметический QM-кодер принадлежит IBM) на практике оно используется редко. В популярную библиотеку libjpeg последних версий включена поддержка арифметического кодирования, но с просмотром сжатых с использованием этого метода изображений могут возникнуть проблемы, поскольку многие программы просмотра не поддерживают их декодирование.

Матрицы, используемые для квантования коэффициентов ДКП, хранятся в заголовочной части JPEG-файла. Обычно они строятся так, что высокочастотные коэффициенты подвергаются более сильному квантованию, чем низкочастотные. Это приводит к огрублению мелких деталей на изображении. Чем выше степень сжатия, тем более сильному квантованию подвергаются все коэффициенты.

При сохранении изображения в JPEG-файле указывается параметр качества, задаваемый в некоторых условных единицах, например, от 1 до 100 или от 1 до 10. Большее число обычно соответствует лучшему качеству (и большему размеру сжатого файла). Однако даже при использовании наивысшего качества (соответствующего матрице квантования, состоящей из одних только единиц) восстановленное изображение не будет в точности совпадать с исходным, что связано как с конечной точностью выполнения ДКП, так и с необходимостью округления значений Y, Cb, Cr и коэффициентов ДКП до ближайшего целого. Режим сжатия Lossless JPEG, не использующий ДКП, обеспечивает точное совпадение восстановленного и исходного изображений, однако его малая эффективность (коэффициент сжатия редко превышает 2) и отсутствие поддержки со стороны разработчиков программного обеспечения не способствовали популярности Lossless JPEG.

Разновидности схем сжатия JPEG

Стандарт JPEG предусматривает два основных способа представления кодируемых данных.

Наиболее распространённым, поддерживаемым большинством доступных кодеков , является последовательное (sequential JPEG) представление данных, предполагающее последовательный обход кодируемого изображения поблочно слева направо, сверху вниз. Над каждым кодируемым блоком изображения осуществляются описанные выше операции, а результаты кодирования помещаются в выходной поток в виде единственного «скана», то есть массива кодированных данных, соответствующего последовательно пройденному («просканированному») изображению. Основной или «базовый» (baseline) режим кодирования допускает только такое представление. Расширенный (extended) режим наряду с последовательным допускает также прогрессивное (progressive JPEG) представление данных.

В случае progressive JPEG сжатые данные записываются в выходной поток в виде набора сканов, каждый из которых описывает изображение полностью с всё большей степенью детализации. Это достигается либо путём записи в каждый скан не полного набора коэффициентов ДКП, а лишь какой-то их части: сначала - низкочастотных, в следующих сканах - высокочастотных (метод «spectral selection» то есть спектральных выборок), либо путём последовательного, от скана к скану, уточнения коэффициентов ДКП (метод «successive approximation», то есть последовательных приближений). Такое прогрессивное представление данных оказывается особенно полезным при передаче сжатых изображений с использованием низкоскоростных каналов связи, поскольку позволяет получить представление обо всём изображении уже после передачи незначительной части JPEG-файла.

Обе описанные схемы (и sequential, и progressive JPEG) базируются на ДКП и принципиально не позволяют получить восстановленное изображение абсолютно идентичным исходному. Однако стандарт допускает также сжатие, не использующее ДКП, а построенное на основе линейного предсказателя (lossless, то есть «без потерь», JPEG), гарантирующее полное, бит-в-бит, совпадение исходного и восстановленного изображений. При этом коэффициент сжатия для фотографических изображений редко достигает 2, но гарантированное отсутствие искажений в некоторых случаях оказывается востребованным. Заметно большие степени сжатия могут быть получены при использовании не имеющего, несмотря на сходство в названиях, непосредственного отношения к стандарту JPEG ISO/IEC 10918-1 (ITU T.81 Recommendation) метода сжатия JPEG-LS , описываемого стандартом ISO/IEC 14495-1 (ITU T.87 Recommendation).

Синтаксис и структура

Файл JPEG содержит последовательность маркеров , каждый из которых начинается с байта 0xFF, свидетельствующего о начале маркера, и байта-идентификатора. Некоторые маркеры состоят только из этой пары байтов, другие же содержат дополнительные данные, состоящие из двухбайтового поля с длиной информационной части маркера (включая длину этого поля, но за вычетом двух байтов начала маркера то есть 0xFF и идентификатора) и собственно данных. Такая структура файла позволяет быстро отыскать маркер с необходимыми данными (например, с длиной строки, числом строк и числом цветовых компонентов сжатого изображения).

Основные маркеры JPEG

Маркер	Байты	Длина	Назначение	Комментарии
SOI	0xFFD8	нет	Начало изображения
SOF0	0xFFC0	переменный размер	Начало фрейма (базовый, ДКП)	Показывает что изображение кодировалось в базовом режиме с использованием ДКП и кода Хаффмана. Маркер содержит число строк и длину строки изображения (двухбайтовые поля со смещением соответственно 5 и 7 относительно начала маркера), количество компонентов (байтовое поле со смещением 8 относительно начала маркера), число бит на компонент (байтовое поле со смещением 4 относительно начала маркера), а также соотношение компонентов (например, 4:2:0).
SOF1	0xFFC1	переменный размер	Начало фрейма (расширенный, ДКП, код Хаффмана)	Показывает что изображение кодировалось в расширенном (extended) режиме с использованием ДКП и кода Хаффмана. Маркер содержит число строк и длину строки изображения, количество компонентов, число бит на компонент, а также соотношение компонентов (например, 4:2:0).
SOF2	0xFFC2	переменный размер	Начало фрейма (прогрессивный, ДКП, код Хаффмана)	Показывает что изображение кодировалось в прогрессивном режиме с использованием ДКП и кода Хаффмана. Маркер содержит число строк и длину строки изображения, количество компонентов, число бит на компонент, а также соотношение компонентов (например, 4:2:0).
DHT	0xFFC4	переменный размер	Содержит таблицы Хаффмана	Задает одну или более таблиц Хаффмана.
DQT	0xFFDB	переменный размер	Содержит таблицы квантования	Задает одну или более таблиц квантования.
DRI	0xFFDD	4 байта	Указывает интервал повторений	Задает интервал между маркерами RST n в макроблоках.
SOS	0xFFDA	переменный размер	Начало сканирования	Начало первого или очередного скана изображения с направлением обхода слева направо сверху вниз. Если использовался базовый режим кодирования, используется один скан. При использовании прогрессивных режимов используется несколько сканов. Маркер SOS является разделяющим между информативной (заголовком) и закодированной (собственно сжатыми данными) частями изображения.
RSTn	0xFFDn	нет	Перезапуск	Вставляется в каждом r макроблоке, где r - интервал перезапуска DRI маркера. Не используется при отсутствии DRI маркера. n , младшие 3 бита маркера кода, циклы от 0 до 7.
APPn	0xFFEn	переменный размер	Задаётся приложением	Например, в EXIF JPEG-файла используется маркер APP1 для хранения метаданных, расположеных в структуре, основанной на TIFF .
COM	0xFFFE	переменный размер	Комментарий	Содержит текст комментария.
EOI	0xFFD9	нет	Конец закодированной части изображения.

Достоинства и недостатки

К недостаткам сжатия по стандарту JPEG следует отнести появление на восстановленных изображениях при высоких степенях сжатия характерных артефактов : изображение рассыпается на блоки размером 8x8 пикселов (этот эффект особенно заметен на областях изображения с плавными изменениями яркости), в областях с высокой пространственной частотой (например, на контрастных контурах и границах изображения) возникают артефакты в виде шумовых ореолов. Следует отметить, что стандарт JPEG (ISO/IEC 10918-1, Annex K, п. K.8) предусматривает использование специальных фильтров для подавления блоковых артефактов, но на практике подобные фильтры, несмотря на их высокую эффективность, практически не используются. Однако, несмотря на недостатки, JPEG получил очень широкое распространение из-за достаточно высокой (относительно существовавших во время его появления альтернатив) степени сжатия, поддержке сжатия полноцветных изображений и относительно невысокой вычислительной сложности .

Производительность сжатия по стандарту JPEG

Для ускорения процесса сжатия по стандарту JPEG традиционно используется распараллеливание вычислений, в частности - при вычислении ДКП. Исторически одна из первых попыток ускорить процесс сжатия с использованием такого подхода описана в опубликованной в 1993 г. статье Касперовича и Бабкина , в которой предлагалась оригинальная аппроксимация ДКП, делающая возможным эффективное распараллеливание вычислений с использованием 32-разрядных регистров общего назначения процессоров Intel 80386. Появившиеся позже более производительные вычислительные схемы использовали SIMD -расширения набора инструкций процессоров архитектуры x86. Значительно лучших результатов позволяют добиться схемы, использующие вычислительные возможности графических ускорителей (технологии NVIDIA CUDA и AMD FireStream) для организации параллельных вычислений не только ДКП, но и других этапов сжатия JPEG (преобразование цветовых пространств, run-level, статистическое кодирование и т.п.), причём для каждого блока 8х8 кодируемого или декодируемого изображения. В статье была впервые [источник? ] представлена реализация распараллеливания всех стадий алгоритма JPEG по технологии CUDA, что значительно ускорило производительность сжатия и декодирования по стандарту JPEG.

В 2010 году ученые из проекта PLANETS поместили инструкции по чтению формата JPEG в специальную капсулу, которую поместили в специальный бункер в швейцарских Альпах. Сделано это было с целью сохранения для потомков информации о популярных в начале XXI века цифровых форматах.

См. также

Примечания

Ссылки

• Спецификация JFIF 1.02 (текстовый файл)
• Оптимизация JPEG. Часть 1 , Часть 2 , Часть 3 .

«Реализация алгоритмов

JPEG и JPEG2000»

Выполнил:

студент группы 819

Угаров Дмитрий

Принципы работы алгоритмов JPEG и JPEG2000

1. Алгоритм JPEG

JPEG (англ. Joint Photographic Experts Group - объединённая группа экспертов в области фотографии) - является широко используемым методом сжатия фотоизображений. Формат файла, который содержит сжатые данные обычно также называют именем JPEG; наиболее распространённые расширения для таких файлов.jpeg, .jfif, .jpg, .JPG, или.JPE. Однако из них.jpg самое популярное расширение на всех платформах.

Алгоритм JPEG является алгоритмом сжатия с потерей качества .

Область применения

Формат является форматом сжатия с потерями, поэтому некорректно считать что JPEG хранит данные как 8 бит на канал (24 бит на пиксель). С другой стороны , так как данные, подвергающиеся компрессии по формату JPEG и декомпрессированные данные обычно представляются в формате 8 бит на канал, иногда используется эта терминология. Поддерживается также сжатие чёрно-белых полутоновых изображений.

При сохранении JPEG-файла можно указать степень качества, а значит и степень сжатия, которую обычно задают в некоторых условных единицах, например, от 1 до 100 или от 1 до 10. Большее число соответствует лучшему качеству, но при этом увеличивается размер файла. Обыкновенно, разница в качестве между 90 и 100 на глаз уже практически не воспринимается. Следует помнить , что побитно восстановленное изображение всегда отличается от оригинала. Распространённым заблуждением является мнение о том, что качество JPEG тождественно доле сохраняемой информации.

Этапы кодирования

Процесс сжатия по схеме JPEG включает ряд этапов:

1. Преобразование изображения в оптимальное цветовое пространство;

В случае применения цветового пространства яркость/цветность (YCbCr) достигается лучшая степень сжатия. На данном этапе кодирования с помощью соответствующих соотношений цветовая модель RGB преобразуется в YCbCr:

Y = 0.299*R + 0.587*G + 0.114*B

Cb = - 0.1687*R – 0.3313*G + 0.5*B

Cr = 0.5*R – 0.4187*G – 0.0813*B.
Во время декодирования можно использовать соответствующее обратное преобразование:
R = Y + 1.402*Cr

G = Y – 0.34414*Cb – 0.71414*Cr

B = Y + 1.772*Cb.
Примечание, связывающее Y,Cb,Cr в человеческой визуальной системе:

Глаз, особенно сетчатка, имеет как визуальные анализаторы два типа ячеек: ячейки для ночного видения, воспринимающие только оттенки серого (от ярко-белого до темно-черного) и ячейки дневного видения, которые воспринимают цветовой оттенок. Первые ячейки , дающие цвет RGB, обнаруживают уровень яркости, подобный величине Y. Другие ячейки, ответственные за восприятие цветового оттенка, - определяют величину, связанную с цветоразностью.

2. Субдискретизация компонентов цветности усреднением групп пикселей;

Большая часть визуальной информации, к которой наиболее чувствительный глаза человека , состоит из высокочастотных, полутоновых компонентов яркости (Y) цветового пространства YCbCr. Две другие составляющие цветности (Cb и Cr) содержат высокочастотную цветовую информацию, к которой глаз человека менее чувствителен. Поэтому определенная ее часть может быть отброшена и, тем самым, можно уменьшить количество учитываемых пикселей для каналов цветности.

1)тип 4:2:0 (когда изображение разбивается на квадраты 2х2 пикселей и в каждом из них все пиксели получают одинаковые значения каналов Cb и Cr, а яркость Y у остается у каждого своя)

2) тип 4:2:2 (объединение по компонентам цветности происходит только по горизонтали в группах по два пикселя).

3)тип 4: 4: 4 подразумевает, что каждому пикселю в каждой строке соответствует собственное уникальное значение компонентов Y, Cb и Cr. (рис.1 а)

4) тип 4:2:2. Выполнив субдискретизацию сигнала цветности с коэффициентом 2 по горизонтали, мы получим из потока 4: 4: 4 YCbCr поток 4: 2: 2 YCbCr. Запись «4: 2: 2» означает , что в отдельно взятой строке на 2 значения цветности приходятся 4 значения яркости (см. рис.1 б). Сигнал 4: 2: 2 YCbCr очень немного проигрывает по качеству изображения сигналу 4: 4: 4 YCbCr, зато требуемая ширина полосы сокращается на 33% от исходной.

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

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

DCT непосредственно применяемый к блоку (в нашем случае 8х8 пикселей) изображения будет выглядеть так:

где х, y - пространственные координаты пикселя (0..7) ,

f(x,y) - значения пикселей исходного макроблока (допустим, яркость)

u,v - координаты пикселя в частотном представлении (0..7)

w(u) =1/SQRT(2) при u=0, в остальных случаях w(u)=1 (SQRT - квадратный корень)

w(v) =1/SQRT(2) при v=0, в остальных случаях w(v)=1

Или в матричной форме:

4. Квантование каждого блока коэффициентов ДКП с применением весовых функций , оптимизированных с учетом визуального восприятия человеком;

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

5. Этап Вторичного Сжатия

Заключительной стадией работы кодера JPEG является кодирование полученной матрицы.

5.1 Зигзагообразная перестановка 64 DCT коэффициентов

Так, после того, как мы выполнили DCT-преобразование над блоком величин 8x8, у нас есть новый блок 8x8. Затем, этот блок 8x8 просматривается по зигзагу подобно этому:

(Числа в блоке 8x8 указывают порядок , в котором мы просматриваем 2-мерную матрицу 8x8)

0, 1, 5, 6,14,15,27,28,

2, 4, 7,13,16,26,29,42,

3, 8,12,17,25,30,41,43,

9,11,18,24,31,40,44,53,

10,19,23,32,39,45,52,54,

20,22,33,38,46,51,55,60,

21,34,37,47,50,56,59,61,

35,36,48,49,57,58,62,63

Как Вы видите, сначала - верхний левый угол (0,0), затем величина в (0,1), затем (1,0), затем (2,0), (1,1), (0,2), (0,3), (1,2), (2,1), (3,0) и т.п.

После того, как мы прошли по зигзагу матрицу 8x8, мы имеем теперь вектор с 64 коэффициентами (0..63) Смысл этого зигзагообразного вектора – в том, что мы просматриваем коэффициенты 8x8 DCT в порядке повышения пространственных частот. Так, мы получаем вектор отсортированный критериями пространственной частоты: первая величина на векторе (индекс 0) соответствует самой низкой частоте в изображении – она обозначается термином DC. С увеличением индекса на векторе, мы получаем величины соответствующие высшим частотам (величина с индексом 63 соответствует амплитуде самой высокой частоте в блоке 8x8). Остальная часть коэффициентов DCT обозначается AC.

5.2 RunLength кодирование нулей (RLE)

Теперь у нас есть вектор с длинной последовательностью нулей. Мы можем использовать это, кодируя последовательные нули. ВАЖНО: Вы увидите позже почему, но здесь мы пропускаем кодировку первого коэффициента вектора (коэффициент DC), который закодирован по-другому. Рассмотрим исходный 64 вектор как 63 вектор (это - 64 вектор без первого коэффициента)

Допустим, мы имеем 57,45,0,0,0,0,23,0,-30,-16,0,0,1,0,0,0,0,0,0, только 0,...,0

Здесь - как RLC JPEG сжатие сделано для этого примера:

(0,57); (0,45); (4,23); (1,-30); (0,-16); (2,1); EOB

Как Вы видите, мы кодируем для каждой величины, отличающейся от 0 количество последовательных ПРЕДШЕСТВУЮЩИХ нулей перед величиной, затем мы добавляем величину. Другое примечание: EOB - короткая форма для Конца Блока , это - специальная кодированная величина (маркер). Если мы достигли в позиции на векторе, от которого мы имеем до конца только нули вектора, мы выделим эту позицию с EOB и завершим сжатие RLC квантованного вектора.

[Заметьте, что если квантованный вектор не оканчивается нулями (имеет последний элемент не 0), мы не будем иметь маркер EOB.]

(0,57); (0,45); (4,23); (1,-30); (0,-16); (2,1); (0,0)

Другая ОСНОВНАЯ вещь: Допустим, где-нибудь на квантованном векторе мы имеем:

57, восемнадцать нулей, 3, 0,0 ,0,0 2, тридцать-три нуля, 895, EOB

Кодирование Хаффмана JPG делает ограничение, по которому число предшествующих нулей должно кодироваться как 4-битовая величина - не может превысить 15.

Так, предшествующий пример должен быть закодирован как:

(0,57); (15,0) (2,3); (4,2); (15,0) (15,0) (1,895), (0,0)

(15,0) - специальная кодированная величина, которая указывает , что там следует за 16 последовательными нулями.

5.3 Конечный шаг - кодирование Хаффмана

Сначала ВАЖНОЕ примечание: Вместо хранения фактической величины, JPEG стандарт определяет, что мы храним минимальный размер в битах, в котором мы можем держать эту величину (это названо категория этой величины) и затем битно кодированное представление этой величины подобно этому:

7,..,-4,4,..,7 3 000,001,010,011,100,101,110,111

15,..,-8,8,..,15 4 0000,..,0111,1000,..,1111

31,..,-16,16,..,31 5 00000,..,01111,10000,..,11111

63,..,-32,32,..,63 6 .

127,..,-64,64,..,127 7 .

255,..,-128,128,..,255 8 .

511,..,-256,256,..,511 9 .

1023,..,-512,512,..,1023 10 .

2047,..,-1024,1024,..,2047 11 .

4095,..,-2048,2048,..,4095 12 .

8191,..,-4096,4096,..,8191 13 .

16383,..,-8192,8192,..,16383 14 .

32767,..,-16384,16384,..,32767 15 .

Впоследствии для предшествующего примера:

(0,57); (0,45); (4,23); (1,-30); (0,-8); (2,1); (0,0)

давайте закодируем только правую величину этих пар, кроме пар, которые являются специальными маркерами подобно (0,0) или (если мы должны иметь) (15,0)

45, аналогично , будет закодирован как (6,101101)

30 -> (5,00001)

И теперь, мы напишем снова строку пар:

(0,6), 111001; (0,6), 101101; (4,5), 10111; (1,5), 00001; (0,4), 0111; (2,1), 1; (0,0)

Пары 2 величин, заключенные в скобки, могут быть представлены в байте, так как фактически каждая из 2 величин может быть представлена в 4-битном кусочке (счетчик предшествующих нулей - всегда меньше, чем 15 и также как и категория [числа закодированные в файле JPG - в области -32767..32767]). В этом байте, старший кусочек представляет число предшествующих нулей, а младший кусочек - категорию новой величины, отличной от 0.

Конечный шаг кодировки состоит в кодировании Хаффмана этого байта, и затем записи в файле JPG , как поток из битов, кода Хаффмана этого байта, сопровождающийся битовым представлением этого числа.

Например, для байта 6 (эквивалент (0,6)) у нас есть код Хаффмана = 111000;

21 = (1,5) - 11111110110

4 = (0,4) - 1011

33 = (2,1) - 11011

0 = EOB= (0,0) - 1010

Конечный поток битов записанных в файле JPG на диск для предшествующего примера 63 коэффициентов (запомните, что мы пропустили первый коэффициент) -

111000 111001 111000 101101 1111111110011001 10111 11111110110 00001

1011 0111 11011 1 1010
Достоинства и недостатки

К недостаткам формата следует отнести то, что при сильных степенях сжатия дает знать о себе блочная структура данных, изображение «дробится на квадратики» (каждый размером 8x8 пикселей). Этот эффект особенно заметен на областях с низкой пространственной частотой (плавные переходы изображения, например, чистое небо). В областях с высокой пространственной частотой (например, контрастные границы изображения), возникают характерные «артефакты» - иррегулярная структура пикселей искаженного цвета и/или яркости. Кроме того, из изображения пропадают мелкие цветные детали. Не стоит также забывать и о том, что данный формат не поддерживает прозрачность.

Однако, несмотря на недостатки, JPEG получил очень широкое распространение из-за высокой степени сжатия, относительно существующих во время его появления альтернатив.

2. Алгоритм JPEG2000

Алгоритм JPEG-2000 разработан той же группой экспертов в области фотографии, что и JPEG. Формирование JPEG как международного стандарта было закончено в 1992 году. В 1997 стало ясно, что необходим новый, более гибкий и мощный стандарт, который и был доработан к зиме 2000 года.

Основные отличия алгоритма в JPEG 2000 от алгоритма в JPEG заключаются в следующем:

1)Лучшее качество изображения при сильной степени сжатия. Или, что то же самое , большая степень сжатия при том же качестве для высоких степеней сжатия. Фактически это означает заметное уменьшение размеров графики "Web-качества", используемой большинством сайтов.

2)Поддержка кодирования отдельных областей с лучшим качеством. Известно, что отдельные области изображения критичны для восприятия человеком (например, глаза на фотографии), в то время как качеством других можно пожертвовать (например, задний план). При "ручной" оптимизации увеличение степени сжатия проводится до тех пор, пока не будет потеряно качество в какой-то важной части изображения. Сейчас появляется возможность задать качество в критичных областях, сжав остальные области сильнее, т.е. мы получаем еще большую окончательную степень сжатия при субъективно равном качестве изображения.

3)Основной алгоритм сжатия заменен на wavelet. Помимо указанного повышения степени сжатия это позволило избавиться от 8-пиксельной блочности, возникающей при повышении степени сжатия. Кроме того, плавное проявление изображения теперь изначально заложено в стандарт (Progressive JPEG, активно применяемый в Интернет, появился много позднее JPEG).

4)Для повышения степени сжатия в алгоритме используется арифметическое сжатие. Изначально в стандарте JPEG также было заложено арифметическое сжатие, однако позднее оно было заменено менее эффективным сжатием по Хаффману, поскольку арифметическое сжатие было защищено патентами. Сейчас срок действия основного патента истек , и появилась возможность улучшить алгоритм.

5)Поддержка сжатия без потерь. Помимо привычного сжатия с потерями новый JPEG теперь будет поддерживать и сжатие без потерь. Таким образом, становится возможным использование JPEG для сжатия медицинских изображений, в полиграфии, при сохранении текста под распознавание OCR системами и т.д.

6)Поддержка сжатия однобитных (2-цветных) изображений. Для сохранения однобитных изображений (рисунки тушью, отсканированный текст и т.п.) ранее повсеместно рекомендовался формат GIF, поскольку сжатие с использованием ДКП весьма неэффективно к изображениям с резкими переходами цветов. В JPEG при сжатии 1-битная картинка приводилась к 8-битной, т.е. увеличивалась в 8 раз, после чего делалась попытка сжимать, нередко менее чем в 8 раз. Сейчас можно рекомендовать JPEG 2000 как универсальный алгоритм.

7)На уровне формата поддерживается прозрачность. Плавно накладывать фон при создании WWW страниц теперь можно будет не только в GIF, но и в JPEG 2000. Кроме того, поддерживается не только 1 бит прозрачности (пиксель прозрачен/непрозрачен), а отдельный канал , что позволит задавать плавный переход от непрозрачного изображения к прозрачному фону.

Кроме того, на уровне формата поддерживаются включение в изображение информации о копирайте, поддержка устойчивости к битовым ошибкам при передаче и широковещании, можно запрашивать для декомпрессии или обработки внешние средства (plug-ins), можно включать в изображение его описание, информацию для поиска и т.д.

Этапы кодирования

Процесс сжатия по схеме JPEG2000 включает ряд этапов:

1. Преобразование изображения в оптимальное цветовое пространство.
На данном этапе кодирования с помощью соответствующих соотношений цветовая модель RGB преобразуется в YUV:

При декомпрессии применяется соответствующее обратное преобразование:

2. Дискретное вейвлет преобразование.

Дискретное wavelet преобразование (DWT) также может быть двух видов - для случая сжатия с потерями и для сжатия без потерь.

Это преобразование в одномерном случае представляет собой скалярное произведение соответствующих коэффициентов на строку значений. Но т.к. многие коэффициенты нулевые, то прямое и обратное вейвлет преобразование можно записать следующими формулами (для преобразования крайних элементов строки используется ее расширение на 2 пикселя в каждую сторону, значения которых симметричны с значениями элементов строки относительно ее крайних пикселей):
y(2*n + 1) = x(2*n + 1) - (int)(x(2*n) + x(2*n + 2)) / 2

y(2*n) = x(2*n) + (int)(y(2*n - 1) + y(2*n + 1) + 2) / 4

и обратное

x(2*n) = y(2*n) - (int)(y(2*n - 1) + y(2*n + 1) + 2) / 4

x(2*n + 1) = y(2*n + 1) + (int)(x(2*n) + x(2*n + 2)) / 2.

3. Квантование коэффициентов.

Так же как и в алгоритме JPEG , при кодировании изображения в формат JPEG2000 используется квантование. Дискретное вейвлет преобразование, так же как и его аналог, сортирует коэффициенты по частотности. Но, в отличие от JPEG, в новом формате матрица квантования одна на все изображение.

4. Этап Вторичного Сжатия

. Как и в JPEG, в новом формате последним этапом алгоритма сжатия является кодирование без потерь. Но, в отличие от предыдущего формата, в JPEG2000 используется алгоритм арифметического сжатия.

Программная реализация

В данной работе реализованы алгоритмы JPEG и JPEG2000. В обоих алгоритмах реализовано прямое и обратное кодирование (отсутствует последний этап вторичного сжатия). Расчет JPEG происходит довольно долго (порядка 30 секунд) в связи «прямым» высчитыванием DCT. Если потребуется увеличить скорость работы , следует изначально вычислить матрицу DCT(изменения производить в классе DCT).

Перейдем к рассмотрению программы:

1. 
   После запуска выводится окно, где

и сможете его сохранить , нажав кнопку (2) и введя желаемое название в диалоговом окне.

При достаточно большом Quality Factor изображение сильно измениться. Если это JPEG алгоритм то будут ярко выражены блоки размера 8x8.(в случае алгоритма JPEG2000, блочного деления не будет)

До:

После:

Алгоритм JPEG

JPEG – достаточно мощный алгоритм. Практически он является стандартом де-факто для полноцветных изображений. Оперирует алгоритм областями 8x8, на которых яркость и цвет меняются сравнительно плавно. Вследствие этого, при разложении матрицы такой области в двойной ряд по косинусам значимыми оказываются только первые коэффициенты. Таким образом, сжатие в JPEG осуществляется за счет плавности изменения цветов в изображении.

Алгоритм разработан группой экспертов в области фотографии специально для сжатия 24-битных изображений. JPEG – Joint Photographic Expert Group – подразделение в рамках ISO – Международной организации по стандартизации. Название алгоритма читается ["jei"peg]. В целом алгоритм основан на дискретном косинусоидальном преобразовании – (Discrete-Cosine Transform – DCT), применяемом к матрице изображения для получения некоторой новой матрицы коэффициентов. Для получения исходного изображения применяется обратное преобразование.

Раскладывает изображение по амплитудам некоторых частот. Таким образом, при преобразовании мы получаем матрицу, в которой многие коэффициенты либо близки, либо равны нулю. Кроме того, благодаря несовершенству человеческого зрения, можно аппроксимировать коэффициенты более грубо без заметной потери качества изображения.

Как работает алгоритм

Рассмотрим алгоритм подробнее. Пусть мы сжимаем 24-битное изображение .

Шаг 1.
Переводим изображение из цветового пространства RGB, с компонентами, отвечающими за красную (Red), зеленую (Green) и синюю (Blue) составляющие цвета точки, в цветовое пространство YCrCb (иногда называют YUV). В нем Y – яркостная составляющая, а Cr, Cb - компоненты, отвечающие за цвет (хроматический красный и хроматический синий). За счет того, что человеческий глаз менее чувствителен к цвету, чем к яркости, появляется возможность архивировать массивы для Cr и Cb компонент с большими потерями и, соответственно, большими степенями сжатия. Подобное преобразование уже давно используется в телевидении. На сигналы, отвечающие за цвет, там выделяется более узкая полоса частот.

Упрощенно перевод из цветового пространства RGB в цветовое пространство YCrCb можно представить с помощью матрицы перехода:

Шаг 2.
Разбиваем исходное изображение на матрицы 8x8. Формируем из каждой три рабочие матрицы – по 8 бит отдельно для каждой компоненты. При больших степенях сжатия этот шаг может выполняться чуть сложнее. Изображение делится по компоненте Y – как и в первом случае, а для компонент Cr и Cb матрицы набираются через строчку и через столбец. Т.е. из исходной матрицы размером 16x16 получается только одна рабочая матрица ДКП. При этом, как нетрудно заметить, мы теряем 3/4 полезной информации о цветовых составляющих изображения и получаем сразу сжатие в два раза. Мы можем поступать так благодаря работе в пространстве YCrCb. На результирующем RGB изображении, как показала практика, это сказывается несильно.

Шаг 3.
В упрощенном виде ДКП при n=8 можно представить так:

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

Шаг 4.
Производим квантование. В принципе, это просто деление рабочей матрицы на матрицу квантования поэлементно. Для каждой компоненты (Y, U и V), в общем случае, задается своя матрица квантования q (далее МК).

На этом шаге осуществляется управление степенью сжатия, и происходят самые большие потери. Понятно, что, задавая МК с большими коэффициентами, мы получим больше нулей и, следовательно, большую степень сжатия. В стандарт JPEG включены рекомендованные МК, построенные опытным путем. Матрицы для большей или меньшей степени сжатия получают путем умножения исходной матрицы на некоторое число gamma.

С квантованием связаны и специфические эффекты алгоритма. При больших значениях коэффициента gamma потери в низких частотах могут быть настолько велики, что изображение распадется на квадраты 8x8. Потери в высоких частотах могут проявиться в так называемом " эффекте Гиббса", когда вокруг контуров с резким переходом цвета образуется своеобразный "нимб".

Шаг 5.
Переводим матрицу 8x8 в 64-элементный вектор при помощи "зигзаг"-сканирования, т.е. берем элементы с индексами (0,0), (0,1), (1,0), (2,0)...

Таким образом, в начале вектора мы получаем коэффициенты матрицы, соответствующие низким частотам, а в конце - высоким.

Шаг 6.
Свертываем вектор с помощью алгоритма группового кодирования . При этом получаем пары типа (пропустить, число), где “пропустить” является счетчиком пропускаемых нулей, а “число” – значение, которое необходимо поставить в следующую ячейку. Так, вектор 42 3000-2 00001... будет свернут в пары (0,42) (0,3) (3,-2) (4,1)...

Шаг 7.
Свертываем получившиеся пары кодированием по Хаффману с фиксированной таблицей.

Процесс восстановления изображения в этом алгоритме полностью симметричен. Метод позволяет сжимать некоторые изображения в 10-15 раз без серьезных потерь.

Конвейер операций, используемый в алгоритме JPEG:

Существенными положительными сторонами алгоритма является то, что:

1. Задается степень сжатия.
2. Выходное цветное изображение может иметь 24 бита на точку.

Отрицательными сторонами алгоритма является то, что:

1. При повышении степени сжатия изображение распадается на отдельные квадраты (8x8). Это связано с тем, что происходят большие потери в низких частотах при квантовании, и восстановить исходные данные становится невозможно.
2. Проявляется эффект Гиббса – ореолы по границам резких переходов цветов.

Не очень приятным свойством JPEG является также то, что нередко горизонтальные и вертикальные полосы на дисплее абсолютно не видны и могут проявиться только при печати в виде муарового узора. Он возникает при наложении наклонного растра печати на горизонтальные и вертикальные полосы изображения. Из-за этих сюрпризов JPEG не рекомендуется активно использовать в полиграфии, задавая высокие коэффициенты матрицы квантования. Однако при архивации изображений, предназначенных для просмотра человеком, он на данный момент незаменим.

Широкое применение JPEG долгое время сдерживалось, пожалуй, лишь тем, что он оперирует 24-битными изображениями. Поэтому для того, чтобы с приемлемым качеством посмотреть картинку на мониторе в 256-цветной палитре, требовалось применение соответствующих алгоритмов и, следовательно, определенное время. В приложениях, ориентированных на придирчивого пользователя, таких, например, как игры, подобные задержки неприемлемы. Кроме того, если имеющиеся изображения, допустим, в 8-битном формате GIF перевести в 24-битный JPEG, а потом обратно в GIF для просмотра, то потеря качества произойдет дважды при обоих преобразованиях. Тем не менее, выигрыш в размерах архивов зачастую настолько велик (в 3-20 раз), а потери качества настолько малы, что хранение изображений в JPEG оказывается очень эффективным.

Стандартизован JPEG в 1991 году. Но уже тогда существовали алгоритмы, сжимающие сильнее при меньших потерях качества. Дело в том, что действия разработчиков стандарта были ограничены мощностью существовавшей на тот момент техники.

Несколько слов необходимо сказать о модификациях этого алгоритма. Хотя JPEG и является стандартом ISO, формат его файлов не был зафиксирован. Пользуясь этим, производители создают свои, несовместимые между собой форматы, и, следовательно, могут изменить алгоритм. Так, внутренние таблицы алгоритма, рекомендованные ISO, заменяются ими на свои собственные. Кроме того, легкая неразбериха присутствует при задании степени потерь. Например, при тестировании выясняется, что "отличное" качество, “100%” и “10 баллов” дают существенно различающиеся картинки. При этом, кстати, "100%" качества не означают сжатие без потерь . Встречаются также варианты JPEG для специфических приложений.

Характеристики алгоритма JPEG:
Степень сжатия 2-200 (Задается пользователем).
Класс изображений: Полноцветные 24 битные изображения или изображения в градациях серого без резких переходов цветов (фотографии).
Симметричность: 1
Характерные особенности: В некоторых случаях, алгоритм создает "ореол" вокруг резких горизонтальных и вертикальных границ в изображении (эффект Гиббса). Кроме того, при высокой степени сжатия изображение распадается на блоки 8x8 пикcелей.

Назад

К cодержанию

Вперёд