Алгоритмы квантования для полутоновых и цветных изображений. Дискретизация и квантование изображений

Привлекает внимание

Например, старый добрый формат GIF использует палитру, максимум на 256 цветов. Если вы захотите сохранить серию своих селфи как gif-анимацию (кому бы это надо было), то первое, что вам, а точнее программе, которую вы будете для этого использовать, надо будет сделать – создать палитру. Можно использовать статическую палитру, например web-safe colors , алгоритм квантизации получиться очень простым и быстрым, но результат будет «не очень». Можно создать оптимальную палитру на основе цветов изображения, что даст результат наиболее визуально похожий на оригинал.

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

Далее вас ждёт скучное и непонятное повествование о методе медианного сечения, алгоритму рассеивания ошибок (шума квантизации) по Флойду-Стейнбергу (и не только), особенностях цветового восприятия человеческого глаза, а так же немного говно кода.

Предыстория

Давным-давно, когда Nokia была тёплой и ламповой главенствовала на рынке смартфонов, а владельцы смартфонов гордо звали себя «смартфонщики», в те стародавние времена писал я простенькие программки на python для series60. На одну из них намедни наткнулся копаясь в архивах. GifTool – программка для создания gif-анимации из набора картинок. В ней я реализовал квантизацию методом медианного сечения, алгоритм сжатия LZW, вся структура файла создавалась самостоятельно, для неизменившихся на следующем слайде пикселей использовалась прозрачность, чтобы уменьшить итоговый размер файла. Захотелось мне освежить свою память, посмотреть – как же она работала. Открыл код и … Это чувство, когда ты не можешь разобраться в своём говнокоде десятилетней давности. Про PEP8 я тогда не знал, поэтому читаемость кода была чуть менее чем никакой (тогда мне нравился минимализм, как и многим начинающим программистам). Прослезился, поплевался, отрефакторил в PyCharm, разобрался как реализовал метод медианного сечения, по быстрому накидал «грязный» скрипт. Работает! Палитра создаётся, изображение на выходе получается сносное. И тут меня закусило – смогу ли я добиться лучших результатов, чтобы картинка визуально была как можно ближе к оригиналу.

Итак – метод медианного сечения. Он прост до безобразия. Первым делом надо из всех уникальных цветов изображения составить RGB куб. Далее рассечь его по самой длинной стороне. Например, диапазон красного у нас от 7 до 231 (длина 231-7=224), зелёного от 32 до 170 (длина 170-32=138), синего от 12 до 250 (длина 250-12=238), значит, будем «резать» куб по синей стороне. Получившиеся сегменты так же рассекаем по длинной стороне и т.д. пока не получим 256 сегментов. Для каждого сегмента высчитать средний цвет – так мы и получим палитру.

Пара картинок почти в тему, для наглядности

Что здесь можно улучшить? Первое, что приходит в голову – вычислять средний цвет не тупо сложив все цвета и разделив на их количество [ sum(color) / count(color) ], а с учётом, сколько раз каждый цвет встречается в изображении. То есть каждый цвет умножаем на количество его вхождений в изображении, полученные значения складываем, результат делим на количество вхождений в изображении всех цветов данного сегмента [ sum(color * total) / sum(total) ]. В результате, наиболее часто встречаемые цвета имеют приоритет при вычислении, но и редкие цвета вносят свои корректировки, поэтому палитра получается лучше, визуальное отклонение цветов меньше. Для лучших результатов желательно ещё учитывать гамму, но я оставил это на потом. Второе не так явно – медианное сечение совсем не учитывает особенности восприятия цвета человеческим глазом. Оттенки зелёного мы воспринимаем гораздо лучше оттенков синего. Я решил исправить это недоразумение и «сплющил» куб – длины сторон помножил на коэффициенты из этой статьи . В результате по зелёной и красной стороне сечений стало больше, по синей меньше. Такого решения я больше нигде не встречал (может плохо искал), но результат на лицо.

