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

РОССИЙСКАЯ АКАДЕШЯ НАУК Сибирское отделение Вычислительный центр

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

Трофимов Олег Евгеньевич

АЛГОРИТМЫ ОБРАЩЕНИЯ И СИНТЕЗА НЕКОТОРЫХ ИНТЕГРАЛЬНЫХ ПРЕОБРАЗОВАНИЙ С ПРИМЕНЕНИЯМИ К ЗАДАЧАМ ТОМОГРАФИИ И ОПТИКИ

Специальность 05.13.16 - применение вычислительной техники и математических методов в научных исследованиях

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

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

46 ОД

- 5 >\ПР '1533

Работа выполнена в Институте автоматики и электрометрии Сибирского отделения Российской Академии Наук.

Социальные оппоненты: академик, доктор (физико-математических

наук, профессор Лаврентьев М. М. , доктор физико-математических наук,

Ведущая организация - Институт электрофизики Уральского

отделения Российской Академии Наук С Екатеринбург}.

Защита состоится " 11 " мая 1993 года в 15 часоь на заседании специализированного совета Д. 002.10.02 по защите диссертаций на соискание * ученой степени доктора наук при Вычислительном центре СО РАН

С 630090, Новосибирск, 90, пр.ак. Лаврентьева, 6).

С диссертацией можно ознакомится в читальном зале библиотеки ВЦ СО РАН по адресу : 630090, Новосибирск,*90, пр.ак. Лаврентьева, 6

профессор Воскобойников Ю. Е.,

доктор физико-математических наук Ерохин Г. Н.

Автореферат разослан

года

Ученый секретарь специализированного совета

ОБЩАЯ ХАРАКТЕРИСТИКА РАбОГЫ

Актуальность темы. Успехи компьютерной томографии хорошо известны. Существуют промышленно выпускаемые томографы различного назначения. Однако проблема снижения дозы облучения при обследованиях по-прежнему остается актуальной. В томографах обычно используется от 100 до 1000 ракурсов и примерно такое же количество отсчетов на каждом направлении. Специфика используемых сигналов и необходимость исследования объектов, меняющих свою форму , например, работающего сердца, делают актуальным построение алгоритмов томографической реконструкции, обеспечивающих заданное разрешение при наименьшем числе проекций. Более полное использование математических фактов, связанных с проблемами томографии, как правило, позволяет повысить качество создаваемых алгоритмов. В работе исследуются некоторые подходы к построению алгоритмов томографической реконструкции, устанавливаются их связи с обычно используемыми, изучается качество получаемых алгоритмов. Компьютерная обработка данных делает актуальной не только разработку алгоритмов нахождения решений соответствующих уравнений, но и разработку эффективных алгоритмов анализа получаемых решений и алгоритмов синтеза объектов с заданными свойствами. Ряд таких алгоритмов, связанных с задачами • выделения связных областей и синтезом киноформных элементов, рассмотрен в работе.

Цель работы. Создание эффективных алгоритмов томографической

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

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

Научная новизна. Новыми являются основные результаты

диссертации: численный алгоритм обращения двумерного преобразования Радона, основанный на свертке с обобщенной'

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

«

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

Апробация работы и публикации. Результаты диссертации докладывались на Всесоюзной конференции "Автоматизация научных исследований на основе применения ЭВМ" (Новосибирск 1979),на 1, 2, 4, 5 Всесоюзных симпозиумах по вычислительной томографии (Новосибирск 1983, Куйбышев 1985, Ташкент 1989, Звенигород Московской области 1991), на международных конференциях по обработке изображений (Лейпциг 1989, Новосибирск 1990), на конференции "Условно-корректные задачи математической физики и анализа" (Новосибирск 1992г).

А

Результаты работы включены в итоговый отчет по теме "Создание дифракционных оптических элементов" Института автоматики и электрометрии С1975-1980). Основные результаты работы опубликованы в [37-57].

Структура работы. Диссертация состоит • иэ введения, четырех глав, объединяющих 11 параграфов, списка литературы.

СОДЕРЖАНИЕ РАБОТЫ

Использование ЭВМ стимулирует постановку и решение задач, в которых исходные данные известны на дискретных решетках. Часто это построение численных алгоритмов для задач, имеющих непрерывные аналоги. Однако во многих случаях числа используются в качестве имен объектов, существенным является лишь то, что исходные данные могут принимать значения иэ некоторого дискретного множества. Математическим аппаратом, естественно возникающим в подобных ситуациях, является аппарат обобщенных функций. В объектах типа ¿-функции хорошо сочетаются свойства непрерывных и дискретных объектов. Однако возможности аппарата обобщенных функций при построении численных алгоритмов не исчерпываются ¿-образными последовательностями. В частности, в задачах томографии существенную роль играет обобщенная функция 1/22

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

Пусть ГС х,у) - плотность объекта в точке Сх,уЭ, а 1С г, 07 -интеграл от ГСх,у) вдоль прямой IX.г,ф), ортогональной лучу, составляющему угол ф с осью ОХ, и отстоящей от начала координат' на расстояние г. В работах [1,39,401 приводятся следующие формулы для восстановления ГСх.у) по 1С г ,ф):

2п

ГСх.уЗ = |эсхсоз0 + У51пф,ф)йф, С13 О

функция 3 определяется следующим образом: ЗСг.фЗ = Кг.фЖ- 1/тгг2) =

■и

Кг-х.фЗ+Нг+х.фЗ-гНг.фЗ ^ ^

х^

здесь * есть свертка обобщенных функций.

В параграфе 1.1 излагается предложенный и обоснованный в работах'[39,40] способ построения алгоритмов восстановления плотности Г(х,уЗ по исходным данным, известным на дискретной решетке, на основании формул типа С13 -и С23. Существо метода заключается в том, что для функции 1Сг,03 при фиксированном ф строится такая аппроксимация, что интеграл по х в С23 сходится, это выполняется, в частности, если аппроксимировать функцию К г, ф) кубическими сплайнами. При этом можно получить явные формулы для аппроксимации функции Я г, ф) [40,41].

Обобщенные функции являются функционалами над пространством

бесконечно дифференцируемых быстро убывающих функций. Однако при

построении аппроксимаций исходных реальных данных по отсчетам,

заданным в дискретных точках, желательно иметь менее жесткие

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

2

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

2 1 1/г заведомо определена и может быть представлена

формулой С 23. Для формулировки соответствующего предложения

нам потребуются обозначения:

Ях,д) =

дСх+у)-д(х-у)-2дСх)±^

2

У

ЗМСх.дЭ = [1^х+у)-д(х-у)-2д(х)1

} у

Теорема 1.1.1. Еоли функция дС х) имеет финитный носитель и такова, что для любой функции р из пространства основных (бесконечно дифференцируемых, быстро убывающих) функция £М

задает функционал С то есть сходится интеграл |зМС х, д) р( х) с1х),

1

то ЗСх,д) = —* дСх).

х

То есть для любой функции р справедливо равенство

1

Ях.дЭрСхМх = С13С х, д) , р( х) ) = С С-* дСх)),рСх)).

х

Условия теоремы выполняются, если дСх) имеет интегрируемый модуль второй производной. Представляет интерес также преобразование Фурье от Б. Так как свертке соответствует умножение в плоскости Фурье, то ЭХ) = |Х|дСХ), здесь ЭСХ) и дС X) - преобразования Фурье от ЗСХ) и д(Х). Бели дСх) удовлетворяет условиям теоремы, то преобразование Фурье существует в обычном смысле и мы приходим к следующему равенству между регулярными функциями

||Х|асХ)е"1Х^х = | дС2^-д(2-у)-2д(2)^

