Вейвлет-сжатие на «пальцах»: практика

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

К чему слова? Давайте напишем программу, сжимающую изображения! (Под катом много картинок!)

Инструментарий

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

Нам потребуется интерпретатор и библиотека PIL (Python Imaging Library). В сети достаточно инструкций по их установке, поэтому не буду на этом останавливаться.

Единственное, хочу порекомендовать для разработки использовать ipython notebook. Замечательная штука, предоставляющая веб-интерфейс, напоминающий пакет Mathematica, но программируемая на Python.

Краткое содержание предыдущей серии

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

Первое рассмотренное преобразование — преобразование Хаара — работает с парами значений яркости. Преобразование заключается в умножении на матрицу

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

Коэффициенты преобразования мы находили решая систему нелинейных уравнений. Вторая строчка — это просто коэффициенты, записанные в обратном порядке с чередованием знаков.

Умножая эту матрицу на вектор-столбец из 4-х значений яркости, получаем только два значения (низкочастотный и высокочастотный коэффициенты). Именно из-за такой несправедливости четвёрки нужно брать не подряд, а «внахлёст», с шагом 2.

Например, если у нас есть значения [1, 2, 3, 4], то нужно взять две четвёрки: [1, 2, 3, 4] (сюрприз, сюрприз!) и [3, 4, 1, 2]. Последняя имеет такой странный вид из-за того, что нам не хватило значений справа и пришлось вернуться по кругу в начало.

К слову, матрица всего преобразования в этом случае будет иметь вид:

Аналогично можно записать матрицы преобразований более высокого порядка: D6, D8 и так далее. (Номера всегда чётные. Подумайте, почему так.)

Начинаем!

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

from math import sqrt CL = [(1 + sqrt(3)) / (4 * sqrt(2)),     (3 + sqrt(3)) / (4 * sqrt(2)),     (3 - sqrt(3)) / (4 * sqrt(2)),     (1 - sqrt(3)) / (4 * sqrt(2))] 

Буква L означает, что это коэффициенты первой строки, относящейся к фильтру низких (L, low) частот.

Это достаточно универсальный подход. Если мы захотим реализовать другое преобразование — просто заменим список.

Коэффициенты высокочастотного фильтра (HPF, high-pass filter) будем находить при помощи отдельной функции, записывающей список в обратном порядке и чередующей знаки.

def hpf_coeffs(CL):     N = len(CL)                    # Количество коэффициентов     CH = [(-1)**k * CL[N - k - 1]  # Коэффициенты в обратном порядке с чередованием знака         for k in xrange(N)]     return CH 

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

Прямое одномерное преобразование

Теперь, когда у нас есть коэффициенты, мы можем выполнить само преобразование. Его результатом будет список, содержащий взвешенные суммы с коэффициентами поочерёдно из списков CL и CH.

Взвешенные суммы (суммы попарных произведений) — это произведения строк матрицы на вектор-столбец. Чётные строки матрицы — это низкочастотый фильтр, а нечётные — высокочастотный. Вот откуда чередование.

Списки из взвешенных суммы, вычисляемых «вдоль» другого списка в математике называются свёртками. Есть эффективные алгоритмы вычисления этих свёрток на основе преобразования Фурье. (Внезапно, правда? Преобразование Фурье — это синусы и косинусы, и вдруг для суммирования используется!) Но мы сторонники простоты и не будем заниматься преждевременными оптимизациями. Посчитаем в лоб, так понятнее. А оптимизации и красивые хаки оставим читателю как упражнение.

Функцию назовём pconv. P — от слова pair (пара), а conv — convolution (свёртка).

def pconv(data, CL, CH, delta = 0):     assert(len(CL) == len(CH))         # Размеры списков коэффициентов должны быть равны     N = len(CL)     M = len(data)     out = []                           # Список с результатом, пока пустой     for k in xrange(0, M, 2):  # Перебираем числа 0, 2, 4…         sL = 0                         # Низкочастотный коэффициент         sH = 0                         # Высокочастотный коэффициент         for i in xrange(N):      # Находим сами взвешенные суммы             sL += data[(k + i - delta) % M] * CL[i]             sH += data[(k + i - delta) % M] * CH[i]         out.append(sL)                 # Добавляем коэффициенты в список         out.append(sH)     return out 

Может быть непонятно, что это за значение delta, равное по умолчанию нулю. Оно позволяет нам начать свёртку не с первого элемента, а с произвольного. Это пригодится нам чуть позже.

Вычисление остатка (выражение «%M») при вычислении взвешенных сумм нужно, чтобы не выходить за границы списка. Если M = 4, а индекс вдруг окажется, равным 5, то мы вернёмся снова к 1-му элементу, так как 5%4 даёт в результате 1.

