Моделирование оптических методов биомедицинской диагностики
Д.
А. ПетровБиблиографическое описание статьи для цитирования:
Петров
Д.
А. Моделирование оптических методов биомедицинской диагностики // Научно-методический электронный журнал «Концепт». –
2016. – Т. 11. – С.
2851–2855. – URL:
http://e-koncept.ru/2016/86601.htm.
Аннотация. Представлен способ Монте-Карло моделирования фотонного транспорта в биологических тканях, основанный на выполнении высокоскоростных вычислений с помощью графических процессоров, который позволяет исследовать особенности современных оптических методов биомедицинской диагностики.
Текст статьи
Петров Денис Алексеевич,Студент 4 курса специальности «Биотехнические системы и технологии» ФГБОУ ВПО «Тамбовский государственный технический университет», г. Тамбовden794@mail.ru
Моделирование оптических методов биомедицинской диагностики
Аннотация.Представлен способ МонтеКарломоделирования фотонного транспорта в биологических тканях, основанный на выполнениивысокоскоростных вычислений с помощью графических процессоров, который позволяет исследовать особенностисовременных оптических методов биомедицинской диагностики.Ключевые слова:метод МонтеКарло, графические процессоры, фотонный транспорт.
Методы медицинской визуализации, основанные на взаимодействии оптического излучения с исследуемой биологической тканью, которые включаютоптическую когерентную томографию, диффузионную оптическую томографию, микроскопию, флуоресцентную ангиографию и т.д., на данный момент являются одними из наиболее используемыхметодов медицинской диагностики. Оптические методы позволяют получить внутреннюю структуру ткани неинвазивно с высоким пространственным разрешением, что при использовании излучения небольшой мощности (~0.1 мВ) позволяет избежать разрушения структуры ткани и полностью минимизировать возможный вред для пациента.Отличительной особенностью большинства методов визуализации, основанных на взаимодействии света с исследуемой биологической тканью, является высокое пространственное разрешение получаемых изображений внутренней структуры, которое, к примеру, для конфокальноймикроскопии составляет 1 микрометр, а для оптической когерентной томографии десятки микрометров. Но, по причине высокого показателя рассеивания излучения для данных методик также характерно крайне низкое пространственное разрешение. Для усовершенствованиясовременных методов оптической диагностики необходимо полное понимание закономерностей, лежащих в основе фотонного транспорта в биологических объектах. Изучение данных закономерностей возможно путем проведения компьютерного моделирования.Как правило, моделирование распространения света в сильнорассеивающих средах проводится с помощью метода, основанного на принципах МонтеКарло. Данный метод реализован с помощью специализированного программного обеспечения MCML, которое находится в открытом доступе с середины девяностых годов прошлого столетия [1]. ДанноеПОпозволяет производить моделирование фотонного транспорта только в полубесконечной слоистой среде, при чем детектирование энергии излученияпроизводится только в точках взаимодействия фотона со средой без учета временных осцилляций оптических характеристик объекта(стационарный режим). Реальные биологические объекты, как правило, имеют гораздо более сложную пространственную структуру, изменяющуюся во времени, поэтому для их исследования необходим способ моделирования, учитывающий эти особенности. При этом стоит также отметить статистический характер метода МонтеКарло, что делает моделирование довольно медленным, в то время как программа MCML, наиболее частоиспользуемаяна данный момент, работает в однопоточном режиме, то есть позволяет использовать только одно ядро центрального процессора(ЦП), в то время как современные ЦП имеют 2128 ядер.Поэтому на данный момент данное ПО является устаревшим, как с точки зрения современных технических систем, так и по причине объективного отсутствия возможностей для исследования реальных биологических объектов и моделирование различных методик медицинской диагностики.Для увеличения скорости проведения моделирования предлагается использовать архитектуру параллельных вычислений CUDA, основанную на высокоскоростных вычислениях с помощью графических процессоров. Первоначально графические процессоры использовались только для вычисления и преобразования компьютерной графики, но с появлением программируемых шейдерных блоков появилась возможность проводить вычисления и в приложениях общего назначения[2]. Моделирование фотонного транспорта предполагает многократные вычисления распространения большого числа фотонов(10100 миллионов), что делает использование графических процессоров наиболее оправданным.Распространение фотонного транспорта можно представить в виде нескольких этапов. Сначала производится запуск фотона, определение дистанции свободного пробега между точками взаимодействия со средой, проверка на пересечение границысреды, внутри которой оптические свойства однородны, рассеивание, которое определяется коэффициентом рассеивания, специфичным для каждой ткани, и поглощения энергии фотонного пучка. Затем эти шаги снова повторяются до тех пора, пока фотон находится внутриисследуемого объекта. При моделировании различных методов медицинской диагностики данные этапы являются аналогичными, а меняется только процесс детектирования фотонов, выходящих из биологической ткани.Для изучения фотонного транспорта в объектах со сложной пространственной структурой предлагается использовать трехмерный массив, каждому элементу которогоприсваивается свой индекс, характеризующий его оптические свойства (коэффициент рассеяния, коэффициент поглощения, показатель преломления и коэффициент анизотропии). Размер и общее количество элементов массиваопределяется исходя из размера объекта исследования и пространственного разрешения методики, то есть для конфокальной микроскопии длина ребра элемента массивабудет небольшой (примерно 1 микрометр), адля методов, основанных на исследовании рассеивания и поглощения объектов, таких как диффузионная томография, размер будет существенно больше (0.10.01 миллиметра). Создание трехмерной среды предлагается проводить путем присваивания группе элементовпо ихиндексамопределенного числа, характеризующего какойлибо слой объекта, оптические свойства которого не изменяются в пространстве, либо путем задания неоднородностей внутри объекта в аналитическом виде (сфера, куб, цилиндр), соответственно те элементы, которые попадают внутрь неоднородности, будут иметь одинаковые оптические свойства. Первый случай используется для создания разреженных в пространстве неоднородностей, к которым можно отнести раковые метастазы, то есть объекты, задание границ которых в аналитическом виде невозможно. Аналитическое задание границ может использоваться при моделировании таких объектов как кровеносные сосуды (сфера, эллипсоид), сляб, фантомы биологических тканей и т.д.
В ПО MCMLрасчет дистанции свободного пробега производится с помощью показателя рассеивания и поглощения среды , что соответствует дискретному(вероятностному) подходу к процессу уменьшения энергии фотонного пучка. При проведении времяразрешенного моделирование гораздо удобнее рассматривать процесс поглощения непрерывным, согласно закону БераБугераЛамбера:
Где W–энергия текущего пакета фотонов, а s–дистанция пролета фотона. Дистанция пробега фотона же рассчитываетсятолько согласно коэффициенту рассеивания:
где
случайное число равномерно распределенное от 0 до 1. Подобный подход также позволяет увеличить скорость моделирования за счет увеличения длины свободного пробега фотона, что связано с исключением показателя поглощения из знаменателя формулы расчеты пробега.По мере распространения частицы внутри объекта фотон будет перемещать из одногоэлемента трехмерного массива в другой. В случае если показатель преломления вэлементах различен, то фотон может отразиться от границы согласно закону Френеля:
Объекты также могут изменять свою структуру во времени, что может быть связано, например, с пульсациями кровеносных сосудов. Для моделирования таких объектов необходимо отслеживать время нахождения фотона в исследуемой среде, котороевычисляется по следующей формуле:
где с –скорость света в вакууме, n –показатель преломления среды, а s–дистанция передвижения.Для каждого слоя/элемента массиваобъекта при этом задают зависимость оптических свойств от времени.
Угол рассеивания получают, используя фазовую функцию ХеньиГринштейна:
где g–коэффициент анизотропии элемента массивасреды, в котором на данный момент находится объект.Полярный угол равномерно распределяется от 0 до . Направления движения фотона при рассеивании определяется по направляющим векторам:
где , и
текущие значения направляющих косинусов, а , и
новыезначения.Проведение вычислений на CUDAпозволяет использовать библиотеки ускоренных функций для деления (__fdividef) и для расчета значений косинусов и синусов(_sinf, _cosf), вычисление которых, как правило, требует большого количества времени[3]. При этом для наиболее высокой скорости вычислений необходимо использовать числа одинарной точности float.Общее количество блоков и нитей (потоков) при вычислении на графических процессорах зависит от характеристик конкретной видеокарты. Каждая нить выполняет вычисление распространения своего фотона, обращаясь к глобальной динамической памяти только при необходимости записи значения энергии фотонной плотности и декрементирования общего числа фотонов. Для того, чтобы избежать конкурирования потоков за запись в однуячейку памяти данныешаги выполняются с помощью библиотеки атомарных функций Cuda, а именно функции atomicAdd. Данная модель реализована с помощью языка программирования CUDAC. Скорость моделирования при использовании полученнойпрограммы в общем случаев 4050 раз превышает скорость моделирования с помощью аналогичного программного кода, полученного с помощью языка C#. По сравнению с программным обеспечением, реализующем стационарную модель фотонного транспорта MCMLскорость моделирования увеличивается в 100120 раз, что позволяет получить в некоторых случаях практически мгновенные результаты моделирования.Процесс моделирования различных методов биомедицинской диагностики основан на использовании описанной выше модели, но при этом имеет свои особенностив зависимости от конкретной методики. Рассмотрим моделирование некоторых из них.Оптическая когерентна томография (ОКТ) –методика визуализации основанная на принципах низкокогерентной интерферометрии. Процедура ОКТ представляет собой последовательное сканирование объекта в поперечном направлении.В случае МонтеКарло производится последовательный запуск группы фотонов в различных точках сред и детектирование отраженных энергии отраженных частиц. После формирования одного Аскана(линии) изображения точка запуска смещается, и процесс повторяется. Интенсивность пикселя изображения в общем случае определяется с помощью следующей формулы:
где W–энергия кванта, 2zи –оптические пути в опорном и объектном плече интерферометра соответственно, а
длина когерентности источника [4].Результат моделирования ОКТ методики исследования среды, содержащей подкожный сосуд, представлен на рисунке 1.
Рис. 1. Результат моделирования подкожного сосудаСтоит отметить, что результат в целом соответствует изображениям, получаемым в экспериментальных ОКТ исследованиях. Сосудистая стенка обладаетсущественно меньшим показателем рассеивания, чем кровь (различие более чем в 200300 раз), поэтому на изображении эта область выглядит затененной, а просвет, соответственно, отличается гораздо более высокой интенсивностью, а максимальной интенсивности соответствует верхний слой кожи stratumcorneum, что связано с высоким различием показателя между ним и воздухом.
Моделирование диффузионной оптической томографии(ДОТ)отличается тем, что данная методика основана на регистрации многократно рассеиваемых фотонов. Результат моделирование ДОТ на основе исследования цилиндрического фантома представлен на рисунке 2. Оптические свойства объекта: коэффициент рассеяния –50 см^1, коэффициент поглощения –0.2 см^1, коэффициент анизотропии –0.9., радиус цилиндра –2 см.
Рис. 2. Плотность распределения энергии пучка в цилиндре, на момент 0.59 наносекунд (справа) и 0.15 наносекунд (слева).
Представленная модель фотонного транспорта позволяет выполнять высокоскоростное моделирование различных медицинских методов визуализации, что имеет большое значение для прикладных исследований. Изменяя оптические свойства среды, детектора, светового пучка при моделировании получают различные результаты и выбирают наиболее оптимальный , что может затем использоваться при реальных исследованиях, например при определении необходимости оптического просветления, настройки используемого лазера, выбора детекторов и т.д.
Ссылки на источники1. L.H. Wang, S. L. Jacques, and L.Q. Zheng MCML—Monte Carlo modeling of light transport in multilayered tissues // Computer Methods and Programs in Biomedicine 47 (2), 1995. Pp. 131–146.2. S. Mittal and J. Vetter, A Survey of CPUGPU Heterogeneous Computing Techniques, // ACM Computing Surveys, 2015.3. Джейсон Сандерс, Эдвард Кэндрот «Технология CUDA в примерах. Введение в программирование графических процессоров» // ДМК Пресс, 2011 г. –232 стр.4. M. Kiriin, I. Meginski, V. Kuzin, E. Sergeev, nd R. Myyä. Siution of otic coherence tomography images by Monte Carlo modeling based on polarization vector approach. Optics Express Vol. 18, Issue 21, 2010 pp. 2171421724
Моделирование оптических методов биомедицинской диагностики
Аннотация.Представлен способ МонтеКарломоделирования фотонного транспорта в биологических тканях, основанный на выполнениивысокоскоростных вычислений с помощью графических процессоров, который позволяет исследовать особенностисовременных оптических методов биомедицинской диагностики.Ключевые слова:метод МонтеКарло, графические процессоры, фотонный транспорт.
Методы медицинской визуализации, основанные на взаимодействии оптического излучения с исследуемой биологической тканью, которые включаютоптическую когерентную томографию, диффузионную оптическую томографию, микроскопию, флуоресцентную ангиографию и т.д., на данный момент являются одними из наиболее используемыхметодов медицинской диагностики. Оптические методы позволяют получить внутреннюю структуру ткани неинвазивно с высоким пространственным разрешением, что при использовании излучения небольшой мощности (~0.1 мВ) позволяет избежать разрушения структуры ткани и полностью минимизировать возможный вред для пациента.Отличительной особенностью большинства методов визуализации, основанных на взаимодействии света с исследуемой биологической тканью, является высокое пространственное разрешение получаемых изображений внутренней структуры, которое, к примеру, для конфокальноймикроскопии составляет 1 микрометр, а для оптической когерентной томографии десятки микрометров. Но, по причине высокого показателя рассеивания излучения для данных методик также характерно крайне низкое пространственное разрешение. Для усовершенствованиясовременных методов оптической диагностики необходимо полное понимание закономерностей, лежащих в основе фотонного транспорта в биологических объектах. Изучение данных закономерностей возможно путем проведения компьютерного моделирования.Как правило, моделирование распространения света в сильнорассеивающих средах проводится с помощью метода, основанного на принципах МонтеКарло. Данный метод реализован с помощью специализированного программного обеспечения MCML, которое находится в открытом доступе с середины девяностых годов прошлого столетия [1]. ДанноеПОпозволяет производить моделирование фотонного транспорта только в полубесконечной слоистой среде, при чем детектирование энергии излученияпроизводится только в точках взаимодействия фотона со средой без учета временных осцилляций оптических характеристик объекта(стационарный режим). Реальные биологические объекты, как правило, имеют гораздо более сложную пространственную структуру, изменяющуюся во времени, поэтому для их исследования необходим способ моделирования, учитывающий эти особенности. При этом стоит также отметить статистический характер метода МонтеКарло, что делает моделирование довольно медленным, в то время как программа MCML, наиболее частоиспользуемаяна данный момент, работает в однопоточном режиме, то есть позволяет использовать только одно ядро центрального процессора(ЦП), в то время как современные ЦП имеют 2128 ядер.Поэтому на данный момент данное ПО является устаревшим, как с точки зрения современных технических систем, так и по причине объективного отсутствия возможностей для исследования реальных биологических объектов и моделирование различных методик медицинской диагностики.Для увеличения скорости проведения моделирования предлагается использовать архитектуру параллельных вычислений CUDA, основанную на высокоскоростных вычислениях с помощью графических процессоров. Первоначально графические процессоры использовались только для вычисления и преобразования компьютерной графики, но с появлением программируемых шейдерных блоков появилась возможность проводить вычисления и в приложениях общего назначения[2]. Моделирование фотонного транспорта предполагает многократные вычисления распространения большого числа фотонов(10100 миллионов), что делает использование графических процессоров наиболее оправданным.Распространение фотонного транспорта можно представить в виде нескольких этапов. Сначала производится запуск фотона, определение дистанции свободного пробега между точками взаимодействия со средой, проверка на пересечение границысреды, внутри которой оптические свойства однородны, рассеивание, которое определяется коэффициентом рассеивания, специфичным для каждой ткани, и поглощения энергии фотонного пучка. Затем эти шаги снова повторяются до тех пора, пока фотон находится внутриисследуемого объекта. При моделировании различных методов медицинской диагностики данные этапы являются аналогичными, а меняется только процесс детектирования фотонов, выходящих из биологической ткани.Для изучения фотонного транспорта в объектах со сложной пространственной структурой предлагается использовать трехмерный массив, каждому элементу которогоприсваивается свой индекс, характеризующий его оптические свойства (коэффициент рассеяния, коэффициент поглощения, показатель преломления и коэффициент анизотропии). Размер и общее количество элементов массиваопределяется исходя из размера объекта исследования и пространственного разрешения методики, то есть для конфокальной микроскопии длина ребра элемента массивабудет небольшой (примерно 1 микрометр), адля методов, основанных на исследовании рассеивания и поглощения объектов, таких как диффузионная томография, размер будет существенно больше (0.10.01 миллиметра). Создание трехмерной среды предлагается проводить путем присваивания группе элементовпо ихиндексамопределенного числа, характеризующего какойлибо слой объекта, оптические свойства которого не изменяются в пространстве, либо путем задания неоднородностей внутри объекта в аналитическом виде (сфера, куб, цилиндр), соответственно те элементы, которые попадают внутрь неоднородности, будут иметь одинаковые оптические свойства. Первый случай используется для создания разреженных в пространстве неоднородностей, к которым можно отнести раковые метастазы, то есть объекты, задание границ которых в аналитическом виде невозможно. Аналитическое задание границ может использоваться при моделировании таких объектов как кровеносные сосуды (сфера, эллипсоид), сляб, фантомы биологических тканей и т.д.
В ПО MCMLрасчет дистанции свободного пробега производится с помощью показателя рассеивания и поглощения среды , что соответствует дискретному(вероятностному) подходу к процессу уменьшения энергии фотонного пучка. При проведении времяразрешенного моделирование гораздо удобнее рассматривать процесс поглощения непрерывным, согласно закону БераБугераЛамбера:
Где W–энергия текущего пакета фотонов, а s–дистанция пролета фотона. Дистанция пробега фотона же рассчитываетсятолько согласно коэффициенту рассеивания:
где
случайное число равномерно распределенное от 0 до 1. Подобный подход также позволяет увеличить скорость моделирования за счет увеличения длины свободного пробега фотона, что связано с исключением показателя поглощения из знаменателя формулы расчеты пробега.По мере распространения частицы внутри объекта фотон будет перемещать из одногоэлемента трехмерного массива в другой. В случае если показатель преломления вэлементах различен, то фотон может отразиться от границы согласно закону Френеля:
Объекты также могут изменять свою структуру во времени, что может быть связано, например, с пульсациями кровеносных сосудов. Для моделирования таких объектов необходимо отслеживать время нахождения фотона в исследуемой среде, котороевычисляется по следующей формуле:
где с –скорость света в вакууме, n –показатель преломления среды, а s–дистанция передвижения.Для каждого слоя/элемента массиваобъекта при этом задают зависимость оптических свойств от времени.
Угол рассеивания получают, используя фазовую функцию ХеньиГринштейна:
где g–коэффициент анизотропии элемента массивасреды, в котором на данный момент находится объект.Полярный угол равномерно распределяется от 0 до . Направления движения фотона при рассеивании определяется по направляющим векторам:
где , и
текущие значения направляющих косинусов, а , и
новыезначения.Проведение вычислений на CUDAпозволяет использовать библиотеки ускоренных функций для деления (__fdividef) и для расчета значений косинусов и синусов(_sinf, _cosf), вычисление которых, как правило, требует большого количества времени[3]. При этом для наиболее высокой скорости вычислений необходимо использовать числа одинарной точности float.Общее количество блоков и нитей (потоков) при вычислении на графических процессорах зависит от характеристик конкретной видеокарты. Каждая нить выполняет вычисление распространения своего фотона, обращаясь к глобальной динамической памяти только при необходимости записи значения энергии фотонной плотности и декрементирования общего числа фотонов. Для того, чтобы избежать конкурирования потоков за запись в однуячейку памяти данныешаги выполняются с помощью библиотеки атомарных функций Cuda, а именно функции atomicAdd. Данная модель реализована с помощью языка программирования CUDAC. Скорость моделирования при использовании полученнойпрограммы в общем случаев 4050 раз превышает скорость моделирования с помощью аналогичного программного кода, полученного с помощью языка C#. По сравнению с программным обеспечением, реализующем стационарную модель фотонного транспорта MCMLскорость моделирования увеличивается в 100120 раз, что позволяет получить в некоторых случаях практически мгновенные результаты моделирования.Процесс моделирования различных методов биомедицинской диагностики основан на использовании описанной выше модели, но при этом имеет свои особенностив зависимости от конкретной методики. Рассмотрим моделирование некоторых из них.Оптическая когерентна томография (ОКТ) –методика визуализации основанная на принципах низкокогерентной интерферометрии. Процедура ОКТ представляет собой последовательное сканирование объекта в поперечном направлении.В случае МонтеКарло производится последовательный запуск группы фотонов в различных точках сред и детектирование отраженных энергии отраженных частиц. После формирования одного Аскана(линии) изображения точка запуска смещается, и процесс повторяется. Интенсивность пикселя изображения в общем случае определяется с помощью следующей формулы:
где W–энергия кванта, 2zи –оптические пути в опорном и объектном плече интерферометра соответственно, а
длина когерентности источника [4].Результат моделирования ОКТ методики исследования среды, содержащей подкожный сосуд, представлен на рисунке 1.
Рис. 1. Результат моделирования подкожного сосудаСтоит отметить, что результат в целом соответствует изображениям, получаемым в экспериментальных ОКТ исследованиях. Сосудистая стенка обладаетсущественно меньшим показателем рассеивания, чем кровь (различие более чем в 200300 раз), поэтому на изображении эта область выглядит затененной, а просвет, соответственно, отличается гораздо более высокой интенсивностью, а максимальной интенсивности соответствует верхний слой кожи stratumcorneum, что связано с высоким различием показателя между ним и воздухом.
Моделирование диффузионной оптической томографии(ДОТ)отличается тем, что данная методика основана на регистрации многократно рассеиваемых фотонов. Результат моделирование ДОТ на основе исследования цилиндрического фантома представлен на рисунке 2. Оптические свойства объекта: коэффициент рассеяния –50 см^1, коэффициент поглощения –0.2 см^1, коэффициент анизотропии –0.9., радиус цилиндра –2 см.
Рис. 2. Плотность распределения энергии пучка в цилиндре, на момент 0.59 наносекунд (справа) и 0.15 наносекунд (слева).
Представленная модель фотонного транспорта позволяет выполнять высокоскоростное моделирование различных медицинских методов визуализации, что имеет большое значение для прикладных исследований. Изменяя оптические свойства среды, детектора, светового пучка при моделировании получают различные результаты и выбирают наиболее оптимальный , что может затем использоваться при реальных исследованиях, например при определении необходимости оптического просветления, настройки используемого лазера, выбора детекторов и т.д.
Ссылки на источники1. L.H. Wang, S. L. Jacques, and L.Q. Zheng MCML—Monte Carlo modeling of light transport in multilayered tissues // Computer Methods and Programs in Biomedicine 47 (2), 1995. Pp. 131–146.2. S. Mittal and J. Vetter, A Survey of CPUGPU Heterogeneous Computing Techniques, // ACM Computing Surveys, 2015.3. Джейсон Сандерс, Эдвард Кэндрот «Технология CUDA в примерах. Введение в программирование графических процессоров» // ДМК Пресс, 2011 г. –232 стр.4. M. Kiriin, I. Meginski, V. Kuzin, E. Sergeev, nd R. Myyä. Siution of otic coherence tomography images by Monte Carlo modeling based on polarization vector approach. Optics Express Vol. 18, Issue 21, 2010 pp. 2171421724
