Neural Network Method for Detecting Blur in Histological Images

Full Text

1. ВВЕДЕНИЕ

Современная гистопатологическая диагностика включает в себя извлечение образцов тканей и проведение гистотехнических процедур для подготовки оцифрованного предметного стекла, называемого полнослайдовым изображением (WSI), для дальнейшего исследования. Неправильное обращение с образцом ткани может привести к появлению артефактов (областей, не имеющих диагностического значения, поскольку внешний вид ткани изменен или поврежден) в полученном WSI [26], что может влиять на системы вычислительной патологии и на поставленный диагноз, если их не исключить или не обработать. Недавние работы по WSI с использованием методов глубокого обучения для диагностического прогнозирования в основном основывались на ручном выборе областей интереса (ROI) без артефактов, имеющих диагностическую значимость [27]. Однако ручной выбор трудоемок и требует много времени. И, поскольку артефакты, включенные во время обучения или вывода, могут ухудшить результаты прогнозирования, то необходимо разрабатывать алгоритмы и программы по автоматическому обнаружению артефактов.

Одним из наиболее часто встречающихся является артефакт размытия изображений, возникающий при смещении фокуса (т. е. профиля сканирования) во время сканирования предметных стекол [28, 29]. Два основных подхода оценки уровня размытия областей гистологических изображений основаны на использовании классических математических методов или методов глубокого обучения. Используя подход, основанный на DCNN, в работе [30] были обучены семь современных архитектур на наборе данных FocusPath¹ [25] и было проведено сравнение результатов с классическими методами. DeepFocus [31] использовал аналогичный подход с CNN для анализа размытых областей в WSI с различными красителями. Другие системы оценки качества фокуса, такие как ConvFocus [32] и FocusLiteNN [33], использовали модели на основе CNN для количественной оценки и локализации размытых областей в WSI. В работе [34] предложен метод, основанный на глубоком ядерном обучении, проанализированный, наряду с оригинальными наборами, также на наборах данных FocusPath [25] и TCGA [33].

При анализе результатов на конкретных наборах данных несомненное преимущество имеют методы глубокого обучения. В то же время хорошо известна проблема плохой устойчивости результатов, возникающая при использовании нейросетевых моделей на отличающихся наборах данных, например, при использовании данных, полученных другими сканерами. В связи с этим несомненный интерес представляют нейросетевые методы, основанные, в том числе, на использовании классических моделей – “нейросетевые методы, основанные на физике” (PINN) [35]. В настоящее время по близкой к PINN тематике в задаче размытия гистологических изображений нам известна только статья [36]. В данной работе нами рассмотрен метод PINN, основанный на использовании нейросетевого оператора Фурье (FNO) и базирующийся на классических методах оценки размытия [37, 38] для задачи бинарной классификации уровня размытия гистологических изображений. Структура статьи включает в себя введение, раздел описывающий моделируемый классический метод, краткое описание нейросетевого оператора Фурье и параметры оператора, выбранные для моделирования, результаты работы алгоритма и заключение.

2. ИСПОЛЬЗУЕМЫЙ КЛАССИЧЕСКИЙ АЛГОРИТМ

Моделируемый классический метод поиска размытых фрагментов на гистологических изображениях состоит из следующих этапов:

1. Построение карты размытия с помощью применения высокочастотного многомасштабного преобразования [37].
2. Оценка степени резкости изображения на основе анализа границ объектов [38].
3. Объединение полученных результатов и детекция размытых областей.

2.1. Высокочастотное многомасштабное преобразование

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

Рис. 1. График отсортированных по возрастанию абсолютных значений коэффициентов ДКП для размытых и четких фрагментов гистологических изображений

Обозначим через I входное изображение в градациях серого размером N₁ × N₂ пикселов. Через $P_{i,j}^M$ обозначим фрагмент (патч) изображения I размером 2M + 1 × 2M + 1 с центром в пикселе (i, j), а через $K_{i,j}^M$ – результат применения к $P_{i,j}^M$ ДКП:

$$\begin{gathered} K_{i,j}^M\left(u,v\right)=\frac{c\left(u\right)c\left(v\right)}{2M+1}\sum _{s=0}^{2M}\sum _{t=0}^{2M}{P}_{i,j}^M\left(s,t\right)S^M\left(s,u,t,v\right); \\ {S}^M\left(s,u,t,v\right)=\cos \left(\frac{\pi (2s+1)u}{2\left(2M+1\right)}\right)\cos \left(\frac{\pi (2s+1)v}{2\left(2M+1\right)}\right); \\ c\left(p\right)=\left\{\begin{array}{l}1, p=0\\ \sqrt{2}, p\ne 0\end{array}\right. \end{gathered}$$

где p – переменная, (1)

В предложенном преобразовании вычисленные коэффициенты ДКП для каждого патча делятся на три группы, а именно, на соответствующие низким, средним и высоким частотам [39]. Рис. 2 иллюстрирует эти три группы коэффициентов на примере патча размером 7 × 7 пикселов.

Рис. 2. Иллюстрация коэффициентов ДКП для блока при их разделении на группы соответствия различным частотам

Обозначим через $H_{i,j}^M$ вектор, состоящий из модулей высокочастотных коэффициентов ДКП $H_{i,j}^M$ патча:

$H_{i,j}^M=\left\{\left|K_{i,j}^M\left(u,v\right)\right|:u+v\ge 2M, 0\le u,v\le 2M\right\}$. (2)

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

$L_{i, j}=sort\left(\underset{r=0}{\overset{m}{\cup }}{H}_{i,j}^{M_r}\right)$, (3)

где $H_{i,j}^{M_r}$ патч размером 2M_r + 1 × 2M_r (M_r = 2^r), а $\underset{r=0}{\overset{m}{\cup }}{H}_{i,j}^{M_r}$ означает объединение модулей высокочастотных коэффициентов ДКП, взятых со всех m + 1 разрешений. Таким образом, L_{i, j} является вектором размерности $\sum _{r=0}^m\frac{{\left(2M_r+1\right)}^2+2M_r+1}{2}$, но исходя из наблюдений достаточно использовать лишь первые $\sum _{r=0}^m\left(2M_r+1\right)$ [37] компонентов для оптимизации времени работы метода.

Пусть L_{i, j ;t} – t-й элемент вектора L_{i, j}, а L_t – t-й слой, образованный путем взятия L_{i, j ;t} элемента во всех векторах L_{i, j} каждого пиксела (i, j) изображения I. Размерность L_t совпадает с размерностью исходного изображения.

Результатом предложенного высокочастотного многомасштабного преобразования являются слои L_t, которые далее min-max нормализуются [37]. По мере увеличения номера слоя, слои начинают состоять из больших по абсолютному значению высокочастотных коэффициентов, тем самым все больше структур и границ объектов исходного изображения появляется на слоях.

2.2. Построение карты размытия

К входному изображению в градациях серого I применяется Гауссов фильтр с параметром σ = 0.5 для подавления высокочастотного шума. Обозначим отфильтрованное изображение I за I_g. Далее для I_g вычисляется карта величин градиентов G следующим образом:

$G=\sqrt{{\left(I_g{h}_x\right)}^2+{\left(I_g{h}_y\right)}^2}$, (4)

где

$h_x=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right] и h_y=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$.

Для карты градиентов G вычисляются слои L_t по формулам (1)–(3). В результате искомая карта размытия изображения рассчитывается как

$D_{i,j}=max\left(\left\{{\tilde{L}}_t\left(i,j\right):t=1, ..., \sum _{r=1}^m{M}_r\right\}\right)$; (5)

$0\le i<N_1, 0\le j<N_2$.

Окончательная карта размытия D сглаживается с помощью фильтра [42]. На рис. 3 приведен пример карты размытия. Чем темнее область на карте размытия, тем сильнее размыта эта область на исходном изображении.

Рис. 3. Визуализация карты размытия для гистологического изображения с отмеченными красным областями размытия

3. ОЦЕНКА РЕЗКОСТИ ИЗОБРАЖЕНИЯ

Границей объекта на изображении (edge) называется область неоднородности или область перехода между двумя различными по интенсивности областями изображения.

Оценка резкости изображения S_I вычисляется на основании оценки размытия вдоль границ объектов на изображении. В оригинальной статье [38] рассматривались лишь горизонтальные и вертикальные границы изображения, однако в данной работе было принято решение рассмотреть, в том числе и диагональные границы.

Более подробно об оценке можно прочитать в оригинальной статье [38], а о ее модификации – в работе [44].

4. ОБЩИЙ АЛГОРИТМ КЛАССИЧЕСКОГО МЕТОДА