Теперь у нас есть оптимальная палитра, не идеальная конечно (я знаю, что её можно ещё улучшить), но достаточно хорошая. Следующий шаг – индексирование цветов изображения. Самый простой вариант – в каком сегменте находится цвет, такой и индекс. Быстро и просто. Но есть одно но, и даже не одно, поэтому к данному шагу мы ещё вернёмся.

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

Вот теперь можно снова вернуться к индексированию. Рассеиванием ошибок мы изменяем цвета пикселей и получаем такие, которых нет в нашем RGB-кубе (напомню, он составлен исключительно из цветов изображения). Теперь нельзя просто посмотреть в каком сегменте находится цвет, чтобы назначить индекс. Решение нашлось сразу – поиск ближайшего цвета в палитре . В данную формулу я подставил всё те же коэффициенты. Сравнивая результаты подбора цвета палитры по индексу сегмента в который входит исходный цвет и результаты поиска ближайшего цвета, наглядно увидел, что ближайший цвет часто оказывается в соседнем сегменте. Если исходный цвет находится ближе к центру сегмента – то индекс сегмента соответствует индексу цвета в палитре, но чем ближе исходный цвет к краям сегмента, тем больше вероятность, что ближайший цвет окажется в соседнем сегменте. В общем, единственный правильный путь индексирования – поиск ближайшего цвета в палитре. Но у поиска есть минус – он медленный, очень медленный. Писать числодробилку на python плохая идея.

Ну вот, хотел объяснить в двух словах, а получилась целая куча непонятной писанины. Надеюсь, код я пишу лучше, чем объясняю, поэтому вот ссылочка на github . Код несколько раз переписывался, сначала совершенствовался алгоритм, пока результат меня не устроил, потом оказалось, что он жрёт слишком много оперативы при обработке фотографий (сначала тестировал на небольших картинках), пришлось перенести RGB-куб, медианное сечение и карту пикселей в базу данных (sqlite). Скрипт работает очень медленно, но результат получается лучше, чем квантизация средствами PIL/Pillow и GIMP’ом (в нём эта операция называется индексирование).

Наглядная демонстрация:

Оригинал

Результат квантизации в GIMP, оптимальная палитра на 256 цветов + размывание цвета по Флойду-Стенбергу (нормальное)

Результат квантизации PIL/Pillow image.convert(mode="P", dither=PIL.Image.FLOYDSTEINBERG, palette=PIL.Image.ADAPTIVE, colors=256)

Результат квантизации моим кодом

На что обратить внимание: рассеивание ошибки у GIMP сильно «шумит», PIL/Pillow создает не очень оптимальную палитру и практически не рассеивает ошибки (резкие переходы между цветами).
Если не видите разницу - посмотрите другие примеры на github .

P.S.: есть замечательная программа Color Quantizer , которая справляется с данной задачей лучше и быстрее, поэтому практического смысла мой скрипт не имеет, сделан исключительно из «спортивного» интереса.
UPD: обновил проект на github . Добавил алгоритм квантизации Octree (октодерево), популярные формулы рассеивания ошибок, поиск ближайшего цвета по среднему значению красного.
При цифровой обработке изображений непрерывный динамический диапазон значений яркости делится на ряд дискретных уровней. Эта процедура называется квантованием. Квантователь преобразует непрерывную переменную в дискретную переменную , принимающую конечное множество значений
. Эти значения называются уровнями квантования. В общем случае преобразование выражается ступенчатой функцией (рис. 8). Если яркость отсчета изображения принадлежит интервалу
(т.е., когда
), то исходный отсчет заменяется на уровень квантования , где
- пороги квантования. При этом полагается, что динамический диапазон значений яркости ограничен и равен
.

