Видимость сквозь турбулентную атмосферу. Компьютерная коррекция изображений удаленных объектов

Авторский пересказ двух публикаций с демонстрационным фильмом.

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


Астрономическое отступление.


Каждый, кто смотрел в телескоп с большим увеличением, например, на Луну, наблюдал дрожащее, расплывчатое и переменчивое изображение, лишь изредка «проясняющееся» в моменты замирания воздушных колебаний. Области пространства с разной плотностью воздуха можно представить себе в виде множества слабых неправильных случайных линз, постоянно возникающих и исчезающих с большой частотой.
В астрономии для исправления изображений, искаженных атмосферной турбулентностью, применяют классические методы адаптивной оптики (датчики волнового фронта и быстродействующие деформируемые зеркала). Это достаточно сложные и дорогие системы, доступные только для больших и сверхбольших телескопов, однако дающие существенный положительный эффект. Вот показательный пример (взят из: Ralf Flicker. Methods of Multi-Conjugate Adaptive Optics for Astronomy, Lund Observatory, Sweden, 2003).



Левая фотка – участок звездного неба, сфотографированный с длинной экспозицией. Случайные искажения, накладываясь друг на друга, приводят к значительному размытию изображения, в результате чего многие звезды просто не видны, а другие сливаются друг с другом. В центре – тот же участок неба, сфотографированный при использовании адаптивного зеркала. Датчик волнового фронта, нацеленный на яркую (опорную) звезду чуть ниже центрального скопления, измеряет искажения волнового фронта в нескольких точках входной апертуры. Компьютер вычисляет и подает управляющие сигналы через высоковольтные усилители на пьезоэлектрические толкатели. Эти толкатели, надавливая на гибкое зеркало сзади в разных местах, искривляют его так, чтобы компенсировать искажения. Такая процедура происходит сотни раз в секунду. Таким образом можно улучшить изображение в некоторой области вблизи опорной звезды. Как видно из средней фотки, на краях изображение не улучшается, потому что в разных частях изображения существуют другие искажения, отличные от тех, которые измеряет датчик в центре. С помощью всяких ухищрений (многосопряженных адаптивных систем с несколькими опорными звездами и зеркалами) можно добиться расширения области улучшенного изображения – фотка справа. Кстати, угловой размер этой области на небе – всего 90 секунд, или 1/20 лунного диска.

Для подзорной трубы, используемой в «земных условиях», разумеется, нет смысла применять столь сложную технику, как адаптивная оптическая система. Не только из-за сложности и дороговизны, но и потому, что отсутствуют яркие опорные точки, по которым можно было бы определить искажения волнового фронта. В отличие от звезд, удаленные приземные объекты, как правило, малоконтрастны и протяженны. И еще потому, что турбулентность в приземном слое существенно выше, чем в высокогорных обсерваториях.
Поэтому вычислительные методы улучшения видимости удаленных объектов в земных условиях, основанные на постдетекторной обработке изображения, имеют, я полагаю, неплохую перспективу.
Сначала рассмотрим, какие бывают искажения.

Упрощенная модель формирования изображения в атмосфере.


В первом приближении можно принять модель формирования изображений объекта, учитывающую прохождение светового потока через случайно-неоднородную среду, описанную в [1]. Она включает в себя три фактора, по-разному искажающих конечное изображение.

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

Второй фактор — мелкомасштабные флуктуации плотности, приводящие к случайному (для каждого кадра) искажению мелких деталей изображения. Умные люди [2,3], основываясь на теории турбулентности атмосферы и статистической оптике, предлагают следующий способ. Поскольку из-за стохастического характера флуктуаций исправить каждый кадр невозможно, надо применить статистику. Усредненные по времени искажения можно теоретически вычислить и получить вполне определенную для данного состояния атмосферы оптическую передаточную функцию. Иначе говоря, усредненное изображение есть результат действия специфического для атмосферы (негауссова) сглаживающего фильтра. Применив к изображению обратный фильтр, теоретически мы должны получить картинку дифракционного качества.

