Development of an automated procedure for solving the problem of reconstructing blurry digital images

Korobeinikov Anatolii Grigor'evich Doctor of Technical Science professor, Pushkov institute of terrestrial magnetism, ionosphere and radio wave propagation of the Russian Academy of Sciences St.-Petersburg Filial 199034, Russia, g. Saint Petersburg, ul. Mendeleevskaya, 1 Korobeynikov_A_G@mail.ru Other publications by this author     Fedosovskii Mikhail Evgen'evich PhD in Technical Science Ph.D. in Technical Sciences, Professor, Head of the Systems and technological safety technology Department, St. Petersburg National Research University of Information Technologies, Mechanics and Optics, General Director of "Diakont" 197101, Russia, Saint Petersburg, Kronverskii prospekt, 49 diakont@diakont.com Other publications by this author

Aleksanin Sergei Andreevich graduate student, St. Petersburg State University of Information Technologies, Mechanics and Optics 197101, Russia, Saint Petersburg, Kronverkskii Prospekt, 49 Aleksanin@diakont.com Other publications by this author

Достаточно часто в реальных условиях на цифровых изображениях присутствует “смаз”, который возникает при взаимном относительном движении оптического сенсора и объекта во время экспозиции. В этом случае наблюдаемое изображение является результатом наложения со смещением последовательности входных изображений. У некоторых исследователей присутствует мнение, что смаз является необратимой операцией. Вследствие чего информация безвозвратно теряется. То есть любой пиксель отображается в пятно и все смешивается. А случае большого радиуса размытия получается однородный цвет на всем изображении. Но это не совсем правильная позиция. Дело в том, что информация перераспределяется по определенному закону и поэтому ее можно восстановить с некоторыми ограничениями. Другими словами можно сказать так – в ходе процесса искажения любой пиксель исходного изображения отображается в случае расфокусировки в пятно, и в случае простого смаза – в отрезок. Но можно сформулировать эту задачу и для искаженного изображения. В этом варианте любой пиксель искаженного изображения есть образ пикселей какого-то множества в исходном изображении. Именно соответствие какому-то закону, по которому собирается один пиксель и называют функцией искажения. Другие названия – синонимы: PSF (Point spread function, т.е. функция распределения точки), ядро искажающего оператора, kernel и другие. Размерность этой функции, обычно меньше размерности самого изображения

Операцию применения PSF–функции к другой функции (в данном случае к изображению) называют сверткой (convolution) и обозначают как «*». То есть в результате операции «*» некоторое множество в исходном изображении сворачивается в один пиксель искаженного изображения. Операцию, являющеюся обратной по отношению к свертке, называют деконволюцией (deconvolution).

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

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

Для решения задачи улучшения изображения (удаления смаза) для начала рассмотрим математическую постановку этой проблемы.

Пусть за время экспозиции τ камера с ПЗС-матрицей совершила равномерное и прямолинейное смещение (сдвиг) со скоростью v= const вдоль некоторого направления на величину Δ = vτ. Зададим неподвижную прямоугольную систему координат ηOξ. Направление оси ξ совместим с направлением сдвига. В этом случае изображение на ПЗС-матрице станет смазанным (смещенным, сдвинутым) по оси ξ. А отсюда возникает задача о восстановлении (реконструкции) истинного (неискаженного изображения), зная смазанное изображение, направление и величину смаза Δ.

Дополнительно введем еще одну подвижную прямоугольную систему координат xOy. Пусть в начальный момент система координат xOyсовпадает с системой координат ηOξ, причем ось x направлена вдоль оси ξ. Далее, пусть на выбранную точку ПЗС-матрицы с координатами (x,y) за время экспози­ции τ проектируется непрерывное множество точек P с абсциссами ξ =x до ξ=x+Δ с разными интенсивностями f(ξ,y). Интегральная интенсивность g(x,y) в точке (x,y) равна сумме (интегралу) интенсивностей:

(1)

Перед интегралом множитель 1/Δ служит для того, чтобы в случае отсутствия смаза (Δ→0) выполнялось g(x,y) → f(x,y). Кроме того, в случае постоянного изображения ( f(x,y) = const) всегда будет выполняться условие g(x,y) = const.

