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

На правах рукописи

Казанцев Иван Гаврилович

Численные и геометрические методы математического моделирования в многомерных задачах томографии и обработки изображений

05.13.18 — математическое моделирование, численные методы и комплексы программ

Автореферат

диссертации на соискание ученой степени доктора физико-математических наук

19 ПАП 2014

Новосибирск - 2014

005549311

Работа выполнена и Федеральном государственном бюджетном учреждении науки Институте вычислительной математики и математической геофизики Сибирского отделения Российской академии наук

Научный консультант:

доктор физико-математических наук Пикалов Валерий Владилгирович Официальные оппоненты: Аксенов Валерий Петрович доктор физико-математических наук,

ФГБУН Институт оптики атмосферы им. В.Е.Зуева СО РАН, г. Томск, ведущий научный сотрудник

Голубятников Владимир) Петрович

доктор физико-математических наук, профессор,

ФГБУН Институт математики имени С.Л. Соболева СО РАН,

г. Новосибирск, главный научный сотрудник

Чеверда Владимир Альбертович

доктор физико-математических наук,

ФГБУН Институт нефтегазовой геологии и геофизики

им. А.А.Трофимука СО РАН, г. Новосибирск, заведующий

лабораторией численных методов обращения геофизических нолей

Ведущая организация: Федеральное государственное бюджет-

ное учреждение науки Институт систем обработки изображений Российской академии наук, г. Самара.

Защита состоится 7 октября 2014 г. в 15.00 часов на заседании диссертационного совета Д 003.061.02 на базе Федерального государственного бюджетного учреждения науки Института вычислительной математики и математической геофизики СО РАН по адресу: 630090, г. Новосибирск, проспект Академика Лаврентьева, 6.

С диссертацией можно ознакомиться в библиотеке и на сайте Института вычислительной математики и математической геофизики СО РАН http://www.sscc.ru

Автореферат разослан « 7 » мая 2014 года. Ученый секретарь

диссертационного совета Д 003.061.02 доктор физико-математических наук

Сорокин

Сергей Борисович

Общая характеристика работы

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

В работе исследуется такой вид дистанционно получаемых данных, как интегральные проекции; особенно актуально их применение в томографии. Хотя математические основы были заложены в начале XX века в работах по интегральной геометрии (Г. Минковский 1905 г., П. Функ 1913 г., И. Радон 1917 г.), для полного решения проблемы реконструкции потребовались современные средства соединения цифровой техники, систем программирования, обработки изображений и численных методов обработки экспериментальных данных в единую систему вычислительной томографии. Математический аппарат многомерных интегральных преобразований лежит также в основе цифровой обработки сигналов и изображений. Известная взаимосвязь преобразований Фурье и Радона приводит к применению (лучевого) преобразования Радона в задачах анализа структур на изображениях и созданию новых быстрых многомерных алгоритмов обработки сигналов. Широкое распространение получил термин "обработка изображений в пространстве Радона". Современный уровень научных исследований требует эффективного математического моделирования и программного обеспечения для исследования влияния множества факторов на качество восстанавливаемых изображений и решения конретных прикладных задач.

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

В соответствии с поставленными целями в диссертационной работе решены следующие (и сопутствующие им) задачи:

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

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

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

4. Дискретизация четырехмерного лучевого сферического преобразования в виде квазирегулярной сетки на единичной сфере в четырехмерном пространстве и создание нового конструктивного некомбинаторного быстрого алгоритма вычисления центральных сечений многомерного куба.

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

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

Научная новизна.

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

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

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

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

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

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

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

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

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

Основные положения, выносимые на защиту

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

2. Созданные двух- и трехмерные модели веерного преобразования в задаче Радона в полосе на основе разложений в базисах веерных плоских волн и характеристических конечных элементов.

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

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

куба.

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

Апробация работы. Основные результаты диссертации докладывались и обсуждались на симпозиумах по вычислительной томографии (Самара 1985, Ташкент 1989, Москва 1991, Новосибирск 1993), международной конференции "Обработка изображений и дистанционные исследования" (Новосибирск 1990), международных конференциях "Обратные и некорректные задачи математической физики" (Новосибирск 2002, 2007, 2012), Российско-немецком семинаре "Распознавание образов и понимание изображений" (Нижний Новгород 2011), XXIII Российской конференции по электронной микроскопии (Черноголовка, 2010), IX Международном научном конгрессе и выставке "ИнтерЭкспо Гео-Сибирь-2013", международном симпозиуме "Вычислительная томография для промышленных приложений" (Берлин 2003), международных симпозиумах "Биомедицинская визуализация'' (IEEE Int. Symp. on Biomedical Imaging), (2008, 2010), международных конференциях "Медицинская визуализация" (IEEE Nuclear Sei. Symp. and Medical Imaging Conf.), (1999, 2000, 2002, 2003, 2004, 2005, 2006, 2008), международных конференциях "Обработка изображений" (Int. Conf. on Image Processing, ICIP), (1996, 1997).

На различных стадиях выполнения работа обсуждалась на семинарах, руководимых ведущими специалистами, в российских и зарубежных институтах и университетах: Институт математики им. С.Л. Соболева СО РАН (Романов В.Г., Аниконов Д.С.), Институт автоматики и электрометрии СО РАН (Потатуркин О.И., Киричук B.C.), Институт физики полупроводников им. A.B. Ржанова СО РАН (Латышев A.B.), Институт прикладной математики, Мюнстерский университет, Германия (F. Natterer), Группа цифровой визуализации, Городской университет Ныо Йорка, США (G.T. Herman), Группа обработки медицинских изображений, Пенсильванский университет, США (R. Lewitt, S. Matej), Факультет информатики и математического моделирования, Датский технический университет, Да1шя (P.C. Hansen).

Личный вклад автора. Основные научные и практические результаты диссертации получены автором лично. Вывод формул в работе [3] принадлежит автору, соавторам И.П. Яровенко и И.В. Прохорову принадлежит численная верификация методом статистического моделирования. Доказательства в работе [13] принадлежит автору, а чис-

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

Публикации. Результаты исследований но теме диссертационной работы опубликованы в 21 печатных работах в изданиях списка ВАК.