Третий фактор — дифракционные искажения, связанные с конечностью приемной апертуры и с которыми ничего сделать нельзя. Идеальное изображение точки всегда будет не точкой, а диском с кольцами (так называемая картина Эйри).

Таким образом, если у нас есть последовательный набор изображений объекта, случайно искаженных турбулентной атмосферой, то для восстановления изображения необходимо последовательное выполнение следующих операций. Измерение случайных смещений изображения, исключение этих смещений, усреднение центрированных изображений, вычисление обратной оптической передаточной функции, фильтрация усредненного изображения. Производя эти операции непрерывно с каждым последующим кадром, постепенно получаем исправленное изображение [4]. Применяя процедуру к различным участкам кадра, можно вычислить оптический поток (вектор смещения каждого пикселя как функцию времени) на всем изображении и расширить корректируемое поле зрения [5].

Алгоритм.


Коррекция смещений. В опорном кадре выбираем квадрат размером в несколько десятков пикселей, внутри которого есть более-менее контрастные детали. В текущем кадре надо найти квадрат такого же размера, лучше всего совпадающий с опорным. Сравнение квадратов производится прямым вычислением взаимно-корреляционной функции S(p,q):

S(p,q) = ∑ mod[In(x+p,y+q) — I0n(x,y)],

где In(x,y), I0n(x,y)– яркости пикселей текущего n-ного кадра и опорного; x,y – координаты пикселей квадрата; p,q – координаты области поиска, примерно равной по величине опорному квадрату; суммирование по всем пикселям квадрата.
Минимум S(p,q), означающий место наилучшего совпадения, находится по известному методу «два шага» поиска экстремумов многомерных функций. Координаты минимума pmin, qmin определяют вектор смещения текущего кадра относительно опорного. Затем квадрат из текущего кадра сдвигается к опорному и суммируется с ним, образуя новый опорный кадр. Так происходит накопление полезного сигнала в опорном квадрате. Суммирование производится рекурсивно с коэффициентом R порядка 0,01–0,05:

I0n+1(x,y) = (1 — R)*I0n(x,y) + R*In(x+pmin,y+qmin)

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

Фильтрация. Или устранение остаточного размытия полученного усредненного изображения, обусловленного влиянием мелкомасштабной турбулентности. Это остаточное размытие, например, для одной яркой точки, представляет собой пятно с колоколообразным распределением яркости, ширина которого зависит от интенсивности турбулентности и в несколько раз больше дифракционного диска. Роль фильтрации состоит в преобразовании этого размытого пятна в дифракционно-ограниченную картину Эйри. Фильтрация изображения и получение конечного результата I(x,y) производится путем двумерного прямого Фурье-преобразования (оператор F) усредненного изображения I0n+1(x,y), умножения полученного спектра на обратную оптическую передаточную функцию H(x,y), вычисленную в [2,3], и обратного Фурье-преобразования (оператор F-1):

I(x,y) = F-1[F(I0n+1(x,y))*H(x,y)],

где
H(x,y)] = [arccos ( r ) – r*(1- r2)1/2]*exp[3.44*rr5/3*(1- r1/3)];
r = (x2+y2)1/2/D;
rr = (x2+y2)1/2/d;
D – величина, зависящая от относительного отверстия оптической системы;
d – величина, зависящая от интенсивности турбулентности (параметр Фрида для атмосферы).

Техническая реализация.


В качестве оптической системы применялись телеобъективы МТО-1000, Варио-Гоир-1Т, Sigma с диаметрами входной апертуры от 70 до 140 мм и фокусными расстояниями от 200 до 5000 мм. Пока у меня не было скоростной видеокамеры, приходилось использовать камеру RT-1000DC фирмы «Raster Tecknolodgy» и наблюдать только неподвижные объекты. Основные результаты в фильме получены с ней. Как видно в эпизоде с номером машины, необходимо ждать несколько секунд, пока изображение «проявится».
Эпизоды с вышкой и зданием – попытка реализации расширения корректируемой области на весь кадр. Пришлось работать с заранее записанными видеофайлами, поскольку объем вычислений очень большой, и вообще там еще много нерешенных проблем.
Последний эпизод фильма («Метель») снят совсем недавно с сетевой камерой RM-6740 фирмы JAI, частота кадров 200 гц, время получения исправленного изображения менее 1 сек. В этом эпизоде наблюдается дополнительный эффект уменьшения шумов (в виде падающего снега) в области коррекции.