Перепишем уравнение (1) в следующем виде:

(2)

где g(x,y) – измеренное изображение, а Δ – величина смаза;

f(x,y) – искомое распределение интенсивности на неискаженном изображе­нии (то есть та интенсивность, которая была бы на изображении в случае отсутствия смаза, т.е. при Δ = 0).

Уравнение (2) – это базовое уравнение при решении задачи реконструкции смазанных изобра­жений ^[1-5]. В этом уравнении ось x направлена вдоль смаза. Необходимо отметить, что данное уравнение является одномерным интегральным уравнением типа Вольтерра I рода относительно f(ξ,y) при любом фиксированном значении y, которое играет роль пара­метра. Другими словами, можно сказать, что (2) является множеством одномерных инте­гральных уравнений. Следовательно, учитывая помехи на изображении, (2) можно записать так ^[5]:

(3)

где δy – помеха.

Следуя терминологии ^[6], данное уравнение будет являться неклассическим уравнени­ем Вольтерра I рода. Это связано с тем, что пределы интегрирования являются переменными. Кроме того, это уравнение можно отнести к нестандартным, так как в нем не содержится в явном виде ядро (но можно считать ядро равным 1/Δ =const.)

Необходимо отметить, что значение Δ часто априори неизвестно и его обычно определя­ют путем подбора на основе визуальной оценки получаемых для ряда значений f(ξ,y) решений ^[7-10].

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

В результате, правильно выбрав направление оси х (вдоль смаза) и величину сма­за Δ, т.е. определив параметры смаза, можно, решив интегральное уравнение (3) (точнее, совокупность уравнений), восстановить в принципе неискаженный снимок (истинное изображение) f(x,y) по искаженному снимку–изображению g(x,y).

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

Достаточно часто оценку величины Δ, а также направле­ние смаза определяют по штрихам на снимке. Этот метод особенно эффективен если хотя бы один из этих штрихов есть результат смазывания яркой точки на изображении.

Изображение, полученное в результате смаза, может быть представлено в виде исходного изображения, к которому применена PSF–функция. В таком случае задача восстановления смазанного изображения будет состоять из двух этапов: нахождение функции импульсного отклика и получение исходного изображения по найденной функции ^[11,12]. Если смаз инвариантен к сдвигу, он может быть смоделирован с помощью конволюции (свертки) скрытого изображения с ядром сдвига, где ядро описывает след от датчика. Тогда удаление смаза может быть произведено при помощи операции деконволюции. Существует несколько подходов к восстановлению изображений. При «неслепой» деконволюции ядро сдвига известно, и задача состоит в том, чтобы восстановить скрытое изображение из искаженного, используя это ядро – классические линейные методы восстановления. Существует много методов деконволюции, применяемых в такой ситуации, таких как обратная фильтрация, фильтр Винера, фильтрация методом наименьших квадратов, рекурсивный фильтр Калмана и принудительные итеративные деконволюционные методы ^[13].

Во многих практических ситуациях функция смаза часто неизвестна, а количество информации о реальном изображении очень мало. Поэтому исходное изображение должно быть восстановлено непосредственно из смазанного с использованием неполной информации о процессе смаза. Такие оценочные подходы, в которых модель искажения предполагается линейной, называют слепой деконволюцией ^[14,15]. Методы слепого восстановления содержат много проблем, не имеющих однозначного решения, поэтому в этой области возможно проведение многочисленных исследований. Слепое восстановление изображений применяется в различных технических областях, таких как астрономические изображения, дистанционное зондирование, рентгенография, оптика, фотография, приложения высоких разрешений, приложений отслеживания движений и многих других ^[16].

