Устранение шумов в изображении и другие многомерные вариационные задачи

Temesgen Kindo 21/09/2018
Share this on Facebook Share this on Twitter Share this on LinkedIn

Ранее мы обсуждали, как решать одномерные вариационные задачи с помощью программного пакета COMSOL Multiphysics® и применять сложные граничные условия с использованием стандартной системы наложения ограничений. Здесь мы продолжим обсуждение этой темы, но теперь уже для многомерных задач, производных высокого порядка и нескольких неизвестных на примере, который, мы надеемся, вам понравится — устранение шумов на изображении. Мы завершаем эту серию статей о вариационных задачах некоторыми рекомендациями по дальнейшим исследованиям.

Вариационные задачи с бо́льшим количеством пространственных измерений

Мы рассмотрели одномерные задачи, определяя функцию u(x), которая минимизирует функционал

E[u(x)] = \int_a^b F(x,u,u^{\prime})dx.

Теперь давайте увеличим в задаче количество пространственных измерений, но ограничимся производными первого порядка. Рассмотрим вариационную задачу минимизации

E[u(x,y,z)] = \int_{\Omega} F(x,y,z,u,u_x,u_y,u_z)d\Omega,

,

определенную в фиксированной области \Omega.

Для близкой функции u+\epsilon\hat{u} получим

E[u+\epsilon\hat{u}] = \int_{\Omega} F\left(x,y,u+\epsilon\hat{u},\frac{\partial}{\partial x}(u+\epsilon\hat{u}),\frac{\partial }{\partial y}(u+\epsilon\hat{u}),\frac{\partial }{\partial z}(u+\epsilon\hat{u})\right)d\Omega.

Как и в линейном случае необходимое условие экстремума первого порядка —

\frac{d}{d\epsilon}\bigg|_{\epsilon=0}E[u+\epsilon \hat u] = 0,

из чего следует

\int_{\Omega} \left[\frac{\partial F}{\partial u}\hat{u}+\frac{\partial F}{\partial u_x}\hat{u}_x+\frac{\partial F}{\partial u_y}\hat{u}_y+\frac{\partial F}{\partial u_z}\hat{u}_z\right]d\Omega = 0, \quad \forall \quad \hat{u}(x,y,z).

Чтобы получить вариационные производные, мы строим, по сути, функцию скалярного параметра \epsilon и используем методы анализа функций с одной переменной. Более быстро тот же результат можно получить, используя формальный подход и записывая вариационную производную, а именно

\delta E = \int_{\Omega} \delta Fd\Omega=\int_{\Omega} \left[\frac{\partial F}{\partial u}\delta u + \frac{\partial F}{\partial u_x}\delta u_x + \frac{\partial F}{\partial u_y}\delta u_y + \frac{\partial F}{\partial u_z}\delta u_z\right]d\Omega

для фиксированных областей.

Здесь \delta u соответствует \hat{u} в нашей предыдущей записи. Стоит отметить, что мы рассматриваем только вариацию функции и ее производных, поскольку область является фиксированной. Если размеры области сами зависят от решения, в вариационной производной появится дополнительное слагаемое, связанное с изменением положения границы области.

Устранение шумов на изображении: многомерная вариационная задача

За последние десятилетия с помощью вариационных методов были разработаны эффективные и строгие алгоритмы обработки изображений, например, для устранения шумов, пятен, ретуширования дефектов и сегментации. Сегодня мы рассмотрим в качестве примера методы устранения шумов. Один из таких методов называется минимизацией полной вариации. Предположим, у вас есть данные в виде изображений u_0, содержащие шумовые дефекты. В изображении имеются черные и белые вкрапления. Вам требуется восстановить оригинальное изображение u с максимальным качеством. Для этого необходимо удалить шумы в изображении.

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

TV(u) = \int_{\Omega}|\nabla u|^2d\Omega = \int_{\Omega}\nabla u \cdot \nabla ud\Omega.

.

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

