\documentclass[a4paper,12pt]{extarticle}
\usepackage{/home/eddy/ed, verbatim}
\title{Практикум \No6: обработка изображений}
\author{}\date{}\nocolon

\long\def\task#1{\par\noindent\leavevmode\refstepcounter{sect}\llap{\textbf{\thesect}\;}\indent\textit{#1}\par}
\def\t#1{{\upshape\ttfamily #1}}

\usepackage{etoolbox}
\makeatletter
\preto{\@verbatim}{\topsep=2pt \partopsep=1pt} % remove big space before verbatim
\AtEndEnvironment{verbatim}{\vspace{5pt}}
\makeatother


\begin{document}
\maketitle

\section{Работа с изображениями в Octave}
\subsection{FITS-файлы}
Наиболее популярной утилитой просмотра FITS-файлов является \t{ds9}.

Для работы с обычными изображениями в Octave используются функции \t{imread}, \t{iminfo},
\t{imshow}. Основные функции обработки изображений находятся в пакете \t{image}.
FITS-файлы~--- пакет \t{fits}. Узнать входящие в состав пакета функции можно так:
\begin{verbatim}
pkg describe -verbose fits
---
Package name:
fits
Version:
1.0.7
Short description:
The Octave-FITS package provides functions for  reading, and writing FITS
(Flexible Image Transport System) files.  This package uses the libcfitsio library.
Status:
Loaded
---
Provides:
Reading and writing FITS files
read_fits_image
save_fits_image
save_fits_image_multi_ext
\end{verbatim}

Для чтения используем \t{read\_fits\_image}:
\begin{verbatim}
[I H] = read_fits_image('stars.fit');
imshow(I,[1000 4000])
\end{verbatim}
Необходимо указывать пороги интенсивностей для \t{imshow}. Более точно можно определить диапазоны
из статистики по изображению:
\begin{verbatim}
Imin = min(I(:))
 Imin =  1321
Imax(I(:))
 Imax =  52826
Imean = mean2(I)
 Imean =  2054.6
Istd = std2(I)
 Istd =  1275.5
Imed = median(I(:))
 Imed =  2009
Low = Imin; High = Imed + Istd;
imshow(I, [Low High]);
% or: imshow(I/High)
\end{verbatim}
Как видите, изображение повернуто вправо на $90\degr$.
Если мы захотим отображать изображения как в оригинале, нужно заранее повернуть их:
\begin{verbatim}
I = rot90(I);
imshow(I, [Low High])
\end{verbatim}
Следует помнить, что в Octave координаты изображения начинаются с верхнего левого угла (и ось~Y
направлена вниз), а в стандарте FITS~--- с левого нижнего (и ось~Y направлена вверх).

Убедимся в этом:
\begin{verbatim}
I2=I;
I2([1:100],[1:100])=65000;
imshow(I2, [Low High]);
\end{verbatim}
Сравним, как отображаются координаты пикселя в ds9.

\subsection{Работа с гистограммой}
Для улучшения отображения попробуем эквализовать гистограмму:
\begin{verbatim}
HE = histeq(I/High, 65536);
imshow(HE);
% compare to HE = histeq(I/High, 256);
\end{verbatim}

Посмотрим, как меняется гистограмма изображения в зависимости от распределения шума.
Сгенерируем сначала пробное изображение:
\begin{verbatim}
I = uint8(zeros(300,300));
I([1:50 250:end], :) = 200;
I([51:249], [1:30 270:end]) = 200;
I([101:199], [100:200]) = 100;
imshow(I)
imhist(I)
Iuniform = I + uint8(50*rand(size(I)));
imshow(Iuniform)
imhist(Iuniform)
Inorm = I + uint8(25*randn(size(I)));
imshow(Inorm)
Inorm = uint8(I + 25*randn(size(I)));
imshow(Inorm)
imhist(Inorm)
Ipoisson = imnoise(I, "poisson");
imshow(Ipoisson)
imhist(Ipoisson)
Isaltpepper = imnoise(I, "salt & pepper", 0.2);
imshow(Isaltpepper)
imhist(Isaltpepper)
Ispeckle = imnoise(I, "speckle", 0.03);
imshow(Ispeckle)
imhist(Ispeckle)
\end{verbatim}