Последнее равенство можно сопоставить с тем, что умножение на X в области Фурье эквивалентно дифференцированию исходных

функций. Преобразование Фурье от X есть 6 , свертка с <5 дает

производную, преобразование Фурье от |х| есть 1/у , и мы приходим к соответствующей свертке.

I При переходе к дискретному варианту в работе предполагается что fCx.y) = 0 вне круга радиуса R с центром в нуле. ' Исходными данными является величины ICr^pj) , здесь г.- отсчеты в интервале [-R.R], 1 < i < М, р^- отсчеты в интервале [0,п], 1 < j < N. Если теперь при заданных значениях функции 1Сг,р) в отсчетах Cr^.fO построить аппроксимацию IС г, р) так, что для SCz,p) вшолняется равенство CID, то, используя CID и С23, можно получить приближение к fCx,y). При изложении существа способа будем предполагать, что отсчеты на осях г и р являются равноотстоящими.

При каждом фиксированном pj функция Krопределяется следующим образом.

1. Функция ICr.pj!) имеет непрерывную первую производную по г.

2. В .узлах решетки аппроксимирующая функция совпадает с заданными отсчетами, а ■ ее производная в этих точках равна

выборочной. То есть, справедливы равенства: ICr^fO = ICr^ ,<fO,

I Cr.,р.) = С 1С г. , ,p.)-iCr. . ,p.)3/2h , где h = 2R/CM-D, г i j 1+1 j «l-i j

Kr0'V =ICrM+rV =0.'i=l....M.

3. На интервале trifгфункция Ifr.pj) есть полином третьей степени от г.

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

М

ICr.pJ = ^Г ICr. ,pj)QCr+R+Cl-i)h3, i=l

где

QCxJ =

C3/2h3)x3 - C5/2h2)x2 + 1, " если О ч< x ч< h , C-l/Sh^x^ СЭ/гЬ^х2 -C4/Wx^+2, если h < хч< 2h , QCx) = QC-x) , QCxD = 0 при |x| > 2h.

Функция ОСхЗ имеет раэрьвы второй производной, но модуль второй производной интегрируем. Из теоремы 1.1.1 следует, что свертка

ЭфСг) = (ХгЖ -1/гг2^) выражается формулой С 2). Непосредственными вычислениями можно получить

=. С -1 /гг2} С 2С Зг+2Ю Сг+Ю1п (+2С Зг-2Ь) С г-Ю 1 п | --9221п | -С Зг+4Ь) С 2+2Ю 1п|г+2Ь|-С Зг-41г) С г-2Ю 1п\z-Zh |. Таким образом, 4

М

ЗСг.р^) = ГСг.рЖ-Х/пг2) ^ГС^ ,р030Сг+1г+С1-1)Ю.

1=1

Г Сх,у) - приближение к функции Пх,у):

Заменяя в С1Э Б на ¡3 и интеграл частной суммой, получаем

нкц N

Ах,у) =С1/2лСН-1)) ^ §Схсогр. +/31 пр., р.). СЗ) 1=1

В работе приведены результаты численных экспериментов по восстановлению функций указанным способом. Построенные таким образом алгоритмы позволяют уверенно определять круговые отверстия при ¿/бг = 0.7 С с± - диаметр отверстия, бг - величина шага дискретизации по г). Сетка по г и р равномерная, количество направлений 18 С<5р.равно 10 градусам).

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

Утверждение 1.1.2. Пусть

1) { х) > - последовательность функций класса С1,

2) { I Сх) > сходится равномерно к 1Сх) в С1,

3) существуют А, 6, такие, что

|1 Сх) -I СО) -х1 Сх) I < Ах2, для всех |х| < 6.

'_п п п 1 11

Сп= 0, 1... ; 1фСх) = I Сх)). Тогда

lim S = lim n

ao

г I СхЗ + I c-x) - 21 СхЗ n__n n ,

-2- dx =

X

,. Г I Cx) + I C-x3 - 21 (x) , lim -g- dx = S.

J x

0

Условие 3 заведомо выполняется, если функции КхЗ и х) имеют ограниченные вторые производные. Из Утверждения 1.1.2 следует Утверждениие 1.1.3. Если

13 € С*, по переменной х при любом <р, п=0,1...,

23 { IпС х,рЭ > сходится равномерно к Кх,р) = I^Сх,<р) , .33 существуют А и 6, не зависящие от z и <р, такие, что

ffljz.p)

|lnCz+x,p3 - In(z,p3 - х—^- | < Ах? п= 0,1... ,

для всех гири для всех |х| < <5, то последовательность обратных преобразований Радона функций I^С х,рЗ сходится поточечно к обратному преобразованию Радона функции Кх,рЗ.

Если функции In(z,p3 и Кг.рЗ финитны, то сходимость будет-равномерной. Сформулированные утверждения позволяют получить достаточные условия сходимости приближенных решений к точному в случае, когда исходные данные известны на дискретной сетке, а шаги сетки по г и <р 'стремятся к нулю. Для этого достаточно получить условия, которым должна удоволетворять функция I(z3, чтобы последовательность аппроксимирующих кусочно-полиномиальных функций удовлетворяла условиям утверждений. В точке z значение функции ICz3 заменяется значением полинома

3 2

PCz,h3 = а, £z,h3xz +a_Cz,h3xz +a_Cz,h3xz +a.Cz,h3. 1 с. 3 4

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

a„Cz,h3=ICz3-IC1)Cz3zKIC2\z3/23z2+CIC33Cz3/23z3+OCh3, 4

с п с ?■) с г

а3Сг,Ы= I ^СгЗ-Г^СгЗг +СГлЧ2)3/2)2 +ОСЫ,

а^.М* 1С23С23/2 -а^СгЗЗ/Эг +ОСМ,

а^.Ь^ 1СЗ)Сг)/2 +ОСМ.

Иэ этих выражений следует, что РСг.М = Кг) +ОСЮ, причем ОСЮ, величины порядка Ь , вообще говоря, зависят от г .Но, при условии ограниченности четвертой производной от 1С г), стремление ОСЮ к нулю, при Ь —> 0, будет равномерным По г.

Боли функция 1Сг,<р) такова, что ■ 4

тах|—-——| < М, то функции, получаемые на основании

д 7.

предложенного алгоритма, сходятся к точному решению при бг и бр стремящихся к нулю , бг и бр - шаги сетки, в которой известны исходные данные. Условия сходимости сформулированы нами в терминах функции Кг.р), то есть преобразования Радона искомой функции ГСх,у). Так как интегрирование не уменьшает порядка гладкости, то условия сходимости будут выполняться, если равномерно ограничена четвертая производная при дифференцировании функции Кх.у) вдоль любой прямой. Численные эксперименты показывают, что при конечных <5г и бр можно получать хорошие приближения даже в случае разрывных искомых функций.

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

5С г) можно сделать следующее замечание. Рассмотрим процесс 7?,

2

являющийся сверткой с 1/г процесса Спектральная плотность

этого линейного преобразования есть |\|. Дня спектральных

плотностей процессов ? и 7? получаем соотношение р

Г СЗО = |\| Дисперсия процесса г? конечна, если