Давайте, испытаем нашу функцию на ненормированном преобразовании Хаара (которое с полусуммами и полуразностями).

>>>C = [0.5, 0.5] >>>pconv([1, 2, 3, 4], C, hpf_coeffs(C)) [1.5, -0.5, 3.5, -0.5] 

Работает!

А обратно?

Как читатель помнит, мы специально так конструировали наши матрицы, чтобы обратное преобразование выполнялось как можно проще. Благодаря ортогональности матриц H и D4, обратные матрицы для них можно найти простым транспонированием.

Рассмотрим ещё раз матрицу преобразования D4, но для вектор-столбца произвольной длины:

Умножая такую матрицу на вектор-столбец со значениями яркостей пикселей мы получим декоррелированные значения. Правда, матрица очень большая, поэтому вместо явного умножения мы воспользовались свёрткой.

Обратную матрицу можно получить транспонированием:

Очень похоже на исходную матрицу, но нет квадратного куска в нижнем левом углу. Хотя если «отрезать» последние два столбца и «приклеить» в начало, то получится матрица такой же формы (хоть и с другими коэффициентами) и мы сможем воспользоваться уже готовой функцией pconv. Стоп! А ведь у нас там есть параметр delta, который задаёт смещение в исходных данных. А ведь это и есть то самое отрезание и приклеивание. Тайна параметра delta раскрыта!

Для D4 нужно «отрезать» 2 последних столбца, для D6 — четыре и т. д.

Дело за малым — определить, порядок коэффициентов в строках обратной матрицы. Для D4 порядок очевиден:

Но как быть в случае D6, D8 и других? Посмотрим на матрицу преобразования D4. Нас интересуют строки транспонированной матрицы, то есть столбцы исходной.

Наши векторы — это третий и четвёртый столбики. Независимо от порядка преобразования смещение строк всегда происходит на два элемента, поэтому закономерность тут такая.

Первый вектор:
[предпоследний элемент первой строки, предпоследний элемент второй, первый элемент первой строки (она совпадает с третьей), первый элемент второй].

Второй вектор:
[последний элемент первой строки, последний элемент второй, второй элемент первой строки (она совпадает с третьей), второй элемент второй].

То есть, коэффициенты поочерёдно выбираются с шагом 2 из CL и CH начиная с предпоследнего для первой строки обратной матрицы или последнего для второй строки.

Это правило будет верно и для D4, и для D6 и даже для преобразования Хаара (которое на самом деле то же, что и D2).

В Python к предпоследнему элементу можно обратиться по индексу -2, а к последнему по индексу -1. Эту полезную особенность мы и использовали в pconv.

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

def icoeffs(CL, CH):     assert(len(CL) == len(CH))         # Размеры списков коэффициентов должны быть равны     iCL = []  # Коэффициенты первой строки     iCH = []  # Коэффициенты второй строки     for k in xrange(0, len(CL), 2):         iCL.extend([CL[k-2], CH[k-2]])         iCH.extend([CL[k-1], CH[k-1]])     return (iCL, iCH) 

Проверим её работу:

>>>C = [0, 1, 2, 3] >>>icoeffs(C, hpf_coeffs(C)) ([2, 1, 0, 3], [3, 0, 1, -2]) 

Похоже, что коэффициенты переупорядочиваются правильно.

Можно протестировать работу функций, преобразовав какой-то список (скажем, [0, 1, 2, 3]), а потом выполнив обратное преобразование.

>>>CL = [(1 + sqrt(3)) / (4 * sqrt(2)), >>>    (3 + sqrt(3)) / (4 * sqrt(2)), >>>    (3 - sqrt(3)) / (4 * sqrt(2)), >>>    (1 - sqrt(3)) / (4 * sqrt(2))] >>>X = [0, 1, 2, 3] >>>CH = hpf_coeffs(CL) >>>iCL, iCH = icoeffs(CL, CH) >>>Y = pconv(X, CL, CH) >>>X2 = pconv(Y, iCL, iCH, len(CL) - 2) [2.5077929123346285e-16, 1.0, 1.9999999999999991, 2.9999999999999991] 

Что за ерунда! Должно было получиться снова [0, 1, 2, 3]! Стоп, ложная тревога. Первое число — это всего-навсего , то есть практически ноль. Остальные числа тоже близки к исходным. Дело в том, что мы работаем с иррациональными числами, а в таких делах неизбежны ошибки округления.

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

Так, преобразования мы запрограммировали. Но преобразования-то эти одномерные! А картинки, которые мы планируем сжимать, очень даже двумерны. Но это легко решается. Двумерное вейвлет-преобразование выполняется как одномерное по всем строчкам, а потом по всем столбикам (или наоборот, по вкусу).

