Вот ссылка, содержащая изображения и скрипты, которые я описал здесь:
Контраст — это различие в цвете или люминесценции изображения, которое делает предметы на нем различимыми. Давайте посмотрим на два набора изображений, одно в оттенках серого и одно в цвете, где изображения в первом наборе имеют низкую контрастность, а изображения в другом — высококонтрастные изображения.
Что мы видим, когда сравниваем эти два набора изображений? На высококонтрастных изображениях мы видим, что цвета намного богаче, свет местами переходит от очень темного к очень светлому. Если мы, в свою очередь, посмотрим на низкоконтрастные изображения, то увидим, что диапазон цветов намного меньше, разные объекты, как правило, имеют очень похожие цвета. То же самое можно заметить, если перевести все изображения в серый цвет, а затем посмотреть на их гистограммы.
Глядя на гистограммы, ясно видно, как пиксели на низкоконтрастных изображениях принимают значения из гораздо меньших наборов значений, из чего следует, что объекты принимают цвета из очень узких наборов значений, что в конечном итоге означает, что разные объекты труднее различить.
Самое простое, что можно сделать для повышения контрастности изображения, — убедиться, что пиксели отображаются на полный набор доступных нам значений (в данном случае 0–255, так как мы работаем с 8-битными изображениями). Мы также хотим сохранить отношения между пикселями — то есть, если пиксель (i,j) темнее пикселя (k,l), мы хотим, чтобы он оставался таким после преобразования.
Посмотрим на гистограмму изображения моста. Мы видим, что диапазон значений пикселей составляет примерно от 70 до 245. Если мы хотим сопоставить их с набором значений от 0 до 255 и сохранить взаимосвязь между ними, нам нужно использовать линейную функцию, которая отображает 70 в 0. , от 245 до 255 и является линейным. Matlab предоставляет функцию, которая выполняет это преобразование, и она называется imadjust.
Давайте посмотрим на среднюю форму функции, на ту, где гамма равна 1. Это именно та функция, которую я описал выше. Давайте применим это к первым двум изображениям и посмотрим на результаты.
bridge_adjusted = imadjust(bridge, [70/255 200/255], [0 1], 1); forest_adjusted = imadjust(forest, [90/255 220/255], [0 1], 1);
Обратите внимание, что эта функция предполагает, что пределы входа и выхода находятся в диапазоне от 0 до 1.
Как и ожидалось, форма гистограммы осталась прежней, просто теперь она растянута.
Вернемся к функции imadjust. Что с ним происходит, когда мы увеличиваем или уменьшаем гамму? Когда для гаммы установлено значение 1, первая половина значений отображается на первую половину нового диапазона. По мере того, как мы уменьшаем гамму, все больше и больше пикселей в конечном итоге отображаются в верхний диапазон новых значений. Как будто мы берем все большую и большую часть, начиная с самых высоких значений исходной гистограммы, и сжимаем их все больше и больше, в то же время растягивая все остальные во все большем и большем диапазоне. Обратное происходит, когда мы увеличиваем гамму — большие значения в исходной гистограмме в конечном итоге все больше и больше растягиваются, в то время как другие сжимаются. Чем больше гамма отдаляется от 1, тем больше удаляется от середины воображаемая линия, отделяющая части гистограммы, растягиваемые от сжимаемой. Когда гамма меньше 1, мы растягиваем меньшие значения, в то время как для гаммы больше 1 происходит обратное.
Давайте посмотрим, что происходит с изображениями, когда мы берем разные значения гаммы.
Здесь происходит что-то очень интересное. Установка гаммы на 1,5 полностью испортила контрастность изображения моста, но с изображением леса мы получили результаты, которые даже лучше, чем при установке гаммы на 1. Так почему же это происходит? Обратите внимание, что когда мы преобразовываем пиксели, мы всегда сопоставляем их с полным набором значений, но мы меняем форму. Таким образом, не только диапазон должен быть полным, чтобы иметь высокую контрастность — форма получаемой гистограммы также очень важна. Взгляните на гистограммы и обратите внимание, чем больше гистограмма преобразованного изображения больше похожа на универсальную функцию, тем больше изображение выглядит более контрастным. То есть, помимо взятия полного диапазона значений, пиксели должны быть как можно более равномерно распределены, чтобы увеличить контрастность.
И вот как мы добираемся до выравнивания гистограммы. Вместо использования функции imadjust и угадывания значений гаммы до тех пор, пока мы не получим удовлетворительные результаты, мы можем использовать кумулятивную сумму (нормализованной) гистограммы в качестве передаточной функции. Как и imadjust, эта функция монотонна и будет равна 0 для значений, которых нет на изображении, и 1 для значений, превышающих максимальное значение пикселей в исходном изображении. За исключением того, что в этом случае гистограмма, которую мы собираемся получить, будет наиболее однородной. Так почему же? Потому что мы по существу используем интегральное преобразование вероятности:
Вот код Matlab, который выполняет это преобразование на 8-битном изображении:
h = imhist(I); hnorm = h./numel(I); % Normalize histogram. cdf = cumsum(hnorm); % Find cumulative sum. T = uint8(floor(cdf*255)); % Create transform function. I_equalized = inlut(I,T); % Equalize histogram.
Обратите внимание, что поскольку наша передаточная функция не является непрерывной, нам пришлось использовать справочную таблицу. Конечно, есть функция Matlab histeq, которая делает это, но это для того, чтобы показать, как работает алгоритм.
Итак, теперь понятно, почему иногда imadjust уменьшает, а иногда увеличивает контраст при использовании одних и тех же параметров на разных изображениях — когда это работает, это потому, что это похоже на кумулятивную функцию распределения. Давайте еще раз посмотрим на общие результаты, чтобы понять, как обстоит дело.
Обратите внимание, что чем более прямолинейно кумулятивное распределение, тем лучше мы получаем контраст.
Рассмотрим компонент Y изображения леса в цветовом пространстве YCbCr. Я проделаю те же операции, но на этот раз только с люминесцентным компонентом. Я буду использовать функцию Matlab stretchlim, чтобы найти значения high_in и low_in для функции imadjust.
forest_ycbcr = rgb2ycbcr(forest); low_high = stretchlim(forest_ycbcr(:,:,1), [0.05]); forest_1(:,:,1) = imadjust(forest_ycbcr(:,:,1), low_high, [0 1], 1); forest_15(:,:,1) = imadjust(forest_ycbcr(:,:,1), low_high, [0 1], 1.5); forest_05(:,:,1) = imadjust(forest_ycbcr(:,:,1), low_high, [0 1], 0.5); forest_eq(:,:,1) = histeq(forest_ycbcr(:,:,1));
Итак, мы нашли функцию, увеличивающую контраст в общем случае. Но это всего лишь улучшенная версия функции imadjust, она имеет те же общие недостатки. Основная проблема с гистограммным эквалайзером заключается в том, что он на самом деле не различает разные области изображения. Обратите внимание, как мы теряем некоторые детали в выровненном изображении слонов — например, некоторые темные части зубов исчезли, но на исходном изображении они были в порядке. Например, пиксели, окруженные темными пикселями, можно сделать еще ярче, чтобы повысить контрастность. То есть нам не нужно делать самую яркую часть всего изображения белой, мы можем сделать самую яркую часть области белой, чтобы еще больше увеличить контрастность. Эта логика используется в адаптивном алгоритме выравнивания гистограммы.
Проходим по каждому пикселю и выполняем выравнивание гистограммы с передаточной функцией, соответствующей блоку, в центре которого находится данный пиксель. Таким образом, значение, которому соответствует пиксельная карта, определяется только близкими к нему пикселями, а не всем изображением. Это означает, что если изображение имеет разный свет в некоторых областях, это не испортит весь контраст. Наконец, есть пиксели, для которых мы не сможем смотреть на окружающий блок, потому что они находятся на краю. Я буду размещать изображение симметрично с каждой стороны, так что даже крайние пиксели будут иметь полные блоки вокруг них.
Вот код, который выполняет это.
% Read image and convert it to YCbCr color space.
elephants = imread('ex3.jpg');
elephants_ycbcr = rgb2ycbcr(elephants);
% Pick block size such that there are approx 32 blocks.
[M N] = size(elephants(:,:,1));
block_size = round(sqrt(M*N/32));
%Pad image on each side by block_size/2. Mirror pixels.
pad_val = round(block_size/2);
I = padarray(elephants_ycbcr(:,:,1),[pad_val pad_val],'symmetric','both');
clip_limit = 0.02;
equalized = I;
% Go through each pixel.
for i=pad_val+1:M+pad_val
for j=pad_val+1:N+pad_val
% Find surrounding block for this pixel.
block = I(i-pad_val:i+pad_val, j-pad_val:j+pad_val);
hist = imhist(block)/numel(block(:));
cdf = uint8(round(255*cumsum(hist)));
% Map value based on cdf.
equalized(i,j) = intlut(I(i,j), cdf);
end
end
%Crop padding and convert to rgb.
elephants_ycbcr(:,:,1) = equalized(pad_val+1:M+pad_val, pad_val+1:N+pad_val);
elephants_eq = ycbcr2rgb(elephants_ycbcr);
imshow(elephants_eq);
С одной стороны, мы увеличили контраст намного сильнее, чем при выравнивании гистограммы, но мы добавили столько шума, что это не оправдывает себя. Шум возникает из-за того, что в частях областей постоянного цвета мы пытаемся отображать значения в весь диапазон, а там это действительно не нужно.
Эту проблему можно решить с помощью так называемого ограничения контраста. Причина этого в том, что если это область постоянного цвета, не меняйте слишком сильно значение исходного пикселя.
Давайте посмотрим на гистограмму для области постоянного цвета.
Как видите, большинство пикселей сгруппированы вокруг одного центрального значения. Это означает, что кумулятивная сумма (которую мы используем в качестве передаточной функции) будет иметь очень крутой наклон и в основном будет равна 0 и 1 не более чем в нижнем и верхнем диапазоне значений пикселей. Когда такое кумулятивное распределение используется в качестве передаточной функции, оно отображает значения, которые очень близки друг к другу, в значения, которые находятся далеко друг от друга. Мы не хотим, чтобы это произошло.
Способ остановить это — взять некоторое значение отсечки гистограммы, и все, что выше этого, равномерно распределяется. Это приведет к тому, что передаточная функция будет иметь наклон, очень близкий к 1, что означает, что значение x тесно связано с x.
Давайте добавим это в код и посмотрим на результат.
% Read image and convert it to YCbCr color space.
elephants = imread('ex3.jpg');
elephants_ycbcr = rgb2ycbcr(elephants);
% Pick block size such that there are approx 32 blocks.
[M N] = size(elephants(:,:,1));
block_size = round(sqrt(M*N/32));
%Pad image on each side by block_size/2. Mirror pixels.
pad_val = round(block_size/2);
I = padarray(elephants_ycbcr(:,:,1),[pad_val pad_val],'symmetric','both');
clip_limit = 0.02;
equalized = I;
% Go through each pixel.
for i=pad_val+1:M+pad_val
for j=pad_val+1:N+pad_val
% Find surrounding block for this pixel.
block = I(i-pad_val:i+pad_val, j-pad_val:j+pad_val);
hist = imhist(block)/numel(block(:));
% Contrast limit.
hist(hist>clip_limit) = clip_limit;
delta = (1 - sum(hist(:)))/256;
hist = hist + delta;
cdf = uint8(round(255*cumsum(hist)));
% Map value based on cdf.
equalized(i,j) = intlut(I(i,j), cdf);
end
end
%Crop padding and convert to rgb.
elephants_ycbcr(:,:,1) = equalized(pad_val+1:M+pad_val, pad_val+1:N+pad_val);
elephants_eq = ycbcr2rgb(elephants_ycbcr);
imshow(elephants_eq);
Итак, мы позаботились и об этом. Но теперь есть другая проблема. Эта программа смехотворно медленная. И неудивительно, выполнять выравнивание гистограммы на каждом пикселе действительно много. Вы могли бы ускорить его, перейдя на C, C++ и написав код более разумно (например, используя скольжение окна и вычисляя кумулятивную сумму от начала или конца и останавливаясь, когда вы добираетесь до пикселя, для которого вам нужно преобразование), но это все равно выиграло. не будет намного быстрее, так как алгоритм действительно сложный.
Другой способ значительно ускорить это — использовать адаптивное выравнивание гистограммы с интерполяцией. Логика в том, что мы делим изображение на блоки и находим кумулятивные распределения для каждого блока. Затем мы просматриваем каждый пиксель и делаем сопоставление путем интерполяции сопоставления на основе окружающих блоков.
Вот код:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This progam is an example of adaptive histogram equalization
% algorithm with interpolation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear all
% Read image and convert it to YCbCr color space.
elephants = imread('ex3.jpg');
elephants_ycbcr = rgb2ycbcr(elephants);
clip_limit = 0.02;
number_of_blocks = 32;
% Pick block size such that there are approx number_of_blocks blocks.
[M N] = size(elephants(:,:,1));
block_size = round(sqrt(M*N/number_of_blocks));
block_size = block_size + ~mod(block_size, 2);
% Split image into blocks. Find cdf with contrast limit for each block.
block_rows = ceil(M/block_size);
block_columns = ceil(N/block_size);
mapings = zeros(block_rows, block_columns, 256);
% Add symmetricaly rows and columns so that all blocks are of the same size.
pad_rows = block_size*block_rows - M;
pad_columns = block_size*block_columns - N;
I = padarray(elephants_ycbcr(:,:,1),[pad_rows pad_columns],'symmetric','post');
equalized = I;
% Find cdf for each block.
for i=0:block_rows-1
for j=0:block_columns-1
block = I(i*block_size+1:(i+1)*block_size, j*block_size+1:(j+1)*block_size);
hist = imhist(block)/numel(block(:));
% Contrast limit.
hist(hist>clip_limit) = clip_limit;
delta = (1 - sum(hist(:)))/256;
hist = hist + delta;
% Save cdf for this block.
mapings(i+1, j+1, :) = uint8(round(255*cumsum(hist)));
end
end
% Set up center pixels for each block.
center_rows = zeros(block_rows, 1);
center_columns = zeros(block_columns,1);
for i=0:block_rows-1
center_rows(i+1) = i*block_size + round(block_size/2);
end
for i=0:block_columns-1
center_columns(i+1) = i*block_size + round(block_size/2);
end
% Go through each pixel and interpolate its maping.
for i=1:M+pad_rows
for j=1:N+pad_columns
% Find block indexes.
block_down = find(center_rows >= i, 1, 'first');
block_up = find(center_rows <= i, 1, 'last');
block_right = find(center_columns >= j, 1, 'first');
block_left = find(center_columns <= j, 1, 'last');
%Edge cases.
if (isempty(block_left)) block_left = 1; end
if (isempty(block_up)) block_up = 1; end
if (isempty(block_right)) block_right = block_columns; end
if (isempty(block_down)) block_down = block_rows; end
% Find neighbour blocks centers for this pixel value.
x_right = center_columns(block_right);
x_left = center_columns(block_left);
y_down = center_rows(block_down);
y_up = center_rows(block_up);
% Find weight coefficents.
if (y_up == y_down)
a = 0.5;
else
a = (i-y_up)/(y_down-y_up);
end
if (x_right == x_left)
b = 0.5;
else
b = (j-x_left)/(x_right-x_left);
end
% Interpolate maping.
m_ul = mapings(block_up, block_left, I(i,j));
m_ur = mapings(block_up, block_right, I(i,j));
m_dl = mapings(block_down, block_left, I(i,j));
m_dr = mapings(block_down, block_right, I(i,j));
m = (1-a)*((1-b)*m_ul+(b)*m_ur)+(a)*((1-b)*m_dl+(b)*m_dr);
% Take that value.
equalized(i,j) = round(m);
end
end
%Crop padding and convert to rgb.
elephants_ycbcr(:,:,1) = equalized(1:M, 1:N);
elephants_eq = ycbcr2rgb(elephants_ycbcr);
figure
imshow(elephants_eq);
Как видите, мы разбили изображение на блоки, и каждый блок имеет свой центральный пиксель. Каждый блок также имеет собственное отображение. Мы берем пиксель и находим ближайшие левый, правый, верхний и нижний центральные пиксели. Если там нет пикселя (например, наш пиксель находится слева от первого левого центрального пикселя), мы снова берем тот, что на краю. По сути, в крайних случаях мы делаем интерполяцию из одного или двух блоков вместо четырех. Наконец, мы берем сопоставленное значение выбранного пикселя из каждого преобразования соседних блоков и берем общее значение в зависимости от того, насколько близко наш пиксель находится к центру блока.
Итак, вы видите, что результат не изменился, но если вы запустили программу, то смогли увидеть, насколько велика разница в скорости.