интегрируема Т ХЗ, то есть процесс ? дифференцируем в среднеквадратическом. ■ Для того, чтобы свертка выражалась формулой С23, на.процесс 2; нужно наложить дополнительные условия, потребовав, например, чтобы выборочные функции с вероятностью единица имели конечную вторую производную. Это условие будет выполняться, если будет конечным шестой момент ^(ХЗ. То есть спектральная плотность шума должна достаточно быстро убывать на бесконечности.

При численном решении задач томографии часто используется

метод свертки и обратного проецирования [11]. Его связь с

изложенным выше методом заключается в следующем. Свертку с 2

обобщенной функцией 1/г можно заменить дифференцированием и сверткой с 1/г (преобразованием ГильбертаЗ. Таким образом, функцию 2 можно представить в виде

ЗСг.фЗ = Кг,.ф)*(.- 1 /кг. 3. Вместо обобщенной функции 1/г обычно рассматривается последовательность регулярных функций рдС г), сходящаяся к 1/г (здесь А - параметр регуйяризацииЗ. Используя интегрирование по частям, дифференцирование переносят на функцию рд(гЗ и таким

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

2 • 2 сходящуюся к 1/г . То есть свертка с обобщенной функцией 1/г

заменяется последовательностью сверток с регулярными функциями

рд(гЗ. Явное использование' обобщенных функций позволяет

целенаправленно строить соответствующие регуляризирующие

последовательности, учитывать априорную информацию, использовать

аналитические формулы, обращения. Последнее особенно важно в

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

большинство имеющихся здесь формул обращения записывается в

терминах обобщенных функций.

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

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

В работе [ 6 ] получена формула обращения для функций, имеющих финитный носитель,' и для кривых, удовлетворяющих определенным условиям. Главным в этих условиях является то, что любая плоскость, пересекающая объект, пересекает кривую, по которой движется источник. Примером кривой, удоволетворяющей условиям работы .[6], является совокупность двух единичных окружностей, лежащих во взаимно перпендикулярных плоскостях. Другим примером является винтовая линия, охватывавшая .объект цилиндрической формы. Однако построение численных алгоритмов непосредствено на основании формулы, приведенной в [6], затруднительно Сем. [11, ст.201]). Дело, в частности, в том, что формула обращения основана на преобразовании Фурье от однородной функции, получаемой из исходных данных. Причем преобразование Фурье понимается в смысле обобщенных функций, а преобразование Фурье в обычном смысле может не существовать.

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

Пусть заданы функция ГСх3 = Ях^х^хд) и точка г = Для интеграла от функции ГСх) вдоль прямой, проходящей через точку г в направлении вектора а, будем использовать обозначения из [73:

оо

В некоторых ситуациях используется функция

С Е^ПСг.оО^СгНеОсИ, О

отличающаяся от С К^ПСг.а) тем, что интегрирование ведется лишь при I > 0. Пусть задана кривая , по которой движется источник, ФСХЗ = СФ^ХЗ ,Ф,,СХЗ .ФдСХЗ,... ,ФСХЗЗ, параметр X пробегает некоторый интервал Л действительной прямой. Для любого а = Са^.е^.Од,... и X € Л определим функции

оо оо

дСа.ХЗ^ГЗСФСХЗ.аЗ = |гС ФС ХЗ ИаЗ сИ = |гС«ХЭ+1-эд-0<И. .

- оо О

оо оо

д+С а, ХЗ =С С ФС ХЗ , аЭ =£С ФС ХЗ НаЗ сИ=-|-|—ФС X) И-щ-З сИ. С 4)

а о

Функция дСа,ХЗ есть интеграл от функции ГСх) вдоль прямой, проходящей через точку ФСХЗ в направлении вектора а.

Для функций, имеющих финитный носитель в трехмерном пространстве, в [61 получена формула обращения: 2 п п/2

г г 1 + ПхЗ = С05Р----& ¿<рйв. С Я

О -п/2 2H.ClCM.fD

Функция в С^?,Х) есть преобразование Фурье от функции д+Са,ХЗ по переменной а, (1=(со2всовр,51пвсо$<р,£1пр), X = ХСх,/В такое, что скалярное произведение С^?,хЭ равно (.р, ФСХЗ) и с/?,ф'сл)3 не равно нулю. Предполагается, что такое X существует для любого х, принадлежащего носителю функции КхЭ, т.е. любая плоскость, пересекающая носитель функции, пересекает кривую ФСХЗ так, что'знаменатель в С5) не обращается в нуль. Примером такой кривой является совокупность двух единичных окружностей, лежащих во взаимно перпендикулярных плоскостях, если носитель лежит в единичном шаре. Построение численных алгоритмов непосредственно на основании формулы, приведенной в [6], затруднительно Сем.[11,

стр. 201 D. Дело в том, что преобразование Фурье от? функции g+(a,X3, понимаемое в обычном смысле, не существует, так как на бесконечности д+СаДЗ имеет порядок 1/|а|. Это связано с переходом от исходных данных, заданных на поверхности |а|=1, к однородной функции, заданной во всем пространстве. В [6] преобразование Фурье понимается в смысле обобщенных функций.

Чтобы использовать формулы типа С 5) для построения алгоритмов, желательно иметь выражения для функций G ,G и f через регулярные функции g и д. В параграфе доказана

Теорема 1.2.1: Пусть G (Ç3 есть преобразование Фурье в смысле обобщенных функций от однородной функции g+Coû, тогда

G+(Ç3 = С-1 /Ап2)

=1

d гр. С 63

Для действительных функций КхЗ в формуле С 53 нужна мнимая часть в+С ? 3.

1пй+С?3= С1/4гй| д+Ср)<5 CC0.?ЗЭd2p. №1=1

Используя обобщенные функции, сосредоточенные на поверхности, [12, стр. 259], получаем из Теоремы 1.2.1 следующее Следствие:

1шв+С?3= С1/4ттЗ I ЬС?,БЗд'С)0ф'.

k дг

)| Щ.БЗд* Зф

Здесь а?3 = <у 6 Б?! С?,=03, Щ,БЗ = ) ?

¡А

производная по направлению Подставляя в С53 функции д С?ДЗ и 6+(?,ХЗ, зависящие от параметра получаем формулу обращения, пригодную для построения численных алгоритмов.

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

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

Замечание. А. С. Денисском независимо и другим методом, без явного использования преобразования Фурье обобщенных функций, получены формулы обращения функции д+ в Rn. Иэ Теоремы 1.2.1, в частности, следует, что при п = 3 формулы A.C. Денисюка и формулы, получаемые изложенным способом иэ формулы Туя, совпадают. В. П. Паламодов и А. С. Денисюк сообщили автору настоящей работы, что при п = 3 аналогичные формулы получены также П. Гранжа С163 и Д. Финчем С результат доложен на школе по томографии в Обервольфахе 1990г). Здесь следует отметить, что' наличие различных видов формул обращения расширяет возможности при построении численных алгоритмов, так как формулы, совпадающие в предельном непрерывном варианте, могут давать существенно различающиеся результаты на дискретных решетках данных. Различные методы получения формул обращения, как правило, дают возможность учитывать различные типы априорной информации о решении и шумах.

Приведенные выше формулы позволяют вычислять преобразование Фурье лучевых данных по их значениям на. единичной сфере.

Представляет также интерес связь функций GCf.AJ и G+C?,X) с

V

исходной функцией fCx) и ее преобразованием Радона fC?,r).

Формулы Восстановления функции fCx) по функции (К^ПСФСХ) ,о0 получены в работе С7], при этом существенным образом используется введеная там функция

CSD?n = [—--^----„dr. х J Сг-Сх,?)) 1

Иэ результатов, полученных в параграфе 1.2, следует, что

(S2f)C?) = CK^fXz,?), то есть функция CS ПС?) является преобразованием Фурье в смысле обобщенных функций исходных лучевых данных.

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

переменных по интегралам вдоль пространственных прямых, если

используется формула С 6). Для этого ее нужно преобразовать к

виду, удобному для регуляризации. Боли ? = СО,0,1), то

2

знаменатель в формуле С 6) равен Csin0). При произвольном £ нужно перейти к сферической системе координат, в которой полюсом является точка Через ,ф,в,Х) будем обозначать функцию д(.ф,в,\) в такой системе координат, а через -2тгОС?,фД) ее интеграл по 9. В указанных обозначениях формула (6) приводится к виду