Литература

1. Оптика атмосферы и океана, 11, 522 (1998).
2. Гудмен Дж. Статистическая оптика (М.: Мир, 1988).
3. Fried D.L. J. Opt. Soc. Am., 56, 1372 (1966).
4. Квантовая электроника, т.40, №5 (2010), с.418-420.
5. Квантовая электроника, т.41, №5 (2011), с.475-478.
Share post
AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 59

    +6
    отличные результаты, респектус
      +2
      Безумно круто. Оно наверно ещё и шумы матрицы скушает неплохо если с большими iso снимать.
        +2
        А не поделитесь реализацией? А то есть у меня камера и 3M-5CA тоже интересный объектив. По итоговой резкости иногда превосходит mto100. Хочется попробовать что-нибудь такое.
          0
          Если камера скоростная, то есть смысл попробовать. Иначе уже просто неинтересно
          0
          Очень даже неплохо убирает шумы, без всяких дополнительных фильтров. Сказывается избыток информации при высокочастотной съемке
          +13
          Впечатлило практическое применение распознавания номера машины аж за 2 километра.
            +3
            Ух ты ж, как круто!!! А камеры с такими «антишоками» когда начнут выпускать? :) Уже хочу!
              –2
              Вы сможете себе позволить камеру за пару миллионов долларов? =D
                0
                Серьезно?!?? Почему такие суммы?
                  0
                  Скоростная светочувствительная камера + активная оптика.
              +1
              На сколько усложняется модель для RGB?

              Результаты — очень впечатлили!
                +1
                Не сильно. Ввиду слабой зависимости турбулентных искажений от длины волны можно все цвета обрабатывать одинаково. Ну нет у меня подходящей цветной камеры. Собственно, из-за воздушной дымки насыщенность цветов сильно падает, поэтому большой разницы между цветной и черно-белой картинкой не будет
                +4
                Автор молодец! Мне периодически приходится фотографироваь через марево — результат как правило печален. Хотя мои объекты довольно крупные, но марево тем не менее дает существенные искажения.

                Не пробовали сравнивать свой результат с AviStack, Registax, DeepSkyStacker?
                  0
                  Плюсую. Использовал регистакс, вроде суть та же. Непонятно чем готовый продукт хуже/лучше того, что в статье.
                    0
                    Чуть ниже я написал, чем лучше и хуже
                    0
                    Регистакс у меня есть, естественно, сравнивал. В принципе, результаты одинаковые, если долго ковыряться в его многочисленных настройках и смотреть только финальный кадр. Достоинство Регистакса в обилии всяких опций (программист все-таки писал), основные недостатки — фильтрация, настраиваемая вслепую, «методом тыка», ненадежное определение смещений (я проверял, если вычислять корреляцию через фурье, будет хуже!) и, разумеется, конечный результат в виде одной фотки, а не видеопоток в реальном времени. Вообще Регистакс, как я понимаю, предназначен для ручной творческой работы, причем специально с астрономическими объектами.
                      0
                      Корреляция через преобразования фурье должна теоретически давать в точности тот же рузультат, что и прямая свёртка. В вашей статье же написано про сумму модулей. Правильно ли я вас понял, что результаты с суммой модулей были лучше, чем с корреляцией?
                        0
                        Да, я почти уверен, что лучше. Для массивов размером менее 100х100 пикселей еще и считается быстрее. Для больших массивов (я пробовал делать стабилизацию участка 512х512), естественно, с фурье считается намного быстрее и вполне точно.
                  • UFO just landed and posted this here
                      +2
                      К плагину для обеспечения нужной эффективности необходима соответствующая оптическая система и сопряженная с ней скоростная камера, установленные на опорно-поворотном устройстве
                        0
                        Насколько скоростная?
                        А вы можете рассказать про «опорную звезду», формируемую лазером?
                          0
                          Для почти полной компенсации атмосферных искажений в активной оптике необходимо отслеживать волновой фронт с частотой не менее ~100Гц (а лучше — больше). Использовать реальные звезды на таких частотах невозможно (слишком мало света), поэтому-то и «зажигают» опорную искусственную звезду. Датчик волнового фронта (чаще всего это — датчик Шака-Гартманна) формирует матрицу углов наклона волнового фронта по зеркалу, по которой генерируется набор управляющих сигналов для пьезокерамических актуаторов активной оптики вторичного зеркала. Так как первый коэффициент Цернике (наклон) почти полностью компенсирует атмосферное дрожание изображения, коррекцию коэффициентов более высоких порядков обычно не проводят (это будет сильно дороже).
                          Примерная стоимость такого счастья (активная оптика с искусственной звездой) не менее ста миллионов рублей (как-то были разговоры о внедрении на БТА, но из-за слишком высокой цены отказались — если начать внедрять, то реально цена наверняка намного больше выйдет).
                            0
                            Ну, скажем, 500 кадров в секунду. Будет видно значительно лучше, чем к меня.

                            Лазерная звезда — это свечение слоя атомов натрия в стратосфере, возбуждаемое лазером. Насколько я понимаю, она выглядит ярким точечным объектом, похожим на настоящую звезду, но зажигать ее можно в любом нужном нам месте (в поле зрения телескопа)
                              0
                              Ну, скажем, 500 кадров в секунду

                              Такое реально лишь для искусственной «звезды».

                              «Лазерная звезда» на самом деле отнюдь не точечный объект: это скорее эдакая «игла». Однако, угловые размеры ее действительно невелики.

                              А «зажигать» ее нужно как можно ближе к реальному изображению объекта, иначе свет от объекта и свет от искусственной звезды пройдут через разные неоднородности и ничего не получится.
                                0
                                На лазерную звезду датчик ВФ смотрит практически вдоль луча лазера, а апертура у него обычно гораздо меньше, чем у телескопа, так что для него область свечения будет почти точка
                              +1
                              Кстати, к вопросу «как работает лазерная звезда». Вот небольшой фильм, имитирующий работу адаптивной системы, снятый в лабораторных условиях.
                              В качестве объекта использовалась спираль от электролампочки, а лазерной звездой был, собственно, лазер. Атмосфера имитировалась нагретым чайником, стоящим перед лампочкой и создающим вертикальные потоки воздуха. Расположив звезду поближе к спирали, можно добиться, чтобы «атмосферные» искажения были примерно одинаковы как для звезды, так и для объекта. Датчика волнового фронта и гибкого зеркала у меня не было, поэтому пришлось искажения для каждого кадра считать программно. Левый верхний квадрат (ЛВ) – это исходное изображение объекта (будем считать, что это «галактика») и опорной лазерной звезды чуть правее и ниже. Правый верхний квадрат (ПВ) – то, что видит «датчик ВФ» (только искаженную опорную звезду, или, иначе говоря, мгновенную функцию рассеяния точки). Правый нижний квадрат (ПН) – идеальное изображение точки для данной оптической системы, то, что должно быть в отсутствие искажений. Эту картинку (ПН) легко вычислить. Дальше надо вычислить ОПФ для каждого кадра как отношение фурье-образов ПВ и ПН. Зная ОПФ, делим фурье-образ ЛВ на ОПФ, берем обратное фурье и получаем левый нижний квадрат (ЛН), т.е. исправляем все точки изображения «галактики» в каждом кадре. Еще надо имитировать длительную экспозицию, как в астрономии, т.е. просто сложить все кадры и сравнить, что получилось:

                              img-fotki.yandex.ru/get/6436/46744491.0/0_8194b_cd5fef4c_XXL.jpeg
                              Видим, что разрешение «галактики» в ЛН по сравнению с ЛВ улучшилось в несколько раз, а опорная звезда, как и следовало ожидать, вообще идеально восстановилась, даже дифракционное кольцо видно.
                              В настоящей адаптивной системе датчик ВФ, сравнивая ПВ и ПН, должен выдавать на гибкое зеркало такие сигналы, чтобы ПВ, отраженный в гибком зеркале, получился такой же, как ПН. Тогда автоматически и ЛВ, отраженный в том же зеркале, преобразуется в ЛН. Т.е. суть одинакова, только реализация механическая, а не программная.
                          0
                          В астрофизике такие методы называются superresolution. На практике не используются почти никогда, т.к. шаг
                          умножение полученного спектра на обратную оптическую передаточную функцию H(x,y), вычисленную в [2,3]

                          приводит к делению на ноль (в худшем случае) или сильным ошибкам вблизи нулей и полюсов передаточной функции (в лучшем).

                          Для задач распознавания объектов (когда сохранение распределения интенсивности излучения — вопрос второстепенный) этот метод еще кое-как можно применять, но и в этом случае основная проблема — нахождение передаточной функции. А найти ее можно только используя специальную калибровочную миру. Причем калибровку нужно делать как можно чаще, т.к. атмосфера довольно быстро меняется. Без калибровки можно лишь немного улучшить разрешение, обрабатывая ряд изображений.
                            0
                            Сверхразрешение — это не совсем то же самое. Если вы заметили, у меня в ОПФ стоит дополнительный множитель с арккосинусом. Он там стоит не просто так, а для того, чтобы не было «сильных ошибок вблизи нулей».

                            Усредненная атмосферная передаточная функция в принципе зависит только от радиуса Фрида, который можно периодически вычислять по дисперсии дрожаний. Я, правда, пока так не делаю, а двигаю ползунок на форме в небольших пределах.
                              0
                              Это хорошо, когда у вас есть несколько кадров. Но в реальности получить за одну ночь несколько кадров невозможно: даже за те десятки минут экспозиции у вас могут измениться атмосферные условия.
                              А вообще, фотометрия не столь интересна, как спектроскопия. Но там уже совсем другие экспозиции: для многих объектов даже ночи не хватит, чтобы нормальный спектр получить.
                                0
                                На астрономические наблюдения я не претендую, хотя и совсем не чужд этому. И спектрометр у меня есть
                                  +1
                                  Ну я и тугодум. Только сейчас до меня дошло, что речь идет не о получении научных данных, а об обычном любительском астрофото! Там — да, приходится сильно ухищряться, т.к. любителю нереально потратить миллионы на телескоп, светоприемник и приличный фотометр (я уже молчу об активной оптике с искусственной звездой).
                                    0
                                    Почему миллионы. Вот Meade 250 мм у меня на работе есть, всего 150 тыс., полный автомат
                                      0
                                      Науки на миде вы не сделаете (кроме исследования атмосферы или наблюдения за метеорами). Нужен как минимум метровый телескоп.
                                      +1
                                      Я — совсем тугодум! Тут, оказывается, еще и замечательное видео было. Пока в хромом не открыл — не замечал (т.к. оно на флеше, а не html5).

                                      Замечательные примеры! Я как раз с другом лет 6 назад подобными вещами пытался заниматься (а он, кстати, даже немного смог заработать, распознав лицо человека сквозь пелену табачного дыма на видеосъемке).
                                      Действительно, ваши бы труды — да в библиотечку!
                              0
                              здорово! молодцы!
                                +3
                                Наверняка кто-то зайдёт в коменты, чтобы найти эту картинку:
                                image
                                  +2
                                  Бредятина, разумеется. В кино можно сделать все, что угодно
                                  +1
                                  Добавляйте в OpenCV
                                    0
                                    Это я не умею
                                      0
                                      Выложите код на github, вам помогут, например я. Если есть желание, конечно.
                                      Если алгоритм актуален, есть математическое описание и пример реализации, думаю сообщество примет. Это тоже спросить можно.
                                    0
                                    вот бы на самом деле реализовать это аппаратно. в камере или телефоне.
                                      +1
                                      Нет смысла. Обычно камеры в телефонах с их микроскопическими объективами просто не замечают никакой турбулентности. Искажения, как правило, имеют величину порядка секунд, или десятков секунд (угловых), что для маленьких объективов недоступно.
                                      0
                                      А нельзя ли формулами как то показать вашу модель искажения/получения данных? А то совершенно непонятно что именно вы компенсируете…

                                      Например, что то вроде этого Y(x,y) = X(x+a,y+b) + c, где Y — наблюдаемое изображение, X — истинное, а a,b,c — случайные величины какого-то распределения.
                                        0
                                        Такой формулой, как вы написали, нельзя. Со всякими интегралами и усреднениями — наверное, можно
                                        0
                                        Это офигенно!
                                          +1
                                          Круто! Думаю, то что-то подобное делает человеческий глаз (постоянно, подсознательно), принимая на себя различные виды вибрации.
                                            0
                                            Именно так. В этом легко убедиться, если сравнить какой-нибудь один нескорректированный стоп-кадр из фильма с исправленным. Разница бывает очень большая. В видеопотоке глаз в некоторой степени подсознательно делает что-то подобное, как мне кажется.
                                              0
                                              А можно такой вопрос: как вы определяли смещение изображения при наклоне (деревья в метель)?
                                              Обычные смещения легко кросс-корреляцией получить, а вот с наклонами дело куда хуже…
                                                0
                                                Ничего особенного не делал. Смещение определялось так, как написано. Очевидно, вклад фона в корреляционную функцию был больше, чем у тех веток, которые наклоняются. Я думаю, из-за большой частоты кадров изображение успевает заново сформироваться и «размазывание» веток при их поворотах просто малозаметно.
                                                  0
                                                  Понятно. Просто я уже давно интересуюсь способом, позволяющим получить значение не только смещения объекта, но и поворота. Увы, пока работает только велосипед, базирующийся на ключевых точках (точность его, понятное дело, крайне невелика).
                                                    0
                                                    Вы можете разбить окно на части, для каждой вычислить смещение с помощью корреляции, например, затем проинтерполировать вектор смещения, и применить это проинтерполированное смещение к исходному изображению. Такие итерации повторяются несколько раз. Такой подход применяется в PIV (particle imaging velocimetry). Другой вариант — считать optical flow.
                                                      0
                                                      Тогда окно придется делать слишком мелким (чтобы поворот не влиял на ККФ), а это приведет к значительным ошибкам отождествления.
                                                      А optical flow тоже не является универсальной штукой.
                                                        0
                                                        Поток, по-моему, еще и считать очень долго
                                                          0
                                                          Зависит от размера изображения и мощности железа. В принципе, если бы данные из ОЗУ компьютера в ОЗУ видеокарты и обратно перемещались бы хотя бы на пару порядков быстрей, то вполне можно было бы на мощной видеокарте обсчитывать 720p с частотой кадров 5 в секунду.
                                                          0
                                                          Для этого и нужны итерации. На первой мы определяем смещения грубо. Но зато устраняем бОльшую часть смещений и поворотов, влагодаря чему на следующей итерации всё отождествляется уже намного точнее.
                                              0
                                              А такая коррекция работает только на «статичных» изображениях или движущиеся объекты тоже нормально выглядят?
                                                0
                                                Если использовать камеру, скажем, 400 кадров/сек, то и движущиеся (не слишком быстро) объекты можно наблюдать примерно с таким же улучшением. Например, идущего человека. Правда, возникнут новые проблемы — надо ускорить обработку, компьютер помощнее, сделать другой канал передачи данных, побыстрее, чем 1-Гбитный Ethernet.
                                                  +1
                                                  У таких скоростных камер главная проблема — интерфейс. Невозможно сделать «full-HD» камеру, которая могла бы снимать видео с частотой ~400 кадров в секунду в течение длительного времени. Те, что есть, имеют буфер (скажем, 64ГБ оперативной памяти DDR3), в который и скидывают накомпленное, т.е. реально они снимают всего лишь несколько секунд.
                                                  0
                                                  Статья — супер!

                                                  Only users with full accounts can post comments. Log in, please.