Еще о HDR, теперь с графиками

[eddy_em]

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

Уровни интенсивности оригинального ("Original") и восстановленного по "фотографиям" с разными экспозициями ("HDR") сигнала.

Итак, смоделируем простой сигнал, имеющий большой разброс по амплитуде + мелкие детали как в области высоких интенсивностей, так и в области низких. Сигнал представлен на рисунке выше. Будем считать, что интенсивность его измеряется в фотонах на пиксель камеры в секунду. Также будем считать, что регистрируются кадры идеальной камерой с квантовой эффективностью, равной единице, у которой отсутствуют все шумы (шумы добавлять лень, т.к. код раздуется еще на медианную фильтрацию). Для простоты изображение берем сразу монохромным, без разбиения на RGB компоненты.
Процесс регистрации за счет ограниченного динамического диапазона камеры приводит к искажению сигнала. Если мы разделим величину в ADU зарегистрированного сигнала на время экспозиции, получим следующее:

Отклонения восстановленного зарегистрированного сигнала от оригинала.

Видно, что из-за переэкспозиции изображения с бóльшей выдержкой имеют бóльшие отклонения от оригинала в области высокой яркости. И наоборот, в области малой яркости отклонения тем больше, чем ниже экспозиция:

Область вблизи малых интенсивностей.

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

Если просто разделить интенсивность в каждом кадре на время экспозиции, из-за большого количества "засвеченных" областей мы получим неправильное восстановление сигнала в ярких областях. Моя моделька не очень хороша, на реальных изображениях эффект должен быть более значимым, а у меня получилось вот так:

Восстановление медианой по всем изображениям. На рисунке — отклонение (в процентах) от оригинала.

Если же мы отбросим данные с высоким или низким уровнем сигнала, то медиана оставшихся величин позволит более точно восстановить сигнал:

Отклонение восстановленного HDR от оригинала (в процентах).

Восстановленный сигнал мы можем тем или иным способом привести к ограниченному диапазону jpg или png изображения. Это можно сделать линейно, или же использовав функциональную зависимость:

Разные виды квантования восстановленного сигнала.

Как я уже говорил в предыдущей записке, логарифм сильно уродует результирующую гистограмму, смещая ее в область высокой яркости; но при этом он еще и акцентирует слабые детали в области низкой яркости. Квадратный корень более оптимально преобразует гистограмму.
Ну и напоследок код, который использовался при написании этой заметки:

function graph()
	X = 1:200;
	I0 = 64; I1 = 192;
	Y = 30*sin(X).^2+150*sin(X/11).^2+1000*cos(X/41).^2;
	plot(X,Y)
	print -dpng Signal.png
	T = 2 .^ [-5:1];
	I = [];
	L = [];
	ii = [1:size(T,2)];
	for i = ii
		I(i, :) = uint8(Y * T(i));
		II(i, :)= I(i, :) / T(i);
		Iplot(i,:) = II(i,:) - Y;
		L = [L; sprintf("T = %.3f", T(i))];
	endfor
	plot(X, Iplot)
	legend(L, "location", "southeast")
	print -dpng ALL.png
	axis([45 75 -20 20])
	legend(L, "location", "northeast")
	print -dpng ALL_crop.png
	clear Iplot
	plot(X, abs(median(II) ./ Y - 1)*100)
	print -dpng ALL_median.png
	OUT=[];
	for x = X
		PT = []; nn = [];
		for i = ii
			g = I(i, x);
			if(g >= I0 && g <= I1)
				PT = [PT g];
				nn = [nn i];
			endif
		endfor
		if(size(PT) == 0)
			if(g < I0)
				colr = I(end,x) / T(end);
			else
				colr = I(1,x) / T(1);
			endif
		else
			PT = PT./T(nn);
			[tmp s] = sort(PT./T(nn));
			S = ceil(size(s,2)/2);
			i = PT(S); S = nn(S);
			colr = I(S,x) / T(S);
		endif
		fflush(1);
		OUT = [OUT colr];
	endfor
	plot(X, abs(OUT ./ Y - 1)*100)
	print -dpng Filtered.png
	plot(X,Y,X,OUT)
	legend(["Original"; "HDR"]);
	print -dpng Signal_filtered.png
	LO = log(OUT+1); SO = sqrt(OUT);
	Lmax = max(LO); Lmin = min(LO);
	Smax = max(SO); Smin = min(SO);
	plot(X,OUT/max(OUT)*255, X,(LO-Lmin)/(Lmax-Lmin)*255, X,(SO-Smin)/(Smax-Smin)*255)
	legend(["out";"log(out)"; "sqrt(out)"])
	print -dpng log_sqrt.png
endfunction