Работать с изображениями мы будем при помощи библиотеки PIL, основанной на numpy. А там есть удобный способ для обращения к элементам двумерных массивов.

Если значения яркостей пикселей нашего изображения хранятся в переменной image, то получить список яркостей 4-й, например, строки, можно так:

image[4, :] 

А список яркостей в 3-м столбце так:

image[:, 3] 

Возьмём для экспериментов какое-то изображение. Важно, чтобы его размеры были степенями двойки! (Вообще, важно, чтобы они были чётными, но мы потом будем делать повторные преобразования.)

Например, можно взять такое (раз уж Лена всем поднадоела):

У него размеры 512×512, так что оно нам подходит.

Загрузим его в память (заодно преобразовав к чёрно-белому виду и сформировав матрицу яркостей):

import PIL.Image as Image image = Image.open('/tmp/boat.png').convert('L') image = array(image) / 255.0 # Диапазон яркостей — [0, 1] imshow(image, cmap=cm.gray)  # Отобразим на экране 

На экране должно появиться что-то вроде:

Итак, переменная image — это массив яркостей.

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

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

def dwt2(image, CL):     CH = hpf_coeffs(CL)   # Вычисляем недостающие коэффициенты     w, h = image.shape    # Размеры изображения     imageT = image.copy() # Копируем исходное изображение для преобразования     for i in xrange(h):   # Обрабатываем строки         imageT[i, :] = pconv(imageT[i, :], CL, CH)     for i in xrange(w):   # Обрабатываем столбцы         imageT[:, i] = pconv(imageT[:, i], CL, CH)      # Переупорядочиваем столбцы и строки     data = imageT.copy()     data[0:h/2, 0:w/2] = imageT[0:h:2, 0:w:2]     data[h/2:h, 0:w/2] = imageT[1:h:2, 0:w:2]     data[0:h/2, w/2:w] = imageT[0:h:2, 1:w:2]     data[h/2:h, w/2:w] = imageT[1:h:2, 1:w:2]     return data 

Проверим работу нашей функции на преобразовании D4, коэффициенты которого мы поместили в CL и CH.

data2 = dwt2(image, CL) imshow(data2, cmap=cm.gray) 

Сразу предупреждаю! Так как мы ничего не оптимизировали, то работать будет медленно. На моём скромном нетбуке расчёт занял полтора десятка секунд.

В итоге должно получиться такая картинка.

Всё, как мы и предполагали! Уменьшенное изображение в углу и небольшие по модулю (о чём говорит чёрный цвет) коэффициенты в остальной части.

А вот график значений вдоль 50-й строки преобразованного изображения.

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

Можно попробовать с вейвлетом Хаара:

C = [1/sqrt(2), 1/sqrt(2)] data3 = dwt2(image, C) imshow(data3, cmap=cm.gray) 

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

Или даже взять коэффициенты для D6 из википедии (не забудьте разделить их на они там нормированы на двойку). А самые смелые могут вообще получить коэффициенты аналитически.

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

Обратное двумерное преобразование

То, что мы умеем делать двумерное преобразование — это здорово, но это умение бесполезно без обратного преобразования.

Его как раз сделать несложно. Выполняем те же шаги, что и с прямым преобразованием, но в обратном порядке. И не забываем, про другие коэффициенты и ненулевое значение delta в pconv.

def idwt2(data, CL):     w, h = data.shape   # Размеры изображения          # Переупорядочиваем столбцы и строки обратно     imageT = data.copy()     imageT[0:h:2, 0:w:2] = data[0:h/2, 0:w/2]     imageT[1:h:2, 0:w:2] = data[h/2:h, 0:w/2]     imageT[0:h:2, 1:w:2] = data[0:h/2, w/2:w]     imageT[1:h:2, 1:w:2] = data[h/2:h, w/2:w]      CH = hpf_coeffs(CL)     iCL, iCH = icoeffs(CL, CH)     image = imageT.copy() # Копируем исходное изображение для преобразования     for i in xrange(w):   # Обрабатывем столбцы         image[:, i] = pconv(image[:, i], iCL, iCH, delta=len(iCL)-2)     for i in xrange(h):   # Обрабатывем строки         image[i, :] = pconv(image[i, :], iCL, iCH, delta=len(iCL)-2)      return image 

Для проверки выполним обратное преобразование полученного ранее изображения image2.

>>>image_rec = idwt2(image2, CL) >>>imshow(image_rec, cmap=cm.gray) 

Получили исходное изображение (ну, или очень-очень похожее), так что у нас всё работает!

Не останавливаемся на достигнутом!

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

