В нескольких статьях нашего блога уже затрагивалась тема различной цифровой обработки изображений, но, в основном в них мы касались вопросов явной трансформации изображений (перевод в оттенки серого, фильтрация, применение фильтров и некоторые другие операции) и никак не обращались к способам изменения избыточности изображений. Под избыточностью, в данном случае, понимается некая абстрактная мера, связанная с количеством информации, которое содержится в некоем изображении. Обычно, когда упоминаются выражения наподобии «способ изменения избыточности» или «метод уменьшения информационной энтропии», то речь идет об алгоритмах сжатия или алгоритмах эффективного представления неких данных.
Конечно, я бы хотел рассмотреть здесь методы сжатия изображений, но больше всего мне бы хотелось продемонстрировать суть подобных приемов, обращаясь непосредственно к внутреннему устройству и самым базовым их принципам, но, тематика сжатия изображений слишком обширна для данного блога.
Поэтому, я хотел бы рассказать сегодня о дискретном вейвлет-преобразовании, а точнее, об одной из самых известных его версий — дискретном преобразовании Хаара.
Дискретное вейвлет-преобразование — это некоторое преобразование дискретных данных (т.е данных, полученных выборкой через определенные постоянные промежутки времени), позволяющее анализировать (получать) различные частотные компоненты. Существуют целые семейства таких преобразований, которые носят общее название — вейвлеты, и которые реализованы не только для дискретных данных.
Однако, в нашем случае, мы рассмотрим дискретное вейвлет-преобразование, которое стало не только исторически самым первым, но и положило начало само теории (и практике) вейвлет-анализа, а называется оно дискретным преобразованием Хаара. Здесь и далее, я буду называть его либо ДВП (дискретное вейвлет-преобразование), либо просто преобразованием Хаара (вейвлет Хаара).
Вейвлет Хаара устроен на редкость очень просто: само преобразование берет вектор данных длиной 2 * N и группирует его элементы попарно, после чего для кажой пары вычисляется полусумма и полуразность. Получившийся вектор (опять длины 2 * N), образованный из векторов, содержащих N полусумм и N полуразностей, и является результатом преобразования. При необходимости, само ДВП можно применить повторно (вызвав само преобразование рекурсивно), получив новый набор данных.
Представим, что в роли вектора (набора) входных данных служит дискретная версия некоего сигнала (т.е по сути, набор отсчетов сигнала снятых через некоторые последовательные промежутки времени). Тогда, после применения преобразования Хаара, мы получим вектор выходных данных с той же самой длиной, что и входной вектор, но состоящий из двух равных по размеру частей с различной смысловой нагрузкой.
Первая часть выходного вектора представляет собой набор полусумм из парных отсчетов входного сигнала и эта часть является огрубленной версией исходного сигнала (исходный сигнал «прорежен» в два раза по частоте). Эту часть принято называть аппроксимацией или приближением исходного сигнала.
Вторая часть выходного вектора представляет собой полуразности парных отсчетов сигнала и является набором дополняющей или детализирующей информации, которая нужна для восстановления исходного сигнала. Эту часть принято именовать детализацией исходного сигнала.
Также, иногда можно встретить и иную интерпретацию: аппроксимацию сигнала называют низкочастотной(низкие частоты), а детализацию — высокочастной компонентой(высокие частоты).
Такая удачная схема построения дискретного вейвлет-преобразования, генерирующая данные компоненты, и привела к тому что подобные преобразования естественным образом образуют схему компактизации (сжатия) данных. Этому способствуют и очевидные свойства вейвлета:
- ДВП позволяет разделить набор данных на некое подобие частотных полос (высокие и низкие частоты), которое можно выполнить за n * log2 (n) операций;
- Из полученных компонентов можно исключить наборы самых высокочастотных представлений, тем самым уменьшив количество данных, которые можно сохранить и использовать для воссоздания исходного сигнала;
- Также преобразование выполняет разбиение исходного сигнала на области с различным масштабом, а также выполняет привязку некоторых ключевых особенностей сигнала во временной области.
На этом с теорией все, перейдем к практике: для этого попробуем реализовать самую простую версию вейвлета Хаара для случая с одномерным входным вектором, т.е реализовав вычисление преобразования Хаара для строки данных.
Реализация одномерного ДВП выглядит следующим образом:
// 1D Haar Transform
auto haarTransform1D(Range)(Range range)
{
alias BaseType = ElementType!Range;
auto halfSum(Range)(Range range)
{
return 0.5 * fold!((a, b) => (a + b))(
range,
to!BaseType(0)
);
}
auto halfSub(Range)(Range range)
{
return 0.5 * fold!((a, b) => (a - b))(
range,
to!BaseType(0)
);
}
auto window = range.chunks(2);
auto approximation = window.map!(a => halfSum(a));
auto detail = window.map!(a => halfSub(a));
return chain(approximation, detail);
}
И здесь реализован самый обобщенный алгоритм для ДВП: во-первых, длина входного вектора значений необязательно должна быть четной (т.е длина входного блока данных необязательно 2 * N); во-вторых, набор значений сигнала может быть абсолютно любой структурой, допускающей последовательный перебор в стиле D-шных диапазонов, не требуя каких-то особенных ограничений на тип входных данных; в-третьих, алгоритму абсолютно безразличен тип элементов входного диапазона, главное, чтобы сам диапазон был согласованным (т.е содержал элементы однотипных числовых, либо приводимых к ним, данных).
Реализация начинается с определения одного очень нужного нам типа, а именно, того типа, элементы которого представлены в диапазоне, который попадет в алгоритм (алгоритм, естественно, в исходно D-шном смысле). Для этого используется трайт (trait) из std.range под названием ElementType, который принимает тип диапазона и выводит из него тип его элементов. Полученный тип, подведенный под псевдоним потребовался для реализации двух функций halfSum и halfSub, которые представляют собой алгоритмы вычисления полусуммы и полуразности для некоторого диапазона.
Эти функции построены по одному и тому же принципу — в них использован мощный шаблон fold, осуществляющий свертку диапазона с помощью некоторой функции двух переменных и начального значения (именно под корректное начальное значение нам и потребовался тип, выраженный псевдонимом BaseType). Самое важное в упомянутых функциях то, что им безразличен размер входного диапазона, что приводит к тому, что разбиение в группы для вычисления полусумм/полуразностей может идти таким образом, что не требует формирования пар значений (при нечетном размере входного сигнала последний элемент не имеет пары для образования полусуммы/полуразности — и это особый случай, в котором нужно выполнить простое деление на два).
После того, как были определены функции для вычисления полусумм/полуразностей, дальнейшие операции просты: разбиваем входной поток данных на пары значений с помощью алгоритма chunks и сохраняем результат в переменную. Данная переменная используется для обеспечения независимого вычисления полусумм/полуразностей с помощью алгоритма map и соответствующей функции из описанных выше, после чего результаты вычислений (в форме диапазонов) направляются в переменные approximation (аппроксимация) и detail (детализация).
Итоговый результат получается «склеиванием» диапазонов с помощью алгоритма chain, который обеспечивает генерацию нового диапазона, который по очереди проходит выданные в chain в качестве аргументов диапазоны, «разряжая» их последовательно.
Применяется полученное одномерное преобразование очень просто: подаем алгоритму на вход некий диапазон данных, получаем на выходе из него новый диапазон (уже состоящий из результатов ДВП), который в дальнейшем можно где-нибудь использовать.
Одномерное преобразование Хаара, конечно, прекрасно, но меня интересовало та его версия, которая безошибочно сработает для любого поданного в алгоритма цифрового изображения. Иными словами, мне нужна была реализация двумерного преобразования Хаара.
К счастью, эту реализацию можно легко получить, имея процедуру вычисления одномерного ДВП. Двумерное ДВП (по определению) является композицией одномерных ДВП, выполненных следующим образом: сначала выполняются одномерные преобразования для каждой строки двумерной матрицы входных данных, после чего выполняются одномерные вейвлет-преобразования для каждого столбца полученной матрицы. А поскольку, цифровое изображение является дискретным двумерным сигналом да еще и представленным в форме двумерной матрицы, то все вышесказанное относительно двумерного ДВП также справедливо и для цифровых изображений.
В описанном выше явно вырисовывается алгоритм для двумерного преобразования Хаара, однако, он требует построения или реализации одномерного ДВП для столбца картинки, или способа превращения столбца в строку. И тот, и другой варианты неудобны для реализации, и даже создают некоторые проблемы: поскольку для работы с изображениями проще всего использовать dlib, то реализация ДВП для изображения упирается в возможности этой библиотеки — в частности, для типа SuperImage из dlib не существует методов, выдающих по запросу диапазон строк или столбцов изображения, а меж тем их построение — это не нетривиальная задача…
Именно, поэтому мне пришла в голову мысль, как несколько творчески переработать алгоритм: выполним сначала преобразование Хаара для каждой из строк исходного изображения, после чего транспонируем его как обычную матрицу и выполним опять преобразование для каждой из строк полученного изображения. И тогда, для получения двумерного преобразования останется только применить еще одно траспонирование…
Хотелось мне было испытать transposed из std.range, но тут внезапно вылезла неприятность там, откуда ее меньше всего ждешь: в готовой реализации (уже сделанной цепочки методов) вылезла ошибка, связанная с дизайном самого одномерного преобразования Хаара, а именно, с реализацией вычисления полусуммы/полуразности. Дело в том, что нельзя привести значение типа int к Color4f с помощью стандартных средств D (cast и to), а значение из свойства init дл я типа (если немного поправить для этого шаблоны функций полусуммы/полуразности) дает провальный результат (так инициализируется вот в это: Colo4f(float.nan, float.nan, float.nan))…
И тут я понял, что замечательная идея с обобщенным кодом здесь ничем не поможет, и нужно делать частную реализацию, но основанную на тех же идеях двумерного вейвлета Хаара…
Двумерное вейвлет-преобразование для случая с цифровым изображением было реализовано следующим образом: сначала реализуем версию преобразования для одной строки, сформировав строку из полусумм и полуразностей так, чтобы полусуммы занимали первую половину строки изображения, а полуразности — вторую; полученное изображение (в котором все строки были преобразованы указанным образом) подвергнем транспонированию; а затем просто повторим эти процедуры еще раз.
Для реализации задуманного необходима процедура транспонирования изображения, отсутствующая в dlib, но которую можно легко реализовать следующим образом:
auto transposeImage(SuperImage source)
{
auto target = image(source.height, source.width);
foreach (x; 0..source.width)
{
foreach (y; 0..source.height)
{
target[y, x] = source[x, y];
}
}
return target;
}
Тогда, реализация преобразования Хаара выглядит так:
auto discreteHaarTransform(SuperImage superImage)
{
auto transformRows(SuperImage source)
{
auto width = source.width;
auto height = source.height;
auto target = image(source.width, source.height);
import std.math : SQRT2;
for (int x = 0; x < width; x++)
{
for (int y = 0; y < height; y += 2)
{
Color4f s, d;
if ((y + 1) > height)
{
s = source[x,y];
d = s;
}
else
{
s = source[x, y] + source[x, y + 1];
d = source[x, y] - source[x, y + 1];
}
target[x, y / 2] = s / SQRT2;
target[x, (height + y) / 2] = d / SQRT2;
}
}
return target;
}
alias transformImage = pipe!(
transformRows,
transposeImage,
transformRows,
transposeImage
);
return transformImage(superImage);
}
Фактически, код реализует то, что описано выше относительно алгоритма, но с некоторыми деталями.
Во-первых, вычисляем мы не полусуммы и не полуразности, а несколько иные значения: делим не на 2, а на корень из двух (он присутствует как константа в модуле std.math). Во-вторых: исходные данные не затираются и требуют обязательного наличия изображения-приемника результата преобразования. В-третьих: версия специфична под конкретный тип (Color4f), но ей безразлично содержит ли строка четное количество элементво или нет (учтен случай с одиночным значением). В-четвертых: реализация самого алгоритма преобразования использует шаблон pipe из std.functional, который последовательно применяет указанные в нем функции к некоторому аргументу (в нашем случае весь алгоритм — это поток преобразований, который внесен под псевдоним).
(И еще: в принципе, можно было убрать повторы преобразований просто свернув одномерное ДВП и транспонирование в отдельную функцию, однако, был важен принцип…)
А теперь возьмем наше любимое изображение Lenna.png и подвергнем его преобразованию:
auto A = load(`Lenna.png`); auto B = A.discreteHaarTransform; B.savePNG(`haar.png`);
В результате чего получаем следующую картинку:
Двумерный вейвлет Хаара дает четыре блока данных (а не два, как одномерный к примеру): первая четверть изображения — это то, что является аппроксимацией; в то время как остальные, являются различными его детализациями — по горизонтали, вертикали и диагонали. Зная тот факт, что одно изображение — огрубленная копия исходной картинки, а все остальное ее приближения (причем, высокочастотные), можно исключить их часть, уменьшив тем самым размер файла с цифровым изображением.
Данное преобразование можно легко обратить, например, так:
auto inverseHaarTransform(SuperImage superImage)
{
auto transformRows(SuperImage source)
{
auto width = source.width;
auto height = source.height;
auto target = image(source.width, source.height);
import std.math : SQRT2;
for (int x = 0; x < width; x++)
{
for (int y = 0; y < height; y += 2)
{
Color4f s = SQRT2 * source[x, y / 2];
Color4f d = SQRT2 * source[x, (height + y) / 2];
if ((y + 1) > height)
{
target[x, y] = s;
}
else
{
target[x, y] = (s + d) / 2;
target[x, y + 1] = (s - d) / 2;
}
}
}
return target;
}
alias transformImage = pipe!(
transformRows,
transposeImage,
transformRows,
transposeImage
);
return transformImage(superImage);
}
И соответственно, можно используя обратное двумерное ДВП, восстановить исходное изображение:
Кроме того, можно провести и более занятные эксперименты: поскольку преобразование Хаара делит изображение на число блоков равное степени двойки, то нетрудно вычислить координаты областей, ответственных за детализацию изображения и «затереть» одну (или более) из них, сделав простую заливку интересующей области черным цветом.
Вот, например, мой эксперимент:
foreach (x; 0..256)
{
foreach (y; 256..512)
{
B[x, y] = Color4f(0.0f, 0.0f, 0.0f);
}
}
B.savePNG("haar_1.png");
B.inverseHaarTransform.savePNG("i_haar_1.png");
Восстановленное изображение:
Кроме того, экспериментируя с данной реализацией ДВП, я заметил, что итоговые файлы имеют меньший размер в файловой системе (процентов на 10 — 30), что явно намекает на то, что можно таким образом сжимать картинки, экономя место на диске. Однако, ценой за такое хранение будет некоторая потеря информации из изображения (приглядитесь внимательно к восстановленной картине Lenna из примера без удаления высокочастотных компонентов и вы меня поймете), что может быть очень критичным для некоторых случаев. Еще мне попадалась информация и о том, что данное преобразования не любит резкие цветовые переходы в изображениях, иногда давая «лестничный эффект» (т.е целые полосы из пикселей разной яркости)…
На этом все, а почитать чуть более подробно про двумерный вейвлет Хаара, вы можете в этой работе.