Сравним с эквализацией гистограммы: \t{imshow(histeq(Inorm))}.

Попробуем применить к изображению гамма-коррекцию:
\begin{verbatim}
imshow(imadjust(Iuniform, [], [], 0.5))
imshow(imadjust(Iuniform, [], [], 5))
\end{verbatim}

Теперь попробуем преобразовать яркость кусочно-линейной функцией
$$\begin{cases}
I_{out} = 0.2 I_{in}, & I_{in} <  70,\\
I_{out} = 14 + 1.5 (I_{in}-70), & 70 \le I_{in} < 130, \\
I_{out} = 209 + 0.368*(I_{in}-130),& 130\le I_{in} < 256.
\end{cases}
$$
\begin{verbatim}
Ilinhyst = zeros(size(Inorm));
idx = find(Inorm < 70);
Ilinhyst(idx) = 0.2*Inorm(idx);
idx = find(Inorm >= 70 & Inorm < 130);
Ilinhyst(idx) = 14 + 1.5*(Inorm(idx)-70);
idx = find(Inorm >= 130);
Ilinhyst(idx) = 209 + 0.368*(Inorm(idx)-130);
Ilinhyst = uint8(Ilinhyst);
imshow(Ilinhyst);
imhist(Ilinhyst)
\end{verbatim}

Попробуем сделать ступенчатое преобразование яркости: на изображении \t{Inorm} заменим все
интенсивности в диапазоне $[50, 150]$ на 100:
\begin{verbatim}
Ilinhyst = Inorm;
Ilinhyst(find(Inorm >= 50 & Inorm <= 150)) = 100;
imshow(Ilinhyst)
imhist(Ilinhyst)
\end{verbatim}

\subsection{Фильтрация изображений}
Попробуем простую медианную фильтрацию на примере пары изображений.
\begin{verbatim}
subplot(2,2,1)
subimage(Isaltpepper)
subplot(2,2,2)
subimage(medfilt2(Isaltpepper, [5 5])))
subplot(2,2,3)
subimage(Inorm)
subplot(2,2,4)
subimage(medfilt2(Inorm, [5 5]))
\end{verbatim}

Видим, что медианный фильтр вполне неплохо справился с шумом <<соль-перец>> и значительно хуже~---
с гауссовым шумом.

Рассмотрим применение различных фильтров. Для начала создадим матрицы фильтров:
\begin{verbatim}
fgauss = fspecial("gaussian", 5, 1);
flog = fspecial("log", 5, 1);
fsobelv = fspecial("sobel");
fsobelh = rot90(fsobelv);
\end{verbatim}

И применим их к оригинальному изображению, \t{Isaltpepper} и \t{Inorm}:
\begin{verbatim}
imshow(imfilter(I, fgauss, "symmetric"));
imshow(imfilter(Isaltpepper, fgauss, "symmetric"));
imshow(imfilter(Inorm, fgauss, "symmetric"));
imshow(histeq(imfilter(I, flog, "symmetric")));
imshow(histeq(imfilter(Isaltpepper, flog, "symmetric")));
imshow(histeq(imfilter(Inorm, flog, "symmetric")));
\end{verbatim}

Как видите, у дифференциального фильтра отсекается отрицательная компонента~--- т.к. изображения
имеют тип \t{uint8\_t}, однако, можно сделать так:
\begin{verbatim}
imshow(histeq(imfilter(double(Isaltpepper), flog, "symmetric")));
imshow(histeq(imfilter(double(Inorm), flog, "symmetric")));
imshow(histeq(imfilter(double(I), flog, "symmetric")));
\end{verbatim}