f(u)=\int_{\Omega}(u-u_0)^2d\Omega.

.

Раз теперь у нас есть несколько потенциально противоречащих целей, попробуем достичь компромисса, введя функционал

E(u) = \int_{\Omega}\left[ \nabla u \cdot \nabla u + \mu (u-u_0)^2\right]d\Omega

в котором параметр регуляризации \mu определяет, насколько важно сохранить детали, а насколько — устранить шум. Этот положительный параметр устанавливается пользователем.

Теперь мы можем вывести условие оптимальности первого порядка, как в предыдущем разделе. Приравнивая вариационную производную нулю, получаем

\delta E = \int_{\Omega}\left[2 \nabla u \cdot \delta (\nabla u) + 2\mu (u-u_0)\delta u\right]d\Omega = 0 , \quad \forall \quad \delta u.

Чтобы продемонстрировать этот процесс, импортируем изображение в COMSOL Multiphysics и добавим в него случайный шум. Вот что получается, если испортить изображение гуся, которое любезно предоставил мой коллега:

Изображение гуся, испорченное случайным шумом.
Тестовое изображение, намеренно испорченное случайным шумом.

Слагаемое в слабой форме, указанное выше, записано в векторной форме. Для дальнейшего использования запишем его в декартовой системе координат и опустим общий множитель 2.

\int_{\Omega}\left[u_x\hat{u}_x + u_y\hat{u}_y + \mu (u-u_0)\hat{u}\right]d\Omega = 0 , \quad \forall \quad \hat{u}.

Это выражение можно ввести в программное обеспечение COMSOL®, как показано на следующем снимке экрана. На всех границах данные сохраняются в первозданном виде с помощью граничного условия Дирихле u=u_0. Используется параметр регуляризации, равный 1e6. В более простых алгоритмах параметр регуляризации определяется путем проб и ошибок. Мы повышаем параметр, если в получающемся изображении отсутствуют важные детали, и понижаем его при большом количестве шума.

Снимок экрана, демонстрирующий постановку вариационной задачи в двух измерениях в COMSOL Multiphysics.
Постановка вариационной задачи в двух измерениях.

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

Изображение гуся без шума.
Оригинал изображения гуся до устранения шума.

Изображение с устраненным шумом (слева) и оригинальное изображение (справа).

Регуляризация по Тихонову слишком интенсивно сглаживает изображение и не сохраняет важные геометрические элементы, например края и углы. Так называемая модель ROF (модель Рудина — Ошера — Фатеми) с функционалом

E(u) = \int_{\Omega}\left[ |\nabla u| + \mu (u-u_0)^2\right]d\Omega

в этом отношении работает лучше.

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

\delta E = \int_{\Omega}\left[ \frac{\nabla u}{|\nabla u|}\cdot \delta (\nabla u) + 2\mu (u-u_0)\delta u\right]d\Omega = 0 \quad \forall \delta u.

.

Обратите внимание, что абсолютное значение в функционале приводит к заметно нелинейному поведению задачи. Кроме того, знаменатель |\nabla u| в численных итерациях может стать нулевым, что вызовет ошибку при попытке деления на 0. Чтобы избежать этого, можно прибавить к нему малое положительное число. В COMSOL Multiphysics мы часто используем параметр относительной точности чисел с плавающей запятой, равный примерно 2.2204 \times 10^{-16}.

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

Вариационные задачи с производными более высокого порядка

Вариационные задачи с производными высокого порядка используются в высокоточной обработке изображений и других областях. Традиционный пример применения — анализ упругих балок и плит. Например, в эйлеровской теории изгиба балок балка с модулем Юнга E и моментом инерции поперечного сечения I, которая находится под воздействием поперечной силы f, изгибается так, что минимизирует полную потенциальную энергию

\int\left[\frac{1}{2}EI(\frac{d^2u}{dx^2})^2-fu\right]dx.

.

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