гг/2

icos^i ОС? ,р,\) -g-dp. С 7)

Csinp)

0 2

Используя каноническую регуляризацию функции 1/Csinp) ,

построенную в С12]С стр 90), получаем гг/2

¡•г cos<£QC?,p,X) QC? ,0,X3

GC?,X) = -g---=-dp. С 8)

JL Csinp) с <pc J

0

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

При выводе формул обращения в работах [5-6] существенным образом используется то, что любая плоскость, пересекающая носитель функции, пересекает кривую, по которой движется источник. Однако, как следует из работ [9,10], для функций, имеющих финитный носитель, это условие не является необходимым.

Боли носитель функции содержится в единичном шаре, то для ее восстановления достаточно, чтобы источник двигался по единичной окружности с центром в начале координат, лежащей в плоскости z = О. В [10] приводятся формулы, которые могут быть использованы для построения численных алгоритмов.' Причем задача разделяется на 2 части:

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

В парарафе 1.3 изложены численные алгоритмы для первой части задачи . Как показано в Г14), задача в целом является существенно неустойчивой, однако для первой части оказалось возможным разработать численные алгоритмы, дающие хорошую точность восстановления в модельных экспериментах [37). Как отмечено в [15), зная интегралы вдоль прямых, проходящих через два взаимно перпендикулярных единичных круга, лежащих в плоскостях 2 = 0 и у = 0, можно восстановить итегралы вдоль прямых, лежащих в любой плоскости, перпендикулярной оси ОХ и пересекающей объект. Этот прием дает еще один класс численных алгоритмов восстановления функции в условиях Кириллова-Туя. Здесь важно то, что задача сводится к хорошо разработанному двумерному случаю. Результаты соответствующих модельных экспериментов приведены в [15,42).

В парарафе 1.4 изложена структура способов построения «

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

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

При переходе от непрерывного варианта задач обработки данных к дискретному и выборе шага дискретизации часто используется теорема Котельникова для функций с ограниченным . спектром [11,18,19,211. Классический вариант теоремы отсчетов относится к равномерной сетке в декартовых координатах. Однако в ряде задач, в частности, в задачах томографии встречаются другие •решетки задания исходных данных [18,19,37,38]. В параграфе

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

Пусть ИЗ - класс функций двух переменных, Фурье-спектр которых содержится в круге с радиусом К и центром в нуле. Можно ли указать шаги дискретизации -1бг,бр) в полярных координатах г и <р так, что любая функция из класса ИЗ может быть восстановлена по своим значениям на соответствующей решетке? Можно построить примеры, показывающие, что ответ на этот вопрос отрицательный 119,37,38]. В [37] соответствующие примеры строятся следующим образом.

Пусть ГС г) - функция двух переменных, здесь 2=ге'1®Р - точка двумерной плоскости. Если Кг) = ГСгЭе*^, где к целое число, то

•V <\,

ГСЮ - двумерное преобразование Фурье от ГС г) - имеет вид ГС и 1к0 ' 1к0

=РСр)е , \ = ре - точка в плоскости Фурье [20, стр. 155]. Функция РСр) есть преобразование Бесселя к-ого порядка от Кг):

оо сх>

РСр) = 2лС-1)к

Кг) ^С2лгр)пЗг = 2гг1к

ГСгП_кСггггр)п±-, С 9)

здесь

2гт

Л"

т

^СУ = С1/2л) е ^¿¡д _ функция Бесселя к-ого

° 1кб порядка. При обратном двумерном преобразовании Фурье от РСр)е

функция ГС г) вычисляется по формулам, аналогичным С9). Из равенства -Г^СI) = О следует, что преобразование Фурье от ГС г) = ГСг)е ^имеет вид РСр)е с той же функцией Пр). Отсюда следует, что функции вида ГСг)с02кр при двумерном преобразовании Фурье переходят в функции РСр)созкр, а функций Иг^пкр - в РСр^пкр. При заданном К выберем РфС р) - функцию

одного переменного, принадлежащую Ь2 , такую, что ^Ср)= 0 при

|р| > К. Рассмотрим радиальную функцию двух переменных =

ре ] - р]. Возьмем от нее двумерное обратное

преобразование Фурье, получим радиальную функцию =

16

^02*"Ге = ^о*" Г">' ФУНКЦ.ИИ виДа ГфСгЭсозкр будут иметь преобразование Фурье, равное РдСр)соэкб, но Р^Ср) =0 при |р| > Б, следовательно, функции Г^СгЗсозкр будут иметь спектр, содержащийся в круге радиуса К. Для любого к функция тЗсоэкр равна нулю на лучах вида кр = л/2+2лп, расстояние между ними по углу р равно Ар = 2гс/к, но функция Г^СгЗсозкр отлична от тождественного нуля. Построенные примеры показывают, что финитности носителя спектра Фурье недостаточно для восстановления функции по ее значениям на равномерной решетке в полярных координатах.

В рассматриваемом параграфе выделены классы функций КШ, для которых справедлива теорема отсчетов в полярных координатах. Функция принадлежит классу КШ, если удовлетворяет условиям следующего Утверждения:

Если существует конечное N такое,.что N

£(г,<р) = У^Сг^е"'"^, и функция ТСгимеет ограниченный носитель двумерного преобразования Фурье, то

1) каждая из функций ^СгЗе*^ имеет ограниченный носитель двумерного преобразования Фурье;

2) каждая из функций Г г) имеет ограниченный носитель одномерного преобразования Фурье;

3) существует равномерная сетка в полярных координатах такая, что функция Т(г,р) может быть восстановлена по своим значениям на этой сетке (Ар = 2тг/(2Ы+1) , Аг = л/Ю. Для действительных функций Ар = 2я/С N+15.

Эта часть результата несколько другими методами получена в работе Старка [34].опубликованной в Журнале Американского Оптического Общества. Отметим, что при к, изменяющихся в

бесконечных пределах, в указанном виде может быть представлена любая функция из КБ 120, стр. 1553.

Построенные выше примеры показывают также, что ограниченности спектра Фурье недостаточно для восстановления функции по значениям на равномерной сетке ее преобразования Радона. Действительно, функция FCzD вида fCrOcoshp при преобразовании Фурье перейдет в функцию вида FCp)coskp. В силу соотношений между преобразованиями Фурье и Радона Cl,стр.19], И FC z),s,р) - преобразование Радона функции FCz)- при каждом фиксированном ç может быть получено с помощью обратного одномерного преобразования Фурье по р от FCp)coskp и, следовательно, будет равно нулю при coskp = О. Выбирая Кр)=0 при р > К, получаем нужные примеры. Функции из классов KKN могут быть восстановлены по значениям преобразования Радона на равномерной сетке. Действительно, преобразование Радона таких функций будет иметь вид

так как оно может быть получено из двумерного преобразования Фурье обратным одномерным на каждом луче. Бзли Др=2гг/( 2N+1 ), то для каждой точки s по величинам PC К z),s,ç^) могут быть найдены функции g^Cs). Таким образом, можно считать известной каждую функцию gfcCs) в дискретной последовательности точек. Спектр Фурье этих функций имеет финитный носитель, следовательно, может быть применена теорема отсчетов. В задачах томографии искомые функции обычно имеют ограниченный носитель, спектр Фурье таких функций не может содержаться в круге конечного радиуса.' Тем не менее, если основная часть спектра сосредоточена в круге радиуса S, то говорят, что эффективная ширина спектра равна" К fil, стр. 66]. Результаты типа теоремы отсчетов используются при этом для получения оценок размеров разрешаемых деталей.

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