Octave неправильно рассчитывает значение фильтров, например:
\begin{verbatim}
sum(flog(:))
ans = -0.024092
sum(fgauss(:))
ans =  1.00000
sum(fsobelv(:))
ans = 0
\end{verbatim}
У дифференциальных фильтров сумма элементов должна быть равна нулю (чтобы константные уровни
превращались в нуль), а у усредняющих~--- единице (чтобы сохранять среднюю интенсивность по
изображению). Расширим  размер матрицы фильтра лапласиана гауссианы, чтобы улучшить приближение к
нулю:
\begin{verbatim}
flog = fspecial("log", 17, 1);
imshow(histeq(imfilter(double(I), flog, "symmetric")));
imshow(histeq(imfilter(double(Isaltpepper), flog, "symmetric")));
imshow(histeq(imfilter(double(Inorm), flog, "symmetric")));
\end{verbatim}
Теперь на примере первого изображения мы видим, что фильтр работает правильно.

Горизонтальная и вертикальная компоненты градиента:
\begin{verbatim}
ih = imfilter(double(Inorm), fsobelh, "symmetric");
iv = imfilter(double(Inorm), fsobelv, "symmetric");
imshow(histeq(iv)));
imshow(histeq(ih));
grad = sqrt(ih.*ih + iv.*iv);
imshow(histeq(grad));
\end{verbatim}

Если мы попытаемся проделать то же самое для \t{I}, увидим, что \t{histeq} в данном случае не
сработает, и отображать придется командой
\begin{verbatim}
imshow(grad,[min(grad(:)) max(grad(:))]);
\end{verbatim}
Здесь \t{histeq} не сработал из-за слишком <<кривой>> гистограммы: подавляющее большинство пикселей
приходится на нулевой уровень.

Для смягчения общего контраста можно  использовать прием, называемый <<нерезким маскированием>>
(применяемый еще в фотографии на фотоэмульсию). Для этого воспользуемся формулой
$$I_{out}=N\cdot I_{in} - F(I_{in}),$$
где $F$~-- некоторый усредняющий фильтр (в оригинале~--- расфокусированное позитивное изображение).
\begin{verbatim}
imshow(2*Inorm-imfilter(Inorm, fgauss, "symmetric"));
\end{verbatim}

Посмотрим, как изменится размытое изображение при нерезком маскировании:
\begin{verbatim}
ig = imfilter(I, fgauss, "symmetric");
subplot(1,2,1)
imshow(ig)
subplot(1,2,2)
imshow(histeq(2*ig-imfilter(ig, fgauss, "symmetric")))
\end{verbatim}

Очистим все и вернемся к FITS-файлам. Загрузим изображение и построим дифференциальные фильтры с
более крупным ядром:
\begin{verbatim}
clear
[I H] = read_fits_image('stars.fit');
imshow(histeq(I,65536));
fgauss = fspecial("gaussian", 30, 5);
flog = fspecial("log", 37, 5);
\end{verbatim}

Попробуем теперь применить к изображению нерезкое маскирование и различные фильтры:
\begin{verbatim}
imshow(histeq(I-imfilter(I, fgauss, "symmetric"), 65536));
imshow(histeq(imfilter(I, fgauss, "symmetric"), 65536));
imshow(histeq(imfilter(I, flog, "symmetric"), 65536));
imshow(histeq(2*I-imfilter(I, fgauss, "symmetric"), 65536));
imshow(histeq(medfilt2(I), 65536));
imshow(histeq(I-medfilt2(I, [15 15]), 65536))
\end{verbatim}
Последний вариант является первым приближением к вычитанию фона.

Посмотрим, как можно найти приближенную маску для поиска звезд:
\begin{verbatim}
Imed = medfilt2(I, [3 3]);
 % [3 3] is default, we can omit it here
