CC BY
41
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
ОБ ОДНОЙ ЗАДАЧЕ ЧИСЛЕННОГО СЕКЦИОНИРОВАНИЯ В ОФТАЛЬМОЛОГИИ
А.В. Разгулин1, Н.Г. Ирошников2, А.В. Ларичев2, С.Д. Павлов1, Т.Е. Романенко1 1Московский государственный университет имени М.В.Ломоносова, факультет вычислительной математики и кибернетики, Москва, Россия,
2Московский государственный университет имени М.В.Ломоносова, физический факультет, Москва, Россия
Аннотация
Рассматривается задача восстановления изображений (секционирования) сечений глазного дна, находящихся на различной глубине, на основе стека изображений, полученных при быстрой перефокусировке изображающей системы и представляющих собой суперпозиции истинных сечений трёхмерного объекта с размытыми изображениями соседних по глубине сечений. Для решения соответствующей задачи трёхмерной деконволюции используется неявный итерационный метод регуляризации в Фурье-плоскости. Приведённые результаты математического моделирования свидетельствуют о возможности эффективного численного секционирования в задачах офтальмологии при наличии помех и шумов различной природы.
Ключевые слова: секционирование, трёхмерная деконволюция, свёртка, глазное дно, итерационный метод регуляризации.
Цитирование: Разгулин, А.В. Об одной задаче численного секционирования в офтальмологии / А.В. Разгулин, Н.Г. Ирошников, А.В. Ларичев, С.Д. Павлов, Т.Е. Романенко // Компьютерная оптика. - 2015. - Т. 39, № 5. - С. 777-786. - DOI: 10.18287/0134-2452-2015-39-5-777-786.
Введение
Важной задачей, возникающей в офтальмологии, является неинвазивное восстановление изображения 3-мерной структуры различных отделов глаза (например, сетчатки) для последующей диагностики их заболеваний. Наиболее распространённым заболеванием сетчатки в мире на сегодня является диабетическая ретинопатия, которая служит основной причиной необратимой слепоты в развитых странах мира [1]. В клинических условиях активность ретинопатии оценивают по параметрам офтальмоскопической картины, толщины сетчатки, сосудистых изменений сетчатки, количеству микроаневризм, экссудатов и кровоизлияний. В большинстве случаев современная диагностика ретинопатий основана на классических способах офтальмоскопии или фоторегистрации изображений сетчатки в видимом спектре, а при дифференциальном диагнозе типа ретинопатии применяют флюоресцентную ангиографию [2]. Метод оптической когерентной томографии [3] также приобретает популярность в качестве альтернативной флюоресцентной ангиографии малоинвазивной методики изучения микроструктуры сетчатки. При этом трёхмерная реконструкция изучаемого участка программным методом реализована лишь в нескольких моделях оптических когерентных томографов сетчатки, а получаемые томограммы не позволяют исследователю в полной мере оценить ангиоархитектонику сетчатки. Кроме того, сосуды сетчатки, лежащие на её поверхности, имеют гиперрефлективные свойства, затрудняя оценку расположенной под ними ткани сетчатки. Также приходится учитывать и соответствующим образом компенсировать артефакты типа разрыва участков изображения, вызванные движением живого глаза [4]. В этой связи быстрая трёхмерная деконволюция сосудов с последующей трёхмерной реконструкцией прилегающей ткани может оказаться превосходящей по информативности клинической методикой оценки экссудативных форм ретинопатии.
Один из перспективных методов восстановления основан на быстрой перефокусировке изображающей системы, например с использованием методов адаптивной оптики [5], для получения стека изображений глазного дна, находящихся на различной глубине, с последующим восстановлением трёхмерной структуры известными методами. Получающиеся на этом пути изображения в каждой фокальной плоскости представляют собой суперпозицию истинного сечения трёхмерного объекта в данной фокальной плоскости с размытыми изображениями соседних по глубине сечений. Кроме того, типичной является ситуация, при которой на истинное изображение в фокальной плоскости накладываются аберрации оптической системы глаза, флуктуации фиксации, а при его регистрации приходится также учитывать искажения светочувствительных сенсоров. Т аким образом, возникает проблема устойчивого к помехам получения «очищенного» от указанных искажений стека изображений сечений глазного дна по глубине для его последующего использования в трёхмерной реконструкции.
Известным аналогом данного подхода является метод цифрового секционирования в биомикроскопии (см., например, [6]; [7], гл. 22; [8], гл. 14). В приближении полностью некогерентной изопланарной оптической системы изображение трёхмерного объекта описывается уравнением
i(x, y, z) = o(x, y, z) *h(x,y, z), (1)
где i(x,y, z) - наблюдаемое изображение (интенсивность), o(x,y, z) - искомый объект, h(x, y, z) - трёхмерная функция точечного источника (point spread function, PSF), * - знак операции трёхмерной свёртки. С помощью известной теоремы о свёртке связь между искомым объектом и его изображением в Фурье-пространстве записывается в виде поточечного произведения
I(u, v, w) = O(u, v, w) ■ H(u, v, w) (2)
Компьютерная оптика, 2015, том 39, № 5
777
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
соответствующих 3D фурье-образов I, O, H от функций i, o, h. Однако даже в идеальном случае, когда наблюдаемое изображение и функция точечного источника известны точно, однозначное нахождение фурье-образа O искомого объекта о из уравнения (2) затруднительно вследствие возможного обращения в ноль (или близости к нулю) мультипликатора H(u, v, w) в некоторых точках. Ситуация ещё более усложняется при работе с реальными данными эксперимента, когда в распоряжении имеются результаты наблюдения в конечном наборе секущих плоскостей zе {z1,z2,...,zN}, при этом как сами измерения, так и функция точечного источника, как правило, известны лишь приближённо. В этой связи особое значение имеет разработка специальных устойчивых методов решения уравнения вида (j), которые бы в наилучшей степени учитывали как специфику физической постановки задачи, связанной с характерными особенностями функции точечного источника рассматриваемой оптической системы глаза, так и требуемое соотношение производительность / точность в условиях реальных искажений, типичных для рассматриваемого класса задач офтальмологии.
1. Постановка задачи и метод решения
Конкретизируем постановку задачи, связанную с отмеченной выше спецификой. Поскольку для обработки доступны сечения трёхмерного изображения глазного дна лишь в конечном наборе плоскостей, задаваемых его сечениями при z e{zj, z2,..., zN }, причём количество таких сечений (10 - 15), как правило, на порядки меньше количества точек в каждом горизонтальном сечении, то имеет смысл перейти от симметричной по координатам (x,y, z) формы записи уравнения (j) к её дискретному по z аналогу следующего вида:
im ° i(x,У, zm ) = 0m * h0 + (3)
+ Z„^m0(x,y,zn)*hn_m(x,y), m = 1,2,.,N,
где * - знак операции двумерной свёртки, hn_m(x,y) = = h(x, y, zn - zm)Dz. Нетрудно видеть, что с помощью уравнения (3) фактически описывается ситуация, в которой искомый достаточно тонкий объект аппроксимируется с помощью стека секущих плоскостей z = zm, расположенных на одинаковом расстоянии Dz друг относительно друга, что приводит к соответствующей аппроксимации интеграла свёртки по переменной z. Уравнение
(3) связывает наблюдаемое в плоскости zm изображение im с входящими в сумму изображением om*h0 искомого объекта в фокальной плоскости и дефокусированными изображениями от соседних плоскостей zn Ф zm. Отметим, что уравнение (3) может использоваться также и вне связи с трёхмерной моделью (1), например, для моделирования трёхмерной изображающей системы со слоистой полупрозрачной структурой [9], причём толщина слоёв может быть различной.
Заметим, что при решении задачи секционирования в биомикроскопии в целях уменьшения количе-
ства ресурсоёмких вычислений обычно предполагается [6], что только ближайшие p (как правило, p = 1), лежащие снизу и сверху, соседние плоскости вносят существенный вклад в формирование изображения im. Кроме того, рассматриваются только эффекты, вызванные дефокусировкой, причём игнорируется её влияние для функции точечного источника в фокальной плоскости, т. е. предполагается, что h0 - дельтафункция Дирака. В результате получается приближённое соотношение
о = i _ У о * h .
m m *—< n n_m
n=m±1,m±2,...,m± p
Поскольку под знаком суммы искомые функции on также неизвестны, то в простейшем случае используются их аппроксимации вида
0m = C2 {im _ С1 У (in * f ) hn_j ^ (4)
у n=m±1,m±2,.,m± p J
с эмпирически подбираемыми коэффициентами c1,2, отвечающими за относительный вклад информации из фокальной и соседних плоскостей в формирование изображения, и высокочастотным фильтром f, который убирает низкочастотную информацию из областей размытия. Для повышения точности вычислений также используются итерационные варианты (3) и
(4), при этом для сокращения вычислений итерации проводятся в фурье-образах. Например, итерационный аналог (4) имеет вид [6]:
om+j) = С2 {Im _ cj
У
oik} ■ h
k = 0,1,.,
n=m±1,m±2,.,m± p
O(0) = I.
(5)
Описанные методы (4), (5) отличаются простотой реализации и хорошо зарекомендовали себя для сравнительно грубой деконволюции в задачах флуоресцентной биологической микроскопии. Однако для задачи секционирования изображений отделов глазного дна, как отмечалось выше, приходится учитывать аберрации оптической системы глаза, флуктуации фиксации, а также учитывать помехи регистрирующих устройств, что приводит к необходимости разработки специальных устойчивых методов решения уравнения вида (3) в полной постановке, что предъявляет жёсткие ограничения на соотношение производительность / точность для рассматриваемого класса задач офтальмологии.
Хорошо известно, что задачи секционирования и деконволюции трёхмерных объектов являются весьма затратными относительно вычислительных ресурсов [6]. Поэтому большую роль играет выбор экономичных алгоритмов, использующих новейшие компьютерные технологии (см., например, [10], [11]). Имея в виду широкие возможности реализации идей распараллеливания алгоритмов на современных многопроцессорных системах (CPU или GPU), рассмотрим вытекающее из (3) соотношение, которое связывает соответствующие двумерные Фурье-образы
778
Компьютерная оптика, 2015, том 39, № 5
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
Im, On,Hn функций i, o, h в плоскости спектральных переменных (u, v):
N
Im (u, v) = Z On (u, v) • Hn-m (u, v), m = 1,2,.,N. (6)
n=1
Таким образом, в каждой точке (u, v) Фурье-плоскости уравнение (6) задаёт систему линейных алгебраических уравнений (СЛАУ) с матрицей S размера N* N относительно вектора коэффициентов Фурье O = O(u, v) = (O1(u, v), O2(u, v),..., ON (u, v)) искомого объекта, которую, опуская параметры (u, v), можно записать в виде
SO = I, I = (Ii, I2,..., In ). (7)
ёв. Если апертура системы имеет характерный размер 2r, фокусное расстояние f то дефокусировка слоя, смещённого на расстояние Dz <<f равна D = Dz/f2, а соответствующая стрелка прогиба на апертуре 2r равна d = (Dz/2)*(r /f)2. При этом степень размытия определяется фазовым сдвигом ps(x,y, Dzm) и зависит от удалённости от плоскости фокусировки:
p (x, y, Dz ) = d
г s V ^ у ’ nm ' n
x2 + y2 2n r2 l
= 0,375n(x2 + y2) In - ml
d Dznmr 2
nm 2 f2
Типичная для рассматриваемых прикладных задач ситуация состоит в том, что, несмотря на относительно небольшую размерность матрицы S, точно обратить её либо нельзя из-за нетривиального ядра, либо нецелесообразно из-за вычислительных ошибок.
Рассмотрим неявный итерационный метод для устойчивого нахождения решения СЛАУ (7) (см., например, [12] гл. 2; [13] гл. 2, п. 5). В обозначениях (7) метод зависит от параметра m > 0:
O(k+1} = (e+mS*s)-1 O(k}+m(E+mS*S)~1 s* I,
—(0) (8)
k = 1,2,., O = (0,0,.,0).
В [12] гл. 2 показано, что (8) задаёт регуляризирую-щий алгоритм (РА). В связи с этим отметим, что первая итерация метода (8) в сочетании с дополнительным правилом выбора параметра m приводит к классическому методу регуляризации А.Н.Тихонова [14], что позволяет делать выводы о роли этого параметра.
Отметим, что метод (8) в терминологии ([12], гл. 5) относится к классу так называемых псевдоитерационных методов. При соответствующей организации вычислительного процесса итерации можно реализовать в виде умножения на соответствующие степени матриц перехода в Фурье-плоскости. Это позволяет существенно экономить машинное время в ситуации, когда входные данные заданы с высокой точностью и, следовательно, получение соответствующего приближённого решения потребовало бы большого числа итераций. Кроме того, поскольку процесс (8) должен быть реализован в каждой точке выбранной области Фурье-плоскости, то в целом описанный подход допускает естественное и эффективное распараллеливание вычислительного алгоритма при использовании многоядерной структуры современных CPU и GPU.
2. Функция точечного источника
Остановимся подробнее на виде используемой PSF. В наиболее простой постановке источником изображения являются N полупрозрачных слоёв, отстоящих друг от друга на расстояние Dz вдоль оптической оси. При точной фокусировке системы на один из слоёв возникает расфокусировка других сло-
где n - номер слоя, вклад которого учитывается при фокусировке на слое с номером m, Dznm - смещение слоя n относительно слоя m. Глубина объекта по координате z составляла 300 микрон. Расстояние между соседними слоями для случая стека из N = 5 секущих плоскостей тестовых изображений составляло 300 / (N- 1) = 75 микрон. Остальные параметры: 1 = 0,5 микрон, f = 20 мм. Т акже в модели учитываются аберрации второго порядка p(x, y) = 0,5n(x2 + y2), не зависящие от номера плоскости. Таким образом, PSF изображающей системы в рассматриваемом случае принимает вид
h(^y, zn - zm) =
I |2
= \F(M(x, y) • exp{i • ps (x,y, Dznm) + i • p(x, y)})| ,
где F - оператор двумерного преобразования Фурье, i 2= -1. Отметим, что наличие зрачка в рассматриваемой оптической системе учитывается с помощью функции зрачка M(x, y), которая задаётся индикаторной функцией круга радиуса r:
M(x, y) = ind (x2 + y2 < r2).
3. Результаты секционирования объекта «Секторы»
Рассмотрим тестовый объект «Секторы», который представляет собой трёхмерную структуру из пяти плоскостей, в каждой из которых расположено изображение секторов, но повёрнутое в каждой плоскости на определённый угол (см. рис. 1а). Поскольку в рамках принятого приближения и в силу симметрии объекта в равноудалённых от центральной плоскости сечениях наблюдаются качественно схожие результаты, далее приводятся изображения только для слоя m = 2. Присутствие пространственно-однородных областей, характеризуемых различными пространственными масштабами, в сочетании с резкими переходами интенсивности приводят к широкому использованию такого рода объектов для тестирования качества различных методов деконволюции, и в частном случае секционирования. Наблюдаемое изображение при фокусировке в каждую плоскость состоит из исходного изображения, искажённого аберрациями второго порядка, и дефокусированных изображений
Компьютерная оптика, 2015, том 39, № 5
779
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
из остальных плоскостей. В результате на рис. 1б хорошо видны искажения, связанные с размытием и наложением изображений, а также частичной инверсией интенсивности.
а)
Рис. 1. Исходное сечение объекта «Секторы» (а) и наблюдаемое изображение сечения (б) в секущей плоскости m = 2
Задача секционирования состоит в восстановлении исходных изображений в каждой плоскости. Рассмотрим результаты работы неявного итерационного метода без шума в зависимости от параметров метода - коэффициента m и количества итераций K. Эксперименты показали, что при m > 0,1 восстановленное изображение получается излишне сглаженным. Рассмотрим полученные результаты при m = 0,1 (см. рис. 2) на примере слоя m = 2. Видно, что увеличение числа итераций приводит к улучшению качества восстановления. Вместе с тем следует отметить эффект «насыщения», когда дальнейшие итерации метода (более 40) уже перестают оказывать существенный вклад в улучшение изображений. Аналогичная картина наблюдалась и при других значениях параметра m.
K = 5
Рис. 2. Зависимость от K (m = 0,1,m = 2)
Уменьшение параметра m в целом приводит к ухудшению качества восстановления вследствие усиления эффекта Гиббса («ringing effect») на контрастных границах объекта и вследствие этого нецелесообразно. Сравнение результатов секционирования при фиксированном числе итераций при разных m показывает, что наилучшие результаты получаются при m = 0,1 и 40 итерациях (см. рис. 3), а увеличение итераций не даёт эффекта, как и выбор других значений р,
Рассмотрим особенности секционирования, когда на вход методу подаются дефокусированные изображения с аддитивным Гауссовым шумом порядка 1 % либо с пуассоновским шумом порядка 2 %, что соответствует порядка 100 000 фотонов единичной интенсивности пикселя.
Характерно, что шум практически не заметен на фоне исходных дефокусированных изображений с соседних плоскостей (рис. 4а, б). Однако наличие шума влияет на качество секционирования.
На примере секционирования второго слоя рассмотрим влияние шума на качество работы алгоритма в зависимости от числа итераций при разных зна-
чениях параметра m. Из анализа рис. 5 видно, что, в отличие от случая без шума, с увеличением числа итераций восстановление может ухудшаться. Количество итераций должно быть согласовано с уровнем шума для получения приемлемых результатов: в рассматриваемом случае достаточно ограничиться 40 итерациями.
m=0,1 m=0,01 m=0,001
Рис. 3. Результаты секционирования слоя m = 2 при K = 40 и различных m
Рис. 4. Дефокусированные изображения слоя m = 2 объекта «Секторы» с аддитивным Гауссовым шумом 1 % (а) и пуассоновским шумом 2 % (б)
K = 5 K = 40 K = 100
Рис. 5. Зависимость от K (m = 0,1, слой m = 2 ) в случае Гауссова шума 1 % (а) и пуассоновского шума 2 % (б) Аналогичные выводы справедливы и для значения параметра m = 0,01. Однако надо иметь в виду, что при уменьшении параметра m качество секционирования ухудшается (сравн. рис. 3 и 6).
m = 0,1 m = 0,01 m = 0,001
Рис. 6. Зависимость от /и (K = 40, слой m = 2 ) в случае Гауссова шума 1 % (а) и пуассоновского шума 2 % (б)
780
Компьютерная оптика, 2015, том 39, № 5
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
Таким образом, возможности выбора параметра m для подстройки работы метода к специфике восстанавливаемых изображений и искажений ограничиваются при малых m эффектом Гиббса, приводящим к ложному оконтуриванию, однако это может быть важно в задаче выделения контуров изображений при априорной информации об однородности изображения между контурами.
4. Секционирование изображений глазного дна
Для работы с изображениями глазного дна рассматривался стек из N = 3 секущих плоскостей объекта (см. рис. 1а) с глубиной 200 микрон по координате z. Заметим, что, в отличие от объекта «Секторы», в рассматриваемом случае существенно больше полутонов и актуальных для диагностики патологий высокочастотных неоднородностей, отражающих
структуру сетчатки, а также кровеносных сосудов.
m = 1 m = 2 m = 3
Рис. 7. Исходные изображения глазного дна (а) и стек наблюдаемых изображений без шума (б)
Соответствующие функции дефокуса между плоскостями и аберрации второго порядка имели вид:
Ps (х, У, Dznm) = 0,5 п( х2 + y2) \п - т\,
Р( х, У) = 0,5п( х2 + у2).
В каждой плоскости фокусировки вместо исходных изображений наблюдаются результаты их конво-люции с изображениями соседних плоскостей в сочетании с общими аберрациями второго порядка. Соответствующие изображения, приведённые на рис. 1б, характеризуются сильным размытием и исчезновением структуры сетчатки.
Для количественной оценки качества восстановления изображений в каждой секущей плоскости фокусировки используется векторный частотный критерий (ВЧК), позволяющий контролировать особенности восстановления в спектральной плоскости и задаваемый функцией
A(r)
1
m(r)
z
(u,v):r2 <u2 +v2 <(r+1)2
F (E (D(o))) F (E (D(o)))
(u, v),
где о - восстановленное изображение, сдвинутое в неотрицательную область, o - исходное изображение, D (X) = X/D(X), E (X) = X/E (X), DQ - дисперсия, E(-) - среднее, m(r) - число точек (u, v), удовлетво-
ряющих соотношению r2 < u2+v2 < (r+1)2. Таким образом, для вычисления ВЧК сравниваемые изображения нормируются на свои средние по апертуре значения и на дисперсию. Далее вычисляются Фурье-образы и находится отношение их модулей. Полученное распределение усредняется по концентрическим кольцам с центром в нуле шириной в несколько пикселей (в наших экспериментах - 2). Составленная из полученных значений результирующая кривая даёт зависимость нормированного отношения спектров от модуля спектральной переменной в относительных единицах и позволяет охарактеризовать результаты восстановления в зависимости от интересующей частотной области.
0 128 256 0 128 256
K = 10 K = 40
Рис. 8. Результаты секционирования слоя изображений глазного дна и ВЧК в зависимости от K (m = 0,1, слой m = 2 без шума)
Как и в случае объекта «Секторы», при обработке реальных изображений глазного дна без шумов при увеличении количества итераций наблюдается эффект «насыщения»: дополнительные (после 40) итерации уже визуально не заметны (см. рис. 8). Однако применение ВЧК позволило выяснить, что на больших итерациях метод отрабатывает в основном высокочастотные составляющие изображений, приближая ВЧК к равномерному распределению.
При уменьшении параметра m до значения 0,01 ухудшается качество восстановления высокочастотных составляющих изображения. Увеличение количества итераций (не менее 100) частично исправляет ситуацию, однако неравномерность векторного критерия выше по сравнению с m = 0,1. В связи с этим дальнейшее уменьшение m нецелесообразно вследствие ухудшения равномерности ВЧК и появления эффекта Гиббса (см. рис. 9).
m=0,1 m=0,01
Рис. 9. Зависимость от m (K =
m = 0,001
100, m = 2, без шума)
Компьютерная оптика, 2015, том 39, № 5
781
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
Таким образом, для восстановления изображений глазного дна без шума диапазон рабочих значений параметра m составляет от 0,01 до 0,1 с соответствующим выбором количества итераций K в зависимости от требуемого распределения ВЧК.
Приведём результаты восстановления изображений глазного дна с аддитивным Г ауссовым шумом порядка 1 %. В качестве исходных данных имеется стек наблюдаемых зашумлённых дефокусированных изображений, приведённых на рис. 10. Заметим, что в реальной ситуации «очищенное» от соседних слоёв изображение также содержит шумы измерения. Поэтому для контроля качества секционирования используется модифицированный ВЧК, при вычислении которого в моделируемой нами ситуации к искомому изображению в слое добавляется шум такой же природы и величины (но, возможно, другой реализации), что и в наблюдаемом дефокусированном изображении.
m = 1 m = 2 m = 3
Рис. 10. Стек дефокусированных изображений глазного дна с аддитивным Гауссовым шумом 1 %
В отличие от тестового объекта «Секторы», содержащего контрастные границы и участки однородного распределения интенсивности, в случае полутоновых изображений глазного дна с шумами для проявления тонких структур (сосуды) и волнистой топологии фона глазного дна, несущей важную офтальмологическую информацию о рецепторах, значение параметра m=0,1 недостаточно для качественного восстановления, причём увеличение числа итераций лишь ухудшает картину. Улучшение качества восстановления наблюдается при уменьшении параметра до значения m = 0,01. При этом метод успешно производит «секционирование» смешанных изображений, однако шум остаётся заметным, несмотря на существенное выравнивание векторного критерия (рис. 11).
Я Я
;’7 А 5,3 -7Г
o,9k~^J \
0 128 256 0 128 256
K = 10 K = 40
Рис. 11. Секционирование изображений глазного дна (m = 0,01, слой m = 2) для Гауссова шума 1 %
Уменьшение значения параметра до m = 0,001 (см. рис. 12) даёт наиболее адекватный результат как с
точки зрения ВЧК, так и с точки зрения качественных особенностей восстановления. Для 100 итераций на среднем восстановленном изображении хорошо проявились как характерная фоновая текстура глазного дна, так и кровеносные сосуды.
0 128 256 0 128 256
K = 40 K = 100
Рис. 12. Секционирование изображений глазного дна (m = 0,001, слой m = 2) для Гауссова шума 1 % Рассмотрим влияние пуассоновского шума около 2 % (100 000 фотонов на единицу интенсивности пикселя) на наблюдаемые изображения (см. рис. 13).
m = 1
m = 2
m = 3
Рис. 13. Стек дефокусированных изображений глазного дна с пуассоновским шумом 2 %
Приведём результаты работы итерационного метода в зависимости от числа итераций при различных значениях параметра m = 0,1, 0,01, 0,001.
Рис. 14. Секционирование изображений глазного дна (m = 0,1, слой m = 2) для пуассоновского шума 2 %
Из рис. 14 видно, что для достижения приемлемой равномерности распределения векторного критерия при m = 0,1 на низких и средних частотах достаточно ограничиться небольшим количеством итераций (не более 10). Увеличение числа итераций приводит к росту искажений в высокочастотной области и соответствующей деградации изображения.
Анализ рис. 15 показывает, что при уменьшении параметра m до значения 0,01, соответствующего предыдущему случаю, результата удаётся добиться ценой увеличения числа итераций до 20 - 40. При
782
Компьютерная оптика, 2015, том 39, № 5
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
этом наблюдается равномерность векторного критерия в диапазоне низких и средних частот. Вместе с тем характерная фоновая текстура глазного дна всё ещё остаётся зашумлённой.
О 128 256
K = 5
Рис. 15. Секционирование стека изображений глазного дна (Ц = 0,01, слой т =2) для пуассоновского шума 2 %
1
1,0 1,0 у 1,0 S A
0,5 0,5 V-Л 0,5 \J\
0 128 256 0 128 256 0 128 256
K = 10 K = 40 K = 100
Рис. 16. Секционирование изображений глазного дна (т = 0,001, слой m = 2) для пуассоновского шума 2 % Уменьшением параметра m до значения 0,001 с одновременным увеличением числа итераций до 100 (см. рис. 16) удаётся добиться хорошей равномерности векторного критерия вплоть до самых высоких частот с проявлением характерной фоновой текстуры глазного дна.
5. О возможности распараллеливания метода
Важной особенностью метода является возможность эффективного распараллеливания вычислительного алгоритма при использовании многоядерной структуры современных CPU и GPU. Об эффективности распараллеливания на многоядерных архитектурах в смысле сильной масштабируемости можно судить по графику ускорения программной реализации алгоритма в зависимости от числа используемых рабочих процессов в вычислительном пуле среды MATLAB (см. рис. 17). Расчёты проводились на персональном компьютере с процессором Intel Core i7-4790K, имеющем 4 процессорных ядра с поддержкой технологии Intel Hyper-Threading, для задачи с размерностью изображений 512^512. Как видно из графика, даже использование стандартных средств распараллеливания, имеющихся в среде MATLAB, позволяет сократить фактическое время расчёта в 4 раза.
Заключение
В работе рассмотрена задача восстановления изображений (секционирования) сечений глазного дна на
основе трёхмерной деконволюции неявным итерационным методом регуляризации в Фурье-плоскости. Полученные результаты численного моделирования свидетельствуют о возможности эффективного секционирования в условиях помех и шумов различной природы. Для типичного в задачах офтальмологии уровня шумов метод достаточно хорошо адаптируется как к пуассоновскому шуму, так и в несколько меньшей степени к Гауссову шуму. Наличие двух управляющих параметров (число m и количество итераций K) позволяет гибко настраивать метод в зависимости от преференций восстановления: быстрое первичное секционирование с наименьшим числом итераций для грубой визуализации слоёв и выделения кровеносных сосудов при относительно больших значениях m либо углубленное восстановление с выявлением фоновой текстуры глазного дна за счёт дополнительных итераций при малых m.
Speed-up
4.0
3.5
3.0
2.5
2.0
1.5
1,0
У
J
у
и
„ Numb гг
roj processes
Рис. 17. График ускорения в зависимости от числа рабочих процессов в пуле MATLAB
Для полноты картины отметим, что после секционирования «очищенные» изображения сечений глазного дна могут пройти дополнительную специализированную цифровую постобработку, нацеленную, например, на визуализацию сосудов, выделение особенностей сетчатки, фильтрацию шумов и др. Аналитический обзор существующих возможностей и потребностей в цифровом анализе изображений глазного дна проведён в [15]. Также отметим примыкающие к рассмотренным вопросам методы постобработки плоских офтальмологических изображений как с использованием некоторой информации о PSF (см., например, [16]), так и без неё ([17], [18]).
Благодарности
Работа выполнена при поддержке гранта РФФИ № 13-02-12219.
Литература
1. Diabetes in America. - 2nd ed. - Bethesda, MD: National Institutes of Health, National Institute of Diabetes and Digestive and Kidney Diseases, 1995. - 782 p.
2. Bursell, S.E. Stereo nonmydriatic digital-video color retinal imaging compared with Early Treatment Diabetic Retinopathy Study seven standard field 35-mm stereo color photos for determining level of diabetic retinopathy / S.E. Bursell, J.D. Cavallerano, A.A. Cavallerano, A.C. Clermont, D. Birkmire-Peters, L.P. Aiello, L.M. Aiello // Ophthalmology. - 2001. - Vol. 108(3). - P. 572-585. - ISSN 01616420. - D0I:10.1016/S0161-6420(00)00604-7.
3. Buabbud, J. Optical coherence tomography imaging for diabetic retinopathy and macular edema / J. Buabbud,
Компьютерная оптика, 2015, том 39, № 5
783
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
M. Al-latayfeh, J. Sun // Current Diabetes Reports. -2010. - Vol. 10(4). - P. 264-269. - ISSN 1534-4827. -DOI: 10.1007/s11892-010-0129-z.
4. Braaf, B. Real-time eye motion correction in phase-resolved OCT angiography with tracking SLO / B. Braaf,
K.V. Vienola, C.K. Sheehy, Q. Yang, K.A. Vermeer, P. Tiruveedhula, D.W. Arathorn, A. Roorda, J.F. de Boer // Biomedical Optics Express. - 2013. - Vol. 4(1). - P. 5165. - ISSN 2156-7085. - DOI: 10.1364/BOE.4.000051.
5. Larichev, A.V. Adaptive system for eye-fundus imaging / A.V. Larichev, P.V. Ivanov, N.G. Iroshnikov, V.I. Shmal-hauzen, L.J. Otten // Quantum Electronics. - 2002. -Vol. 32. - P. 902-908. - ISSN 1063-7818.
6. Agard, D. Optical sectioning microscopy: cellular architecture in three dimensions / D. Agard // Annual Review of Biophysics Bioengineering. - 1984. - Vol. 13. - P. 191219. - DOI: 10.1146/annurev. bb.13.060184.001203.
7. Castleman, K. Digital Image Processing / K. Castleman. -Pearson Education, 2007. - 686 p.
8. Wu, Q. Microscope image processing. / Q. Wu, F. Merchant,
K. Castleman. - Academic Press, 2008. - 576 p.
9. Матвиенко, А.Н. Метод обработки изображений полупрозрачных слоистых сред / А.Н. Матвиенко, Т.Н. Новикова // Вестник Московского университета. Серия 15. Вычислительная математика и кибернетика. - 1988. -№ 4. - С. 33-37. - ISSN: 0137-0782.
10. Bruce, M. Real-time GPU-based 3D deconvolution /
M. Bruce, M. Butte // Optics Express. - 2013. - Vol. 21(4). - P. 4766-4773. - ISSN: 1094-4087. - DOI:
10.1364/OE.21.004766.
11. Zanella, R. Towards real-time image deconvolution: application to confocal and STED microscopy / R. Zanella, G. Zanghirati, R. Cavicchioli, L. Zanni , P. Boccacci , M. Bertero, G. Vicidomini // Scientific Reports. - 2013. -Vol. 3(2523). - DOI: 10.1038/srep02523.
12. Итеративные методы решения некорректных задач /
A. Б. Бакушинский, А.В. Гончарский. - М.: Наука, 1989. - 128 с.
13. Введение в теорию обратных задач / А.М. Денисов. - М.: Издательство Московского университета, 1994. - 208 с.
14. Методы решения некорректных задач / А.Н. Тихонов,
B. Я. Арсенин. - М.: Наука, 1986. - 288 с.
15. Ильясова, Н.Ю. Методы цифрового анализа сосудистой системы человека. Обзор литературы / Н.Ю. Ильясова // Компьютерная оптика. - 2013. - Т. 37, № 4. - С. 511-535.
16. Iroshnikov, N. A modified bispectral image reconstruction method in ophthalmology / N. Iroshnikov, A. Larichev, A. Razgulin, A. Starostin // Computational Mathematics and Modeling. - 2015. - Vol. 26(4). - P. 534-545. - DOI: 10.1007/s10598-015-9290-1.
17. Крылов, А.С. Компьютерный анализ изображений
глазного дна / А.С. Крылов, А.В. Насонов, А.С. Семашко, А.А. Черноморец, В.В. Сергеев, В.С. Акопян,
А.С. Родин, Н.С. Семенова // СПб.: Труды VIII Российско-Баварской конференции по биомедицинской инженерии, 2012. - С. 129-133.
18. Насонова, А.А. Выделение сосудов на изображениях глазного дна и его оценка качества / А. А. Насонова, А.С. Крылов // Биотехносфера. - 2014. - Т. 2. - С. 24-25.
References
[1] Diabetes in America. 2nd ed. Bethesda, MD: National Institutes of Health, National Institute of Diabetes and Digestive and Kidney Diseases; 1995.
[2] Bursell SE, Cavallerano JD, Cavallerano AA, Clermont AC, Birkmire-Peters D, Aiello LP, Aiello LM. Stereo nonmydriatic digital-video color retinal imaging compared with Early Treatment Diabetic Retinopathy Study seven standard field 35-mm stereo color photos for determining level of diabetic retinopathy. Ophthalmology 2001; 108(3): 572-5. DOI:10.1016/S0161 -6420(00)00604-7.
[3] Buabbud J, Al-latayfeh M, Sun J. Optical coherence tomography imaging for diabetic retinopathy and macular edema. Current Diabetes Reports 2010; 10(4): 264-9. DOI: 10.1007/s11892-010-0129-z.
[4] Braaf B, Vienola KV, Sheehy CK, Yang Q, Vermeer KA, Tiruveedhula P, Arathorn DW, Roorda A, de Boer JF. Real-time eye motion correction in phase-resolved OCT angiography with tracking SLO. Biomedical Optics Express 2013; 4(1): 51-65. DOI: 10.1364/BOE.4.000051.
[5] Larichev AV, Ivanov PV, Iroshnikov NG, Shmalhauzen VI, Otten LJ. Adaptive system for eye-fundus imaging. Quantum Electron [In Russian]. 2002; 32: 902-8.
[6] Agard D. Optical sectioning microscopy: cellular architecture in three dimensions. Ann. Rev. Biophys. Bioeng. 1984; 13: 191-219. DOI: 10.1146/annurev.bb.13.060184.001203.
[7] Castleman K. Digital Image Processing. Pearson Education; 2007.
[8] Wu Q, Merchant F, Castleman K. Microscope image processing. Academic Press; 2008.
[9] Matvienko AN, Novikova TN. The translucent stratified media image processing method [In Russian]. Vestnik Moskovskogo Universiteta. Seriya 15. Vychislitel'naya Matematika i Kibernetika 1988; 4: 33-7.
[10] Bruce M, Butte M. Real-time GPU-based 3D deconvolution. Optics Express 2013; 21(4): 4766-73. DOI: 10.1364/OE.21.004766.
[11] Zanella R, Zanghirati G, Cavicchioli R, Zanni L, Boc-
cacci P, Bertero M, Vicidomini G. Towards real-time image deconvolution: application to confocal and STED microscopy. Scientific Reports 2013; 3(2523).
DOI:10.1038/srep02523.
[12] Bakushinskiy AB, Goncharskiy AV. Iterative methods of solving ill-posed problems [In Russian]. Moscow: “Nau-ka” Publisher; 1989.
[13] Denisov AM. Introduction to the inverse problems theory [In Russian]. Moscow: Publisher of Moscow University; 1994.
[14] Tikhonov AN, Arsenin VY. Ill-posed problems solving methods [in Russian]. Moscow: “Nauka” Publisher; 1986.
[15] Ilyasova NYu. Methods for digital analysis of human vascular system. Literature review. Computer Optics 2013; 37(4) 511-35.
[16] Iroshnikov N, Larichev A, Razgulin A, Starostin A. A modified bispectral image reconstruction method in ophthalmology. Computational Mathematics and Modeling 2015; 26(4): 534-45. DOI: 10.1007/s10598-015-9290-1.
[17] Krylov AS, Nasonov AV, Semashko AS, Chernomorets AA, Sergeev VV,. Akopyan VS, Rodin AS, Semenova NS. Computer analysis of fundus images [In Russian]. SPb: Proceedings of the VIII Russian-Bavarian conference on biomedical engineering; 2012: 129-33.
[18] Nasonova AA, Krylov AS. Allocation of blood vessels on the fundus image and its quality evaluation [In Russian]. Biotekhnosphera 2014; 2: 24-5.
784
Компьютерная оптика, 2015, том 39, № 5
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
ON A PROBLEM OF NUMERICAL SECTIONING IN OPHTHALMOLOGY
A.V. Razgulin, N.G. Iroshnikov2, A.V. Larichev2, S.D. Pavlov1, T.E. Romanenko1 1Lomonosov Moscow State University, Computational Mathematics and Cybernetics Faculty, Moscow, Russia, 2Lomonosov Moscow State University, Physics Faculty, Moscow, Russia
Abstract
We consider a problem of eye-fundus image reconstruction from different-depth fundus sections based on an image stack, with the images obtained via imaging system's fast refocusing as superposition of 3D object sections and blurred images of adjacent-depth sections. An implicit iterative regularization method in the Fourier plane is used for 3D deconvolution. The results of mathematical modeling have demonstrated that the numerical sectioning shows promise when processing ophthalmological images distorted by various types of noise.
Keywords: sectioning, 3D deconvolution, convolution, fundus, iterative regularization method.
Citation: Razgulin AV, Iroshnikov NG, Larichev AV, Pavlov SD, Romanenko TE. On a problem of numerical sectioning in ophthalmology. Computer Optics 2015; 39(5): 777-86. DOI: 10.18287/0134-2452-2015-39-5- 777-786.
Acknowledgements: The work was funded by the Russian Foundation of Basic Research grant No 13-02-12219.
Сведения об авторах
Разгулин Александр Витальевич, 1963 года рождения, в 1985 году окончил Московский государственный университет имени М.В. Ломоносова по специальности «Прикладная математика», кандидат физикоматематических наук с 1988 года, доктор физико-математических наук с 2007 г., работает профессором на факультете вычислительной математики и кибернетики МГУ имени М.В. Ломоносова. Область научных интересов: функционально-дифференциальные уравнения, оптимальное управление и аппроксимация задач адаптивной и нелинейной оптики, физиологическая оптика ,математическое моделирование.
E-mail: razgulin@cs. msu. ru .
Alexander Vitalievich Razgulin (b. 1963) graduated from Lomonosov Moscow State University in 1985, majoring in Applied Mathematics, got Ph.D. in Physics and Mathematics in 1988, Doctor in Physics and Mathematics in 2007. Currently he works as the professor at Lomonosov Moscow State University, Computational Mathematics and Cybernetics faculty. Research interests are functional-differential equations, optimal control and approximation to the problems of adaptive and nonlinear optics, physiological optics, mathematical modeling.
Ирошников Никита Г еоргиевич, 1965 года рождения, в 1988 году окончил Московский государственный университет имени М.В. Ломоносова по специальности «Физика», работает доцентом в МГУ имени М.В. Ломоносова. Область научных интересов: нелинейная оптика, медицинская физика, программирование.
E-mail: nikita@optics. ru .
Nikita Georgievich Iroshnikov (b. 1965) graduated from Lomonosov Moscow State University in 1988, majoring in Physics. Currently he works as the associate professor at Lomonosov Moscow State University. Research interests are nonlinear optics, medical physics, programming.
Ларичев Андрей Викторович, 1965 года рождения, в 1989 году окончил Московский государственный университет имени М.В. Ломоносова по специальности «Физика», кандидат физико-математических наук с 1995 года, работает доцентом в МГУ имени М.В. Ломоносова. Область научных интересов: адаптивная оптика, физиологическая оптика, нелинейная оптика, оптическая обработка информации, нелинейная динамика, цифровая обработка изображений.
E-mail: larichev@,optics.ru.
Andrey Victorovich Larichev (b. 1965) graduated from Lomonosov Moscow State University in 1989, majoring in Physics, got Ph.D. in Physics and Mathematics in 1995. Currently he works as the associate professor at Lomonosov Moscow State University. Research interests are adaptive optics, physiological optics, nonlinear optics, optical information processing, nonlinear dynamics and digital image processing.
Павлов Станислав Дмитриевич, 1992 года рождения, в 2014 году окончил Московский государственный университет имени М.В. Ломоносова по специальности «Прикладная математика и информатика», проходит обучение в аспирантуре МГУ имени М. В. Ломоносова. Область научных интересов: нелинейные функциональнодифференциальные уравнения диффузии, самоорганизация в моделях нелинейной оптики, программирование.
E-mail: spavlov@cs.msu.ru .
Компьютерная оптика, 2015, том 39, № 5
785
Об одной задаче численного секционирования...
Разгулин А.В., Ирошников Н.Г., Ларичев А.В., Павлов С.Д., Романенко Т.Е.
Stanislav Dmitrievich Pavlov (b. 1992) graduated from Lomonosov Moscow State University in 2014, majoring in Applied Mathematics and Informatics. He is a PhD student at Lomonosov Moscow State University. Research interests: nonlinear diffusion functional differential equations, self-organization in nonlinear optical models, programming.
Романенко Татьяна Евгеньевна, 1988 года рождения, в 2010 году окончила Московский государственный университет имени М.В. Ломоносова по специальности «Прикладная математика и информатика», работает ассистентом в МГУ имени М.В. Ломоносова. Область научных интересов: параболические функциональнодифференциальные уравнения с запаздыванием, самокомпенсация искажений в нелинейных оптических системах, обработка графических изображений, программирование.
E-mail: romanenko@cs.msu.ru .
Tatiana Evgenievna Romanenko (b. 1988) graduated from Lomonosov Moscow State University in 2010, majoring in Applied Mathematics and Informatics. Currently she works as the assistant professor at Lomonosov Moscow State University. Research interests are parabolic functional differential equations with time delay, distortion self-compensation in nonlinear optical systems, computer graphics processing, and programming.
Поступила в редакцию 15 октября 2015 г. Окончательный вариант - 1 декабря 2015 г.
786
Компьютерная оптика, 2015, том 39, № 5
