Как мы исследовали светоприемник Apogee Alta U16M

[eddy_em]

Решил я (в основном - для себя, чтобы не забыть) описать результаты тестирования светоприемника Apogee Alta U16M с матрицей Kodak KAF-16803 (судя по стоимости светоприемника — немногим более 10 килобаксов, матрицу эту подобрали на помойке). Начальство поручило нам с Тимуром проанализировать, насколько поганый этот светоприемник и можно ли его использовать для ненаучных задач (например, для оценки экстинкции). В конце марта-начале апреля мы провели нужные измерения, а потом долго их обрабатывали.

Итак, мы наделали в общей сложности порядка 30ГБ снимков: «плоское поле» для определения общих параметров матрицы, а также изображения с узкополосного источника (монохроматора) для определения спектралки.
Довести температуру камеры до -35°С, как было написано в отчетах по тесту светоприемника самой фирмой, мы не смогли: камера постоянно отключалась из-за перегрева холодильника, поэтому работали на -30°С. Еще для исследования поведения темнового тока, мы делали темновые снимки на разных температурах (от -5°C до -30°C).
Стоит сказать, что я пока так и не дал ума библиотечке для работы с этой камерой, поэтому мы пользовались демо-версией MaximDL, которая работает только под мастдаем и имеет ужаснейший интерфейс (я вообще не понимаю, как за такое г. можно просить деньги). К счастью, оказалось, что к этому г. можно писать сценарии на (!) visual-basic! Это жесть!!! Но все-таки, при помощи гугла и чьей-то матери мы справились и написали сценарии для получения комплектов снимков — очень уж не хотелось вручную получать эти сотни изображений. Темновые и плоские мы снимали в течение ночи: вечером запускали сценарий, и он щелкал нам потихонечку кадры.
Основные результаты измерений обрабатывал Тимур (т.к. мне просто было некогда — на меня навалилось еще несколько срочных дел), я же попробовал сделать попиксельный анализ светоприемника (ниже расскажу). Итак, для начала результаты, полученные Тимуром.

Общие данные