ImedLG = imfilter(Imed, flog, "symmetric");
imshow(ImedLG, [-0.1, 0.1]);
mask = logical(zeros(size(I)));
\end{verbatim}

Теперь нужно подобрать подходящий пороговый уровень:
\begin{verbatim}
mask = logical(zeros(size(I))); mask(find(ImedLG < 0)) = 1;
imshow(mask)
mask = logical(zeros(size(I))); mask(find(ImedLG < -0.005)) = 1;
imshow(mask)
\end{verbatim}
Уже такое грубое приближение позволило обнаружить некоторое количество звезд на изображении. Если
использовать более точный (<<робастный>>) метод выделения фона (а здесь фон достаточно перекошен),
получатся результаты куда лучше.

Отобразим оригинал, заменив нулями все точки, где маска истинна, а затем~--- наоборот:
\begin{verbatim}
J=I; J(mask)=0; % hide stars
imshow(histeq(J, 65536));
J=I; J(~mask)=0; % hide background
imshow(J);
\end{verbatim}
Наша грубая маска определила все звезды, а также некоторое количество шумов, но при помощи
морфологических операций мы можем удалить мелкие объекты (правда, не факт, что не потеряем таким
образом слабые звезды; кроме того, мы видим на первом изображении, что особо яркие звезды
маскированы не полностью, наш способ достаточно плох и годится лишь для очень  грубой астрометрии).

\subsection{Морфологические операции, поиск связных областей}
Попробуем подчистить немного нашу маску.
\begin{verbatim}
imshow(imopen(mask, strel('disk', 1, 0)));
imshow(imopen(mask, strel('disk', 2, 0)));
imshow(imopen(mask, strel('disk', 3, 0))); % use this
mask = imopen(mask, strel('disk', 3, 0));
\end{verbatim}

Теперь можем аналогично просмотреть маскированные изображения:
\begin{verbatim}
J=I; J(mask)=0; % hide stars
imshow(histeq(J, 65536));
J=I; J(~mask)=0; % hide background
imshow(J);
\end{verbatim}

Но нам интересней выделить все звезды из изображения. Используем поиск связных областей:
\begin{verbatim}
[L NUM]= bwlabel(mask, 8);
NUM
NUM =  200
imshow(L, [1 NUM])
\end{verbatim}
Как видим, нумерация объектов идет слева-направо, сверху-вниз. Чтобы выделить конкретный объект, мы
можем применить маску такого типа:
\begin{verbatim}
J = zeros(size(I)); idx = find(L == 100); J(idx) = I(idx);
imshow(J);
\end{verbatim}
Мы выделили одну звезду. Можем выделить тонкую полоску фона вокруг нее:
\begin{verbatim}
mask100 = logical(zeros(size(I)));
mask100(idx) = 1;
mask100 = xor(imdilate(mask100, strel('disk',2,0)), mask100);
imshow(mask100);
\end{verbatim}

И посчитать по ней статистику:
\begin{verbatim}
bg100 = I(mask100);
mean(bg100)
ans =  2063.1
std(bg100)
ans =  61.711
\end{verbatim}

Но для вычисления положений центроидов удобней использовать другую функцию поиска связных областей:
\t{bwconncomp}.
\begin{verbatim}
K = bwconncomp(mask, 8);
p = regionprops(K);
\end{verbatim}
По умолчанию в \t{p} (это~--- массив структур) содержатся поля \t{Area} (площадь объекта),
\t{BoundingBox} (координаты бокса, в котором содержится объект: левый верхний угол и размер) и
\t{Centroid} (геометрический центр объекта). По первым двум параметрам мы можем отфильтровать
слишком мелкие и слишком крупные объекты, а также объекты, у которых отношение ширины к высоте
сильно отличается от единицы (т.е.  там внутри~--- явно не круг). Третий параметр позволяет судить
о примерном расположении центра объекта (более точно, конечно, мы можем вычислить его лишь взяв
данные с оригинального изображения, вычтя фон и посчитав точный центроид).

