Заметки по поводу фильтрации FITS-ов
Mar. 30th, 2012 09:45 am![[personal profile]](https://www.dreamwidth.org/img/silk/identity/user.png){kind=small}
eddy_emВ связи с тем, что я уже пятый день мучаю тех. задание на систему регистрации ИК-спектрофотометра, у меня возникли кое-какие соображения касательно фильтрации шумов на изображении.
Итак, принцип работы КМОП-мультиплексоров таков, что позволяет производить в процессе накопления недеструктивные считывания сигнала. Пошукав по литературе и пообсуждав этот факт, было решено добавить в систему регистрации фильтрацию изображения посредством линейной аппроксимации сигнала в каждом пикселе изображения. А вот как эту аппроксимацию сделать - уже особый разговор.
Допустим, что мы делаем экспозицию в 1 час с квантом времени в 5с на считывание промежуточных изображений. В итоге получаем 3600/5=720 16-битных картинок размером 1024х1024, на хранение которых понадобится 1024*1024*2*720=1.5ГБ (плюс еще какое-то место на «шапки»)! Теоретически, все их можно «запихать» в оперативку, для этого рассчитаем предельный размер оперативной памяти: если экспозиция будет длиться всю ночь (самая длинная ночь у нас — где-то 15 часов, вычтя время на сумерки и всякие подготовки, получим в пределе часов 13), квант времени будет равен 1с, то максимальное количество файлов составит 13*3600=46800 и займут эти файлы как минимум 93ГБ! Есть материнки, поддерживающие 128ГБ оперативки, т.е. теоретически такой компьютер собрать можно, практически же стоить он будет слишком дорого. Поэтому идем другим путем.
Другой путь — складывать все временные файлы на диске (лучше всего - на SSD, для увеличения скорости доступа). Т.е. каждое промежуточное изображение с соответствующей шапкой мы сохраняем на диск (всю шапку заполнять не надо - только меняющиеся данные, вроде температур и т.п. параметров, а также - времени от предыдущего считывания, т.к. выдержать требуемый квант с точностью в микросекунды вряд ли получится; постоянные параметры достаточно занести в шапку только первого файла). В итоге по окончании экспозиции у нас на диске получается n-е количество FITS-файлов, которые можно в фоновом режиме, не мешая наблюдениям, обработать и отдать под утро наблюдателю все его файлы в «чистеньком» виде.
Для обработки нам надо для каждого пикселя проанализировать кривую зависимости уровня сигнала от времени экспозиции (в основном - ее линейную часть, аппроксимирующуюся выражением I=I₀+ΔIt), выдернув оттуда такие параметры, как bias (I₀), произведение gain на освещенность в данном пикселе (ΔI), предел линейности (тот участок графика, где наша кривая начинает «загибаться» вправо) и уровень насыщения (предельный уровень сигнала в пикселе, не изменяющийся при повышении времени экспозиции). Нелинейные параметры, в принципе, наблюдателю неинтересны, поэтому расчет их можно сделать опциональным. Помимо этих характеристик обработка позволит выявить «горячие» пиксели (находящиеся в насыщении вне зависимости от уровня сигнала), «битые» пиксели (вроде «горячих», но с низким уровнем сигнала) и «неправильные» пиксели (с нелинейной характеристикой или прочими проблемами).
Теория, конечно, хороша, но опять вернемся к нашим баранам: современные компьютеры пока еще не позволяют за разумные деньги сделать такие вычисления без ухищрений. Во-первых, cfitsio не позволяет открыть такое количество файлов (а если бы библиотека и позволяла, ее особенности в том, что она, по сути, mmap'ит файлы в память - а то и просто копирует их туда, так что, оперативки не хватит), следовательно, если захотеть открывать эти файлы «разом» и поочередно считывать, придется писать свой FITS-парсер (чего делать не хочется — это во-вторых). В-третьих, можно придумать еще уйму «отмазок», чтобы закончить эти размышления и послать все нафиг. Но работать-то надо!
Итак, единственным более-менее приличным (и легко реализуемым) вариантом предварительной обработки данных мне предоставляется следующий:
До обработки создаем итоговый FITS и заполняем его шапку: заносим всю статику, обрабатываем динамику и заносим туда же экстремальные, средние и медианные значения, а также стандартное отклонение. Далее начинаем заполнять данные изображения:
Да, забыл добавить: благодаря тому, что этот способ очистки компенсирует «перекоп» от слишком ярких объектов (т.к. анализирует лишь линейный участок, а затем умножает крутизну на время экспозиции), данные надо будет сохранять уже как минимум в формате 32-битных целых, а то и 32-битных числах с плавающей точкой. FITS-файлы с дополнительной информацией позволят упростить предварительную обработку изображений, да и всякие калибровки будут уже значительно точнее.
Итак, принцип работы КМОП-мультиплексоров таков, что позволяет производить в процессе накопления недеструктивные считывания сигнала. Пошукав по литературе и пообсуждав этот факт, было решено добавить в систему регистрации фильтрацию изображения посредством линейной аппроксимации сигнала в каждом пикселе изображения. А вот как эту аппроксимацию сделать - уже особый разговор.
Допустим, что мы делаем экспозицию в 1 час с квантом времени в 5с на считывание промежуточных изображений. В итоге получаем 3600/5=720 16-битных картинок размером 1024х1024, на хранение которых понадобится 1024*1024*2*720=1.5ГБ (плюс еще какое-то место на «шапки»)! Теоретически, все их можно «запихать» в оперативку, для этого рассчитаем предельный размер оперативной памяти: если экспозиция будет длиться всю ночь (самая длинная ночь у нас — где-то 15 часов, вычтя время на сумерки и всякие подготовки, получим в пределе часов 13), квант времени будет равен 1с, то максимальное количество файлов составит 13*3600=46800 и займут эти файлы как минимум 93ГБ! Есть материнки, поддерживающие 128ГБ оперативки, т.е. теоретически такой компьютер собрать можно, практически же стоить он будет слишком дорого. Поэтому идем другим путем.
Другой путь — складывать все временные файлы на диске (лучше всего - на SSD, для увеличения скорости доступа). Т.е. каждое промежуточное изображение с соответствующей шапкой мы сохраняем на диск (всю шапку заполнять не надо - только меняющиеся данные, вроде температур и т.п. параметров, а также - времени от предыдущего считывания, т.к. выдержать требуемый квант с точностью в микросекунды вряд ли получится; постоянные параметры достаточно занести в шапку только первого файла). В итоге по окончании экспозиции у нас на диске получается n-е количество FITS-файлов, которые можно в фоновом режиме, не мешая наблюдениям, обработать и отдать под утро наблюдателю все его файлы в «чистеньком» виде.
Для обработки нам надо для каждого пикселя проанализировать кривую зависимости уровня сигнала от времени экспозиции (в основном - ее линейную часть, аппроксимирующуюся выражением I=I₀+ΔIt), выдернув оттуда такие параметры, как bias (I₀), произведение gain на освещенность в данном пикселе (ΔI), предел линейности (тот участок графика, где наша кривая начинает «загибаться» вправо) и уровень насыщения (предельный уровень сигнала в пикселе, не изменяющийся при повышении времени экспозиции). Нелинейные параметры, в принципе, наблюдателю неинтересны, поэтому расчет их можно сделать опциональным. Помимо этих характеристик обработка позволит выявить «горячие» пиксели (находящиеся в насыщении вне зависимости от уровня сигнала), «битые» пиксели (вроде «горячих», но с низким уровнем сигнала) и «неправильные» пиксели (с нелинейной характеристикой или прочими проблемами).
Теория, конечно, хороша, но опять вернемся к нашим баранам: современные компьютеры пока еще не позволяют за разумные деньги сделать такие вычисления без ухищрений. Во-первых, cfitsio не позволяет открыть такое количество файлов (а если бы библиотека и позволяла, ее особенности в том, что она, по сути, mmap'ит файлы в память - а то и просто копирует их туда, так что, оперативки не хватит), следовательно, если захотеть открывать эти файлы «разом» и поочередно считывать, придется писать свой FITS-парсер (чего делать не хочется — это во-вторых). В-третьих, можно придумать еще уйму «отмазок», чтобы закончить эти размышления и послать все нафиг. Но работать-то надо!
Итак, единственным более-менее приличным (и легко реализуемым) вариантом предварительной обработки данных мне предоставляется следующий:
- создать во временной директории (назовем ее пока /tmp) 1024 (в общем, в зависимости от количества строк) поддиректории, в каждой из которых завести по 1024 (в зависимости от количества столбцов) файла;
- открыть очередной FITS-файл;
- прочитать «шапку» и сохранить в оперативке динамические данные (статику возьмем из нулевого изображения);
- перейти к сегменту данных изображения и начать их попиксельное считывание;
- значения интенсивности в пикселе с координатами (X, Y) дописать в файл /tmp/Y/X;
- считать так все пиксели изображения;
- перейти к следующему изображению;
До обработки создаем итоговый FITS и заполняем его шапку: заносим всю статику, обрабатываем динамику и заносим туда же экстремальные, средние и медианные значения, а также стандартное отклонение. Далее начинаем заполнять данные изображения:
- открываем очередной файл /tmp/Y/X;
- считываем все данные в оперативку (или просто mmap'им его);
- получаем графики первой и второй производных;
- из анализа графиков производных находим точки перегиба (определяя таким образом участок «разгона», линейный участок, переход в насыщение и собственно насыщение);
- если мы не имеем в пикселе линейного участка — определяем, «горячий» он или «битый», увеличивая соответствующий счетчик и занося в соответствующий пиксель 65535 или 0 (плюс можно создавать дополнительные изображения-маски, в которых помечать все «косяки») и переходим к следующему пикселю;
- аппроксимируем линейный участок, вычисляем требуемые параметры (I₀, ΔI и т.п.) и заносим в очередной пиксель итогового изображения значение, равное ΔI·T (T — время экспозиции), производя таким образом сразу же очистку от bias'а; сопутствующую информацию можно записывать в дополнительные файлы;
- переходим к следующему пикселю.
Да, забыл добавить: благодаря тому, что этот способ очистки компенсирует «перекоп» от слишком ярких объектов (т.к. анализирует лишь линейный участок, а затем умножает крутизну на время экспозиции), данные надо будет сохранять уже как минимум в формате 32-битных целых, а то и 32-битных числах с плавающей точкой. FITS-файлы с дополнительной информацией позволят упростить предварительную обработку изображений, да и всякие калибровки будут уже значительно точнее.
![[identity profile]](https://www.dreamwidth.org/img/silk/identity/openid.png){kind=small}
no subject
Date: 2012-03-31 04:42 am (UTC)Для таких (и больших) объемов памяти применяют серверные платы, серверную память (с буферизацией и коррекцией ошибок) и серверные процессоры. Стоить будет порядка первых сотен тысяч рублей.
Возможно есть смысл приспособить для вычислений видеокарту (CUDA, OpenCL)? Тогда скорость обработки будет ограничена практически скоростью работа диска.
Во всяком случае, билатеральный фильтр для картинки ~4000х4000 на карточке у меня работает быстрее раз в 20-30.
Карточка старая, Geforce GTX 260, процессор Core Quad 3.2Ггц, память 8Гб.
no subject
Date: 2012-03-31 09:33 am (UTC)Для такой обработки использовать CUDA бессмысленно: копирование между ОЗУ и оперативкой видеокарты будет происходить дольше, чем обработка на CPU. Для более сложной обработки, естественно, я пользуюсь CUDA. Правда, приходится иногда вместо GPU использовать CPU, т.к. оперативки на видеокарте слишком мало (с 512МБ оперативки на видеокарте я даже Фурье не могу сделать для изображения 3Кх3К).
no subject
Date: 2012-04-02 02:33 am (UTC)Телескоп направлен всю ночь в одну точку ?
no subject
Date: 2012-04-02 05:14 am (UTC)А как еще вы построите для каждого пикселя сотни-другой изображений аппроксимацию? Либо открыть все файлы одновременно и поочередно считывать значение каждого пикселя из каждого файла, вычислять, записывать это значение в выходное изображение и переходить к следующему пикселю (но это не позволяет библиотека cfitsio, а свою, как я уже говорил, писать лень); либо поочередно открывать каждый файл с заполнением отдельных файлов для каждого пикселя и последующей обработкой.
> Телескоп направлен всю ночь в одну точку ?
Нет, за ночь обычно наблюдают несколько объектов, но очень слабые требуют длительной экспозиции (вплоть до длительности ночи). Правда, т.к. с ростом времени экспозиции из-за шумов светоприемника ухудшается S/N, приходится делать набор коротких экспозиций с последующим суммированием оных.
Аппроксимация промежуточных недеструктивных считываний в принципе может улучшить S/N (но это еще бабушка надвое нагадала - надо будет сначала изучить светоприемник).