Рис.8.Функция, описывающая квантование
Задача построения квантователя состоит в определении значений порогов и уровней . Простейший способ решения этой задачи состоит в разбиении динамического диапазона на одинаковые интервалы. Однако такое решение не является наилучшим. Если значения яркости большинства отсчетов изображения сгруппированы, например, в «темной» области и число уровней ограничено, то целесообразно квантовать неравномерно. В «темной» области следует квантовать чаще, а в «светлой» реже. Это позволит уменьшить ошибку квантования .

В реальных системах в основном используются два типа квантования – линейное гамма-корректированное. В последнем случае аналоговый сигнал перед квантованием подвергается нелинейному преобразованию x’=x 1 /  . Такая функция реализована практически во всех промышленно выпускаемых ПЗС камерах. Стандартное значение  равно 1.4.

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

Оптимальный выбор числа уровней дискретизации в значительной степени зависит от характеристик приемного устройства (ПЗС-камеры, например). ПЗС камеры общего назначения редко имеют соотношение сигнал-шум больше чем 46дБ. Соотношение сигнал- шум определяется следующим выражением:
, где
-максимальная амплитуда полезного сигнала,
-среднеквадратичная амплитуда шума. Соответственно, при отношении сигнал шум равном 46дБ полезное число уровней квантования составляет 200, что говорит о целесообразности использования восьмибитного квантователя.

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

Рис.1.5.Функция, описывающая квантование

Задача построения квантователя состоит в определении значений порогов и уровней . Простейший способ решения этой задачи состоит в разбиении динамического диапазона на одинаковые интервалы. Однако такое решение не является наилучшим. Если значения яркости большинства отсчетов изображения сгруппированы, например, в «темной» области и число уровней ограничено, то целесообразно квантовать неравномерно. В «темной» области следует квантовать чаще, а в «светлой» реже. Это позволит уменьшить ошибку квантования .

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

(1.12)

в предположении, что яркость - случайная величина с известной плотностью вероятности .

Cреднеквадратическая ошибка квантования (1.12) равна

. (1.13)

Дифференцируя (1.13) по переменным , и приравнивая производные нулю, получаем нелинейных уравнений

.

Следует отметить, что крайние пороги и определяются динамическим диапазоном яркости. Уравнения (1.14) нетрудно привести к виду

.

Из (1.15) следует, что пороги должны располагаться по середине между двумя соседними уровнями и . Решение этих уравнений можно найти итеративным способом. Оптимальный квантователь, удовлетворяющий критерию (1.12), называется квантователем Ллойда-Макса , а среднеквадратическая ошибка для такого квантователя равна

(1.16)

При равномерном распределении яркости нелинейные уравнения (1.15) можно представить в виде

,

а среднеквадратическая ошибка равна .

В системах цифровой обработки изображений стремятся уменьшить число уровней и порогов квантования, т.к. от их количества зависит длина двоичного кодового слова, которым представляются проквантованные отсчеты в ЭВМ. Однако при относительно небольшом числе уровней на проквантованном изображении появляются ложные контуры. Они возникают вследствие скачкообразного изменения яркости проквантованного изображения (рис.1.6) и особенно заметны на пологих участках ее изменения.

Ложные контуры значительно ухудшают визуальное качество изображения, т.к. зрение человека особенно чувствительно именно к контурам. При равномерном квантовании типичных изображений требуется не менее 64 уровней. На рис.1.7.а и 1.7.б приведены результаты равномерного квантования изображения «Портрет» соответственно на 256 и 14 уровней квантования.

Рис.1.6. К механизму возникновения ложных контуров

Рис.1.7. Результаты равномерного квантования
Рис.1.8. Результат неравномерного квантования	Рис.1.9. Гистограмма изображения “Портрет”

В темных частях изображения на рис. 1.7.б заметны ложные контуры. Использование квантователя Ллойда-Макса позволяет существенно снизить их уровень (см. рис. 1.8, где число уровней квантования также равно 14). На рис. 1.9 приведена гистограмма яркости изображения «Портрет» при 256 уровнях квантования и отмечены пороги при . Из рисунка следует, что чаще квантуются те области динамического диапазона, в которых сгруппированы значения яркости отсчетов.

