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

  • Tutorial

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

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



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



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

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

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

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



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

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


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


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

Умножая эту матрицу на вектор-столбец из 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 готов! ;)

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

Similar posts

Comments 31

    +8
    Ну у тебя и «пальцы» :D. А вообще — спасибо! Статья очень содержательная!
      0
      спасибо, всё круто!
      недавно проходили в универе вейвлет функции Хаара, теперь есть возможность практики :D
        +2
        Добавьте python в тэги.
        Т.к. все таки используете python
          +3
          Язык программирования это всего лишь инструмент для достижения поставленной задачи. Так что тег тут будет лишним, ведь сия статья вовсе не о языке.
            +2
            Почему же. Человек искал что-нибудь про питон, а найдет он реализацию несложного алгоритма сжатия картинок на питоне. Стоит того, чтобы быть найденным, мне кажется.
              +4
              Добавил тег. Думаю, не повредит.
          +1
          А почему нет ни слова про PSNR?
            +5
            О, тут про очень многое нет ни слова! (О чём я жалею, но, увы, времени всё описывать нет.)

            Например, про то же квантование упомянул лишь мельком, но не упомянуть было сложно. А PSNR тут особой роли не играет, потому решил не вводить лишние сущности. Если всё упоминать, то получится не статейка-tutorial, а книжка Гонсалеса и Вудса. :)
            +3
            После прочтения подобных статей, смачно сдобренных «как читатель понял», «как мы помним» — я себя чувствую деревом :(
            • UFO just landed and posted this here
              –19
              Раньше на BASIC примеры приводили, во всяких журналах особенно, типа того же «Радиолюбителя». А теперь вот его место занял Python, который также пихают во все щели по поводу и без повода. Ничего не скажешь, «достойная» замена. Нет бы чего нормальное выучили.
                +14
                \Ну, так назовите язык, который лучше других подошёл бы.

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

                Я мог написать программу на Си (или даже ассемблере). И погряз бы в бесчисленных деталях вроде выделения/освобождения памяти.

                Я мог написать программу на Java/C#. Но тогда ради двух команд пришлось бы городить целый класс.

                Можно было воспользоваться F#. Но ведь есть же люди, для которых функциональное программирование непривычно.

                Хотя, знаете, хорошо подошёл бы MatLab. Там вообще есть toolbox для вейвлетов. Можно было всё, что тут написано реализовать парой строк. Но к чему тогда вообще была бы нужна вся эта писанина?

                Да что бы я не выбрал, всё равно у этого выбора нашлись бы недостатки. Об этом я в самом начале поста написал.

                А с другой стороны. Вам, допустим, не нравится Python. Но Вы поняли, что именно делают приведённые в посте куски кода? Если поняли, значит Python со своей задачей — донести суть алгоритма — справился.

                Если не согласны, то предложите свою альтернативу! Я с радостью принимаю любые конструктивные замечания.
                  –6
                  Ну, если вы хотели показать применимость алгоритма — то почему бы действительно и не MatLab или тот же Maple? Никто ж не заставляет использовать поставляемые «из коробки» фичи платформы и не мешает написать своё. А городить целый класс на Java из единственного метода public static void main(String[] args) { } (утрируя) — это проблема, да.

                  Я понимаю, что все большая часть ЯП плюс-минус универсальна, но всё же существуют более пригодные и менее пригодные для решения конкретных задач языки. И вот это запихивание Python в каждую дырку приводит к тому, что люди перестают осознавать разницу. От MatLab или Maple никто не ждёт чудес производительности, зато там есть адекватные для решения подобных задач синтаксис и семантика. Это делает их подходящими для демонстрации решения задачи сжатия изображений — возвращаясь к теме поста. На Python это тоже можно сделать. Но так как Python менее специфичный, кто-нибудь это потом прочитает и скажет «а мы будем жать изображения на Python!». Опять же, утрируя, но результат всё равно понятен.

                  Мне кажется, что имеет смысл не только обучать людей решению каких-либо задач, но и обучать их делать это правильным образом. Так-то можно гвозди и микроскопом забивать. Можно даже зубочисткой, при некоторой сноровке и большом запасе кинетической энергии. Но лучше это всё-таки делать молотком.
                    +12
                    Про инструменты согласен. Я и сам вычисления по работе вовсе не в Python делаю.

                    Но сами посмотрите: в научных статьях везде MatLab (не видел, чтобы Maple для обработки изображений активно использовали), в промышленности те же самые алгоритмы на Си. Цели-то разные.

                    Моя цель была — показать алгоритм и позволить максимальному количеству людей воспроизвести результаты. Мне кажется, что Python всё же есть у большей доли целевой аудитории, нежели MatLab/Octave.

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

                    Если я напишу формулу с использованием обращения матриц, вычисления интегралов и чего-то такого в Maple, но наверняка у кого-то возникнет вопрос, как это запрограммировать в его любимом C#/Java/C++/C/Ruby и т.д., чтобы самому попробовать поковыряться в алгоритме.

                    Да и раз уж заговорили о numpy. Скажите, много ли принципиальных отличий между Python+numpy/pylab и MatLab? Чем Python не молоток в данном случае? Очень многие используют его для прототипирования при исследованиях, свяязанных с обработкой изображений, машинным обучением, статистикой. Есть неплохие библиотеки вроде PIL, numpy, pylab, scikits. Есть биндинги к OpenCV. Если какая-то часть требует быстродействия — можно написать её на C.

                    Думаю, авторы статей из «Радиолюбителя» тоже не от хорошей жизни все программы на BASIC приводили. Просто в своё время это был один из самых распространённых языков на ПК. Зачем писать программу на Fortran, если у людей BASIC был вместе с операционкой уже.

                    Впрочем, этот спор уже принимает религиозный оттенок.
                +3
                Кто увидел отличия — тот молодец! Я вот с трудом нашёл.

                Отличия конечно есть, на оригинальном изображении отчетливо виден «теплый ламповый звук» зерен пленки, на втором изображении их нет. Также отличия видны на наклонных линиях канатов, их контрастность заметно меньше. Понятно, что чудес не бывает и любое сжатие приводит к потерям и вот тут как раз оценка по PSNR дала бы количественный результат, что лучше, а что хуже. Также для наглядности было бы здорово привести сравнение с тем же jpeg и png. Опять же в цифрах.
                В целом статья отличная. Спасибо!
                  0
                  Да, PSNR даёт количественные оценки, но они не всегда совпадают с мнением человека об отличиях. PSNR распространён только потому, что он очень легко вычисляется. Есть более сложные метрики. Например, SSIM или PSNR-HVS/-HVS-M/-HVA. Они хорошо коррелируют с субъективными оценками качества, но формулы нетривиальные.

                  Хотя, конечно же, Вы правы. Стоило бы рассказать и о метриках, сравнить с тем же JPEG. Но уж очень большая статья бы получилась, а свободного времени у меня не так много. Самому обидно, но про очень многое приходится не упоминать. Может, в следующей статье расскажу о метриках.
                  0
                  Скажите, а вы поделитесь опытом в том, как можно оптимизировать скорость работы?
                    +3
                    Переписать на C с SSE intrinsic или ассемблерными вставками и долго гонять через профайлер. Желательно ещё оптимизировать структуры данных с точки зрения минимизации cache misses и unaligned access.
                      0
                      Такие оптимизации, фактически, делаются ведь под конкретный камень и мамку, нет?
                        +1
                        Нет. То есть если надо совсем идеально — то да, но можно сделать и универсально, при этом лишь едва потеряв в скорости.
                      0
                      Оптимизация — непростое дело. У приведённого кода много недостатков. Если не отказываться от Python, то можно увеличить быстродействие, решив следующие проблемы: многократное копирование массивов, неэффективный алгоритм вычисления свёртки, невекторизованные вычисления (можно использовать numpy).
                        0
                        Вот как раз интересует «Есть эффективные алгоритмы вычисления этих свёрток на основе преобразования Фурье.» — можете подсказать, где о них почитать?
                        И еще один вопрос: можно ли как-то обойти ограничение на то, что размеры изображения должны быть кратны двум?
                          0
                          В матанализе есть теорема о свёртке (так и называется, можно гуглить), согласно которой при преобразовании Фурье свёртка переходит в простое умножение.

                          Таким образом свёртку последовательностей A и B можно выполнить так:
                          1. Выполняем дискретное преобразование Фурье (для него есть эффективный алгоритм — быстрое преобразование Фурье, FFT):
                          A' = FFT(A)
                          B' = FFT(B)
                          2. Перемножаем последовательности поэлементно:
                          C'[i] = A'[i] * B'[i]
                          3. Находим обратное преобразование Фурье:
                          C = IFFT(C')
                          Хоть действий и больше, вычислительная сложность алгоритма получается ниже. Вместо квадратичной сложности получаем линейно-логарифмическую.

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

                          Ограничение на размеры можно обойти по-разному. Самый простой способ — разбить изображение на блоки, у которых размеры являются степенями двойки. При необходимости можно изображение дополнить фиктивными строками и столбиками, которые не отображаются при декодировании.
                          0
                          Если не отказываться от Python, то хотя-бы отказаться от floating point — не знаю там, таблицы юзать или fixed point или все вместе сразу.
                            +1
                            Это усложнило бы восприятие программы. Я же не ставил перед собой цель написать промышленное решение. Это просто иллюстрация к предыдущей статье, а там числа вещественные.

                            Если уж оптимизировать, то лучше сделать быструю свёртку. Это даст более заметный прирост в быстродействии, мне кажется.
                              +2
                              Обработка изображений и видео очень ресурсоемкая задача, поэтому даже в образовательных целях лучше делать все на СИ. Ну ни кто в мире не пишет софт для обработки изображений на Python или Java.
                              Основные принципы работы с изображениями:
                              1. Самая медленная часть современных компьютеров система памяти. Речь идет о нормальных изображения 5-10 мегапикселей, которые не влезают целиком в кэш. Поэтому в зависимости от алгоритма, если
                              время на обработку изображения меньше чтения его из памяти в процессор то нет смысла подключать
                              оптимизацию, если больше то можно уже применять всякие SSE.
                              2. Всячески избегать обработки картинки «против шерсти», то есть использовать последовательный доступ, так как она лежит в памяти построчно с лево на право. Любые переупорядочивания подсаживают скорость в разы. Все из-за то, что контролер памяти может читать только блоками по 8-16 и больше слов и получается очень не эффективная работа контролера, когда читаете картинку сверху в низ.
                              3. Стараться использовать в каждом алгоритме свой внутренний минимальных буфер, допустим если используется оператор свертки 3х3 пикселя, то минимально необходимо держать в буфере 3 строки изображения. Это нужно для того чтобы было проще перенести алгоритм в железо, с си на verilog или VHDL. А без железной реализации алгоритмы не получат широкого использования. Потому JPEG такой популярный, что есть аппаратная реализация
                              в каждой железяке.
                              По поводу вейвлетов, статья очень понятная и полезная, но если речь зашла о сжатии, то надо было упомянуть про энтропию. Без этого понятия тяжело объяснить почему вейвлеты могут сжать картинку.
                                0
                                Обработка изображений и видео очень ресурсоемкая задача, поэтому даже в образовательных целях лучше делать все на СИ. Ну ни кто в мире не пишет софт для обработки изображений на Python или Java.

                                Пишут, пишут! Тот же ImageJ очень даже на Java написан. Всё от задач зависит. Кодеки, конечно, лучше писать на си.

                                Советы Вы даёте действительно дельные. Но я в своём посте не пишу промышленную систему, а иллюстрирую идею. Если бы я думал о быстродействии, я бы вообще другие алгоритмы выбрал. Язык тут даже не на первом месте.

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

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

                                Про энтропию стоило бы рассказать, но естественный её ввод потребовал бы много объяснений. Так как из самого понятия энтропии не следует, что вейвлеты решают задачу сжатия. Я же старался вывести идею вейвлетов из минимального числа базовых идей.
                      • UFO just landed and posted this here
                          +1
                          Можно немного критики после восторгов?

                          Всё это мило, но у стандартное БПФ в JPEG даёт большее разрешение по частоте.

                          Если бы была возможность отсекать не половину частот, а скажем только в корень из двух на каждом этапе, объяснение такой математики было бы интереснее. А с полусуммами — полуразностями и так полно объяснений в сети.
                            0
                            Критика только приветствуется!

                            DCT в JPEG и DWT довольно сильно отличаются, чтобы их так просто сравнивать, хоть оба и являются линейными преобразованиями.

                            К преимуществам DWT можно отнести, например то, что искажения после DWT с сильным квантованием будут менее «неприятными» по сравнению с DCT. В первом случае мы увидим размытую и замыленную картинку с полиномиальными перепадами яркости, а во втором — сильный «звон», обусловленный эффектом Гиббса.

                            Про отсечение корня из двух частот не совсем понял.

                            Объяснения с полусуммами и полуразностями я и сам видел, но там либо дальше полусумм и полуразностей не идут, останавливаясь на преобразовании Хаара, либо сразу начинают с интегралов, пугая читателей. (Про биортогональные вейвлеты вообще ) Моя цель — вывести преобразования максимально естественным образом из общих соображений, не вводя сущности без нужды.

                            Но был бы рад увидеть ссылки на хорошие, по Вашему мнению, статьи. Может, какие-то интересные идеи почерпну.
                              0
                              Это не совсем про сжатие, но вы про DT-CWT читали?

                              Статей по scolar.google.com много, а на русском информации почти нет.

                              Если я верно понимаю, то в каком-нибудь neat image используется подобная методика.

                          Only users with full accounts can post comments. Log in, please.