В этой статье мы расскажем вам о том, как можно своими руками реализовать один из методов улучшения качества изображений с помощью усиления контраста отдельных областей. Данный метод называется эквализацией (выравниванием) гистограммы и его реализацию мы покажем далее.
Эквализация гистограммы — это один самых простых методов улучшения качества изображения, однако, в литературе и интернете сложно найти пример какой-либо простой реализации без погружения в смутные детали. Эквализация, это просто процедура выравнивания гистограммы изображения, путем воздействия (т.е. коррекции) яркости отдельных пикселей.
Дело в том, что гистограмма произвольного изображения представляет собой график, отображающий в виде пиков, количество пикселей в изображении с определенной яркостью, и, как правило, для некоторого изображения она представляет собой множество пиков неравномерно распределенных по графику.
Неравномерное распределение, в данном случае, означает что в гистограмме есть области в которых сосредоточена большая часть пиков или наоборот меньшая, т.е. в пространстве графика присутствуют участки, где плотность значений выше/ниже, чем в среднем по графику. Процедура эквализации сделает распределение значений в гистограмме более равномерным, а это значит что в ней почти не останется пропусков и областей с чрезмерно высоким количеством (и высотой) пиков.
Возможно, подробное объяснение непонятно, однако, я предлагаю вам взглянуть на код примера, сделанный с использованием библиотеки ppmformats:
import std.conv;
import std.stdio;
import ppmformats;
auto histogram(PixMapFile pmf)
{
float[256] histogram = 0;
foreach (e; pmf.array)
{
auto I = cast(int) e.luminance;
histogram[I]++;
}
return histogram;
}
auto equalize(PixMapFile pmf)
{
import std.algorithm : filter, map, min, reduce, sum;
import std.range : iota;
auto histogram = pmf.histogram;
auto cdf(float x)
{
return histogram[0..x.to!int].sum;
}
auto m = iota(0, 256)
.map!(a => cdf(a))
.filter!(a => (a > 0))
.reduce!min;
auto S = pmf.width * pmf.height;
foreach (i,e; pmf.array)
{
auto I = e.luminance;
auto X = cdf(I);
float L = ((X - m) / (S - 1)) * 255;
float NI = L / I;
pmf[i] = e * NI;
}
}
auto histogramToPPM(PixMapFile pmf, int width = 256, int height = 128)
{
import std.algorithm : map;
import std.range : array;
auto S = pmf.width * pmf.height;
auto hist = pmf
.histogram[]
.map!(a => a / S)
.array;
auto img = new P6Image(width, height, new RGBColor(255, 255, 255));
foreach (i, e; hist)
{
auto H = cast(int) (100 * e * (height >> 2));
foreach (j; 0..H)
{
img[i, height-j] = new RGBColor(0,0,0);
}
}
return img;
}
void main()
{
auto pmf = new P6Image;
pmf.load(`/home/aquareji/D/Lenna.ppm`);
pmf.histogramToPPM.save(`Lenna_before.ppm`);
pmf.equalize;
pmf.histogramToPPM.save(`Lenna_after.ppm`);
pmf.save(`/home/aquareji/D/Lenna_equalize.ppm`);
}Результат применения скрипта к уже всем приевшемуся изображению Лены:
А вот оригинальное изображение для сравнения:
Спросите, а причем тут гистограмма?
Для ответа на этот вопрос была создана вспомогательная функция histogramToPPM, которая генерирует из некоторого изображения рисунок его гистограммы (в виде PPM P6) и позволит увидеть разницу между двумя изображениями. До процедуры эквализации в изображении, содержащем график яркостей (т.е. гистограмму) будут присутствовать множественные разрывы и большие области, в которых сосредоточена большая часть ярких участков изображения, вот гистограмма оригинального изображения:
А вот для сравнения выровненная гистограмма, которая получена после процедуры эквализации:
Видно, что график стал более равномерным, а локальные области максимальной плотности ярких участков исчезли и гистограмма стала более равномерной. В свою очередь, это привело к тому, что контраст некоторых областей стал сильнее и динамический диапазон изображения стал более насыщенным, а это означает, что качество изображения стало несколько выше.
Теперь давайте приступим к разбору реализации алгоритма эквализации. Прежде всего, большое спасибо автору статьи о реализации данного алгоритма на F#. Объяснения этого автора оказались очень качественными, а некоторые идеи мы применили в своей реализации алгоритма.
Но, возвращаемся к коду.
Процедура equalize использует функцию histogram, которая возвращает массив float с фиксированной размерностью (а, именно — 256, поскольку яркости пикселей находятся в диапазоне от 0 до 255, и их всего 256) и принимает в качестве аргумента универсальный тип изображений в ppmformats под названием PixMapFile. Функция histogram в самом начале своей работы присваивает всем элементам одноименного массива 0 (поскольку иначе не сработает дальнейшая схема подсчета яркостей), а затем проходит по всем пикселям изображения, вычисляет яркость каждого, приводит к int и использует приведенное значение в качестве индекса массива, увеличивая таким образом значение в массиве яркостей на 1. Весь цикл при этом приведет к тому, что из-за ограниченности целого значения яркостей, каждый пиксель внесет свое значение яркости в копилку массива гистограммы.
Примечание автора. Если вначале функции histogram не присвоить 0 фиксированному массиву с тем же именем, то он будет инициализирован значением float.nan, и прибавление каждый раз 1 к этом значению приведет к тому, что весь массив окажется заполнен только значениями float.nan. За подробными объяснениями почему это так, лучше обратиться к спецификации IEEE 754 относительно формата float.
Процедура equalize сначала получает гистограмму и сохраняет ее в переменную histogram.
Примечание автора.Обращаем ваше внимание на то, что здесь точно также как и в функции histogram, в качестве типа входного аргумента использован тип PixMapFile и, соответственно, процедура не возвращает никакого значения, а лишь модифицирует входной аргумент.
Эта переменная нужна для того, чтобы рассчитать функцию распределения для изображения, которая для некоторой яркости I представляет собой сумму всех яркостей в гистограмме histogram от нулевой до I. Соответственно, вычисление функции распределения реализовано по указанной схеме с помощью внутренней функции cdf, принимающей в качестве аргумента некоторую яркость x, приводимую к int. Процедура sum, использованная для суммирования всего среза в cdf, взята из std.algorithm, как указано в импортах.
После определения функции распределения необходимо вычислить ее минимально возможное значение отличное от нуля. Для этого была определена переменная m, в которую помещается вычисленное значение ненулевого минимума.
Само вычисление реализовано следующим образом: поскольку область определения функции распределения лежит в интервале от 0 до 256 (не включая последнее значение), то с помощью iota генерируется весь диапазон входных значений, к которому далее применяется map совместно с cdf (это обеспечивает генерацию всех возможных значений cdf). Данный диапазон значений фильтруется с помощью filter на предмет ненулевых значений и таким образом избавленный от нулей поступает на вход функции reduce c шаблонным параметром в виде функции min, находящей минимум для двух значений. Поскольку reduce реализует поэлементную свертку диапазона (в данном случае, попарную) с помощью некоторой функции, то применение min здесь дает минимальное значение из всей выборки, и именно оно сохраняется в переменную m.
Дальнейшая часть алгоритма предельно простая и для ее работы необходима переменная S, которая содержит количество всех пикселей в изображении (длина, умноженная на ширину). После вычисления площади изображения, необходимо скорректировать яркость каждого пикселя в отдельности для чего вычисляется его яркость I, значение функции распределения для данной яркости X, а также новая яркость L, вычисление которой представляет реализацию формулы из статьи с HabrHabr, так сказать, «в лоб».
На вычислении трех этих значений сходство с упоминающейся статьей заканчивается, и начинается зыбкая область смутных догадок: для коррекции пикселя потребуется вычислить во сколько раз отличаются его новая (L) и исходная яркости (I) и соответственно увеличить в это количество яркость текущего пикселя, что реализовано банальным умножением пикселя на получившийся коэффициент коррекции (NI).
Примечание автора.Данная идея работает отчасти потому, что умножение/деление пикселя на некое число приводит к изменению его яркости: умножение ее увеличивает, а деление — уменьшает. Однако, мы используем умножение, т.к. значение NI большее 1 эквивалентно увеличению яркости (работает так, будто у нас обычное умножение с повышением яркости), а меньшее 1 — делению (обычный математический факт для чисел).
Теперь, мы перейдем к описанию вспомогательной функции histogramToPPm, которая принимает в качестве аргументов общий тип изображений PixMapFile, длину и ширину графического представления гистограммы, которое будет представлено типом P6Image.
Сначала, мы рассчитываем площадь изображения (т.е. общее количество в нем точек-пикселей), затем получаем и несколько обрабатываем гистограмму для того, чтобы вместо количества точек с конкретной яркостью, у нас было процентное отношение количества точек с данной яркостью ко всему изображению: для этого используется map с обычным делением значения на общее количество точек и array, для превращение диапазона в обычный массив.
Примечание автора.Здесь применен небольшой лайфхак, про который вам нигде и никогда не расскажут. Заключается он в том, что к функции, дающей некоторый массив (особенно это касается массивов с фиксированной размерностью) добавлены квадратные скобки, говорящие о том, что берется срез размером с весь массив и его тип — массив того же типа, что и был, но уже с динамической размерностью! Дело в том, что некоторые функции стандартной библиотеки не любят разного рода массивы и предпочитают вместо них диапазоны, а срез, по мнению D, также является входным диапазоном. Таким образом, если вдруг функция из стандартной библиотеки не принимает массив, то просто добавьте квадратные скобки после аргумента — и все будет как надо.
После предварительной подготовки массива-гистограммы, создается новое изображение формата P6 с длиной width, шириной height и в качестве заполнителя используется белый цвет (по умолчанию, в ppmformats используется черный). С полученным изображением производиться отрисовка столбцов гистограммы, для чего происходит проход по всем значениям гистограммы: с помощью текущего ее значения e вычисляется высота столбца H для отрисовки пиков гистограммы.
Высота H считается из того факта, что значения в массиве представлены процентным отношением, выраженном в долях единицы, и предположения о том, что высота занимает какой-то процент от половины ширины изображения (половина выражена сдвигом на 2, что позволяет избавиться от операции деления). Внутренний цикл с переменной j выполняет отрисовку отдельного столбца, просто заполняя пиксели, соответствующие высоте H на картинке черным цветом, рисуя тем самым одну черную линию. Таким образом и получается отрисовка всех столбцов гистограммы.
Примечание редактора.Хотя длина и ширина изображения гистограммы настраиваются, но все же лучше будет, если длина гистограммы останется равной 256 (по количеству возможных яркостей). Кроме того, лишний раз напоминаем о том, что вам потребуется сменить пути в примере, если захотите тестировать пример на своих файлах.
Ради интереса (а также для того, чтобы поместить результирующее изображение и графики гистограмм в статью), мы переписали пример с использованием dlib вместо ppmformats:
#!/usr/bin/env dub
/+ dub.sdl:
name "exp"
dependency "dlib" version="~>0.17.0"
+/
//
import dlib.image;
auto histogram(SuperImage simg)
{
float[256] histogram = 0;
foreach (i; 0..simg.width)
{
foreach (j; 0..simg.height)
{
auto I = cast(int) (simg[i, j].luminance * 255.0);
histogram[I]++;
}
}
return histogram;
}
auto histogramToPNG(SuperImage simg, int width = 256, int height = 128)
{
import std.algorithm : map;
import std.range : array;
auto S = simg.width * simg.height;
auto hist = simg
.histogram[]
.map!(a => a / S)
.array;
auto img = image(width, height);
foreach (i; 0..width)
{
foreach (j; 0..height)
{
img[i, j] = Color4f(1.0f, 1.0f, 1.0f);
}
}
foreach (i, e; hist)
{
auto H = cast(int) (100 * e * (height >> 2));
foreach (j; 0..H)
{
img[cast(int) i, height-j] = Color4f(0.0f, 0.0f, 0.0f);
}
}
return img;
}
auto equalize(SuperImage simg)
{
import std.algorithm : filter, map, min, reduce, sum;
import std.conv : to;
import std.range : iota;
auto histogram = simg.histogram;
auto cdf(float x)
{
return histogram[0..x.to!int].sum;
}
auto m = iota(0, 256)
.map!(a => cdf(a))
.filter!(a => (a > 0))
.reduce!min;
auto S = simg.width * simg.height;
foreach (i; 0..simg.width)
{
foreach (j; 0..simg.height)
{
auto I = simg[i, j].luminance;
auto X = cdf(I * 255);
float L = ((X - m) / (S - 1));
float NI = L / I;
simg[i, j] = simg[i, j] * NI;
}
}
}
void main()
{
auto img = load(`/home/aquareji/Downloads/Lenna.png`);
img.histogramToPNG.savePNG(`Lenna_before.png`);
img.equalize;
img.histogramToPNG.save(`Lenna_after.png`);
img.savePNG(`Lenna_equalized.png`);
}Как видите, код почти такой же, за исключением пары моментов: в функции получения гистограмм histogram, а также внутри функции equalize получаемые значения из luminance умножаются на 255 (поскольку, в dlib яркость может принимать значения от 0.0f до 1.0f); функция histogramToPPM заменена на аналог histogramToPNG, а тип цвета сменился с RGBColor на Color4f. Кроме того, в предыдущем случае, мы запускали приложение с помощью rdmd (при этом неявно предполагается, что файл примера находиться в той же папке, что и файл ppmformats.d, взятый с гитхаба разработчика), но в этом случае используется dub с аргументом —single в опции run, как описывалось в одной из наших статей.
В заключение остается только сказать, что теперь в ваших руках оказался один из инструментов улучшения изображений и мы от имени блога LightHouse Software желаем вам успеха в использовании и улучшении данного инструмента.