Слепая деконволюция одиночного объекта является некорректной задачей, так как число неизвестных превышает число известных данных. Проблемой существующих методов является неустойчивость решения, поэтому для получения точного решения задачи восстановления необходимо использовать априорную информацию об изображении, такую как неотрицательность изображения и параметрическую форму функции импульсного отклика ^[11]. Чтобы уменьшить некорректность слепой деконволюции, накладываются такие ограничения на смаз, как линейность параметризированного движения. Также смаз предполагается однородным по всему изображению ^[14]. На данный момент не существует универсальных алгоритмов решения проблемы. Все алгоритмы или дороги в вычислительном отношении или требуют большой предварительной обработки. Некоторые из них полагаются на экспертные знания пользователя или нуждаются в существенных взаимодействиях с ним для алгоритмической калибровки ^[17]. Из-за плохо изложенной природы частой проблемой алгоритмов являются нехватка стабильности, надежности и сходимости ^[13,18].

Некорректная природа проблемы предполагает, что для измеренного сигнала или должны быть введены дополнительные предположения ^[15]. В ходе исследований в качестве дополнительного условия предполагалось, что Δ неотрицательно и его величина мала по сравнению с размером изображения.

Блок-схема автоматизированной процедуры выбора метода обработки смаза представлена на Рисунке 1.

Для решения (3) применим метод называемый «Тихоновской регуляризацией». Эквивалентные названия: «фильтрация по Тихонову», «сглаживающая фильтрация методом наименьших квадратов со связью».

В связи с тем, что задача восстановле­ния изображений является некорректной, то нельзя получить точного решения основ­ного интегрального уравнения (3), которое устойчиво к малым возмущениям во входных данных. Поэтому возникает задача поиска некоторого при­ближенного решения. Это решение можно получить, ограни­чив тем или иным образом полосу частот входного сигнала и выбирая какое-нибудь единственное приемлемое решение из всего множества возмож­ных решений. Но при произвольном выборе приближенного решения достаточно просто можно получить ненужное решение. Задача избавления зависимости решения от произвольных факторов и определения общих методов решения некорректных задач сделало актуальной разработку новых решений и алгоритмов. Фунда­ментальное значение при разработке новых решений имеют математические понятия регуляризации решения и регуляризующего оператора, которые ввел еще в 1943 г. А. Н. Тихонов ^[19].

Первое понятие, а именно регуляризация решения заключается в построении семейст­ва обратных операторов, которые зависят от некоторого число­вого параметра α, который называется параметром регуляризации. Любой из операторов этого семейства позволяет найти решение корректной за­дачи. Необходимо отметить, что при согласованном стремлении к нулю пара­метра α и ошибок во входных данных решение корректной задачи будет сходится к истинному решению соответствующей некорректной задачи.

Другими словами, если в некорректной задаче Ax=y вме­сто точной правой части мы имеем элемент y_прибл`in` G для ко­торого расстояние r(y,y_прибл) ≤ g , то элемент x_g `in` X можно определить с помощью оператора, зависящего от параметра α, значе­ния которых надо брать согласованными с погрешностью g входных данных y_прибл.

Эта согласованность должна быть такой, что при при­ближении правой части y_прибл к точному значению y (при g → 0) приближенное решение x_g стремилось бы к искомому точному решению x (в метрике пространства X). Таким образом, идея регуляризации сводится к поиску такого оператора R(y,α) , где α – параметр регуляризации, который, действуя на правую часть уравнения Ax=y приводит к решению которое не слишком сильно отличается от точного решения. Таким образом, основной смысл регуляризации заклю­чается в “хорошем” вводе неоднозначности в решение, которое за­ведомо не превосходит заданных границ.

Главной особенностью метода регуляризации А. Н. Ти­хонова, отличающей его от других общих методов решения некорректно поставленных задач, является его применимость в ситуации, когда класс X возможных реше­ний уравнения Ax=y не является компактом. Эта ситуация как раз характерна для задач восстановления изображе­ний, в которых обратный оператор основного интегрально­го уравнения (3) заведомо не является непрерывным.

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

Само решение записывается следующим образом:

(4)

где α – параметр регуляризации;

A(u,v) – Фурье-преобразование оператора Лапласа;

H(u,v) – искажающий оператор (функция);

G(u,v) – Фурье–образ искаженного изображения;

|H(u,v)|² =H^*(u,v)H(u,v) ;

H^*(u,v) – комплексно-сопряженная функция H(u,v);

