Добрый день.
О методах деконволюции изображений уже сказано так много, что добавить, кажется, больше нечего.
Однако всегда есть алгоритм, который лучше и новее предыдущих.
Не так давно был описан итерационный алгоритм, который имеет линейную скорость сходимости при небольших затратах памяти, стабилен и хорошо распараллеливается.
А через некоторое время она была улучшена до квадратичной сходимости.
Встречайте: (быстрый) итеративный алгоритм порогового усадки.
Постановка задачи
Довольно давно была опубликована целая серия статей об алгоритмах деконволюции изображений.В цикле рассмотрен классический метод решения систем линейных алгебраических уравнений и их регуляризация.
Описанный метод с точки зрения постановки задачи практически не отличается от классического, но дьявол кроется в деталях.
Давайте начнем нашу магию математики.
Вот как выглядит классическое выражение для полученного испорченного и зашумленного изображения: 
Матрица 
представляет собой матрицу свертки и с физической точки зрения показывает взаимосвязь между точками исходных и поврежденных данных.
В случае изображений эту матрицу можно получить путем преобразования ядра свертки.
Вектор 
И 
представляют собой исходное и поврежденное изображение в векторизованной форме.
Другими словами, все столбцы пикселей изображения объединяются один за другим.
Первая проблема при такой постановке задачи — это объем выделяемой памяти.
Если изображение имеет размеры 640х480, то задача хранит два вектора изображения размерами 307200х1 и матрицу размером 307200х307200 элементов.
Матрица свертки разрежена, но ее размеры даже в разреженном виде будут большими и занимать много памяти.
Произведение с матрицей свертки можно заменить двумерной сверткой с ядром, что не требует хранения большой матрицы в памяти, это решит проблему с хранением матрицы
.
Для поиска изображения, максимально близкого к исходному, запишем оптимизируемую функцию в классической форме второй нормы несоответствия поврежденного и искомого изображений.
Отличие от классических методов
Данная постановка задачи является одной из самых популярных, но мы дополним ее для повышения производительности.Используемый подход называется «Майоризация-максимизация» и заключается в замене исходной задачи другой, очень похожей.
Новая задача имеет два свойства:
- Задача гораздо проще
- Во всех точках, кроме текущей, она имеет расхождение большее, чем исходная задача.
Звучит довольно сложно, намного сложнее, чем кажется на самом деле.
Окончательная функция минимизации для итерации 
написано следующим образом: 
Матрица правого члена положительно определена.
Это достигается за счет 
.
Функция минимизации состоит из суммы двух величин, каждая из которых неотрицательна.
В то же время в точке 
второе слагаемое равно нулю, благодаря чему выполняется второе условие из списка.
Поиск решения
Начнем поиск решения классическим методом: раскроем скобки, выпишем градиент и приравняем его нулю.
Важно отметить, что это градиент
- эта итерация.
Внимание, много математики! 
В результате производная будет иметь следующий вид: 
Следующий классический шаг в этой ситуации — приравнять градиент нулю и из полученного выражения выразить искомый вектор 
.
Записано
давайте определим это как 
или как решение в следующей итерации.
Полученное выражение называется «итерацией Ландвебера».
Это основная формула между итерациями.
Используя его, вы можете гарантировать, что расхождение уменьшается от итерации к итерации с линейной скоростью.
Базовый алгоритм содержит дополнительный шаг, называемый «мягкий порог», который предполагает разреженность решения.
Этот шаг приравнивает к нулю все значения искомого вектора по модулю меньше заданного значения.
Это снижает влияние шума на результат реконструкции.
Как это выглядит 
Как видно из решения, у нас есть два параметра, которые мы корректируем по своему вкусу.
отвечает за точность аппроксимации, ограничивая сходимость порогу.
отвечает за стабильность работы и скорость сходимости алгоритма.
Модификация алгоритма
Впоследствии алгоритм был усовершенствован дополнительным термином, аналогичным по идее методу сопряженных градиентов.Добавляя на каждом шаге разницу между результатом двух предыдущих итераций, мы увеличиваем сходимость алгоритма к квадратичной.
Итоговая процедура алгоритма следующая:
- Установить параметры
- Набор
- Повторяем базовый алгоритм
- Обновить временные переменные
- Добавить сопряженный компонент к переменной
Цикл расчета порога на каждой итерации может быть выражен в виде поэлементных произведений (произведение Адамара) и операций сравнения, которые включены как в OpenCL, так и в CUDA, поэтому их можно легко распараллелить.
Алгоритм я реализовал в Matlab в двух вариантах: для расчета на CPU и на GPU. Позже я стал использовать только GPU-версию кода, так как она более удобна.
Код представлен ниже.
Код здесь
Давайте ближе к практике и проверим работу алгоритма в зависимости от размера картинки.function [x] = fista_gpu(y,H,Ht,lambda,alpha,Nit) x=Ht(y); y_k=x; t_1=1; T=lambda/alpha; for k = 1:Nit x_old=x; x=soft_gpu((y_k+(1/alpha)*Ht(y-H(y_k))), T); t_0=t_1-1; t_1=0.5+sqrt(0.25+t_1^2); y_k=x+(t_0/t_1)*(x-x_old); end end function [y] = soft_gpu(x,T) eq=ge(abs(x),abs(T)); y=eq.*sign(x).*(abs(x)-abs(T)); end
Конфигурация моего ноутбука Intel ядро I3-4005U NVIDIA GeForce 820M Слева график производительности алгоритма за 100 итераций в зависимости от количества пикселей для CPU и GPU.
- Как видно из результатов, алгоритм GPU практически не зависит от размера задачи и для действительно больших изображений мы ограничены только памятью GPU. Полсекунды для 2000х2000 (без учёта выгрузки из памяти) — достойный результат, вам не кажется?
- Справа – зависимость времени, затраченного алгоритмом, от количества итераций.
По мере увеличения количества итераций алгоритм на графическом процессоре начинает реагировать увеличением времени своей работы.
Скорее всего, это связано с моими ошибками в кодировании, либо с вынужденным снижением рабочей частоты графического процессора.
В большинстве случаев достаточно 100 итераций.
На следующих графиках показано сравнение базового и улучшенного алгоритмов, а также зависимость сходимости алгоритма от
И
параметры.
- Как видно из графика слева, улучшенный алгоритм сходится значительно быстрее базового и после ста итераций практически не показывает прироста точности.
Базовый алгоритм показывает приемлемый результат только к 500-й итерации.
- По центральному графику видно, в зависимости от
порог сходимости алгоритма изменяется, что ограничивает производительность.Это необходимый компромисс, если сигнал очень зашумлен.
- На правом графике видно, что с увеличением
алгоритм сходится гораздо медленнее.Как и в предыдущем случае, существует компромисс между «качеством» матрицы свертки и скоростью сходимости.
И напоследок пример алгоритма работы с FHD картинкой, Исходное изображение
Испорченная и шумная картинка
Восстановленное изображение
Как видите, результат вполне хороший, особенно если это занимает всего 15 секунд, из которых 1,5 — работа алгоритма и 13,5 — выгрузка изображения из памяти графического процессора.
Весь код, используемый для алгоритма, размещен в GitHub .
Рекомендации
- И.
Селесник.
(2009) О проекте: Восстановление разреженного сигнала.
- А.
Бек и М.
Тебулль Алгоритмы на основе быстрых градиентов для решения задач шумоподавления и устранения размытия изображений с ограниченной совокупной вариацией.
IEEE ТРАНЗАКЦИИ ПО ОБРАБОТКЕ ИЗОБРАЖЕНИЙ, VOL. 18, НЕТ.
11 НОЯБРЯ 2009 ГОДА
- П.
Л.
Комбетт и Ж.
-К.
Методы проксимального разделения Песке при обработке сигналов
-
Менделеев Дмитрий Иванович.
19 Oct, 24 -
Google Принял Решение О Рекламе На Youtube
19 Oct, 24 -
Трекер Thepiratebay Закрывается...
19 Oct, 24 -
Твиттер+Веб-Камера
19 Oct, 24