Дано входное изображение I. Для инвариантности к средней интенсивности изображений это изображение нормируется по параметру value цветового пространства HSV способом, описанным в [45]. Далее это изображение переводится в изображение в градациях серого, и для него строится карта размытия способом, описанным в разделе 2.2. Детекция размытия осуществляется с помощью скользящего окна размером 64 × 64 пиксела с шагом в 32 пиксела как по вертикальной, так и по горизонтальной оси. Тем самым задача детекции размытия переходит в задачу классификации патчей на два класса – резкие и размытые. Патч считается размытым при одновременном выполнении трех условий.

1. Патч принадлежит темной зоне на карте размытия. Если же среднее значение патча (обозначим за V) на карте размытия D превышает 80, то патч считается четким (карта размытия рассматривается в градациях серого).
2. Патч не принадлежит фону. Патчи, средняя интенсивность которых превышает 170, классифицируются как фон.
3. Для всех остальных патчей анализируется пара значений (V, S_I). Патчи из обучающей выборки были сгруппированы по значению V, для каждой из групп путем максимизации метрики F₁ было определено пороговое значение S_I(V), ниже которого патч классифицируется как размытый. Полученный график S_I(V) приведен на рис. 4.

Рис. 4. График пороговых значений S_I(V)

Таким образом, для каждого изображения все его патчи делятся на резкие и размытые.

В пункте 2 был приведен порог интенсивности фоновых патчей равный 170, он был получен следующим образом: на базе набора данных PATH-DT-MSU² WSS2 был создан набор фоновых патчей размера 64 × 64 с предварительной нормализацией параметра value в цветовом пространстве HSV [45] и конвертацией патчей в изображения в градациях серого. Для каждого такого патча считалась средняя интенсивность, а затем для массива полученных значений был высчитан 0.1 квантиль, оказавшийся равным 170. Тем самым патчи, интенсивность которых была больше 170, исключались из рассмотрения.

Общая схема работы предложенного метода приведена на рис. 5.

Рис. 5. Общий алгоритм работы классического метода

5. МЕТОД НА ОСНОВЕ НЕЙРОННОГО ОПЕРАТОРА ФУРЬЕ

Fourier Neural Operator (FNO) – это метод, основанный на данных, для решения задач, описываемых классическими моделями, в том числе уравнениями в частных производных, способный сопоставлять параметры классических моделей с его решением. Его ключевыми особенностями по сравнению с классическими методами решения являются быстродействие и инвариантность к дискретизации сетки. В данной работе рассматривается применение этого метода к задаче аппроксимации классического метода обнаружения размытых регионов на гистологических изображениях.

5.1. FNO

Архитектура сети FNO представлена на рис. 6. Схема преобразований входного изображения, используемая в данной работе:

Рис. 6. Архитектура нейронного оператора Фурье

1. Входное изображение подается в блок P, в котором с помощью линейных локальных операторов число каналов входа приводится к числу каналов внутреннего представления сети hidden_channels = 32.
2. Далее выход блока P подается на вход нескольким слоям Фурье. Внутри этих слоев размерность входов и выходов сохраняется.
3. Выход пункта 2 подается в блок Q, в котором также с помощью линейных локальных операторов число каналов приводится к числу каналов выходного представления out_channels = 3.

Внутри Фурье слоя:

1. С помощью быстрого преобразования Фурье (F) для входа вычисляются коэффициенты Фурье, а также параллельно к входу применяется локальное линейное преобразование W.
2. Производится линейное преобразование R над низкочастотными коэффициентами Фурье и исключение из рассмотрения части высокочастотных коэффициентов.
3. Обратное преобразование Фурье.
4. Поэлементное суммирование двух преобразованных входов.
5. Применение нелинейной функции активации к выходу из слоя.

В блоке P находится два сверточных слоя с ядрами размера 1 × 1. В Фурье-слоях в качестве линейного преобразования над низкочастотными коэффициентами Фурье R используется преобразование из [46], а в качестве линейного локального преобразования W – свертка с ядром размера 1 × 1. Нелинейной функцией активации на выходе из Фурье-слоя является сигмоида. Всего в архитектуре данной работы восемь Фурье-слоев. В блоке Q также находятся два сверточных слоя с ядрами размерами 1 × 1 .

5.2. Общий алгоритм и данные нейросетевого метода

Для обучения нейронной сети было использовано ≈ 4.300.000 патчей размерами 64 × 64. Из них ≈ 4.000.000 резких и ≈ 300.000 размытых.