Можно нарисовать несколько боксов:
\begin{verbatim}
rectangle('Position', p(1).BoundingBox, 'EdgeColor','g','LineWidth',2)
rectangle('Position', p(3).BoundingBox, 'EdgeColor','g','LineWidth',2)
rectangle('Position', p(100).BoundingBox, 'EdgeColor','g','LineWidth',2)
\end{verbatim}

\begin{verbatim}
centroids = cat(1, p.Centroid);
plot(centroids(:,1), centroids(:,2), 'ko');
\end{verbatim}

Теперь сравним координаты центроида изображения звезды с координатами центроида маски. Обратите
внимание, что \t{BoundingBox} имеет координаты <<с половинками>>:
\begin{verbatim}
p(100).BoundingBox
502.500   710.500    14.000    14.000
imshow(mask100);
rectangle('Position', p(100).BoundingBox, 'EdgeColor','g','LineWidth',2)
\end{verbatim}
Это связано с тем, что при отображении изображений целым координатам соответствует центр пикселя.
Контур, описывающий нашу область интереса, не должен пересекать пиксели, из-за этого координаты его
левого угла на 0.5\,пикселя меньше, а ширина~--- на 1\,пиксель больше.
Следовательно, для того, чтобы получить часть изображения, содержащую только пиксели зоны интереса,
нам нужно будет сделать преобразование:
\begin{verbatim}
Bstart = p(100).BoundingBox(1:2) + 0.5;
Bsize = p(100).BoundingBox(3:4) - 1;
rectangle('Position', [Bstart Bsize], 'EdgeColor','r','LineWidth',2);
\end{verbatim}
Теперь наш бокс включает все пиксели изображения звезды. Вычислим суммарный поток и центроид:
\begin{verbatim}
bg = mean(bg100); % mean background
[y x] = ind2sub(size(I), idx);
% coordinates of pixels in mask of 100th star
M = sum(I(idx)) - numel(idx)*bg
% energy in circle
M =  55148.54545
xc=0; yc=0; for n=1:numel(x); i=I(idx(n))-bg; xc+=i*x(n); yc+=i*y(n);
printf("%d(%d,%d)\n",i,x(n),y(n)); endfor; xc/=M; yc/=M;
xc
 % X centroid position
xc =  509.36
yc  % Y centroid position
yc =  718.41
p(100).Centroid
 % centroid by mask
ans =
509.63   717.57
imshow(histeq(I, 65536)); % take a look at real image
% real centroid:
rectangle('Position', [xc-1 yc-1 2 2], 'EdgeColor','g','LineWidth',2,
    "Curvature", 1)
% centroid of mask:
rectangle('Position', [p(100).Centroid-1 2 2], 'EdgeColor','r','LineWidth',2,
    "Curvature", 1)
\end{verbatim}

Мы видим, что настоящий центр тяжести изображения звезды достаточно сильно отличается от
геометрического центра маски. Особенно сильно это проявляется в случае деформаций изображения
(вследствие сдвига изображения, комы или прочих аберраций).

Вычислять большое количество характеристик удобней иначе. Например, так:
\begin{verbatim}
regioninfo = regionprops(L,I,'MajorAxisLength','MinorAxisLength','PixelIdxList',
    'PixelValues');
\end{verbatim}
мы получим значения большой и малой полуосей эллипса с идентичными нашей фигурой моментами (что
позволит судить о <<круглости>> каждого объекта), список индексов пикселей и список интенсивностей
пикселей на изображении. Если бы фон у нашего объекта был идеально гладким, можно было бы добавить
еще свойство \t{'WeightedCentroid'}, чтобы сразу получить координаты центра каждой звезды, но, увы,
для каждой звезды придется считать фон аналогично вышесказанному~--- через кольца. Для тесных полей
потребуется убедиться, что зона вокруг звезды свободна от другой звезды (при помощи все тех же
морфологических операций).

