Про HDR

[eddy_em]

Вот посмотрел я сегодня еще раз те жалкие картинки, которые мне сделал hugin вместо HDR, и подумал, что надо бы все-таки реализовать свою HDR-собиралку (теоретически, т.к. hugin — всего лишь обертка, запускающая консольные утилиты, можно ей вместо hugin_hdrmerge свое подсунуть).
Теорию я уже приводил на ЛОРе, да и частично реализовывал при анализе камеры Apogee. Т.к. здесь работа ведется с png-файлами, все сильно упрощается (правда, надо где-то сохранять сведения об экспозиции каждого кадра, иначе ничего не получится): открываем одновременно все файлы, затем построчно считываем их, строка за строкой заполняя результирующее HDR-изображение. Обработка элементарная: отбрасываем значения ЧБ-интенсивности, близкие к 0 или 255; далее делим оставшиеся интенсивности на экспозиции; вычисляем медиану; сохраняем логарифм полученного значения во временный файл (лучше, пожалуй, double использовать), параллельно вычисляя максимальную и минимальную величину логарифма. В другой или этот же файл сохраняем данные о цвете, которые берем из того изображения, где уровень серого в данном пикселе лежит где-нибудь так в диапазоне 176..230 (надо экспериментально подобрать).
По окончании обработки всего изображения преобразуем наш бинарник в нормальный png-файл.
Такой способ обработки все шумы максимально уменьшит, кроме шумов на уж очень темных местах (т.е. на темных участках изображения с максимальной экспозицией). Их придется восстанавливать, предварительно пройдясь медианным фильтром.

В общем, если не забуду об этой идее, как-нибудь от скуки реализую.
Для демонстрации я набросал простенький код. Функция proc_images(im) открывает изображение im, генерирует на его основе "кадры с разными экспозициями", а затем получает нечто, вроде HDR (еще допиливать надо):

function [II, BIN] = proc_images(Imname)
	II = imread(Imname);
	I0 = 64; I1 = 192;
	W = size(II, 2);
	H = size(II, 1);
	I = []; G = [];
	N = [1:7];
	T = 2.^(N-N(ceil(size(N,2)/2)));
	for n = N
		[I_t G_t] = mk_im(II, T(n));
		I_t(:,:,1) = medfilt2(I_t(:,:,1), "replicate");
		I_t(:,:,2) = medfilt2(I_t(:,:,2), "replicate");
		I_t(:,:,3) = medfilt2(I_t(:,:,3), "replicate");
		G_t(find(G_t < I0)) = 0;
		G_t(find(G_t > I1)) = 255;
		I(:,:,:,n) = I_t;
		G(:,:,n) = G_t;
		clear I_t G_t
	endfor
	clear II;
	BIN = [];
	Imax = 0;
	Imin = 1e10;
	for y = 1:H
		fprintf(1, "Y = %d\n", y);
		fflush(1);
		for x = 1:W
			PT = [];
			nn = [];
			for n = N
				g = G(y,x,n);
				if(g >= I0 && g <= I1)
					PT = [PT g];
					nn = [nn n];
				endif
			endfor
			if(size(PT) == 0)
				if(g < I0)
					%colr = [0 0 0];
					colr = I(y,x,:,end) / T(end);
				else
					%i = 255 / T(1);
					%colr = [i i i];
					colr = I(y,x,:,1) / T(1);
				endif
			else
				%i = median(PT./T(nn));
				PT = PT./T(nn);
				[ii s] = sort(PT./T(nn));
				S = ceil(size(s,2)/2);
				i = PT(S); S = nn(S);
				colr = I(y,x,:,S) / T(S);
				
			endif
			BIN(y,x,:) = colr;
			im = max(colr);
			if(Imax < im) Imax = im; endif
			im = min(colr);
			if(Imin > im) Imin = im; endif
		endfor
	endfor
	II1 = uint8((BIN  - Imin) / (Imax - Imin) * 255);
	Imin = log(Imin + 1); Imax = log(Imax + 1);
	II = uint8((log(BIN + 1) - Imin) / (Imax - Imin) * 255);
	imwrite(II, "result_log.png");
	imwrite(II1, "result_lin.png");
endfunction

Вспомогательная фукнция mk_im(i, t) из изображения i (экспозиция которого считается равной 1с) получает изображение с экспозицией t:

function [I, G] = mk_im(I0, T)
	I = uint8(I0 * T);
	name = sprintf("exp_%.3f.png", T);
	imwrite(I, name);
	G = rgb2gray(I);
	name = sprintf("exp_%.3f_gray.png", T);
	imwrite(G, name);
endfunction

И вот что получается:

Оригинал — кусок фотографии для тренировок.

Кадры с разной "экспозицией".

Соответствующие уровни серого.

Результат.

Из-за медианного усреднения результат несколько размыт, но это не беда.
С логарифмом у меня вышел "косячок": т.к. после логарифмирования "внизу" получается очень мало значений, итоговая гистограмма съезжает в яркую область. Надо подумать, как покрасивше делать.

Естественно, с настоящими фото интересней поработать, но это надо на сях писать обработчик, т.к. октава уж очень тормозит...