\delta \int\left[\frac{1}{2}EI(\frac{d^2u}{dx^2})^2-fu\right]dx = \int\left[EI\frac{d^2u}{dx^2}\delta(\frac{d^2u}{dx^2})-f\delta u\right]dx = 0 \quad \forall \delta u .

Обратите внимание, что мы не дифференцируем модуль Юнга в вариационной производной, поскольку рассматриваем линейно упругий материал, свойства которого не зависят от деформации. Для нелинейных материалов необходимо добавить в вариационную производную слагаемое, связанное со свойствами материала. Эта функция встроена в модуль Нелинейные конструкционные материалы, который является расширением COMSOL Multiphysics.

Включение производных более высокого порядка в функционал принципиально не меняет того, как мы будем определять критерий оптимальности первого порядка, но все же приводит к некоторым вычислительным сложностям. Наша интерполяционная функция конечного элемента, выбранная в разделе Discretization (Дискретизация) в окне настроек узла Weak Form PDE (Дифференциальное уравнение в частных производных в слабой форме), должна иметь степень многочлена не ниже высшего порядка пространственных производных в вариационной форме. Например, в задаче о балке имеется производная второго порядка в вариационной форме, поэтому мы не можем использовать линейные пробные функции. Если бы мы это сделали, все производные второго порядка в нашем уравнении стали бы нулевыми. Потребуется использовать квадратичные функции или функции более высокого порядка.

Вариационные задачи с несколькими неизвестными

До сих пор мы рассматривали минимизацию функционалов, зависящих только от одной неизвестной u. Часто функционалы содержат более одной неизвестной. Предположим, что у нас имеется вторая неизвестная функция v и функционал

E(u,v)=\int_{\Omega}F(x,y,z,u,v,u_x,u_y,u_z,\dotsc)d\Omega.

.

Первая вариация этого функционала —

\delta E=\int_{\Omega}\left[\frac{\partial F}{\partial u}\delta u+\frac{\partial F}{\partial v}\delta v + \frac{\partial F}{\partial u_x}\delta u_x+\frac{\partial F}{\partial u_y}\delta u_y+\frac{\partial F}{\partial u_z}\delta u_z + \frac{\partial F}{\partial v_x}\delta v_x+\frac{\partial F}{\partial v_y}\delta v_y+\frac{\partial F}{\partial v_z}\delta v_z+…\right]d\Omega.

Здесь мы предполагаем, что два поля u и v не зависят друг от друга, так что мы можем считать вариации по двум переменным независимыми. Иногда это не так. Переменные могут быть связаны общим ограничением.

В простейшем случае одна из переменных задается в явном виде относительно другой, например, v = g(u). В таких случаях можно устранить v из задачи, рассмотрев функционал

\mathcal{E}(u)=E(u,g(u)).

В общем случае ограничения заданы в форме g(u,v,\dotsc)=0, и связь между переменными нельзя записать в виде явной функции. В таких случаях мы используем методы, рассмотренные в третьей части серии статей, посвященной применению ограничений. Например, при использовании метода множителей Лагранжа имеется расширенный функционал

\mathcal{L}(u,v,\lambda) = \int_{\Omega}\left[F+\lambda g\right]d\Omega,

в котором \lambda = \lambda(x,y,z) в cлучае распределенного ограничения, или

\mathcal{L}(u,v,\lambda) = \int_{\Omega}Fd\Omega\quad +\quad \lambda g,

где \lambda является константой в случае глобального ограничения.

Для постановки задачи с несколькими полями мы используем тот же интерфейс Weak Form PDE (Дифференциальное уравнение в частных производных в слабой форме). Вопрос в том, использовать ли один интерфейс или по одному интерфейсу на каждое неизвестное? Зависит от соотношения между неизвестными. С одной стороны, если u и v являются разными компонентами одного и того же физического вектора, например, смещения или скорости, можно использовать один интерфейс и указать количество зависимых переменных в разделе Dependent Variables (Зависимые переменные). С другой стороны, если u обозначает температуру, а v — электрический потенциал, то можно использовать разные интерфейсы дифференциальных уравнений в частных производных. Если по какой-то причине приходится использовать разный порядок дискретизации или разные масштабы для компонентов одного и того же неизвестного вектора, тогда тоже можно добавить разные интерфейсы. При наличии нескольких неизвестных необходимо тщательно выбирать интерполяционные функции для каждого поля (как было сказано в предыдущей статье о порядке элементов в мультифизических моделях).