Чтобы избежать неравномерного квантования, которое не может быть выполнено с помощью стандартного АЦП, используют нелинейные преобразования (рис.1.10). Отсчет исходного изображения подвергается нелинейному преобразованию, чтобы плотность распределения вероятностей преобразованных отсчетов была равномерной, т.е. выполняется процедура эквализации, которая подробно описана в главе 2. Затем отсчеты квантуются с равномерным шагом и подвергаются обратному нелинейному преобразованию.

Рис.1.10. Квантование с предварительным нелинейным преобразованием

Для разрушения ложных контуров Робертс предложил перед равномерным квантованием к отсчетам яркости добавлять шум с равномерной плотностью распределения вероятностей. Добавленный шум переводит одни отсчеты изображения на уровень выше, а другие на уровень ниже. Тем самым разрушаются ложные контуры. Дисперсия добавляемого шума должна быть небольшой, чтобы не привести к искажениям, воспринимаемым как «снег» на изображении, и в то же время достаточной для разрушения ложных контуров. Обычно используют равномерно распределенный шум на интервале . Результаты равномерного квантования на 14 и 8 уровней изображения «Портрет» с предварительным добавлением шума приведены на рис.1.11.а и 1.11.б. При 8-ми уровнях квантования добавляемый шум становится слишком заметным, однако ложные контуры разрушены практически полностью.

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

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

.

Эта операция повторяется для всех блоков. Получаемое при этом изображение квантуется на два уровня. На рис. 1.12.а приведено полутоновое изображение «Портрет» с добавленным возмущающим сигналом. На рис. 1.12.б,в приведены результаты бинарного квантования изображения «Портрет» с добавленным возмущающим сигналом (рис.1.12.б) и без него (рис.1.12.в).

Рис.1.12.Растрирование изображений

Бинарное растровое изображение обеспечивает значительно лучшее зрительное впечатление, чем обычное бинарное изображение. Передача шкалы яркости при растрировании достигается благодаря изменению геометрических размеров белого пятна, наблюдаемого на черном фоне. Если в блоке сгруппировались «светлые» отсчеты, то геометрические размеры белого пятна максимальны и равны размеру блока. При уменьшении яркости его геометрические размеры также уменьшаются. Глаз человека выполняет локальное усреднение, создавая иллюзию наблюдения полутонового изображения. Процедура растрирования особенно эффективна при печати изображений с высоким разрешением, когда одиночное пятно едва различимо глазом.

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

Дискретизация изображений

Дискретизация – это преобразование непрерывного сигнала в последовательность чисел (отсчетов), то есть представление этого сигнала по какому-либо конечномерному базису. Это представление состоит в проектировании сигнала на данный базис.

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

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

Для случая равномерной дискретизации справедлива теорема Котельникова , опубликованная в 1933 году в работе “О пропускной способности эфира и проволоки в электросвязи”. Она гласит: если непрерывный сигнал имеет спектр, ограниченный частотой , то он может быть полностью и однозначно восстановлен по его дискретным отсчетам, взятым с периодом , т.е. с частотой .

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

.

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

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

Приведенное соотношение означает возможность точного восстановления сигналов с ограниченным спектром по последовательности их отсчетов. Сигналы с ограниченным спектром – это сигналы, спектр Фурье которых отличен от нуля только в пределах ограниченного участка области определения. Оптические сигналы можно отнести к ним, т.к. спектр Фурье изображений, получаемых в оптических системах, ограничен из-за ограниченности размеров их элементов. Частоту называют частотой Найквиста . Это предельная частота, выше которой во входном сигнале не должно быть спектральных компонентов.

Квантование изображений

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

Рис. 1. Функция, описывающая квантование

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

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