Кристаллография, 2020, T. 65, № 6, стр. 845-850

ВВЕДЕНИЕ

В [1–3] предложен последовательный подход к решению обратной задачи дифракционной рентгеновской топо-томографии на основе полукинематического решения уравнений Такаги–Топена для амплитуды дифрагированной σ-поляризованной волны на примере точечного дефекта кулоновского типа в кристаллической пластине Si(111) в условиях симметричной дифракции Лауэ и набора наклонных двумерных топографических проекций, отвечающих вращению плоскопараллельного кристалла-образца вокруг вектора дифракции [2$\bar {2}$0]. При некоторых a priori ограничениях, наложенных на класс функций поля смещений, минимизационный алгоритм моделирования квазиньютоновского градиентного спуска показал хорошую сходимость решений к глобальному минимуму по данным одной или нескольких 2D-проекций, отвечающих изображению точечного дефекта на рентгеновской топограмме в случае отсутствия шумовой составляющей в данных. Однако при наличии шума достичь глобального минимума в процессе компьютерного моделирования трехмерной функции поля смещений точечного дефекта кулоновского типа в кристалле становится невозможным. Возникает задача фильтрации (предварительной очистки двумерных томографических проекционных изображений от шумов) для того, чтобы улучшить качество восстановления при решении обратной задачи дифракционной томографии (ДТ).

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

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

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

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

Использование такого подхода позволило выбрать метод фильтрации, наиболее подходящий для целей последующей 3D-реконструкции изображения дефекта.

ОЧИСТКА ПРОЕКЦИОННЫХ ИЗОБРАЖЕНИЙ ОТ ШУМА

Для изучения возможностей фильтрации были подготовлены тестовые наборы проекционных изображений, полученных на основе расчетных изображений точечного дефекта (рис. 1а). Особенностью данных изображений является практически постоянное значение интенсивности вне небольшой окрестности дефекта. Изображение дефекта представляет собой два минимума интенсивности, разделенных только одной точкой пространственной сетки пикселей. Размер области перехода от минимума интенсивности к постоянному значению составляет несколько пикселей.

К тестовым изображениям добавлен трехпроцентный аддитивный гауссовский шум с нулевым средним значением (рис. 1б). Добавление шума приводит к случайным отклонениям интенсивности от постоянного значения на плато и непредсказуемым статистически независимым изменениям значений в области дефекта. Указанные особенности хорошо видны на рис. 2, где приведены профили интенсивностей исходного теоретического и зашумленного изображений вдоль прямых, проведенных через область дефекта (рис. 2а) и на удалении от него в области постоянных значений интенсивности (рис. 2б).

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

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

где ${{f}_{{ij}}}$ и $\widetilde {{{f}_{{ij}}}}$ – интенсивности исходного и выбранного для сравнения с ним изображений в точке, стоящей на пересечении i-й строки и j-го столбца матрицы изображения, сумма берется по всем N точкам изображения. Величина среднеквадратичного отклонения для зашумленного изображения представлена в табл. 1.

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

С этой целью использована величина отклонения (1), вычисленная по прямоугольной симметричной окрестности дефекта шириной 11 и высотой 9 пикселей, покрывающей область перехода от минимальных значений интенсивности к постоянному значению вне дефекта (рис. 2а). Качество фильтрации периферической части изображения оценивается по той же формуле (1), но в этом случае сумма берется по точкам, лежащим вне окрестности дефекта. Полученные величины приведены в табл. 1. Для зашумленного изображения они различаются незначительно.

Для сглаживания отвечающих случайному шуму резких изменений интенсивности часто используют пространственные усредняющие фильтры [4], в которых новое значение интенсивности изображения в пикселе, расположенном в i-й строке и  j-м столбце матрицы изображения, представляет собой взвешенное среднее значение, определенное по его окрестности размером n × m пикселей (2). Указанную окрестность называют окном фильтра.