Аргументом \t{regionprops} может выступать и структура, полученная при помощи \t{bwconncomp}.

Теперь найдем все объекты более-менее круглой формы:
\begin{verbatim}
roundness = [regioninfo.MajorAxisLength] ./ [regioninfo.MinorAxisLength];
staridx = find(roundness < 1.1 & roundness > 0.9);
 % indexes of "stars"
imshow(ismember(L, staridx)); % show filtered image
\end{verbatim}

Теперь, если пройтись итеративно по каждому объекту с номером \t{staridx}, можно посчитать для него
индивидуальную маску фона (как было показано выше~--- через дилатацию), определить таким образом
уровень фона и использовать его для вычисления прочих параметров: полной интенсивности внутри маски
и координат центроида.

Если мы добавим в \t{regionprops} свойство \t{'Area'}, у нас появится возможность до получения
номеров объектов, возможно являющихся звездами, отсортировать их по убыванию площади (скорей
всего~--- и по убыванию яркости):
\begin{verbatim}
[~, i] = sort([regioninfo(staridx).Area], 'descend'); % get sorted indexes
staridxsorted = staridx(i);
\end{verbatim}
И можем отобразить окрестности первых 10 по площади маски:
\begin{verbatim}
imshow(histeq(I, 65536));
for i=1:10;
c = hsv2rgb([i/10. 1 1]);
rectangle('Position', p(staridxsorted(i)).BoundingBox, 'EdgeColor', c,
    'LineWidth',2);
endfor
for i = 1:10; ltxt{i} = sprintf("star #%d", i); endfor
legend(ltxt)
\end{verbatim}
Видим, что у ярких звезд маска потеряла приличную часть крыльев, а также одна достаточно яркая
звезда была пропущена.
Более надежным будет отсортировать по значению интенсивности.

\section{Первичная обработка снимков}
Наиболее распространенными светоприемниками в астрофизике оптического диапазона являются
ПЗС-матрицы. Среди прочих светоприемников этого диапазона они имеют наилучшие характеристики.
Однако, есть у них и недостатки: в тонких подложках возможна интерференция (особенно ближе к
ИК-диапазону), приводящая к возникновения <<фрингов>>; а наличие хотя бы одного очень яркого
источника приводит к <<растеканию заряда>>, усугубляющемуся при считывании изображения.
КМОП-светоприемники в значительной мере избавлены от проблем растекания заряда, но в них есть
другая проблема: каждый пиксель по сути снабжен собственным усилителем, в результате точная
фотометрия с КМОП превращается в сплошной ад!

Процесс получения изображений на ПЗС можно в упрощенном виде свести к формуле:
$$ ADU_{x,y}= \gamma F_{x,y}(e_{x,y} + D_{x,y})T +  B_{x,y} + R,$$
(здесь опущен процесс дискретизации, растекания заряда и т.п.), где $ADU$~-- отсчеты, считанные с
пикселя с координатами $(x,y)$; $\gamma$~-- <<gain>>, коэффициент преобразования фотоэлектронов в
$ADU$ (зависит от длины волны излучения); $F$~-- поглощение излучения вследствие неидеальности
оптической системы (компенсируется делением на <<плоские поля>>); $e$~-- количество фотоэлектронов,
возникших при попадании кванта в данную ячейку (зависит от длины волны), оно зашумляется множеством
факторов: фотонным шумом, турбулентностью атмосферы и т.д., и т.п.; $D$~--
<<темновой ток>>, определяющийся электронами, возникающими в результате тепловых процессов в теле
полупроводника; $B$~-- <<bias>>, начальное смещение (даже при нулевой экспозиции~--- собственно,
так кадры <<bias>> и получают) для устранения нелинейности; $T$~-- время накопления сигнала
(экспозиция); $R$~-- шум считывания.