В данной постановке необходима информация о параметре регуляризации α. Поэтому в ходе исследований были проведены численные эксперименты для определения параметра регуляризации в (4).

В первом эксперименте было взято изображение, представленное на Рисунок 2.

Рисунок 2. Исходное изображение

Исходное изображение было обработано процедурой смаза с длиной в 10 пикселов. К результату был добавлен однопроцентный гауссовский шум (Рисунок 3).

Рисунок 3. Зашумленное изображение

Восстановление изображения происходило при помощи метода квадратур с регуляризацией Тихонова: 1 – с размытыми краями и переопределенной системой линейных алгебраических уравнений (СЛАУ), 2 – с усеченными краями и СЛАУ.

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

(5)

где f(i,j) – исходное изображение;

f_α(i,j) – восстановленное изображение при заданном α.

На рисунке Рисунок 4 представлена в логарифмическом масштабе зависимость σ от параметра регуляризации α. Метод 1 представлен кривой 1, и, соответственно, метод 2 представлен кривой 2. В результате оптимальные значения α (то есть для минимального σ в случае конкретного метода) таковы:

- для метода 1: α =0.05011 `Sigma`_rel(α) =0.15958

- для метода 2: α =0.02511 `Sigma`_rel(α) = 0.08582.

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

Рисунок 4. Зависимость σ от параметра регуляризации α.

На Рисунок 5 представлены результаты восстановления зашумленного изображения методом 1 и 2.

Рисунок 5. Результаты восстановления зашумленного изображения методом 1 и 2

Во втором эксперименте определялись параметры регуляризации для различных видов изображений. На Рисунок 6 представлена фотография с трещиной на графите.

Рисунок 6. Исходное изображение

Трещина

На Рисунок 7 представлена фотография, показывающая расстояние между графитовыми стержнями.

Рисунок 7. Исходное изображение.

Расстояние между стержнями

На Рисунок 8 представлена фотография цветка (Роза). И, наконец, на Рисунок 9 представлена фотография акватории Невы с Метеором (Кораблик).

Рисунок 8. Исходное изображение Рисунок 9. Исходное изображение

Роза Кораблик

Все изображения были обработаны процедурой смаза с длиной в 10 пикселов. К результатам был добавлен полуторапроцентный гауссовский шум. Результаты представлены на Рисунок 10 – 13.

Рисунок 10. Смазанное изображение Рисунок 11. Смазанное изображение

Трещина Расстояние между стержнями

Рисунок 12. Смазанное изображение Рисунок 13. Смазанное изображение

Роза Кораблик

Далее для всех изображений методом 2 рассчитывалось `Sigma` _rel(α) при различных α. Результаты расчетов представлены в логарифмическом масштабе на Рисунок 14.

Рисунок 14. Результаты расчетов `Sigma` _rel(α) в зависимости от параметра регуляризации для различных видов изображений

Полученные в результате расчетов оптимальные значения α_опт и для разных видов изображений представлены в таблице.

Таблица.

Оптимальные значения α_опт и `Sigma`_rel(α_опт) для разных видов изображений

Восстановленные изображения представлены на Рисунок 15 – 18.

Рисунок 15. Восстановленное изображение Рисунок 16. Восстановленное изображение

Роза Кораблик

Рисунок 17. Восстановленное изображение Рисунок 18. Восстановленное изображение

Трещина Расстояние между стержнями

На восстановленных изображениях появились эффекты на краях (эффект Гибса). С ним можно бороться различными методами, например, регуляризирующие множители Ланцоша и Фейера ^[20]. Но, так как нас интересует в основном внутренняя область, то для решения задачи обнаружения дефектов, эти эффекты не имеют существенного значения.

В статье представлена разработанная автоматизированная процедура выбора метода обработки смаза. Все методы базируются на решении неклассического уравнения Вольтерра I рода. Используя разработанную процедуру, была произведена оценка параметра регуляризации при восстановлении “смазанных” изображений для различных видов изображений используя две модификации метода квадратур. Был сделан вывод о выборе предпочтительного метода.

Все программное обеспечение было реализовано при помощи системы компьютерной алгебры ^[21].