Структура и объем работы.

Диссертация состоит из введения, четырех глав, заключения, списка литературы из 190 наименований и двух приложений. Содержание основного текста диссертации изложено на 238 страницах, содержит 60 иллюстраций и 8 таблиц.

Основное содержание работы

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

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

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

В разделе 1.1 вводится лучевое преобразование (Рис. 1 (a)) Raf, отображающее функцию / на ее интегралы вдоль пучка параллельных прямых, параметризуемых углом а € [0,7г) и расстоянием s G [— 1,1] от начала координат

Raf(s)=pa(s)= / _f(s cosa — ísina, ssina + t cosa) dt. (1)

J —x/l — s7

Известна формула обращения (Й. Радон, 1917)

= - ¿УТЛ г г Ра(*Т ^ <2>

¿к /о J -1 (s- (—isma + j/cosa)J

Если множество направлений проекций представлено n-вектором ш = (соi,...,ojn) € [0,7r)n, то соответствующий набор "полных" (известных для всех s € [—1,1]) проекций обозначаем как Buf = (ВШ1/,..., BWrif). Среди аналитических методов восстановления наиболее известен метод фильтрованных обратных проекций (ФОП), являющийся следствием формулы (2), существенно использующий направления и>, равномерно распределенные на окружности. В случае произвольного набора направлений и) предпочтительнее подход (В. Logan, L. Shepp, 1975), основанный на существовании и единственности аппроксимации /* функции /, согласованной с проекционными данными B.uf, в виде суперпозиции (вообще говоря, не единственных) плоских волн заданных направлений /* = Н^(х,у) = ^2i_1hi(xcosuJi + у sinw,). Плоские волны hi(x cos u>i + у sin Uli) представляют собой обратно проецированные версии функций hi одной переменной. В отечественной литературе имеются работы (Б.А. Вострецов, М.А. Крейнес, 1961) с теоремами существования приближения методом СПВ непрерывных функций в общем многомерном случае. Для Нш справедливо равенство ||/ — Нш ||2 = ||/||2 — ||ЯШ||2. Из этого следует, что знание нормы ||ЯШ||2, называемой в диссертации информативностью множества проекций Buf, полезно в оценивании точности восстановления по данному набору проекций.

В разделе 1.1 получен новый конструктивный алгоритм вычисления плоских волн, в сумме составляющих на основе сингулярного разложения лучевого преобразования (Теорема 1.4). Показано , что функции hi,..., hn имеют вид

i ОО 72 Г1

hM = ;ЕЕ n^Uk-i(s) / RvjWk-im, (3)

fc=ij=i J

Рис. 2. Томограммы головного мозга 160 х 160 пикселов, полученные из 256 проекций. (а) Реконструкция по методу ФОП. (б) Реконструкция по методу плоских волн, (в) Центральные столбцы изображений (а) и (б) - вверху и внизу, соответственно.

где и к ~ полиномы Чебышева второго рода, щ- - элементы матрицы, обратной к матрице Л*, = {А^}, г.] = 1 ,п, А-^ = кт(/с(шг — и)у))/(А; зт(ц — ш^)). Впервые получена формула обращения, в которой I информация о геометрии просвечивания сконцентрирована в матрицах А/с- В случае равномерно распределенных углов Шг алгоритм становится сравнимым по быстродействию с ФОП в силу аналитической формулы обращения матриц Л/;. Доказывается (раздел 1.1.1, Теорема 1.6), что

к2

для к < п, обобщенная обратная Л^ имеет вид = —^Л^. В разделе 1.1.2 получены аналогичные новые формулы в случае проекций с отсчетами в виде интегралов по полосам с шириной, определяемой апертурой детекторов. В разделе 1.1.3 полученные формулы апробированы на реальных данных магнитно-резонансного томографа (Рис. 2). '

В разделе 1.2 дискретная версия лучевого преобразования рассматривается в алгебраической постановке в виде СЛАУ I

А/ = р, (4) |

где / = е М-7, р = {рг} еШ1, А = {аи} е М/„7 есть матрица с I строками и Л столбцами. Вектор / представляет цифровое изображение с прямоугольной 3\ х сеткой пикселов, вытянутое по лексикографи-чекому правилу в вектор размером ./ = ./] х ./г- Вектором данных р является набор лучевых (полосовых, конусных) интегралов, А - большая разреженная бесструктурная матрица с неотрицательными весами, определяющими вклад пиксела изображения в формирование проекционных данных р : простейшая модель предполагает, что агз есть длина (площадь) части г-то луча (полосы, конуса, призмы), проходящего через ,7~й пиксел (Рис. 1 (б)). Матрица А называется проекционной матрицей.

Таблица 1. d-меры по 4 проекциям (без шума).

Итерация w = 0 w = 1 w — 1.5 w = 2 w = 3 w = 4

1 0.062 0.059 0.089 0.119 0.150 0.212

3 0.050 0.048 0.045 0.045 0.047 0.059

6 0.047 0.046 0.044 0.043 0.042 0.045

8 0.046 0.045 0.043 0.042 0.041 0.044

Wloc

Методы решения систем (4) в томографии условно объединены в группу ART ("Algebraic Reconstruction Techniques").

В разделе 1.2 исследуется возмущенная система Awf = pw, где параметр w есть ширина полосы интегрирования. Доказывается Утверждение 1.10 о существовании малого w > 0 такого, что оценка /* = A^pw является локально лучшей, что впервые обосновывает варьирование ширины полосы как параметра контроля за точностью реконструкции. В подразделе 1.2.1 проведены вычислительные эксперименты по проверке доказанного утверждения. Точность восстановления оценивается в виде нормализованной среднеквадратичной ошибки, определяемой

как d = Е;=1 (/,(fc) - /,) 7 Ej=1 if} - Л2- Здесь /(*) = {/}*> } есть восстановленная томограмма после к итераций, {fjj - тестовое изображение, / - его среднее значение. Новая оценка/(к+11 определяется из f(k) после к итераций коррекционной процедурой ART (С. Качмаж, 1937)

/№+i) = /№) + А ft - у.а<) { = k{modI) + lf (5)

где Л - параметр релаксации, 0 < Л < 2, аг - г-я строка матрицы А, /° з 0. Проекции от тестового объекта в виде суперпозиции нескольких гауссиан генерировались при различных уровнях шума. Оценивались реконструкции, полученные для различных значений ширины w полосы интегрирования. Локальная нелинейная регуляризация линейной системы с помощью варьирования параметра носителя (ширина полосы, угол раствора конуса) дает улучшение точности восстановления при ширине полос интегрирования и>[ос (Таблица 1).

В разделе 1.3 исследуется информативность Q(f,ui) = \\Нш{х, у)\\2 набора проекций Rwf. Информативность отдельной проекции рш вычисляется как Q(f.u)j) = 1/2 f (1 -s2)~1/2p2 (s)ds. Наиболее инфор-

7—1 J

мативный набор направлений и>о находится глобальным поиском

ш0 = argmax Q(f, cj). (6)

о»

Рис. 3. (а) тестовая текстура "Орнамент", 128 X 128, ||/||2 = 85473; (б) лучевое преобразование по горизонтали отложены углы а; (в) информативность Q

отдельных 180 проекций ; (г) реконструкция Ншо по 4 проекциям в направлениях ш0 = (0°, 48°, 90°, 131°), ||ЯЫ0 ||2 = 83757, ||ЛГШо ||2/||/||2 = 0.98.

Фунционал информативности Q впервые используется для моделирования текстурных изображений и оптимального планирования облучения в радиотерапии. Регулярная текстура неформально определяется как изображение, воспринимаемое как повторяющиеся агрегаты образов (текстоны). Регулярные (в отличие от изотропных стохастических) текстуры обладают несколькими выделенными направлениями цикличности текстонов и предполагается, что такие текстуры / хорошо аппроксимируются суммой небольшого числа (п = 2,3,4) плоских волн, т.е. информативные проекции содержат информацию об одномерных регулярных структурах, суммарно составляющих целую текстуру.

В разделе 1.3.1 проведены численные эксперименты по декомпозиции текстур в информативные плоские волны (Рис. 3). Численно сгенерированы 180 проекций в дипазоне [0, 7г). Наилучшими четырьмя направлениями оказались и>о = (0°, 48°, 90°, 131°); соответствующая декомпозиция по этим направлениям показана на Рис. 3 (г).

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

д(х,у) = }{х,у) + Ъ(х1у), (7)

где д есть наблюдаемое изображение; / - искомое полезное изображение, как правило, с множеством деталей; b - некоторая мешающая текстура, о которой известно, что она хорошо представима суммой НШ0Ь = ^ /гШ; небольшого числа плоских волн в направлениях u>o = (wi,..., и>п), т.е. b » НШоЪ. Задача - восстановить / по данным д и ojq. В силу аддитивности модели (7) и линейности лучевого преобразования (Радона), разложение функций д, f,b в п плоских волн в направлениях и о приводит к равенств}^

(а) (б) (в) (г)

Рис. 4. (а) Тестовое изображение Г - "Фотограф"; (б) сумма g = { + Ь тестового изображения Г и текстуры Ь "Орнамент"; (в) реконструкция изображения g но 4 проекциям в информативных для текстуры Ь направлениях шд = (0°, 48°, 90°, 131°); (г) разница § — Н^оё между изображениями (б) и (в).

НШод = Нша/ + НШоЪ. (8)

Поскольку текстура Ь хорошо аппроксимируется плоскими волнами в направлениях и>о, вычтем уравнение (8) из (7) и получим

д - НШод « / - НШ0/. (9)

Идея метода состоит в том, что функция Нша / является очень размытой версией искомой функции /, поскольку п мало для качественной реконструкции функции /. Поэтому разница (/(ж, у) — НШо/(х,у)) в правой части (9) приближенно равна / со смещенным средним значением. В результате синтезируется изображение, существенно избавленное от текстуры Ъ. Численные эксперименты представлены на Рис. 4.

В разделе 1.3.2 рассмотрена проблема планирования доставки определенной радиационной дозы к выявленной ранее при томографическом сканировании критичной области организма пациента, с целью уничтожения там злокачественных клеток, минимизируя при этом радиационное воздействие на окружающие здоровые ткани и жизненно важные органы. Задача вычисления п (обычно 5-9) профилей радиационной интенсивности облучения по заданной плановой дозе впервые сводится к задаче декомпозиции дозы в сумму плоских волн. Для этого плановая доза рассматривается как известное изображение /, затем вычисляется достаточно большое число (90-180) виртуальных проекций, среди которых перебором находится набор п информативных направлений, из которых производится облучение в соответствии с амплитудами вычисленных плоских волн. Для численного эксперимента выбрана, плановая доза средней сложности (Рис. 5 (а)), которая имеет Гамма-образную форму с органом, требующим щадящего воздействия в правом нижнем углу. Информативными направлениями оказались

(а) (б) (в) (г)

Рис. 5. (а) Тестовое изображение "Гамма"; (6) реконструкция методом ФОП по 5 равноотстоящим проекциям (0°,36°,72°,108°,144°); (в) реконструкция методом ФОП по 5 информативным направлениям (и>о = 0°,26°,64°,90°, 135°) ; (г) реконструкция методом СПВ по 5 информативным проекциям.

Wo = (0°, 26°, 64°, 90°, 135°), найденные по формуле (6). Проведенные численные эксперименты по сравнению метода суперпозиции плоских волн (СПВ) и метода фильтрованных обратных проекций (ФОП) показали преимущество метода СПВ.

При рассмотрении задач главы 1 исследовались три основных метода - ФОП, СПВ и ART. В формулах восстановления но методу СПВ имеется алгебраическая составляющая и этот метод имеет некоторые общие характеристики и с методами ФОП. и с методами ART. Это неслучайно, поскольку метод СПВ является непрерывным аналогом более общего дискретного подхода, называемого методом характеристических конечных элементов (ХКЭ) и кратко введенным в разделе 1.4. В этом методе на изображении / геометрия сбора данных задает I областей интегрирования Di и их характеристические функции ipi (Рис. 1 (в)). Их общее число I определяет общее количество измерений Pi = J J f{x, y)ipi(x,y)dxdy. Оценка функции / производится обрат-

I

ным проецированием и суперпозицией в виде /*(ж, у) — У^ у).

3=1

Это приводит к алгебраической системе Гс = р относительно скалярных коэффициентов с, где Г = {7- матрица Грама с элементами

7ij = J J xjji(x,y)ipj(x,y)dxdy. Этот метод применим к произвольным многомерным нерегулярным геометриям просвечивания и используется в главе 2.

В разделе 1.5 представлены выводы по главе 1.

В Главе 2 рассматривается лучевое преобразование Радона в полосе в двух- и трехмерных постановках, в дискретном и континуальном случаях. Данными являются интегралы по линиям, соединяющим точки ("источники" и "приемники"), находящиеся па противоположных

(а) (б)

Рис. 6. (а) Лучевое преобразование Радона в полосе между линейками А и В. (б) Дискретная версия в трехмерном случае. Данные интерпретируются как интегралы по наклонным призмам (характеристическим конечным элементам) ¡/»¿, ipj . Элементы матрицы Грама -yij — объем пересечения элементов и tpj.

сторонах А и В прямоугольной области, или "полосы" (Рис. 6 (а)). Лучевое преобразование Радона в полосе (C.B. Гольдин, 1996) определятся как набор интегралов от функции / с носителем в прямоугольнике [—а, а] х [—6,6], вдоль отрезков C(s,q) длины L = L(s,q), соединяющих точки (s,—b) и (q,b), для s,g G [—а,а]:

ь

d{s,q) = f f(x,y)dl = If((s,-b) + j{q-s,2b)yi. (10)

C(s,q) 0

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

В разделе 2.1.1 приведены матрицы Грама для различных геометрий сбора данных. Вычислительные модельные эксперименты с матрицами Грама небольших размеров проведены в разделе 2.1.2.

Если панели А и В (Рис. 6 (б)) покрыты сеткой N х N детекторов каждая, то размер матрицы Грама этой геометрии есть N4 х jV4. Большие размеры матриц Грама в практических задачах стали мотивацией для разработки в разделе 2.2 новой аналитической математической модели прямого преобразования Радона в полосе. Данные (10) интерпретируются как веерные проекции d(s,q) с "источником" в точке (s, —Ь) и "детекторами" в точках (q,b). Обобщая понятие плоской параллельной

о

(б)

В 3

(в)

Рис. 7. (а) Двумерный случай двух пар линеек {А, В} и {С, £>}. (б) Принцип динамического охвата объекта / подвижными панелями детекторов, (в) Двумерная схема рассеяния. Фотоны на линии MN могли бы разделить точки С/ и V, однако они не регистрируются линейкой детекторов В. Интегралы по дугам вРд (пунктирные кривые) регистрируются комптоновской томографией, если фотоны попали в точки 5 и д, претерпев рассеяние в точке Р ломаной линии вРд.

волны, определим для функции г>(5,д) = У8(д), € [—а,а], ее веерную обратную проекцию, или веерную плоскую волну у£(х,у), (х,у) £ [-а, а] х [~Ь,Ь]

Ув(х,У) = (х,у) € ¿(в,?). (И)

Ищем аппроксимацию / в виде суперпозиции веерных плоских волн у^

а

Г(х,у) = IуГ(х,у)<1з, (х,у) е [-а,а] х [-6,6], (12)

— а

согласованной с данными ¿(в, q). В подразделе 2.2.1 доказывается (Теорема 2.1), что являются решениями уравнения

—а —а

Следствием уравнения (13) является новый быстрый алгоритм

г; « -¿ВД, (14)

где И есть преобразование Гильберта по одной (любой) из переменных в, д правой части Б. Любая плоскость, пересекающая панели А и В в трехмерной задаче (Рис. 6 (б)), индуцирует множество двумерных задач Радона в полосе, поэтому быстрое решение набора плоских задач является актуальным. Подобное уравнение выведено также для случая

Рис. 8. (а) Восстановление по линейкам А и В. (б) Восстановление по линейкам С и V. (в) Восстановление по объединению данных {А, В}и{С, И}, (г) Восстановление (7-я итерация) алгоритма Качмажа по данным {Л, В} и {С, И}.

(раздел 2.2.3), когда детекторы расположены на дугах окружности. Полученные алгоритмы апробированы в сравнении с методом ART (разделы 2.2.2, 2.2.4).

Рассмотрена также задача с данными, расширенными до отрезков интегрирования, соединящих две другие противоположные стороны прямоугольного носителя функции /, т.е. данные регистрируются двумя парами противолежащих линеек детекторов {А, В} и {С, D} (Рис. 7 (а)). Эта перспективная конфигурация сканера позволяет сдвигать панели {А, В} и {С, D} до оптимального прилегания к объекту просвечивания / (Рис. 7 (б)). Проведены вычислительные эксперименты (Рис. 8) для геометрии сбора данных в виде пары линеек {А, В} (Рис. 8 (а)) и {C,D} (Рис. 8 (б)), двух пар линеек {А, В} U {C,D} (Рис. 8 (в)). Показано,что быстрый алгоритм реконструкции дает результаты, сопоставимые с итерационным алгоритмом Качмажа (Рис. 8 (г)).

Задачи в полосе возникают при межскважинном просвечивании в геофизике, где противоположные стороны полосы моделируют скважины с источниками и приемниками. В разделе 2.3 проведен модельный эксперимент (Рис. 7 (в)), показывающий, что добавление к имеющимся веерным проекциям в задаче Радона в полосе значений интегралов по дугам определенных окружностей, проходящих через точки детекторов А и В , существенно улучшает разрешение. Интегралы по дугам являются данными в модели эмиссионной томографии на однократном комптоновском рассеянии, рассматриваемой в главе 3.

В разделе 2.4 представлены выводы по главе 2.

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

Vf

-Pf^f

f -»■ Pkf

rkf

C2 (17/1) : Коррекция

Сз(|г/|) : Восстановление

VBK~xVkj

BK~lVkf

В

К;

K-Zlvkj

Рис. 9. Диаграмма обработки данных. Объект / отображается в набор расфокусированных проекций Алгоритм ОППФ состоит из трех этапов — послойной фильтрации КГ1, последующем обратном проецировании В и трехмерной декон-волюции Сз с ядром . В качестве альтернативы показана коррекция, состоящая из репроекции V промежуточного результата ВК^1Т3>! I с последующей двумерной деконволюцией С2 с ядром и дающая в результате проекции "Р/ рентгеновской томографии.

виде набора проекций

оо оо со

Kpf^y) = 9аАх,у) = / / J fa'4^,y',z)k(x-x',y-y',z)dx'dy'dz,

— ОО—ОО—ОО

(15)

где fa,,e - вращательный вариант функции /пае [0,2тг), /3 G [0,7г), — сферические координаты, определяющие поворот в системе координат Oxyz. Модель формирования данных в частотной области имеет вид

оо

Gaip(rjx,T]y) = J Fa'p(t]x,i]y, z)K(j]x,T]y,z)dz, (16)

—00

где Ga,p - двумерное преобразование Фурье проекции да,в\ Fa'P и К - преобразования Фурье сечений плоскостью г = const вращательной версии объекта /а"3 и ядра к соответственно. Передаточная функция К имеет форму, зависящую от применямого вида томографии.

В случае классического лучевого преобразования V (Рис. 9), интегральное изображение как результат суммирования по сфере S2 всех трехмерных плоских волн, т.е. обратно проецированных данных ра,р, символически можно записать в виде b(x, y,z) = J pajjdn d[i. Метод

s2

обратного проецирования с послойной фильтрацией (ОППФ), разработанный в разделе 3.1, использует обобщенную версию плоских волн в виде Wat/3(x,y,z) = T21{Ga^(r]x,r]y)K~1(r]x,r]y,z)}, т.е. компенсирующей фильтрации регуляризированным оператором К(К+), произведенной в каждом слое z. Суммирование по сфере всех обобщенных

плоских волн и>а%р дает в случае Тк (Рис. 9) интегральное изображение метода ОППФ (аналогично классическому лучевому преобразованию) в виде и>(х,у,г) = JПрименяя метод стационарной фазы,

доказано (Теорема 3.1), что при наличии полных данных, преобразования Фурье интегральных изображений Ъ и ы асимптотически равны на высоких частотах. Этот новый результат впервые сводит задачу с послойными искажениями к классическому лучевому преобразованию, открывая этим возможность использования всего существующего арсенала методов восстановления изображений по лучевым проекциям.

В разделе 3.2 впервые показано, что модель (15) формирования изображений с послойными искажениями применима и к задаче восстановления распределения активности изотопов в организме по данным комп-тоновского рассеяния позитронной эмиссионной томографии (ПЭТ). Рассматривается традиционная модель формирования данных ПЭТ на первичных (не рассеянных) фотонах в виде:

РАВ = ехр

ff(x,y,z)dl, (17)

где f(x,y,z) - распределение источников активности изотопа внутри среды с линейным коэффициентом ослабления ц(х,у,г), А и В - точки детектирования (Рис. 10).

Задача состоит в реконструкции f при известном ^ по данным РАВ, регистрируемым большим множеством пар детекторов (А, В). Физическая модель эмиссионной томографии основана на использовании пар гамма-квантов, разлетающихся в противоположных направлениях из точек носителя D(f) функции активности /. Данные детекторов РАВ интерпретируются как интегралы от / по прямым линиям распространения первичных фотонов с энергией Е = 511 кэВ. Однако в детекторы попадают и фотоны с меньшими энергиями Е' < Е, претерпевшие комптоновское рассеяние (Рис. 11 (а)) , которые в вычислениях не используются. В зависимости от разрешения детектора по энергии, используется спектральное окно чувствительности W = [f, 511] с порогом t и фотоны с энергией вне окна W в алгоритмах не учитываются. Различными методами в сканерах оценивается отношение числа первичных к рассеянным фотонам и

Рис. 10. ПЭТ на первичных фотонах (и, V) с энергией Е = 511 кэВ, С £ £>(/) - точка аннигиляции.

Рис. 11. (а) НЭТ на однократном комптоновском рассеянии; (и, v) - napa, анни-гиляционных фотонов, S 6 D(fi) - точка рассеяния, v' - фотон v энергии Е1, рассеянный с углом в. (б) Поверхность Sg (тело вращения дуги AS В)- геометрическое место точек S рассеяния с углом в и сферическими координатами (i¡>,<p, |-4S|).

затем данные масштабируются в соответствии с этим отношением.

С улучшением спектрального разрешения детекторов появляется возможность настраивания сканеров на регистрацию фотонов, рассеянных под определенным углом в, соответствующим энергии Е' с последующим использованием информации в этих данных в дополнение к первичным фотонам. В разделе 3.2.1 доказана (Рис. 11 (б)) Теорема 3.2. Мгновенное количество фотонов, рассеянных с углом в, регистрируемых детектором В при условии, что в А попали первичные фотоны, описывается формулой:

в 2тг

САВ f, cosycos(<¿> - в) Г , ц(-ф,у,\А5'\) дас s,

о о

- | JHdl+Jn'dl |AS| е V^ s > f f{4>,<p,r)dr, (18)

o

где - дифференциальное сечение комптоновского рассеяния, т.очка dfi

рассея.ния S определяется в сферических координатах как (ip,ip, |-<4>5|), расстояние от точки А до S есть |Л5| = | А В | Sm^n g ^ , м' - коэффициент ослабления после рассеяния.

В разделе 3.2.2 полученная новая формула (18) исследована в сравнении с проекциями, генерированными статистическим методом Монте-Карло. Были проведены численные эксперименты с фантомом, изображенным на Рис. 12. Детекторы линейки А регистрируют первичные фотоны, а детекторы В — рассеянные. Каждая линейка состоит из 100 детекторов сечения 5 = 0.16 см. Проекции (формула (18))и

Рис. 12. Томографический фантом, сечение уОг (слева); общий вид (справа). Цилиндр диаметром 16 см заполнен водой, внутри помещены четыре стационарно излучающих сферических источника с диаметрами 1.8 см, заполненные фармпрепаратом с единичной активностью.

Мв (статистическое моделирование) вычисляются для 100 пар приемников {А, В}. Во всех расчетах при моделировании траектории фотона отслеживалось до 10 актов взаимодействия со средой. Моделировалось 1011 траекторий (время расчета на компьютере Pentium 3.2 GHz около 5 суток). При расчете профилей сигнала по формуле (18) ко-эфициент ослабления в воде принимался одинаковым для всех энергий f_i — ¡j! = 0.096 см-1. Проекции £в и Мд масштабировались умножением на скаляр по принципу совмещения точек максимума и приводились к условным единицам в интервале [0,10]. Графики на рис. 13 демонстрируют поведение сравниваемых величин в случае, когда отслеживается рассеяние на угол в = 7г/3. Наблюдается хорошее соответствие между проекциями, полученными по формуле (18) и статистическим моделированием.

В разделе 3.3 выводится частный вид нового преобразования (18) в

Рис. 13. Сравнительные результаты моделирования. Прерывистые линии представляют проекции Мд, насчитанные методом Монте-Карло, сплошные - проекции (а), (б), (в), (г)- Проекции от сфер №1, №2, №3 и №4 соответственно.

Рис. 14. Фантом "Молекула" размером 1283. Ряд (а)- Сечения исходного фантома. Ряд (б)- Сечения реконструкции методом ОППФ. Ряд (б)- Сечения реконструкции методом ОППФ с использованием центрального сечения передаточной функции для всех слоев. Ряд (г)- Реконструкции методом ФОП без учета искажений.

случае ц = ц' = const, что является допустимым при моделировании эмиссионной томографии человека. Показывается, что в этом случае модель (18) сводится к лучевому преобразованию с послойной сверткой (15).

В разделе 3.4 проведены численные эксперименты с математической моделью трансмиссионной электронной микроскопии с передаточной функцией в виде

K(V,z)= Sm(27rA42(-^ + ^^))e-"a«oa(c.AV-^WA4)a> (19)

где г] = <Jvx + Vy ~ пространственная частота, Л = 0.033487 А - длина

волны электрона, q0 = 0.00746558 A"1, Cs = 22, 000, 000 А - коэффициент сферической аберрации линзы микроскопа, ф(г) - значение расфокусировки микроскопа в плоскости z, ф{г) £ [1000,3000] А.

Для обращения использован регуляризирующий фильтр Тихонова

К(v,z) = ТТ K}?'Z\-о • (20)

где а = 0.01. Вычислена трехмерная реконструкция методом обратного проецирования с послойной фильтрацией (Рис. 14) по 5,000 проекциям, выбранным в случайных направлениях. Проведены численные эксперименты с добавлением шумов в проекционные данные.

В разделе 3.5 приведены выводы по главе 3.

В главе 4 рассмотрено четырехмерное лучевое сферическое преобразование и вопросы его дискретизации с целью обращения алгебраическими итерационными методами. Сферическое лучевое преобразование (СЛП) определяется (Г. Минковский, П. Функ, С. Хелгасон, В. П. Пала-модов) в виде интеграла скалярной неотрицательной четной функции /

2тг

х[С(РШ/) = ^/с(Р1Я)Пх)*х = ±1 /(«(¿))А (21)

по большим окружностям

•V

C(P2,Q2)

к

C(P,,Q¡)

S3 R4

■ 0.2"

С(Р, Q) = { x(t) I x(t) = Р eos(t) + Q sin(í), í e [0,2тг) },

заданным парой ортогональных единичных векторов P,Q G S3.

Лучевое сферическое преобразование (21) возникает при решении обратной задачи кристаллофизики, а именно задачи восстановления функций / распределения ориентаций (ФРО) кристаллической решетки в образце кристаллического вещества конечного объема по данным дифракционного эксперимента. В дифракционном эксперименте (Рис. 15) функция / недоступна для прямого измерения, однако интенсивности яркостей (Д и - отсчеты оцифрованных интенсив-ностей) на дифракционных картинах пропорциональны интегралам от / по некоторым большим кругам (C(Pi,Qi) и C(P2,Q2)) на сфере S3, лежащим в плоскостях, положение которых определяется из координат отсчетов яркостей и геометрии регистрации данных. Четыре малые сферы (искомая текстура) (Рис. 15) иллюстрируют функцию / большого разрешения, используемую при численном моделировании.

Рис. 15. Томографическая схема дифракционного эксперимента

(а) (б)

Рис. 16. (а) Большая окружность С сферы §3 схематически; I - интеграл от функции распределения ориентаций по окружности, (б) Половина развертки границы гиперкуба, изоморфная половине сферы Б3; набор отрезков - гномо-нический образ окружности С.

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

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

В разделе 4.1 дано описание геометрии границы четырехмерного куба (гиперкуба) [—1,1]4, состоящей из восьми граничных трехмерных кубов Ж-ьь (к = 1,2,3,4), или граней гиперкуба

¥±к = {а; = (х1,х2,х3,х4) е!4 | хк = ±1; -1 < х^ < 1, г ф к}. В общем случае, если <3 = (<?ъ 12, ■ ■ ■, Чы) 6 В^-1, то отображение

^ : Я ^ (22)

где

?тв*==1Шж{|в1|)Ы>...)Ы},

(23)

называется гномоническнм проецированием сферы §Af_1 на границу F = U^jF-tfc куба К = [—1,1]дг. Гномонический проектор устанавливает взаимноодназначное соответствие между §ЛГ_1 и границей F куба К. Любой большой круг на §3 гномонически отображается в некоторое (не более восьми) множество отрезков связной замкнутой ломаной внутри граничных кубов F±fc. Так как / - четная функция, то для ее изоморфного представления достаточно только четыре трехмерных куба из восьми (Рис. 16). При дискретизации граничных кубов в регулярную сетку равных элементарных конечных элементов, мы получаем соответствующую (гномонически) квази-регулярную сетку на сфере §3 в форме элементарных криволинейных сферических вокселов (раздел 4.2). Это определяет дискретизацию, или разбиение большой окружности С в набор элементарных дуг, что дает конечную сумму - аппроксимацию интегралов в лучевом сферическом преобразовании. Веса в интегральной сумме являются длинами дуг большой окружности, попавших в конкретный сферический воксел. Такой подход сводит проблему трассировки окружности на §3 к задачам:

• трассировка граней: какая из граней F±fc пересекается в данный момент t плоскостью, содержащей окружность С — Р eos t+Q sin í;

• трассировка отрезка линии внутри посещаемой грани.

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

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

Шаг 1. Введем систему координат (р, q) и отобразим на ней 2N точек Т±к = ±{pk,qk), к = 1,..., Л7, составленных из компонент ортонорми-рованных векторов Р = {рк}, Q = {Чк}-

Шаг 2. Вычислим выпуклую оболочку точек Т±к-

S = conv({T±k}). (24)

Вычислим число вершин М = |S|.

Шаг 3. Упорядочим соседние точки Ткт в многоугольнике S против часовой стрелки, начиная с некоторой (любой) . Переобозначим их в Sm, то = 1, М + 1, с условием Sm+i = Si = .

(а)

(б)

(в)

Рис. 17. Малоразмерный вычислительный эксперимент с тестовым массивом размера 4 х 643 — 1. 048, 576. (а) .'5Г) изображение размером 64 х 64 двумерного сечения одного из граничных 3-кубов фантома, содержащего гладкий параболический фон и 4 малых сферы (Рис. 15). (б) ЗЭ изображение того же сечения реконструкции по I = 6 х 106 отсчетам, (в) Реконструкция по I = 24 X 10е отсчетам.

Шаг 4• Вычислим угловые координаты тт точек Ет ^ 1:

таких, что ОЕт _L SmSm+i. Иными словами, точки Ет служат основаниями перпендикуляров, опущенных из начала координат на отрезки SmSm+1 (или их продолжения).

Шаг 5. Вычислим точки Dm = PcosTm + Qsinrm.

Шаг 6. Вычислим Vm = Q[Dm\ - гномонические образы точек Drn. Это и есть вершины искомого многоугольника, т.е. сечения куба. Основанием для Шага 2 алгоритма является доказанная

Теорема 4.4. Область сферы §N~l , чей гномонический образ есть грань F-tfc для некоторого к, пересекается окружностью С(Р, Q) тогда и только тогда, когда точка Тпринадлежит выпуклой оболочке S.

В разделе 4.4 проведены вычислительные эксперименты, показавшие эффективность разработанного алгоритма. Шаги 1-6 выполняются достаточно быстро, и это ключевой момент во всей задаче ускорения трассировки большой окружности. При трассировке большой окружности С{, соответствующей отсчету pi,i = 1,1, вычисляется г-я строка проекционной матрицы А — {a,ij},j = 1, J с весами интегрирования дискретного изображения {fj} вдоль C¿. Для решения системы Af = р используется итерационный метод Качмажа (5).

В разделе 4.4.1 проведены маломасштабные вычислительные эксперименты по восстановлению функции распределения ориентаций высокого разрешения с использованием алгебраических итерационных методов - на сетке четырех сферических кубов J = 4 х 643 = 1,048,576 (Рис. 17), в разделе 4.4.2 - крупномасштабные на сетке J = 4 х 7003 = 1,372,000,000 и 103 х 103 х 103 элементов, и в разделе 4.4.3 - экспери-

(а) (б) (в)

Рис. 18. Численный эксперимент, (а) Сечение грани Fi со сферическими объектами фантома, (б) Восстановление на сетке постоянного размера 4 X 2563. (в) Восстановление с использованием последовательности уменьшающихся се ток.

менты с использованием последовательности уменьшающихся сеток (с удвоением разрешения на последующем уровне) на сетках размерами до J = 4 х 2563 = 67,108, 864 элементов (Рис. 18). Вследствие большого числа вокселов многие узлы сетки на уровне самого высокого разрешения остаются нулевыми с момента начала итераций и не участвуют в итерационных обновлениях. Это проявляется в зернистом шуме на реконструкциях (Рис. 18 (б)). Применение уменьшающихся сеток помогает в решении этой проблемы на первых шагах удвоения сетки (Рис. 18 (в)).

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

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

В приложении I описана прикладная программа для моделирования и визуализации численных экспериментов по реконструкции и обработке изображений. Пакет программ объединен в интерактивное диалоговое приложение, созданное на языке МАТЛАБ. В приложении II собраны использованные сокращения сочетаний терминов.

В Заключении диссертации сформулированы основные результаты.

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

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

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

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

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

Публикации в журналах списка ВАК

[1] Kazantsev I.G., Schmidt S. A spherical x-ray transform and hypercube sections // Journal of Inverse and Ill-posed Problems. — 2014. — V. 22, No. 4.— P. 471-483.

[2] Schmidt S., Gade-Nielsen N.F., Hostergaard M., Dammann В., Kazantsev I.G. High resolution orientation distribution function // Materials Science Forum. — 2012, — V. 702-703, No. 1, —P. 536-539.

[3] Казанцев И.Г., Яровенко И.П., Прохоров И.В. Моделирование процесса измерения комптоновского рассеяния в позитронной эмиссионной томографии /'/ Вычислительные технологии. — 2011. —Т. 16, №. 6. —С. 27-37.

[4] Kazantsev I.G., Klukowska J., Herman G.T., Cernetic L. Fully three-dimensional defocus-gradient corrected backprojection in cryoelectron microscopy 11 Ultramicroscopy. —2010. —V. 110, No. 9, —P. 1128-1142.

[5] Kazantsev I.G., Schmidt S., Poulsen H.F. A discrete spherical x-ray transform of orientation distribution functions using bounding cubes /'/ Inverse Problems.— 2009. —V. 25, No. 10.—P. 105009.

[6] Matej S., Kazantsev I.G. Fourier-based Reconstruction for Fully 3-D PET: Optimization of Interpolation Parameters // IEEE Transactions on Medical Imaging.— 2006.—V. 25, No. 7, —P. 845-854.

[7] Kazantsev I.G., Matej S., Lewitt R.M. Inversion of 2D Planogram Data for Finite-length Detectors // IEEE Transactions on Nuclear Science. — 2006.— V. 53, No. 1. —P. 160-166.

[8] Matej S., Fessler J.A., Kazantsev I.G. Iterative Tomographic Image Reconstruction Using Fourier-Based Forward and Back-Projectors // IEEE Transactions on Medical Imaging. — 2004. — V. 23, No. 4,—P. 101-412.

[9] Kazantsev I.G., Matej S., Lewitt R.M. System and Gram Matrices of 3D Planogram Data 11 IEEE Transactions on Nuclear Science. — 2004. — V. 51, No. 5, —P. 2579-2587.

[10] Kazantsev I. G., Lemahieu I., Salov G.I., Denys R. Statistical detection of defects in radiographic images гп nondestructive testing // Signal Processing.— 2002, —V. 82, —P. 791-801.

[11] Kazantsev I.G., Lemahieu I. Reconstruction of elongated structures using ridge functions and natural pixels // Inverse Problems. — 2000. — V. 16, No. 6.— P. 505-517.

[12] Kazantsev I.G., Van de Walle R., Lemahieu I. Ridge Functions, Natural Pixels and Minimal Norm Reconstruction // IEEE Transactions on Nuclear Science—2000,—V. 47, No. 3, —P. 1118-1122.

[13] Kazantsev I.G., Pickalov V.V. On the accuracy of line-, strip- and fan-based algebraic reconstruction from few projections // Signal Processing.—1999.— V. 78, No. 1, —P. 117-126.

[14] Казанцев И. Г. О формуле обращения преобразования Радона для конечного числа проекций // Доклады Академии Наук.—1999.— Т. 364, № 4.—

С. 447-448.

[15] Kazantsev I.G. Tomographic reconstruction from arbitrary directions using ridge functions // Inverse Problems.—1998.- V. 14, No. 3. —P. 635-645.

[16] Kazantsev I.G., Pyatkin V.P., Salov G.I. The Tomographic and Statistical Approach to Detecting Anomalous Structures in Aerospace Images // Pattern Recognition and Image Analysis.—1996,—V. 6, No. 4, —P. 682-686.

[17] Казанцев И.Г. Выделение структур цифровых изображений с полю-щью преобразования Радона // Распознавание образов и анализ изображений.—1992.—Т. 2, № 2,—С. 208-210.

[18] Kazantsev I.G. The criterion of the informative projection choice in computed tomography // Computers and Artificial Intelligence.—1991. — V. 10, No. 6. —P. 581-587.

[19] Kazantsev I.G. Information Content of Projections // Inverse Problems.— 1991, — V. 7, No. 6, —P. 887-898.

[20] Kazantsev I.G. Radon-Space Straight Edge Detection in Digital Images // Computers and Artificial Intelligence.—1989.—V. 8, No. 2.—P. 189-196.

[21] Алексеев А.С., Казанцев И.Г., Пяткин В.П. Томографический подход к выделению лгтеамснтов на аэрокосмических изображениях // Исследование Земли из космоса. —1988. — № 5. —С. 99-103.

Другие публикации автора по теме диссертации

[1] Казанцев И.Г. Преобразование Радона с послойной сверткой // Тезисы докладов международной конференции "Методы создания, исследования и идентификации математических моделей".— Новосибирск, 2013.— С. 42.

[2] Казанцев И.Г., Пяткин В.П. Использование преобразования Радона в полосе для реконструкции структуры грязевого вулкана // Труды IX Международного научного конгресса и выставки "ИнтерЭкспо Гео-Сибирь-2013" , Т. 1. —Новосибирск, 2013.—С. 187-189.

[3] Казанцев И.Г. Использование комптоновского рассеяния в эмиссионной томографии // Тез. докл. международной конференции "Обратные и некорректные задачи математической физики". — Новосибирск, 2012.— С. 202-203.

[4] Kazantsev I.G., Schmidt S. A fast method for computing a central section of hypercube using gnomonic projection // Российско-немецкий семинар "Распознав. образов и понимание изобр." —Нижн. Новгород, 2011.—С. 111-114.

[5] Яровенко И.П., Казанцев И.Г. Комптоновское рассеяние в рентгеновской и позитронно-эмиссионной томографии // Сибирские электронные математические известия. — 2011. — Т. 8.— С. 172-181.

[6] Казанцев И.Г. Восстановление изображений трехмерной структуры макромолекул по проекциям в криозлектронной микроскопии // Тез. докл. XXIII Российск. конф. по электронной микроскопии. — Черноголовка, Московская область, 2010.— С. 363.

[7] Kazantsev I.G., Herman G.T., Cernetic L. Backprojection-based Reconstruction and Correction for Distance-dependent Defocus in Cryoelectron Microscopy // Proc. of the IEEE International Symp. on Biomed. Imaging. — Paris, France,

[8] Kazantsev I.G., Matej S., Lewitt R.M. Optimal Ordering of Projections using Permutation Matrices and Angles between Projection Subspaces // Electronic Notes in Discrete Mathematics. —2005, —V. 20,—P. 205-216.

[9] Kazantsev I.G., Lemahieu I. Inverse planning for radiotherapy using ridge functions // IEEE Nuclear Science Symposium and Medical Imaging Conference. -Lyon, France, —2000,—V. 3, —P. 1942-1946.

[10] Kazantsev I.G. A New Formula of the Radon Transform Inversion // Proceedings of the IEEE International Conference on Image Processing (ICIP'97).— Santa Barbara, USA.—1997, V. 1. —P. 189-191.

[11] Kazantsev I.G. Tomographic Artefacts Suppression via Backprojection Operator Optimization // Proc. of the IEEE International Conference on Image Processing (ICIP'96).—Lausanne, Switzerland.—1996, V. 1,—P. 749-751.

[12] Kazantsev I.G. An Algebraic Approach to Projection Data Informativity in Computerized Tomography /'/ Zeitschrift für Angewandte Mathematik und Mechanik.—1996, —V. 76, Suppl 3. —P. 469-470.

[13] Казанцев И. Г. Критерий выбора информативной проекции в вычислительной томографии // Препринт ВЦ СО АН СССР №656. — Новосибирск, 1986.

Подписано в печать 06.05.2014 Формат 60x84 1\ 16 Усл. печ. л. 2 Объем 32 стр. Тираж 130 экз. Заказ № 102

Отпечатано Омега Принт 630090, г. Новосибирск, пр. Ак.Лаврентьева,6 email: omegap@yandex.ru