Для вейвлета D4 предельным размером «уголка» являются размеры 4×4, так как изображение меньшего размера преобразовать уже не удастся — не хватит значений, чтобы умножить на 4 коэффициента матрицы.

Например, для преобразования D4, подобное рекурсивное преобразование можно запрограммировать следующим образом:

data = image.copy() w, h = data.shape  while w >= len(CL) && h >= len(CL):     data[0:w, 0:h] = dwt2(data[0:w, 0:h], CL)     w /= 2     h /= 2 imshow(data, cmap=cm.gray) 

Картинка при этом получается не особо интересная.

Посмотрим, в каком диапазоне меняются значения. Вот график значений из нулевой строки.

Первые несколько значений на самом деле зашкаливают. Зато остальные малы.

Удаляем лишнее

А теперь давайте проделаем то, ради чего, всё это задумывалось — получим много нулей!

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

В качестве такого порогового значения возьмём 0,05. Чем больше, тем больше нулей, но и тем больше потери!

from numpy import abs threshold = 0.05 data[abs(data)<threshold] = 0 

После этого график станет заметно ровнее.

Подсчитаем количество нулей:

>>>sum(data == 0) 224000 

Неплохо! Учитывая, что в изображении всего-то 512×512=262144 пикселя, выходит, что мы просто отбросили 85,4% имеющихся коэффициентов!!!

Сильно ли это повредило изображению?

Выполним обратное преобразование. (Как именно — упражнение для читателя!)

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

Кто увидел отличия — тот молодец! Я вот с трудом нашёл.

А если взять 0,1 (обнуляется 92,9% коэффициентов)?

Тоже неплохо!

Возьмём теперь порог, равный 0,2 (96,9% отбрасываемых коэффициентов).

Уже заметны серьёзные искажения.

Ну и наконец, попробуем отбросить всё, что меньше по модулю, чем 0,5 (99,2% отброшено).

Жуткое мыло! Но квадратики не такие контрастные, как в пережатых JPEG. Мы использовали D4, поэтому мыло «квадратичное», с плавными градиентами внутри.

Для сравнения попробуем с этим же порогом (то есть, 0,5) сжать изображение, но уже с использованием преобразования Хаара.

Ужас-ужас! Ещё и линейные искажения. И отброшено всего-то 99,2% коэффициентов. Хуже, чем D4.

Выводы

Итак.

1. После преобразования D4 мы можем смело отбросить 90% коэффициентов простым обнулением и не потеряем в качестве. А если применить продвинутые методики квантования, то результат можно ещё улучшить.

2. D4 лучше, чем D2 (преобразование Хаара). D6 будет лучше, чем D4 и так далее. Но тут есть проблема. Чем выше порядок, тем раньше нам придётся остановить процесс рекурсивного преобразования верхнего левого уголка. Так что во всём нужно знать меру.

Что дальше?

На самом деле мы ещё не закончили. Да, мы отбросили 90% коэффициентов. Но это вещественные коэффициенты, занимающие 8 байт! Так что сжатие у нас будет очень и очень небольшое. Поэтому при сохранении в файл коэффициенты округляют и используют для хранения ограниченное количество бит. Кроме того, для лучшей упаковки коэффициенты можно переупорядочить особым образом. Об этом можно написать не одну статью.

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

Сохранение можно осуществить так.

from numpy import save save("boat", data.astype(float16)) 

Несмотря на значительную потерю точности из-за преобразования float64→float16 и сжатие в лоб без каких-либо специальных алгоритмов, получаем довольно неплохой результат:

При этом получили фактор сжатия, равный (512×512)/43604 ≈ 6 или 1,33 бит/пиксель. Не бог весть что, но мы не сильно-то и старались сжать. Правильное квантование и хороший алгоритм сжатия могут очень существенно улучшить этот результат! Так что тут есть куда расти! Но это уже выходит за рамки нашего «проекта на вечер». Может, как-нибудь в другой раз. 😉

Домашнее задание

1. Попробуйте объединить приведённые куски кода в один скрипт. Пусть он принимает в качестве параметров имя графического файла, величину порога и тип преобразования, сжимает указанное изображения, и сохраняет результаты в файл. Предусмотрите два режима работы: сжатие и распаковка.
2. Поэкспериментируйте с кодированием и декодированием разных изображений.
3. Попробуйте другие перобразования.
4. Подумайте, как более экономно сжимать коэффициенты.
5. Попробуйте ускорить работу программы.
6. Реализуйте сжатие цветных изображений.
7. Снимите ограничение на размеры (пока что они должны быть степенями двойки).
8. ???
9. PROFIT!!!

Очередной убийца JPEG готов! 😉

Всем спасибо! Надеюсь, было интересно!