N

k=-N

71/2

С2/гО ^кцпб) с!б = Кх) СЮ) О

наэьвается уравнением Шлемильха [22, стр: 323].Заменой переменных Х2з.пб = г оно может быть приведено к виду

х

С 2/гО I—----dz =Ях). сю')

J Гг "г~ о у х " 2

К решении уравнения СЮ) и близких к нему уравнений

•сводится задача идентификации ряда нелинейных элементов

радиоэлектронных, оптических и биологических систем [23-25,37].

Уравнение СЮ) используется также в задачах интегральной

геометрии, в частности, в [17, стр. 22] при доказательстве

теоремы о носителе. Для уравнения Шлемильха известно следующее

утверждение [22, стр. 323]: если функция Ях) непрерывно

дифференцируема, то уравнение СЮ) имеет не более чем одно

непрерывно дифференцируемое решение, которое задается формулой

п/2

= «0)+^

у/Сх) = «0)+^ СХ51П0^0. СИ)

О

Уравнение СЮ) желательно рассмотреть в более широком классе функций, чем непрерывно дифферецируемые. Бзли мы ограничивемся непрерывно дифференцируемыми искомыми функциями, то не можем идентифицировать нелинейные зависимости, например, вида

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

Использовав замену Х51П0=г, от формулы СИ) можно перейти к формуле

х

вСх) = «0) + х —--------dz СИ')

J /~2 2~

О

Для существования интеграла в СИ') необходима

интегрируемость f (хЗ. Чтобы при вычислении по формуле СИ') разным f(x3 отвечали разные ytxl, функция f(x3 не должна содержать составляющей, которая исчезает при дифференцировании.

Естественным классом правых частей, удоволетворяющих перечисленным условиям, является класс абсолютно непрерывных функций.

В параграфе 2.2 доказана Теорема II.2.2. Бэли f(x3 абсолютно непрерывна на [О, х], то функция уКх), задаваемая правой частью формулы (11'3, определена для почти всех х и является решением уравнения (103.

Доказано также

Предложение II. 2.2. Боли f(x3 задается формулой (103, в которой V<x3 является суммой абсолютно непрерывной функции и функции с конечным числом скачков, то для всех х

п/2

уСх-03 = xjf'(xsin63dS + f(+03, О

где уСх-03 - предел слева функции уКх) в точке х, f(+03 предел справа функции f(x3 в точке х = О.

Отметим, что если не имеет скачка в точке х=0, то f(x3 абсолютно непрерывна и fС +03 = f(03.

В формулы обращения (ИЗ и (11'3 входит производная правой части. При неточна заданных исходных данных это может приводить к неустойчивости решения при малых абсолютных значениях погрешности. В подобны-; ситуациях применяются обычно методы регуляризации [2,3]. Один из способов применения этих методов к уравнению (103 заключается в следующем. Рассмотрим интеграл от искомой функции

X

«хЗ = Jy/(t3dt;

О

учитывая формулу обращения (11'3, получаем

■в

,т, , , zf(z)dz ФСх) = - = х

и/г

f(xsin0)sin0dS.

j .г 2 о у 1 "z о

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

В [25] задача идентификации операторов

I

ГСр) = V(p] + QC jpCridr) + -^-SCp) О

и I TCf) = fjCkj* f2Ck2*p)

сводится к уравнениям несколько более общего вида:

л/2

fCx) = С2/тг) jyCxsir>0)cosCпб)d6 ; n = 0,1,2; С12) О и/2

f(x) = С2/я)Jy<xsiпб)sinCпб)d0 ; n = 1,2,3. (13). 0

Формулы обращения для уравнений (12) и (13) получаются на основании соответствующих формул для 'уравнения (10) [43].

Отметим также, что в [26] в классе абсолютно непрерывных функций исследуется уравнение Абеля, близкое к уравнению Шлемильха.

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

2 Ч

требований к точности выводов, возникает необходимость

использования ЭВМ не только на этапе получения решения, но и на

этапе анализа решения. Так в задачах томографии часто инородныые

включения имеют плотность, сильно отличающуюся от плотности

основного объекта. Необходимо уметь выделять эти включения.

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

два значения, например, 0 и 1. С учетом дискретизации по

пространственным координатам, приходим к задаче выделения

компонент связности на матрице, состоящей из нулей и единиц.

Понятие выделения объекта в разных задачах уточняется

по-разному. В частности, может быть необходимым подсчитать

количество объектов. Задача определения связности объектов С один

или не один) представляет и теоретический интерес как пример

задачи, не являющейся локальной. В работе 1271 подробно

исследуется вопрос о сложности этой задачи при реализации на

персептронах, в работе [28] изучается вопрос сложности логических

схем, необходимых для подсчета объектов. Существует ряд

алгоритмов, подсчитывающих объекты при различных ограничениях

(выпуклость, отсутствие дыр, и т. д.) [29].

В параграфе исследуется задача минимизации памяти при

подсчете числа плоских объектов произвольной формы, заданных на

матрице из п строк и m столбцов, состоящей из нулей и единиц.

Предполагается, что матрица предъявляется поэлементно, причем

сначала подаются слева направо элементы первой строки , затем

второй и т. д. . После предъявления последнего элемента матрицы

в счетчике будет находиться число, равное числу объектов.

Введем некоторые определения. Пусть задана матрица,

элементы которой принимают два значения Oui. Будем говорить,

что элементы матрицы а. ,= а(к,1) и а = aCp.q) принадлежат К, 1 р, CJ

одному объекту, если существует последовательность индексов (ij.jj),... £is.js). такая, что i1= k, ^=1, р, js=q,

ivW+|vW <"2 C1- к - S3' aCVV = !"

Последовательность элементов aCi^.j^), называется путем,

соединяющим точки а, , и а . Последовательность вида k,l p,q

а. ,,... ,а.. а. , . ,...а. . будем называть сечением Б. , .

.Ьк , к J-l.ni .ьк

В параграфе построен алгоритм с объемом памяти

2т+2109^111+109^+1. Другие известные алгоритмы подсчета объектов произвольной формы используют память порядка пйод^К, где К -число объектов [291. Минимизация памяти в предлагаемом алгоритме достигается благодаря специальному способу кодирования связей, суть которого заключается в следующем. Имется алфавит из четырех символов: "С , } ,1 , 0". Отрезкам в том порядке, в каком они появляется в сечении, присваивается один из перечисленных символов по следующему правилу: отрезку, не имеющему связей с другими в предъявленной части матрицы, присваивается символ "о", первому в группе связанных отрезков присваивается символ "(", последнему -")", остальным отрезкам присваивается символ "1". Получаемое при этом слово, характеризующие текущие связи, напоминает алгебраическое выражение и позволяет экономно хранить необходимую промежуточную информацию. Работа алгоритма состоит в перестраивании слова, характеризующего связность объектов, при поступлении очередного элемента матрицы.

Предложенный способ кодирования может быть использован при вычислении других топологических характеристик изображения. Например, таких, как наличие "дыр" в объектах, степень вложеннности объектов друг в друга односвязность или

многосвязность соответствующих областей.

Глава 3 посвящена методам расчета киноформных элементов.

Использование ЭВМ для управления лазерным лучом позволяет создавать устройства для получения объектов, которые с точки зрения геометрии являются ступенчатыми функциями с размерами высоты ступеньки порядка долей длины волны света в видимом диапазоне. Если подобные объекты изготовлены из прозрачного материала, то они могут быть использованы в качестве оптических элементов. Оптические элементы такого рода принято называть киноформными оптическими элементами,или киноформами [30,31]. Термин "киноформ" впервые появился в [30] применительно к одному классу подобных объектов, о котором более подробно будет сказано ниже. Однако в настоящее время этот термин используется для

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

Одна из причин интереса к киноформным элементам заключается в том, что при прохождении оптических лучей размеры, кратные длине волны, не влияют на формируемое изображение , что позволяет сделать оптические элеметы практически плоскими. Очень интересной является также возможность обеспечивать высокие точности по координатам, лежащим в плоскости, перпендикулярной оптической оси элемента. Использование ЭВМ позволяет создавать элементы, толщина которых z сложным образом зависит от координат х и у. Расчет и создание киноформных элементов является сложной вычислительной и технологической задачей. Учет специфики подобных элементов и небходимость изучения их возможностей делают целесообразной постановку • и исследование ряда математических задач [31,47-49]. Некоторые такие задачи рассмотрены в работе.

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

В [30] изложен способ создания объекта, названного там киноформом. Объект представляет собой прозрачную пластину, дающую при освещении монохроматическим светом изображение, близкое к заданному в зоне наблюдения, для которой справедливо приближение Френеля. Разность между толщиной пластины в различных точках может не превосходить длины волны падающего на пластину света. Изображением, использовавшимся в [30], было слово " kinoform" . С математической точки эрения способ создания объекта заключается в следующем [31,49].

Пусть 1С ¡1,17} - заданная неотрицательная функция, изображение которой мы хотим получить. Рассмотрим функцию