Избавиться от шума считывания невозможно, однако, как мы помним, если усреднить $N$~кадров, то
$R$~уменьшится в $\sqrt{N}$~раз. Кроме того, чем больше экспозиция, тем большим будет вклад прочих
шумов. Однако, при малых экспозициях может проявляться <<эффект затвора>>, вызванный тем, что время
открывания\slash закрывания затвора конечно.

Чтобы скорректировать интенсивность на уровень смещения, мы можем сделать как можно больше кадров с
нулевой экспозицией (<<bias>>), получить медиану стопки изображений и вычесть из изображения. При
этом шум bias'а также будет уменьшен в $\sqrt{N}$~раз.

Для коррекции на тепловой шум нам необходимо накопить как можно больше <<дарков>> (в идеале их
экспозиция должна быть такой же, как экспозиция кадра данных, а условия должны быть абсолютно
такими же: температура чипа, окружающей среды и т.д., и т.п.). При закрытом затворе мы получим:
$$ ADU_{x,y}= \gamma D_{x,y}T +  B_{x,y} + R,$$
т.е. <<автоматически>> сюда же включен уровень смещения. Следовательно, если наши <<дарки>> имеют
ту же самую экспозицию, что и основной научный кадр, то для удаления тепловых шумов достаточно
вычесть медианное усреднение стопки <<дарков>> из нашего кадра. Однако, не всегда это возможно:
скажем, при получении спектра высокого разрешения экспозиция может длиться несколько часов!
Естественно, мы не можем получить даже десятка темновых кадров с такой экспозицией (т.к. и условия
поменяться успеют, и начнется новая ночь). В подобных случаях прибегают к хитрости: темновой ток
пропорционален времени экспозиции. Это означает, что можно сделать несколько серий коротких
экспозиций (скажем, одна, пять, десять и двадцать минут) в течение дня, а затем экстраполировать их.
Вот в этом случае для получения такого <<искусственного дарка>> необходимо будет, помимо всех
остальных кадров, снимать еще и bias'ы. Для каждой экспозиции~$T$ мы вычисляем темновой ток как
$D_T=M(D_{Ti})-M(B_i)$, где $M$~-- медиана. При помощи линейной экстраполяции мы получаем
значение~$D$ для заданной экспозиции. Но т.к. это значение уже избавлено от смещения, для коррекции
нам необходимо будет из изображения вычесть не только~$D$, но и~$B$!
Это~--- единственный случай, когда нам необходимо будет получить и стопку кадров смещения (как
говорилось выше, если <<темновые>> сняты с той же экспозицией, что и основной кадр, то смещение
автоматически корректируется при вычитании <<темновых>> из <<объекта>>).

Значительно более сложной является коррекция на неоднородности оптической системы (вплоть до
пылинок на фильтрах и прочих оптических поверхностях), $F$. Термин <<плоскость>> в данном случае
носит двоякий смысл: для прямых снимков это~--- неоднородность освещенности по кадру. Для
спектральных~--- еще и неоднородность кривой спектральной чувствительности. Отвлечемся пока от
спектроскопии. Основными факторами, искажающими интенсивность объектов в кадре, являются
виньетирование (экранирование боковых пучков апертурой оптического тракта) и прочие неоднородные
свойства оптических систем. В идеальной ситуации для компенсирования~$F$ нам необходимо получить
снимок абсолютно однородно засвеченного протяженного объекта с угловыми размерами, значительно
превосходящими поле зрения нашей оптической системы. В случае телескопов с малым полем зрения для
данных целей хватает вечерних или утренних снимков неба при максимально рассеянном солнечном свете
во время гражданских сумерек. При большом поле зрения необходимо использовать специальные источники
света с как можно более равномерным освещением. После получения <<плоских>> кадров следует
вычислить медиану из их стопки, вычесть <<темновые>>, затем нормировать на максимум
интенсивности, а далее~--- разделить подготовленный (с вычетом <<темновых>>) научный кадр на
полученное изображение. При этом мы
получаем:
$$I_{x,y} = \frac{\gamma F_{x,y}e_{x,y}T + R}{\gamma F_{x,y}T + R}\approx  e_{x,y}T.$$
Конечно, <<плоские>> кадры мы можем снимать с любой подходящей экспозицией. Но для коррекции на
темновые токи, если экспозиция <<плоских>> отличается от экспозиции <<научных>> кадров, нам
необходимо будет сделать еще один набор <<темновых>>.

