Простой конвейер для обработки FITS-файлов

[eddy_em]

Для обработки FITS-файлов обычно используются пакеты MIDAS или IRAF. Однако, они довольно тяжеловатые для сравнительно простых задач (которые возникают у меня при технических наблюдениях), а также тяжеловаты для написания своих собственных методик обработки изображений. Поэтому я решил написать простой конвейер, позволяющий удалять шумы с изображений, вычислять и вычитать фон, выделять объекты, строить изофоты и т.д., и т.п.
В результате получилась вот такая штука. Пока еще это — только пре-пре-пре-альфа версия, добавлять нужно еще очень много разных вещей. Но уже кое-что оно умеет (хотя и не исключены сегфолты на малейший чих, тестов пока маловато проведено).

Итак, для образца возьмем длиннощелевой спектр:

и будем всячески его преобразовывать.

Немного отвлекусь: сама утилитка на данный момент имеет следующие параметры командной строки:

./fitsread -h

Usage: fitsread [args] [outfile prefix]

	Where args are:

  -h, --help           отобразить эту справку
  -i, --infile=arg     входной файл
  -o, --outfile=arg    output file
  -c, --conveyer       set pipeline parameters, arg: type=type:[help]:...
		type - transformation type (help for list)
		help - list of available parameters for given 'type'
  --rewrite            перезаписать выходной файл, если он существует (только с опцией -i)
  -v, --verbose        уровень подробностей вывода (каждый -v увеличивает его)
  -s, --stat           отобразить статистические параметры входного и выходного изображения

Обязательный аргумент -i позволяет задать входной файл (пока один, дальше посмотрим); выходной файл можно задать точным именем (флаг -o, с ним же можно использовать флаг --rewrite, чтобы перезаписать выходной файл, если он уже есть) либо (как в apogee_control) префиксом (без флагов, просто свободный параметр); -c задает параметры конвейера, этот флаг может повторяться сколько угодно раз, про него чуть дальше; -v — для вывода всяких сообщений, пока не используется (этот флаг инкрементальный, как в некоторых утилитах, т.е. количество -v определяет уровень информативности); -s показывает статистику по изображению: до обработки и после нее.
Обработка параметров командной строки выполняется моим велосипедом на основе getopt_long (допиленным для возможности "многоразового" использования ключей). Еще есть желание сделать возможным конвейерную обработку вообще всех ключей (чтобы, скажем, можно было на определенном этапе изображения добавить к нему что-нибудь или, скажем, умножить на что-нибудь и т.п.), но пока я не придумал, как это осуществить. Возможно, обойдусь и без этого (используя промежуточные файлы, правда, тогда нужно будет в промежуточных сохранять изображение в виде double, чтобы не было ошибок округления).
Русификация выполнена как обычно — через gettext. Правда, пока еще не совсем до конца. Кстати, неудачное название fitsread надо будет поменять на что-нибудь поприличней.

Итак, ключ -c задает параметры конвейера. Обязательный параметр type определяет вид преобразования. Чтобы получить информацию о всех доступных видах, пишем type=help:

./fitsread -c type=help
Pipeline parameters:
	median       median filter
	adpmed       simple adaptive median filter
	lapgauss     laplassian of gaussian
	gauss        gaussian
	sobelh       horizontal Sobel
	sobelv       vertical Sobel
	simplegrad   simple gradient (by Sobel)
	prewitth     Prewitt horizontal - simple derivative
	prewittv     Prewitt vertical
	scharrh      Scharr (modified Sobel) horizontal
	scharrv      Scharr vertical
	step         posterisation

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

./fitsread -c type=lapgauss:help
Conversion lapgauss  parameters:
sx,sy	sigma by axes x & y
w,h	non-zero window width & height

и видим: преобразование "лапласиан гауссианы" имеет в качестве параметров горизонтальную и вертикальную полуширину (можно задавать "овал" для обработки данных с характерной асимметрией) — обязательные параметры, а также размер окна фильтра (это ненулевая область, которая собственно и заполняется аппаратной функцией фильтра перед его Фурье-преобразованием — свертка у меня выполняется через Фурье, используя библиотеку fftw3_threads) — если этот параметр опустить, окно будет всего-то 5х5 пикселей, что при больших полуширинах фильтра приведет к совсем неожиданным действиям (вместо ожидаемого получится элементарное усреднение).

Вернемся теперь к картинкам.
Медианный фильтр по изображению крайне редко используется в астрофизике, т.к. искажает данные (само собой, медианная фильтрация "стопки" изображений — наоборот, довольно распространенное явление). Но в случае, скажем, с датчиком Шака-Гартманна мы можем и медиану применить. Вот что получится при обработке оригинального спектра медианным фильтром 3х3 пикселя (в параметрах задается радиус фильтра, размер окна вычисляется как 2*r+1):

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

Как видим, кое-где данные черт-те во что превращены, шумы тоже не везде сглажены. Есть еще над чем поработать.
Иногда для последующей обработки требуется размыть изображение гауссианой. Вот как выглядит спектр, пропущенный через гауссов фильтр с полушириной 5 пикселей по обеим осям и размером окна 50х50 пикселей:

А вот для выделения объектов зачастую необходимо применить лапласиан гауссианы. Вот так выглядит лапласиан гауссианы нашего спектра с полушириной 3 пикселя по обеим осям и размером окна 40х40 пикселей:

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

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

Понятно, что нужно знать величины этих самых уровней постеризации. Пока они просто выводятся в терминал (в дальнейшем, видимо, буду загонять их в таблицу внутри итогового изображения). Вот так были рассчитаны уровни для этого случая:

    0: -7.22045
    1: 4.89835
    2: 53.5787
    3: 249.124
    4: 1034.62
    5: 4189.89
    6: 16864.4
    7: 67776.9

В данном случае у изображения есть и отрицательные значения, поэтому чтобы обойти проблему, скажем, логарифмирования, я весь диапазон перед "постеризацией" преобразую к виду 1..(max-min+1), а затем уже его преобразую. Экспоненциальная постеризация по сути преобразует эти новые уровни так, чтобы каждый последующий в N раз был больше предыдущего (т.е. если прологарифмировать полученную шкалу, получим равноотстоящие значения).
Обратная — логарифмическая — постеризация в данном случае показывает сплошную черноту, т.к. уровни "липнут" к верхней части (равноотстоящие значения получаются при вычислении экспоненты от уровней):

    0: 21374.2
    1: 33883.3
    2: 42758.7
    3: 49643
    4: 55267.8
    5: 60023.6
    6: 64143.2
    7: 67776.9

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

Степенная зависимость вида I² тоже подчеркивает меньшие уровни яркости (хоть и не так хорошо, как экспоненциальная):

Но все это — лишь одиночные преобразования. Интересней же конвейер в действии проверить. Пожалуйста. Вот пример двухпроходного конвейера, сначала была выполнена медианная фильтрация по квадрату 7х7 пикселей (для лучшей чистки шумов), а затем вычисление градиента:

Или вот, скажем, сначала медиана по окну 11х11 пикселей, а затем горизонтальный фильтр Прюитта:

Как подначитаюсь Гонсалеса-Вудса, добавлю, наверное, еще какие-нибудь полезности.

Ну и напоследок два примера преобразований моих многострадальных гартманнограмм. Сначала изображение отфильтровано медианой по квадрату 7х7 пикселей, а затем постеризовано. Здесь, в отличие от длиннощелевого спектра, динамический диапазон полезных данных значительно более широкий, поэтому в обоих случаях что-то полезное остается. Постеризация экспонентой:

И постеризация логарифмом:

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