ЦСи,у)=С1АсО || 11/2(?,Г1)ехрС1(пА(и[Си-532-Ку-т))2+рС5,т)3])с1?с1л,

здесь X - длина волны, <3 - расстояние до зоны наблюдения Френеля, рС^.тр - случайная функция, удовлетворяющая условию

МСехрС1рС?1,7]1)-1рС?2,772)))= Л?1-?2,г?1-»?2). Пусть £5Си,у) - фаза функции 1Ки,у) по модулю 2п. Езли взять одну из реализаций случайной функции ?, тр, мы получим одну из реализаций случайной функции ВС и,v). Теперь можно изготовить фазосдвигаюшую пластинку, соответствующую функции 0Си,у), за счет изменения толщины пластинки в пределах длины, волны. В [30] утверждается, что если такую пластинку осветить монохроматическим светом с длинной волны X, то на расстоянии с! будет поле, квадрат амплитуды которого близок к 1С?,г?). Там же обсуждаются результаты подобных экспериментов для изображения, являющегося словом "1с1ПоГогш" .

В [31] приводится обоснование этого факта, которое заключается в следующем. Рассматривается случайная функция УСх,у), являющаяся преобразованием Френеля от случайной функции 0Си,у), то есть

УС х, у) = С1 /ХсЮ || ехрС 1С лЛсЮ [ С и-х) 2+С у-у) 2+6С и, у) ]) с!ибу.

• к —

Пусть 1Сх,у) = МСУСх.у) V Сх.уЗ), то есть Кх,у)

-математическое. ожидание квадрата модуля случайной функции

\Кх,у). В [31] показано, что функция 1Сх,у) может быть

представлена в виде ряда

оа

Кх,у) = £ скГСх,у), к=1

где ^(х,у) = 1Сх,у), а коэффициент с^ составляет 78% от суммы коэф$ицинтов с^. Этот факт является обоснованием изложенного выше метода построения кинофэрмных элементов. Таким образом, с помощью метода, предложенного в [30], можно получать хорошие приближения к заданным изображениям, но он обладает систематической погрешностью.

В общем виде задача может быть поставлена следующим

образом.

Пусть задана неотрицательная.функция Их,у). Нужно найти фазу бСи,у) такую, что квадрат модуля фунции УСх,у), где УСх,'у) =

•Я-

С1/Лс1) II ехр( 1С я Ас!) Г С и-х) 2+( у-у) 2+0С и, у)]) бис!у,

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

дСи,у ) =

е , при а^ и < а2< у < Ь2>

О в противном случае.

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

Преобразования Френеля и Фурье входят в класс преобразований, представимых в виде

ГСхЭ = | ?СйОдСы,х)бы;

Б

для преобразования Фурье дСсо.х) = е1ых, для преобразования

2 2

Френеля д(ш,х) = МЛсО ехрС IС п/Хд) Г С ы^ -х^ ) -Кы^-у^) ]), здесь х =С х^, х,,), ш = С , ш^). В практических ситуациях областью интегрирования В является некоторое ограниченное множество, которое для'дальнейшего, не ограничивая, общности,можно считать квадратом. Модуль функции ?Сш) в общем случае отличен от единицы. Нам нужно найти другую функцию с модулем, равным

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

Суть предлагаемого способа удобно изложить на примере преобразования Фурье в непрерывном случае , а затем рассмотреть

^

дискретный вариант в общем случае. Пусть

ГСх) = [еКХ'х:> ?СХт=]е1а'х3 |?СХ) |е^с1Х,

В Б

и пусть (ЗСх) - преобразование Фурье от фазовой функции вида

где рСХ) случайная функция, то есть

(ЗСх) = |ехрСКСХ,х)+рСХ)))с1Х,

Б ■ ■

где X = СХ^.Х^, с1Х = (ЗХ^сЛ^, х = Сх^х,,), (Х,х) = (.Х^+Х^х^),

В = <Х: 0 < Х1 <Х0, 0 < Х2 <Х0>.

Для математического ожидания случайной функции (ЗС х)

получаем

МйСх) = [е1СХ'х3Ме^Х^Х.

О

Если теперь выбрать случайную величину рСХ) такую, что

?СХ), то получим равенство МОСх) = ГСх). Боли выполняется условие О < |ЯХ)| <1, то соответствующую величину всегда можно подобрать. Действительно, Ме^^ есть значение характеристической функции случайной величины рСХ) в точке и=1, Характеристическая функция дСи) гауссовской случайной величины т? имеет вид

дСи) = ехрС 1 аи-С 1 /2)сг^и2), 2

где а - математическое ожидание г;, а - ее дисперсия. Выбирая

2

случайную величину рСХ) так, что НрСХ) = ц<\") истСрСХ)) = = - 21 пС |?СХ) |), получаем нужное равенство. Здесь ?СХ) = |ГСХ)|е*^Х\ Распределение случайной величины р(Х) может быть и другим, например, равномерным.

Изложенный выше прием обеспечивает равенство математического ожидания случайной функции вС х) и заданной функции ГСх). Для того, чтобы это равенство приближенно выполнялось на одной реализации, необходима малость дисперсии величины |(ЗСх)|. Этого можно добиться, если выбрать полё рСХ)

слабо коррелированным.

Пусть рСХ) - гауссовское поле с математическим ожиданием ФСХ) и корреляционной функцией ИХ,«) = МрСХ)рСы). После

о

прербразований для среднего значения |бСх)| получаем М|(ЗСх) |2 =

= |ехрПСХ,х)-С<о,х)+ФСХ)-ФСи)) х '

■ В1 х ехр - |( с2СХ) + ст2Си)- 2СКХ,ш) + §СХ)«о)))бХс1и,

здесь а2СХ) = КСХ,Х) - Ф2СХ).

Бели поле некоррелированное, то есть ИХ,а) - ФСХ)#Сы) = О,

то

М|(ЗСх) |2 = ||ехрШСХ,х)+ФСХ))) х ехрС- |С(72СХ))ёХ|2.

Езли теперь ГСх) - детерминированная действительная функция такая, что

ГСх) = |е1Хх ?СХ)с1Х =|еШ |ГСХ)| е^сЛ. ' Б О

то генерируя случайное поле рСХ), такое, что МрСХ)=$СХ) =ц<. X) и а2СрСХ)) = -21п |?СХ) |, получаем М|(ЗСх) |2 ' = Г2Сх). С Предполагается, что 0< |ГСХ)| < 1).