Здесь $f({{x}_{j}},~{{y}_{i}})$ – значение интенсивности изображения в точке с координатами $({{x}_{{j~}}},~{{y}_{i}})$, w(s, t) – нормированные на единицу веса, составляющие маску фильтра – матрицу размером n × m, n = 2a + + 1, m = 2b + 1, n, m – нечетные числа. Суммирование проводится по пикселям окна фильтра, s и t – координаты точек окрестности, отсчитанные от ее центра.

Из-за описанных выше особенностей тестовых проекционных изображений результаты применения к ним сглаживающих фильтров оказываются непригодными для реконструкции 3D-изображения даже в случае минимально возможного размера окна фильтра (3 × 3 пикселя). Это хорошо видно на примере фильтра Гаусса. Выбор минимально возможного размера окна фильтра обусловлен стремлением минимизировать его влияние на область дефекта, имеющую ключевое значение для 3D-реконструкции.

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

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

Другим популярным классом фильтров, применяемых для очистки изображений от шума, являются пространственные фильтры, основанные на порядковых статистиках [4]. Это нелинейные не усредняющие фильтры. Новое значение интенсивности в пикселе, расположенном в центре окна фильтра, выбирается после ранжирования интенсивностей в точках окна, исходя из места в этой упорядоченной последовательности. Фильтр максимума использует для этого максимальное значение, фильтр минимума – минимальное, а наиболее популярный медианный фильтр – значение, расположенное посередине. Медианный фильтр отлично подавляет шумы и при этом меньше “размазывает” изображение, чем сглаживающие фильтры с таким же размером окна.

Оценить результат применения медианного фильтра с окном 3 × 3 пикселя к зашумленному проекционному изображению позволяют данные табл. 1. Как и в случае гауссовского фильтра, результат фильтрации непригоден для использования в 3D-реконструкции изображения дефекта. Несмотря на уменьшение значения среднеквадратичного отклонения всего фильтрованного изображения от исходного, отклонение в окрестности дефекта существенно выросло. В то же время фильтр отлично сгладил периферийную область, почти вдвое уменьшив рассчитанное для нее отклонение. Изменения, внесенные в изображение медианным фильтром в области дефекта, исказили его значительно больше, чем шум, который необходимо удалить. Причина этого в малом размере области дефекта. Ее площадь оказалась сравнимой с минимально возможной площадью окна фильтра. Известно, что при использовании медианного фильтра наибольшим изменениям подвергаются области изображения, площадь которых меньше половины площади окна фильтра [4].

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

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

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

Популярным фильтром, сохраняющим границы объектов, является билатеральный фильтр [5, 6]. Для формирования значений интенсивности фильтрованного изображения он использует выражение вида (2), но входящие в него веса w(s, t) зависят не только от взаимного пространственного расположения пикселей в окне фильтра, но и от того, как сильно различаются их интенсивности. Максимальный вес отвечает соседним пикселям, имеющим близкую интенсивность.

Такой выбор весов позволяет сохранить на изображении границы объектов, вблизи которых соседние пиксели существенно различаются по интенсивности и, как следствие, практически не влияют друг на друга при фильтрации. К недостаткам билатерального фильтра относят возникновение артефактов в местах резкого изменения интенсивности (gradient reversal artifacts) в случае, когда вблизи этих точек есть группа пикселей с близкими значениями интенсивности [7].

Указанного недостатка лишен так называемый управляемый фильтр (Guided Image Filter, GIF) [7], успешно использующийся в задачах сглаживания и выделения краев объектов, удаления шумов и артефактов [8], [9]. Кроме того, этот метод фильтрации реализует более быстрый алгоритм по сравнению с билатеральным фильтром.

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

Пусть ${{I}_{i}}$, ${{f}_{i}}$ и ${{g}_{i}}$ – значения интенсивностей опорного, исходного и фильтрованного изображений в пикселе с номером i. Для упрощения формул использована одномерная индексация пикселей изображений. Интенсивность фильтрованного изображения в i-м пикселе представим в виде линейной функции интенсивности опорного изображения, определенной на квадратном окне ${{{\omega }}_{k}}$ вокруг этого пикселя:

где индекс k нумерует все квадратные окрестности заданного радиуса r, покрывающие i-й пиксель (рис. 3). Предполагается, что параметры ${{a}_{k}}$ и ${{b}_{k}}$ постоянны внутри окна ${{{\omega }}_{k}}$.

Параметры линейной зависимости определяют из требования наилучшего соответствия полученного значения ${{g}_{i}}$ интенсивности исходного изображения ${{f}_{i}}$, минимизируя по этим параметрам целевую функцию

Здесь ε – регуляризационный параметр, введенный для предотвращения появления слишком больших значений ${{a}_{k}}$.

Стандартная линейная регрессия дает решение в следующем виде:

где ${{{\mu }}_{k}}$ и ${{{\sigma }}_{k}}$ – среднее значение и дисперсия интенсивности опорного изображения в окне ${{{\omega }}_{k}}$; $\left| {\omega } \right|$ – количество пикселей в окне, равное ${{(2r + 1)}^{2}}$; ${{\bar {f}}_{k}} = \frac{1}{{\left| {\omega } \right|}}\sum\nolimits_{i \in {{{\omega }}_{k}}}^{} {{{f}_{i}}} $ – среднее значение интенсивности исходного изображения в окне ${{{\omega }}_{k}}$.

Поскольку значения ${{g}_{i}}$ различаются для различных окон, включающих в себя i-й пиксель, значение интенсивности фильтрованного изображения в этой точке определяют как среднее по всем окнам ${{{\omega }}_{k}}$. С учетом симметрии окон относительно i-го пикселя это значение выражается через средние значения коэффициентов по окну ${{{\omega }}_{i}}$, в центре которого расположен i-й пиксель:

Выражения (3), (4) позволяют реализовать управляемый фильтр, параметрами которого являются радиус окна r и регуляризационный параметр ε. Согласно (3), если дисперсия интенсивности в границах окна много больше ε, результирующее значение будет мало отличаться от исходного. В противном случае результирующее значение будет равно среднему значению интенсивности в окне. Таким образом, в случае больших изменений интенсивности в границах окна управляемый фильтр почти не изменяет значение в пикселе, а в случае малых изменений выходное значение равно среднему по окну. Величина изменения интенсивности определяется сравнением с величиной регуляризационного параметра.

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

При таком выборе параметров удалось добиться значительного уменьшения среднеквадратичного отклонения фильтрованного изображения от исходного незашумленного как в окрестности дефекта, так и вне ее (табл. 1, рис. 1в).

Полученный в результате применения управляемого фильтра набор очищенных от шумов проекционных изображений был с успехом использован для реконструкции 3D-поля смещений точечного дефекта кулоновского типа в кристалле Si (111). В табл. 2 приведены результаты компьютерного моделирования данных ДТ. Как видно из табл. 2, качество 3D-реконструкции поля смещений дефекта (контрольный параметр – CP) улучшается более чем в 2 раза после применения управляемого фильтра (со значения 0.0296, получаемого без применения фильтрации, до значения 0.0145). Отметим, что полученный эффект выигрыша в качестве 3D-реконструкции данных ДТ является важным и существенным для практического применения разработанных ранее теоретических алгоритмов на основе полукинематического решения уравнений Такаги–Топена при анализе экспериментальных данных ДТ (как в случае рентгеновских измерений, так и при использовании метода электронной микроскопии).

ЗАКЛЮЧЕНИЕ

Проведено сравнение эффективности алгоритмов фильтрации шумовой составляющей 2D-проекций дифракционной томографии и выбран оптимальный алгоритм – управляемый фильтр (Guided Image Filter), позволяющий существенно улучшить качество восстановления 3D-поля смещений точечного дефекта кулоновского типа при решении обратной задачи дифракционной томографии. Результаты компьютерного моделирования выглядят достаточно многообещающими для дальнейшего использования алгоритма фильтрации шумов при анализе экспериментальных данных дифракционной электронной и рентгеновской томографии.

Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН в части по применению томографических алгоритмов и Российского фонда фундаментальных исследований (проект № 19-02-00556 А) в части развития подходов по решению проблемы фильтрации шумов в данных дифракционной томографии.