Заключительные замечания

Вариационные методы представляют собой единый подход к моделированию множества научных задач. Интерфейс Weak Form PDE (Дифференциальное уравнение в частных производных в слабой форме), включенный в состав программного пакета COMSOL Multiphysics, расширяет возможности COMSOL® и позволяет вам описывать и решать собственные вариационные задачи. В этой серии статей мы продемонстрировали эти функции, решив задачи о нахождении формы мыльных пленок, определении пути туриста вокруг озера и восстановлении поврежденных изображений.

Имейте в виду, что некоторые важные дифференциальные уравнения с частными производными не возникают в результате минимизации функционала. Один из примеров — уравнение Навье — Стокса. При этом интерфейс Weak Form PDE (Дифференциальное уравнение в частных производных в слабой форме) можно использовать для таких задач, если записать их в слабой форме.

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

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

  • Классическая литература по вариационному анализу:
    • I.M. Gelfand and S.V. Fomin (English translation by R.A. Silverman), Calculus of Variations, Dover Publications, Inc., 1963. Русское издание: Гельфанд И.М., Фомин С.В. Вариационное исчисление. М.: Физматлит, 1961.
  • Более актуальная литература, в которой рассматривается несколько прикладных задач:
    • K.W. Cassel, Variational Methods with Applications in Science and Engineering (Применение вариационных методов в науке и технике), Cambridge University Press, 2013.
    • J.W. Brewer, Engineering Analysis in Applied Mechanics (Инженерные расчеты в прикладной механике), Taylor & Francis, 2002.
  • Для методов применения ограничений:
    • D.G. Luenberger, Introduction to Linear and Nonlinear Programming (Введение в линейное и нелинейное программирование), Addison-Wesley Publishing Company, 1973.
    • S.S. Rao, Engineering Optimization: Theory and Practice (Теория и практика прикладной оптимизации), John Wiley & Sons Inc., 2009.

В этой серии статей мы постарались описать в общих чертах методы моделирования вариационных задач с помощью COMSOL Multiphysics, особенно в случаях, когда объект моделирования еще не встроен в программный пакет.

Желаем успешной работы, а если у вас появятся вопросы, свяжитесь с нами, нажав кнопку ниже:

Изучите другие статьи из серии о вариационных задачах и условиях


Загрузка комментариев...

Темы публикаций


Теги

3D печать Cерия "Гибридное моделирование" Введение в среду разработки приложений Видео Волновые электромагнитные процессы Глазами пользователя Графен Интернет вещей Кластеры Моделирование высокочастотных электромагнитных явлений на различных пространственных масштабах Модуль AC/DC Модуль MEMS Модуль Акустика Модуль Волновая оптика Модуль Вычислительная гидродинамика Модуль Геометрическая оптика Модуль Динамика многих тел Модуль Композитные материалы Модуль Коррозия Модуль Механика конструкций Модуль Миксер Модуль Нелинейные конструкционные материалы Модуль Оптимизация Модуль Плазма Модуль Полупроводники Модуль Радиочастоты Модуль Роторная динамика Модуль Теплопередача Модуль Течение в трубопроводах Модуль Химические реакции Модуль Электрохимия Модуль аккумуляторов и топливных элементов Охлаждение испарением Пищевые технологии Рубрика Решатели Серия "Геотермальная энергия" Серия "Конструкционные материалы" Серия "Электрические машины" Серия “Моделирование зубчатых передач” Сертифицированные консультанты Технический контент Указания по применению физика спорта