Для дальнейшего нам потребуется следующее

2 2

Утверждение: Пусть М? = а, М|? | = а , 1ша = 0, тогда | = а.

Те) есть дисперсия случайной величины г? = | равна нулю.

Возвращаясь к функции (ЗСх), получаем, что если случайная

1 рС X) ^

функция рСХ) удовлетворяет условиям Ме г = ГСХ) и МрСХ)рСш) = = МрСХ)МрСсо) при X * ы , то с вероятностью единиц?

|ГСх) | = |(ЗСх)| = 1|ехраССХ,х)+р(Х)))с1Х|. II

В реальной ситуации, генерируемое случайное поле не будет некоррелированным, поэтому будет выполняться равенство М|(ЗСх)|2 = Г 2С х)< +а2Сх),

Кх) = ГСойдСи.хЗбш,

,будет тем меньше, чем быстрее убывает : корреляционная функция ИХ).

Гфи создании киноформных элементов с помощью ЭВМ .естественно рассматривать дискретный вариант задачи. В этом случае область Б разбивается на ячейки, внутри каждой ячейки случайная функция постоянна и удовлетворяет условию Ме^^

. Л.

= ПХф), где Х.ф - некоторая средняя точка ячейки, значения поля 1>Ш в различных ячейках независимы При достаточно малых размерах ячеек с большой вероятностью будет приближенно выполняться необходимое равенство. В параграфе 3.1 доказано следующее Предложение-. Пусть

1

Б

причем |?Сса) | < 1. Рассмотрим последовательность случайных функций

п-1 п-1

ГпСх) = £ ^ | е кг дСы.хЗбы, « г-0 ^

здесь х= Сх^.х^З, * =(,,{ - последовательность независимых случайных величин таких, что

1ркг ~ Ме = ГСа1+кА,а2+гЛ),

Д = а/п, а - длина стороны квадрата Б, Са^а^ - его- левая

нижняя граничная точка,

0кг= { о) : а+кД < ^ < а+Ск+ПД, а+гД < (¿г < а+Сг+13Д>. Тогда для всех х при п —> ю

М^СхЗ —> ГСхЗ ; М|ГпСхЭ- КхЭ |2 ~> О; М|ГпСхЗ | —> |ГСхЗ| ;

М||ГпСхЗ|- 112 —> О.

Отметим, что. если |?СыЗ | ¿1, то для любого п существует

последовательность случайных величин С <р У таких, что Ме =

» кг

-V

= ГСа^+кД^+гД). Действительно, пусть (Ка^кД.а^гДЗ | = А^ ,

~ КГ Л

и ГСа^+кД.а^+гДЗ = А^е Пусть находится из условия

51п1кг/1^г= А^г, тогда в качестве р^можно взять, случайную

величину «|СГ+?)СГ. гДе равномерно распределена в. интервале

[-1 Д.]. к к

Напомним, что для преобразования Френеля, используемого при

синтезе оптических кинофэрмньгх элементов, дСа.х) = С1/Ххйк -

2 2 ' ехрС 1С пЛ<13 С С ы^ -х^) ■Кь^-у^З ]).■

Сказанное выше относилось к монохроматическому свету. Теоретический и практический интерес представляет синтез фазовых преобразователей, обеспечивавших заданные характеристики для волн различной длины [32,33,48,49].

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

АШ = 1X2 |ехрС-|^-1СУ'х2+у2+22 + СпСХ)-1Э1(х,уЗЗЗ&сс1у. Е

Здесь X - длина волны, Е - круг радиуса Е с центром в начале координат, пСХ) - показатель преломления, 1Сх,уЗ - толщина фаэосдвигающей пластинки в точке Сх,у,03. Для дальнейшего удобно ввести обозначение .КХЗ = 1ХгАСХЗ, кСХЗ = пСХЗ-1.

Изучается задача нахождения следующего

экстремума шах СшпС |.КА 3 \f\KKp I2). Показано, что если

1Сх,уЗ выполняется условие

кСХ, )Х_

12 т

Х^СХ^- Х^СХ^

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

пипС |.КХ 3 |?|ЛСХ23 I2). > созСг/пЗСтгЙ.

Таким образом, при достаточно больших п модули амплитуд для

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

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

В этой схеме сканирования объект расположен между источником излучения и приемником, представляющим собой шкалу из NDTK детекторов, раскрытую на угол ALDTK. Расстояние от источника до приемника - EIST, радиус круга, в котором расположен объект - RMAX, начало координат находится в точке 0. Томограф состоит из устройства вращения объекта и устройства параллельного перемещения системы "источни^приемник" относительно объекта. В начале работы томографа источник находится в точке -КМАХ, производится снятие данных. При одном снятии получаем NDTK интегралов. Затем система "источник-приемник" перемещается на шаг Н. Движение системы до крайней правой точки называется траверсом. После окончания траверса система возвращается в исходное положение, объект поворачивается на угол ALDTK по часовой стрелке, и начинается новый траверс. Бели рассматривать данные в плоскости Радона с координатами С г,р), то после каждого траверса получаем данные на дискретной сетке в полосе, ширина которой по р равна ALTO. Число траверсов равно NPOV=rc/ALTD.

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

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

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

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

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

ЛИТЕРАТУРА

1. Гельфанд U.M.,Граев М.Н.,Виленкин Н. Я. Обобщенные функции.

Вып. 5: Интегральная геометрия и связанные с ней вопросы теории представлений. М.: Фиэматгиэ, 1962.

2. Лаврентьев М. М. О некоторых некорректных задачах математической физики. Новосибирск: Наука, 1962.

3.. Тихонов А. Н. , Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1974.

4. Денисгж A.C. Исследования по интегральной геометрии в

вещественном пространстве. Дисс. канд. фиэ.- мат. наук. М., МГУ, 1991 г. мех. -мат. факультет.

5. Кириллов А. А., Об одной задаче И. М. Гельфанда. Докл. АН СССР, 1961, 137, N 2, 276-277.

6. Tuy Н. К., An inversion formula for cone-beam reconstruction. SI AM. J. APPL. MATH., 1983, vol.43, No 3,546-552.

7. Гельфанд И. M., Гончаров А. Б. , Восстановление финитной функции, исходя из ее интегралов по прямым, пересекающим данное множество точек в пространстве. Докл. АН СССР, ■ 1986, 290, N 5, 1037-1040.

8. Орлов С. С. Теория трехмерной реконструкции. Оператор восстановления. Кристаллография, 1975, 20, N 4, 701-709.

9. Аниконов Ю. Е. Некоторые методы исследования многомерных обратных задач для дифференциальных уравнений. Новосибирск: Наука, 1978..

10. Благовещенский A.C. О восстановлении функции по известным интегралам от нее, взятым вдоль линейных многообразий. Математические, заметки, 1986, т.39, N 6, с. 841-849.

И. Наттерер Ф. Математические аспекты компьютерной томографии. М.: Мир, 1990.

12. Гельфанд И. М. , Шилов Г. Е. Обобщенные функции. Вып. 1:

Обобщенные функции и действия кад ними. М.: Фиэматгиэ, 1959.

13. Семянистый В. И. Однородные функции и некоторые задачи

интегральной геометрии в пространствах постоянной кривизны. Докл. АН СССР, 1961, 136, N 2, 288-291.

14. Finch D. V. Cone beam reconstruction with sources on a curve. SI AM. J. APPL. MATH. 1985, vol.45, No 4,665-671.

15. Ашова H. Б., Голубятников В. П. Алгоритмы решения обратной задачи рентгеновской томографии. 4-й Всесоюзный симпозиум по вычислительной томографии, Ташкент С тезисы докладов). Новосибирск, 1989г.