Среди четких патчей:

• ≈ 3.600.000 получены из набора данных PATH-DT-MSU WSS2 классическим методом из раздела 4;
• ≈ 400.000 получены из тренировочной выборки набора данных FocusPath. Патчи были выбраны из изображений с уровнем размытия [–2, 2].

Среди размытых патчей:

• ≈ 10.000 получены из набора данных PATH-DT-MSU WSS2 классическим методом из раздела 4;
• ≈ 180.000 получены путем применения синтетического размытия к четким патчам из набора данных PATH-DT-MSU WSS2;
• ≈ 110.000 получены из тренировочной выборки набора данных FocusPath. Патчи были выбраны из изображений с уровнем размытия либо ≥ 3, либо ≤ 3.

В качестве классификатора используется нейронная сеть на базе архитектуры FNO с добавлением к ней classification head, включающую в себя average pooling слой и два полносвязных слоя с функцией активации GeLU между ними.

6. ТЕСТИРОВАНИЕ НЕЙРОСЕТЕВОГО МЕТОДА

Тестирование метода проводилось на двух тестовых выборках из наборов данных PATH-DT-MSU WSS1 и FocusPath. Выборка из WSS1 содержит 100 изображений размерами ≈ 1500 ×1500, к которым было применено синтетическое Гауссовское размытие. Выборка из FocusPath содержит 192 размеченных изображения размерами 1024 × 1024.

В табл. 1 и 2 продемонстрировано сравнение классического и предложенного нейросетевого метода на этих двух тестовых выборках. Для определения порогов классического метода использовался набор данных PATH-DT-MSU WSS2, изображения и размытия которого похожи на присутствующие в наборе WSS1. Соответственно нейросетевой метод, обученный, в том числе, на наборе данных, собранных классическим методом из WSS2, показывает на WSS1 немного более низкие результаты по сравнению с классическим методом. Однако классический метод [37] не способен строить карты размытия для целиком размытых изображений, отмечая все области карты размытия для таких изображений как четкие. Набор данных FocusPath состоит из целиком размытых и целиком четких изображений. В то же время нейросетевой метод лишен такого недостатка и демонстрирует удовлетворительные результаты на тестовой выборке из FocusPath. На рис. 7 зеленым отмечены результаты детекции реального размытия на изображении из набора данных WSS1.

Таблица 1. Значения метрик методов на наборе данных WSS1

Таблица 2. Значения метрик методов на наборе данных FocusPath

Рис. 7. Визуализация детекций размытия

7. ЗАКЛЮЧЕНИЕ

Предложен нейросетевой метод обнаружения размытых фрагментов на гистологических изображениях. Алгоритм FNO был обучен на результатах детектирования размытых областей классическим алгоритмом и дополнительных данных. Он продемонстрировал близкие по качеству работы результаты на изображениях набора PATH-DT-MSU, на части которого обучался классический алгоритм. Однако, метод FNO позволил достаточно хорошо детектировать размытие и для набора FocusPath.

8. БЛАГОДАРНОСТИ

Исследования выполнены при финансовой поддержке Российского научного фонда в рамках научного проекта № 22-41-02002.

¹ https://sites.google.com/view/focuspathuoft

² https://imaging.cs.msu.ru/en/research/histology/path-dt-msu

About the authors

G. S. Nazarenko

Author for correspondence.
Email: s02190303@gse.cs.msu.ru
Russian Federation, Moscow

A. S. Krylov

Email: kryl@cs.msu.ru
Russian Federation, Moscow

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Graph of the absolute values of DDC coefficients sorted in ascending order for blurred and clear fragments of histological images

Download (187KB)

Indexing metadata

3. Fig. 2. Illustration of DDC coefficients for a block when they are divided into groups of correspondence to different frequencies

Download (114KB)

Indexing metadata

4. Fig. 3. Visualisation of the blur map for the histological image with blurred areas marked in red

Download (668KB)

Indexing metadata

5. Fig. 4. Graph of SI(V) threshold values

Download (80KB)

Indexing metadata

6. Fig. 5. General algorithm of the classical method

Download (110KB)

Indexing metadata

7. Fig. 6. Architecture of the neural Fourier operator

Download (165KB)

Indexing metadata

8. Fig. 7. Visualisation of blur detections

Download (173KB)

Indexing metadata