Как пример медианного усреднения стопки изображений, сгенерируем зашумленные изображения (можно
использовать и \t{awgn}):
\begin{verbatim}
for N=1:10; II(:,:,N) = I+100*randn(size(I)); endfor
imshow(histeq(II(:,:,2),65536))
\end{verbatim}

Подобным образом формируются на ПЗС реальные изображения (разве что шумы мы несколько
преувеличили). Именно из-за шумов для получения качественных кадров dark и flat необходимо сделать
их как можно больше! После чего стопка изображений усредняется:
\begin{verbatim}
Imed = median(II, 3);
subplot(2,2,1)
imshow(histeq(II(:,:,2),65536))
subplot(2,2,2)
imshow(histeq(II(:,:,5),65536))
subplot(2,2,3)
imshow(histeq(Imed,65536))
subplot(2,2,4)
imshow(histeq(I,65536))
\end{verbatim}

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

И для изображений объекта мы можем выполнять подобные манипуляции для снижения шума каждого
$e_{x,y}$. Но, в этом случае процесс
восстановления будет намного сложней, т.к., во-первых, телескоп имеет неидеальное сопровождение (и
объекты не только немного <<размазываются>>, но еще и смещаются от экспозиции к экспозиции);
во-вторых, в процессе получения изображений могут изменяться условия (seeing, температура и пр.), в
результате чего PSF изображений также будет меняться (следовательно, придется сделать некоторую
деконволюцию и приведение изображения к единой координатной системе, а это чревато ошибками
интерполяции).

\section{Задание для самостоятельного выполнения}
В данных к заданию содержатся кадры object, dark и flat. Выполните следующее.
\begin{enumerate}
\item Из <<научных>> кадров выберите один с наилучшим качеством изображений.
\item Вычислите медианные значения темновых и плоских, получите нормированное изображение.
\item Подберите подходящий фильтр, который позволит получить маску, выделяющую изображения звезд.
\item Подготовьте m-файлы для выполнения рутинных операций: вычисления маски кольца фоновых
пикселей, определения уровня фона, освещенности от звезды и координат ее центроида.
\item Отсортируйте маски по уменьшению площади, для дальнейших вычислений используйте первые
пятьсот.
\item Вычислите освещенности от звезд в пределах полученных масок. Отсортируйте по уменьшению
освещенности. Выберите первые 20.
\item Выполните астрометрию для этих самых ярких звезд кадра, получив координаты~$(x,y)$ центроидов.
\item Проверьте полученные величины (преобразуйте ADU в звездные величины и проведите нормировку на
стандартные объекты из каталогов; для координат центроидов учтите, что изображения повернуты на
$90\degr$ и ось~$Y$ направлена вниз) по каталогам. Определите погрешность своих расчетов (звездных
величин и координат).
\end{enumerate}
Для определения каталожных величин воспользуйтесь aladin, xephem, stellarium (и т.п.) или on-line
сервисами каталогов. Преобразовать координаты~$(x,y)$ на изображении в $(\alpha, \delta)$ можно при
помощи утилиты \t{xy2sky} из пакета \t{wcstools} (либо при помощи \t{ds9}). Не забывайте о
собственных движениях звезд: координаты из каталога (в т.ч. те, которыми мы можем воспользоваться в
\t{ds9}) привязаны к эпохе J2000, а с тех пор прошло уже два десятилетия.
\end{document}