16. Grangeat P. Analyse d'une système d'imagerie 3D par reconstruction à partir de radiographies X en géométrie conique. These de doctorat, Grenoble, 1987.

17. Хелгасон С. Преобразование Радона. M. : Мир 1983.

18. Айзенберг Л. Д., Кравцов Б. А., Шаимкулов Б. А. Об интерполяции сигналов с финитным спектром Фурье. Вычислительный эксперимент. Автометрия, 1989, N 4, стр. 23-29.

19. Сизых О. Н. Об экстраполяции функций* классов Винера и Харди в произведении полуплоскостей для случая решетки в полярных координатах. Автометрия, 1990, N 6,- стр. 91-95.

20. Огейн И, Вейс Г. Введение в гармонический анализ на евклидовых пространствах.- М. : Мир 1974.

21. Ярославский Л. П. Цифровая обработка сигналов в оптике и голографии. М. : Радио и связь, 1987.

22. Уиттекер Э. Т. ,, Ватсон Дд. Н. Курс современного анализа 1. , М. :Ф.-М. Л. 1963.

23. Малюженец Г. А. О зависимости между амплитудной характеристикой нелинейного элемента. Докл. АН СССР, 1946, т. 4, N 6.

24. Zodeh L. A. On'the Identification Problem. - IRE Tufns. On Cirevity Theory. 1956, v. 3, N 4, pp. 277,-281.

25. Аграновский M. Л., Баглай Pi Д., Смирнов К. К., К идентификации одного класса нелинейных- операторов. - Журнал

вычислительной математики и математической физики, 1978, т. 18, N 2, с. 284- 293.

26. Tonelli L. Su un- problema di Abel. Matematische Annalen, Bd. 99, Hft. 1/2, 1928.

27. Минский M., Пейперт С. Персептрон. M.: Мир. 1971.

28. Башкиров 0. А. Рудаметова С. Б. , Чудинович Б. М. Выделение и счет связных областей на дискретном изображении.

Автоматика и телемеханика, N 11, 1972, стр. 84-91.

29. Роэенфельд А. Распознавание и обработка изображений. Мир, 1972.

30. Lesen' L. В. , Hirsch Р. М. , Jordan J. A. The kinoform: a new wave form reconstruction device. - IBM J. Res. Develop. ,

1969, v. 13, p. 150-154; Кинофэрм, Зарубежная

радиоэлектроника, 1969, N 12, стр. 41-50.

31. Kermisch D. Image reconstruction from phase information only. - J. Opt. Sbc. Am., 1970, v 60, N 1, p 15-17.

32.Слюсарев Г. Г. Оптические элементы с фазовыми слоями. Докл. АН СССР, 1957, т. ИЗ, N 4, стр. 780-783.

33. Тудоровский А. И. Объектив с фазовой пластинкой. - Опт. и спектр, 1959, т 6, выл 2, стр. 68-73.

34. Stark Н. Sampling theorems in polar coordinates. J. Opt. Soc. Am., 1979, -v §9, N 11, p. 1519-1525.

35. Корольков В. П., Коронкевич В. П. и др. Киноформы: Технологии, новые элементы и оптические системы. Автометрия, N 3,4, 1989.

36. Гудмен Дж. Введение в Фурье-оптику. Мир, М. 1970.

Работы автора по теме диссертации

37. Трофимов 0. Е. О теореме Котельникова в полярных координатах. Автометрия, 1985, N стр. 110-113.

38. Trofimov О.Е. On sampling theorem in polar coordinates. Proceedings of third International conference on Automatic

Image Processing, v 2, pp 120- 122, Leipzig, 1989.

39. Смирнов К. К. , Трофимов О. Е. Восстановление функции по ее интегралам на прямых. Всесоюзная конференция "Автоматизация научных исследований на основе применения ЭВМ", 1979 (тезисы докладов), Новосибирск, ИАиЭ СО АН СССР.

40. Жирнов В. Т., Смирнов К, К, Трофимов О.Е. О численных методах решения задач томографии. В сб.: Методы и средства обработки изображений. Новосибирск: ИАиЭ СО АН СССР, 1982, стр. 105-114.

41. Тро<|ммов 0. Е. Тюренкова Л. В. Об одном способе численного восстановления изображения по многоракурсной томограмме.

Препринт Института автоматики и электрометрии СО АН СССР

N 440. Новосибирск, 1989.

42. Трофимов 0. Е. Численное восстановление интегралов от функции трех переменных вдоль прямых, проходящих через круг, по интегралам от той же функции вдоль прямых, проходящих через окружность. Препринт Института автоматики и электрометрии СО АН СССР N 427. Новосибирск,1989.

43. Баглай Р. Д., Смирнов К. К., Трофимов O.E. Об одном классе интегральных уравнений со слабой особенностью и их приложении к задачам идентификации нелинейных систем. Препринт

Института автоматики и электрометрии СО АН СССР N 182. Новосибирск,1982.

44. Старков М. А., Трофимов 0. Е., Фризен Д. Г. Автомат для подсчета числа плоских объектов. Автоматика и телемеханика, N 6, 1976, стр. 134-140.

45. Трофимов 0. Ё. Блок-схема автомата для счета объектов.

Автометрия, 1978, N 2, стр. 68-73. ■

46. Жирнов В. Т., Трофимов 0. Е. Шстрые линейные преобразования плоскости. Журнал вычислительной математики и математической физики. 1982, N 4, стр. 986-989.

47. Трофимов 0. Е. О задаче синтеза фазовых . преобразователей

волновых сигналов. Автометрия, 1979, N 2, стр. 21-25.

48. Трофимов O.E. Расчет киноформов. В сб.: Киноформные оптические элементы. Новосибирск : ИАиЭ СО АН СССР, 1981. стр. 84-92.

49. Трофимов 0. Е. Об одном способе синтеза киноформов. Автометрия, 1978, N 3, стр. 68-71.

50. Трофимов 0. Е. Численное восстановление интегралов от функции 3-х переменных вдоль прямых, проходящих через круг, по интегралам от той же функции вдоль прямых, проходящих через окружность. 4-й Всесоюзный симпозиум по вычислительной томографии, Ташкент С тезисы докладов). Новосибирск, 1989.

51. Трофимов 0. Е. Соотношение между дискретным и непрерывным преобразованиями Фурье. Препринт Института автоматики и электрометрии СО АН СССР N 424. Новосибирск, 1989.

52. Трофимов 0. Е. О восстановлении функции с финитным носителем по ее преобразованию Радона. Препринт Института автоматики и электрометрии СО АН СССР N 427. Новосибирск,1989.

53. Трофимов 0. Е. К задаче восстановления функции трех переменных по ее интегралам на прямых, пересекающих заданную кривую. Автометрия, 1991, N5 , стр. 57-64.

54. Трофимов' O.E. О восстановлении функции по лучевым данным.. 5-й Всесоюзный симпозиум по вычислительной томографии, Звенигород С тезисы-докладов). Москва, 1991.

55. Трофимов O.E. О преобразовании Фурье в смысле обобщенных функций лучевых данных в конусе. Всеросийская конференция: Условно-корректные задачи математической физики и анализа, (тезисы докладов). Новосибирск, 1992.

56 Трофимов 0. Е. К задаче восстановления функции по ее интегралам вдоль прямых, пересекающих заданную кривую. В сб.: Методы решения условно-корректных задач. Новосибирск Институт математики СО РАН , 1991,. стр. 80-94.

57. Трофимов Е. О численных алгоритмах трехмерной томографической реконструкции. Автометрия, 1992, N3, стр. 81-86.