gain = 1.43±0.04 ē/ADU
Bias (median) = 1239 ADU
Bias (RMS) = 10.3 ē
Dark = 0.015 ē/pix/s
Линейность ±0.05% в диапазоне 2000÷40000 ADU, ±1% в диапазоне 0÷45000ADU (с вычтенным Bias'ом)
Очень плохой столбец X = 598
QE - от 24% (380нм) до 60% (550нм)

Методы

Шум считывания определяется как шум Bias'а (лучше использовать для вычисления шума попарные разности кадров, чтобы исключить влияние "плохих" и "горячих" пикселей), деленный на корень из двух: σ_RON=σ_Bias/√2.
Линейность определялась по изображениям «плоских полей» по одному и тому же куску изображения с более-менее равномерной засветкой. Из каждого фрейма вычитался Bias, затем фрейм усреднялся. По зависимости средней интенсивности в фрейме для каждой экспозиции от времени экспозиции и вычислялась линейность (с аппроксимацией прямой I=at). Надо сказать, что кривая имела "загибы" вверх в левой части и вниз в правой.
Коэффициент усиления (gain) вычисляется обычно из зависимости дисперсии интенсивности от величины этой интенсивности. Дисперсия, опять-таки, вычислялась из попарных разностей изображений с одинаковым временем экспозиции. Общие соображения здесь просты. Общий шум складывается из шума считывания, фотонного шума и прочего аддитивного шума (тепловой шум, дробовой и т.п.). Дисперсия фотонного шума, т.к. этот шум имеет статистику Пуассона, равна среднему значению интенсивности сигнала. Шум считывания — некая постоянная величина, а прочий аддитивный шум зависит от различных факторов (температуры, фазы Луны, количества принятого спиртного и т.п.). Так как для перевода в отчеты ADU нам нужно умножить все величины в электронах на g (gain), получим, что общий шум (в ADU²) является суммой I/g (I - средняя интенсивность) и дисперсий шума считывания и добавки. Таким образом, если бы не эта неизвестная добавка, мы получили однозначную линейную зависимость для определения g. Если бы эта добавка не зависела от уровня сигнала, g тоже легко было бы получить (как котангенс наклона прямой σ(I)). Вычисляя общий шум по разности двух изображений, можно попытаться снизить влияние аддитивного нелинейного члена. После построения графика зависимости общего шума от интенсивности, выберем аниболее линейный его участок и по нему определим приблизительное значение g.
Темновой ток измерялся довольно просто: как средняя интенсивность разницы «темновой кадр» - bias. При изменении температуры чипа на 30°C темновой ток изменялся почти на порядок.
Квантовую эффективность чипа мы вычисляли при помощи монохроматора и калиброванного светоприемника (фотоэлемент с известной кривой чувствительности). Начиная от 380нм с шагом в 10нм мы делали по 5 (чтобы было хоть что-то, что можно усреднять) изображений выходной щели монохроматора + 5 «темновых» (для вычитания рассеянного света + bias'а) . Во время экспозиции темновых мы следили за значением тока, текущего через фотоэлемент и записывали среднее показание и ошибку. При приближении к 700нм мы вышли на «полку» и уже решили было, что наш фотоэлемент брешет, поэтому сделали пробные измерения на 725нм - оказалось, что это просто такая кривая QE у матрицы. Дальше мы сделали еще пару измерений на 750нм и 775нм (находиться в душной темной комнате, в которой еще и работал криотайгер для охлаждения другой матрицы, было невыносимо - мы и так каждые полчаса курить бегали, чтобы хоть немного продышаться). По результатам измерений во фрейме с изображением щели находилась суммарная интенсивность (попиксельная сумма значений разности медианы пяти изображений и медианы пяти темновых), по градуировочной таблице и значению фототока определялась освещенность, исходя из длины волны вычислялось количество попавших на светоприемник фотонов и с учетом известного уже g вычислялась QE.

Попиксельная неоднородность

Для измерения попиксельной неоднородности я долго придумывал алгоритм (памяти-то у компьютера не бесконечное количество - все файлы в нее не поместишь; библиотека cfitsio не умеет открывать одновременно много файлов, а писать свою мне не хотелось). В итоге получился такой велосипед: для каждой строки матрицы создавалась директория (4096 директорий), в которую помещалось по одному файлу на каждый пиксель (4116 файлов). В каждый файл заносились значения интенсивности в пикселе для времен экспозиции от 30 до 900 секунд с шагом в 30 секунд. Это выполнялось функцией

function proc_by_pixels(savepath)
% записывает в отдельный файл для каждого пикселя зависимость интенсивности
% от времени экспозиции
% savepath - директория, в которой будет нагромождено директорий/файлов отчета

	W = 4116; H = 4096; % размер изображений
	path = "../FITS/FLAT/Flat-"; % путь к плоским полям с началом имени файла
	printf("Creating file structure\n"); fflush(1);
	for h = 1:H % высота
		dirname = sprintf("%s/%04d", savepath, h);
		[ s e m ] = stat(dirname);
		if(e != 0) % если такой директории нет - создаем
			mkdir(dirname);
		endif
		for w = 1:W % заполняем директории файлами
			F = fopen(sprintf("%s/%04d", dirname, w), "w");
			fclose(F);
		endfor
		if(floor(h/100)*100 == h)
			printf(" %d\n", h); fflush(1);
		elseif(floor(h/10)*10 == h)
			printf("."); fflush(1);
		endif
	endfor
	for exptime = 30:30:900 % пробегаемся по всем экспозициям
		printf("\n\nExp time = %d seconds...\n", exptime); fflush(1);
		% получаем медиану плоского поля для данной экспозиции
		Img = get_median(sprintf("%s%d-", path, exptime));
		printf("Fill pixel structure...\n"); fflush(1);
		for h = 1:H
			hpath = sprintf("%s/%04d", savepath, h);
			for w = 1:W
				fname = sprintf("%s/%04d", hpath, w);
				F = fopen(fname, "a");
				if(F < 0)
					printf("Some error occured, x=%d, y=%d\n", w, h); fflush(1);
					return;
				endif
				fprintf(F, "%d\n", Img(h,w));
				fclose(F);
			endfor
		if(floor(h/100)*100 == h)
			printf(" %d\n", h); fflush(1);
		elseif(floor(h/10)*10 == h)
			printf("."); fflush(1);
		endif
		endfor
	endfor
endfunction

Вспомогательная функция для медианной фильтрации изображений:

function M = get_median(f_prefix)
% медианная фильтрация пяти файлов с префиксом f_prefix
% имя файла образуется как f_prefixN.fit, где N=1:5

	Imgs = []; % массив для накопления изображений
	for i = 1:5
		fname = sprintf("%s%d.fit", f_prefix, i);
		[I H] = fits_read(fname);
		Imgs(:,:,i) = I;
		clear I H;
	endfor
	printf("calculate median...\n"); fflush(1);
	M = median(Imgs,3);
	clear Imgs;
endfunction

Никаких bias'ов здесь делать не нужно было, т.к. по идее в попиксельной зависимости они "уже есть". А вот о чем я не подумал - так это о снятии темновых вместе с "плоскими полями", т.к. в этом случае можно было бы получить попиксельную картину неоднородности g, а так - получаем сумму с темновым током.
Где-то около шести суток работала эта программка. Оно и понятно: слишком частые дисковые операции. Возможно, построчная обработка, несмотря на то, что все изображения пришлось бы перечитывать 4 тысячи раз, была бы быстрей.
Для вычисления кое-какой статистики я запустил еще одну функцию:

function save_results(path, Xstart)
%
% пробегаемся по директории path, считывая данные по столбцам
% и сохраняя их в соответствующие файлы в текущей директории
% считать файлы можно так:
% fid = fopen(filename, "r"); M = fread(fid, [4096 4116], "float");
% fclose(fid)
%

	X0 = 1;
	if(nargin == 2)
		X0 = Xstart;
	endif
	W = 4116; H = 4096; % размер изображений
	fnames = ["frac.dat"; "frac1.dat"; "Imin.dat"; "Imax.dat"; ...
			"med.dat"; "koeffa.dat"; "koeffb.dat"; ...
			"Stda.dat"; "Stdb.dat"; "Stdc.dat"];
	N = size(fnames,1);
	fids = zeros(N);
	for i = 1:N
		fids(i) = fopen(strtrim(fnames(i,:)), "a");
		if(fids(i) == -1)
			printf("error opening file %s\n", fnames(i,:));
			return;
		endif
	endfor
	for x = X0:W
		for y = 1:H
			fields = single(parse_pixel(path, x, y));
			for i = 1:N
				fwrite(fids(i), fields(i), "float");
			endfor
		endfor
		for i = 1:N
			fflush(fids(i)); % синхронизируем буферы по окончании записи столбца
		endfor
		if(floor(x/100)*100 == x)
			printf(" %d\n", x); fflush(1);
		elseif(floor(x/10)*10 == x)
			printf("."); fflush(1);
		endif
	endfor
	printf("\nReady!!!\n\n");
	for i = 1:N
		fclose(fids(i));
	endfor
endfunction

Она считала около недели. Здесь используется вспомогательная функция

function DATA = parse_pixel(path, x, y)
%
% обрабатываем результаты для пикселя (x,y)
% по пути path (читаем файл path/y/x)
% на выходе:
% frac, frac1 - доля не отброшенных точек с критерием в 2 сигмы и сигму
% Imin, Imax - минимальное и максимальное значение интенсивности
%			на линейном участке с отбраковкой в 1 сигму
% med - медианное значение всех интенсивностей
% koeff - коэффициенты [a b] для линейной апроксимации I = at + b,
%		если после отбраковки остается лишь 3 точки, a = 0, b = -1
% std - стандартные отклонения данных от линейной аппроксимации
%		1 - для всех данных, 2 - с отбраковкой 2сигмы, 3 - с отбраковкой 1 сигма
%		если данны слишком мало, std = [-1, -1, -1];
%
	fname = sprintf("%s/%04d/%04d", path, y, x);
	T = [30:30:900]';
	I = textread(fname);
	L = size(I,1);
	D = diff(diff(I));
	koeff = [0.; -1.]; Stda = [-1 -1 -1];
	% отбраковываем выбросы больше 2std второй производной I
	idx = throwBadIdx(L, [ 1; 1 + find(abs(D) < 2*std(D))]);
	% доля неотбракованных пикселей
	frac = (2 + size(idx, 1)) / size(I, 1);
	% и отбраковка по 1 сигме
	idx1 = throwBadIdx(L, [ 1; 1 + find(abs(D) < std(D))]);
	frac1 = (2 + size(idx1, 1)) / size(I, 1);
	Imin = min(I(idx1)); Imax = max(I(idx1));
	med = median(I);
	if(size(idx,1) > 3)
		if(frac1 > 0.5)
			T0 = T(idx1); I0 = I(idx1);
		else
		% если отбраковано больше половины пикселей (overscan/bad?),
		% используем весь ряд данных
			if(frac > 0.5)
				T0 = T(idx); I0 = I(idx);
			else
				T0 = T; I0 = I;
			endif
		endif
		TT = [ T0 ones(size(T0)) ];
		koeff = TT \ I0;
		I1 = T*koeff(1)+koeff(2);
		Std = [ std(I - I1)				...
				std(I(idx) - I1(idx))	...
				std(I(idx1) - I1(idx1)) ];
	endif;
	DATA = [frac frac1 Imin Imax med koeff' Std];
	%plot(T,[I I1], T0, I0, 'o');
endfunction

function newidx = throwBadIdx(L, idx)
% отбраковка выбросов
	badIdx = find(diff(idx) > 1);
	if(~isempty(badIdx))
	% не обращаем внимания на выбросы внутри середины данных
	%	if(size(badIdx, 1) > 1) % выбросов больше одного
			% нижний край - выброс с максимальным номером в
			% пределах первой четверти данных
			low = max(badIdx(find(idx(badIdx) < L*0.25)));
			% верхний край - выброс с минимальным номером
			% в пределах последней четверти данных
			hig = min(badIdx(find(idx(badIdx) > L*0.75)));
			if(~isempty(hig))
				idx = idx(1:hig);
			endif;
			if(~isempty(low)) idx = idx(low+1:end); endif;
	%	endif
	endif;
	newidx = idx;
endfunction

Принцип обработки данных для каждого пикселя таков:

1. вычисляем вторую производную (через разделяемые разности) зависимости I от t и ее стандартное отклонение от нуля;
2. выбросы второй производной, превышающие 2σ, отбраковываем (т.к. они явно уже не лежат на прямой линии):
   1. на выбросы в центре ряда данных не обращаем особого внимания;
   2. обрабатываем на выбросы лишь первую и последнюю четверти ряда данных;
   3. в пределах обрабатываемых частей находим номера ближайших к центру ряда данных выбросов и считаем валидными лишь элементы, располагающиеся между этими выбросами;
   вычисляем долю оставшихся неотбракованными точек;
3. повторяем предыдущий пункт с отбраковкой по 1σ и вычисляем долю неотбракованных точек, далее работаем с точками, оставшимися после отбраковке по σ;
4. вычисляем нижнюю и верхнюю границы интервала интенсивностей;
5. аппроксимируем оставшиеся точки прямой (методом наименьших квадратов) и находим наклон и точку пересечения с ординатой, здесь мы используем оставшиеся точки, если их не менее половины от первоначального количества, иначе используем все точки (т.е. это - явно либо bad, либо hot);
6. вычисляем стандартные отклонения от найденной прямой: всех точек зависимости, отбракованных по 2σ, отбракованных по 1σ.

Вот кое-какая статистика по полученным данным:

frac.dat - доля неотбракованных для линейной аппроксимации пикселей в пределах
	двух сигм

min = 0.5
max = 1.0
mean = 0.96154
std(data-mean) = 0.027642

Количество пикселей (N) с data < I
I		N
0.9		292247
0.8		10156
0.7		1744
0.6		140
0.5		0


================================================================================

frac1.dat - доля неотбракованных для линейной аппроксимации пикселей в пределах
	одной сигмы

min = 0.2
max = 1.0
mean = 0.87400
std(data-mean) = 0.092842

Количество пикселей (N) с data < I
I		N
0.9		10729451
0.8		1835013
0.7		1008734
0.6		353171
0.5		118556
0.4		14560
0.3		194
0.2		0


================================================================================

Imax.dat - максимальная интенсивность на участке с линейной интерполяцией с
	точностью до сигмы

min = 1213
max = 65535
mean = 60521
std(data-mean) = 5012.6

Количество пикселей (N) с data < I
I		N
6553.6		73726
13107.2		73728
19660.8		73728
26214.4		73925
32768		75653
39321.6		79721
45875.2		80460
52428.8		347725
58982.4		4.11578e+06


================================================================================

Imax.dat - минимальная интенсивность на участке с линейной интерполяцией с
	точностью до сигмы

min = 1198
max = 65535
mean = 4371.42
std(data-mean) = 3631.89

Количество пикселей (N) с data < I
I		N
3276.8		77116
6553.6		16035620
13107.2		16079312
19660.8		16365708
26214.4		16845985
32768		16858891
39321.6		16859120
45875.2		16859130
52428.8		16859133
58982.4		16859133


================================================================================

Stda.dat    - стандартное отклонение всех точек по интенсивности от линейной
	аппроксимации

min = 1.60561
max = 13684.4
mean = 1078.75
std(data-mean) = 392.584


Количество пикселей (N) с data < I
I		N
750		3719601
1500		14231936
3000		16851736
4500		16856908
6000		16859075
7500		16859105
9000		16859114
10500		16859128
12000		16859132
13500		16859135


================================================================================

Stdb.dat.fit    - стандартное отклонение точек в пределах двух сигм от линейной аппроксимации

min = 0
max = 14108.1
mean = 249.268
std(data-mean) = 82.5815

Количество пикселей (N) с data < I
I		N
50		73729
100		74200
200		407597
300		16083234
400		16817342
500		16848989
600		16851597
700		16852211
800		16852303
900		16852338


================================================================================

Stdc.dat.fit    - стандартное отклонение точек в пределах одной сигмы от линейной аппроксимации

min = 0
max = 8949.48
mean = 218.044
std(data-mean) = 79.4199

Количество пикселей (N) с data < I
I		N
25		73722
50		74782
100		171754
150		1054747
200		4549563
250		14185746
300		16827374
350		16834100
400		16834699
450		16834980

================================================================================

koeffa.dat.fit  - Коэффициент a линейной аппроксимации (I = a*T + b)
	(пропорционален освещенности в пикселе)

min = -0.00970588
max = 90.5223
mean = 76.2195
std(data-mean) = 5.43172

Количество пикселей (N) с data < I
I		N
70		106235
72		486162
74		1973313
76		5807382
78		12253130
80		16843472

================================================================================

koeffb.dat.fit  - Коэффициент b линейной аппроксимации (I = a*T + b)
	(нечто вроде вычисленного bias'а)

min = -7689.82
max = 65535
mean = 1626.94
std(data-mean) = 197.202

Количество пикселей (N) с data < I
I		N
0		3348
1000		3361
1500		1285936
1700		15928839
2000		16362533


================================================================================

med.dat    - медиана интенсивности (для всех экспозиций)

min = 1210
max = 65535
mean = 37297.8
std(data-mean) = 2575.64

Количество пикселей (N) с data < I
I		N

25000		77079
30000		77099
34000		82209
36000		1362629
38000		11301119
42000		16858959
45000		16859033

Картинки

Stdc - стандартное отклонение от линейной аппроксимации в пределах одной сигмы. Видно, что помимо основного «очень плохого» столбца (слева) есть еще один плохой столбец покороче и «плоховатый» столбец (штрих справа вверху):

koeffa - величина, пропорциональная сумме произведения QE на освещенность пикселя и темнового тока в пикселе. Явственно прослеживаются два плохих столбца и область оверскана справа:

koeffb - эдакий «вычисленный» bias. Из-за всяких нелинейностей и влияния темнового тока он получился выше реального:

Ну и напоследок - как это выглядит в 3D.
bias:

Плоское поле:

Медиана пяти плоских полей:

Шумы:

P.S. К сожалению, у ЖЖшки отвалилась функция добавления картинок в сам ЖЖ, так что пришлось все кидать на imageshack.
P.P.S. Сегодня (10 мая) я таки умудрился домучить libapogee2.2 (которую нашел где-то на просторах интернета), чтобы она считывала мне кадр с ПЗС. Однако, несмотря на то, что экспозиция протекает нормально, температура устанавливается какая надо, все пиксели полученных изображений содержат 65535.
Буду доламывать дальше: все-таки, эта камера нам нужна будет для работы и надо обязательно научиться полностью управлять ею.

Comments

[sevasat.livejournal.com]

Блин, а мускуль и постгрес для слабаков? Иначе зачем было городить этот огород со 100500 каталогами и файлами?

[eddy-em.livejournal.com]

Это было бы значительно дольше! Да я и не представляю себе, как можно организовать БД, в которой бы хранилось 100500 рядов данных!
Да и просто не умею работать с БД - не нужно мне это совершенно!

[sevasat.livejournal.com]

Так мускуль простой как топор. Ты ему срешь и говоришь куда положить, он ложит. Общее устройство можно вечером за чаем раскурить, еще пол-дня на то как с ним работать. А потом всяко пригодится.

[eddy-em.livejournal.com]

Я для простых нужд sqlite использовал когда-то. Но это было давно, и все функции работы с БД запускаются из оберток.
Логи же всякие хранить я в бинарном виде предпочитаю — значительная экономия места, а скорость доступа тоже вполне нормальная, т.к. каждые N записей я делаю в отдельном файле отметку о смещении.