2. В узлах решетки аппроксимирующая функция совпадает с заданными отсчетами, а ее производная в этих точках равна выборочной. То есть справедливы равенства: [pic], [pic], здесь h = 2R/(M-1), I(r0,?j) =

I(rM+1, ?j) = 0, i = 1, -, M.

3. На интервале [ri, ri+1] функция [pic]есть полином третьей степени от r.
Перечисленные условия позволяют в явном виде получить коэффициенты соответствующего сплайна. Непосредственными вычислениями можно получить, что
[pic], где
[pic]
Q(x) = Q(-x), Q(x) = 0 при |x|> 2h, h=ri+1-ri.
Функция Q(x) имеет разрывы второй производной, но модуль второй производной интегрируем, используя это обстоятельство можно показать, что свертка S0(z) = Q(x) [pic](-1/?z2) выражается формулой (1.5.1).
Непосредственными вычислениями получаем
[pic]
Графики функций Q(x) и S0(z) для различных значений h представлены на рис. 1 и рис. 2.
[pic][pic]
Таким образом,
[pic].
Заменяя в (1.5.2) S на [pic]и интеграл частной суммой, получаем f*(x, y)
- приближение к функции f(x, y),
[pic]. (1.5.3)
Как уже отмечалось выше, обычно в компьютерной томографии используется метод свертки и обратного проецирования. Рассмотрим соотношение между этим методом и методом, изложенным в настоящем параграфе. Используя интегрирование по частям, свертку с обобщенной функцией 1/z2 можно заменить дифференцированием и сверткой с 1/z (преобразованием Гильберта).
То есть функцию
S(z, ?) = I(z, ?) [pic]1/z2 можно представить в виде
S(z, ?) = Iz/(z, ?) [pic]1/z
При построении численных алгоритмов вместо обобщенной функции 1/z или, что то же самое, интеграла в смысле главного значения, в методе свертки и обратного проецирования используют некоторую последовательность регулярных функций pА(z), сходящуюся к 1/z (в смысле обобщенных функций) при A стремящемся к бесконечности. Используя интегрирование по частям, дифференцирование переносят на функции pА(z) и таким образом получают регулярные функции, сходящиеся к 1/z2, то есть свертка с обобщенной функцией 1/z2 заменяется последовательностью сверток с регулярными функциями p/А(z).
Таким образом, шаг свертки в классическом методе можно интерпретировать следующим образом: исходные данные аппроксимируются ступенчатой функцией и осуществляется свертка с регулярной функцией, являющейся приближением к обобщенной функцией 1/z2.
В методе настоящего параграфа исходные данные аппроксимируются более гладкими функциями - сплайнами 3-го порядка. Это позволяет точно вычислить свертку с обобщенной функцией 1/z2, причем в явном виде.
Шаг обратного проецирования соответствующий интегрированию свертки в обоих алгоритмах одинаков.
При использовании алгоритмов в реальных ситуациях важно уметь оценивать влияние шумов на точность получаемых приближений. Наличие явного выражения для аппроксимирующей функции позволяет вычислить дисперсию ошибки в любой точке при фиксированных ?r, ?? ? известных статистических характеристиках шума. Для случая независимого, аддитивного, стационарного шума ? (z) можно сделать следующее замечание. Рассмотрим процесс ?, являющийся сверткой с
1/z2 процесса ?. Спектральная плотность этого линейного преобразования есть
|?|. Для спектральных плотностей процессов ? и ? получаем соотношение f ?
(?) = |?|2f? (?). ?исперсия процесса ? конечна, если интегрируема f? (?),
?о есть процесс ? дифференцируем в среднеквадратическом. Для того, чтобы свертка выражалась формулой (1.5.1), на процесс ? нужно наложить дополнительные условия, потребовав, например, чтобы выборочные функции с вероятностью единица имели конечную вторую производную.
Численное моделирование и восстановление плотности реальных объектов с использование метода, изложенного в настоящем параграфе, показало высокую точность метода, особенно при исследовании объектов и дефектов, имеющих сложную конфигурацию и участки с резкими границами.
Примеры восстановления, с использованием методов, изложенных в настоящем параграфе, приведены на рисунке 3. Тестовый объект состоит из 10 частиц.
Рис.3(слева) соответствует 10 поворотам и Рис.3 (справа) соответствует 32 поворотам. [pic]

4.2. Приведение формул обращения томографической реконструкции в конусе лучей к виду, позволяющему строить численные алгоритмы.
В компьютерной рентгеновской томографии трехмерный объект представляется обычно в виде набора тонких срезов. Для восстановления плотности среза решается задача обращения двумерного преобразования Радона. Для исследования ряда объектов более естественной является другая схема, когда источник излучения движется по некоторой пространственной кривой. Каждой точке кривой соответствует конус лучей, проходящих через эту точку.
Исходными данными являются данные об ослаблении излучения при прохождении через объект. Математически задача ставится как задача восстановления функции трех переменных по интегралам вдоль прямых, проходящих через заданную кривую. Была получена формула обращения для функций, имеющих финитный носитель, и для кривых, удовлетворяющих определенным условиям.
Главным в этих условиях является то, что любая плоскость, пересекающая объект пересекает кривую, по которой движется источник. Примером кривой, удовлетворяющей условиям, является совокупность двух единичных окружностей, лежащих во взаимно перпендикулярных плоскостях. Однако построение численных алгоритмов непосредственно на основании этой формулы, затруднительно. Дело, в частности, в том, что формула обращения основана на преобразовании Фурье от однородной функции, получаемой из исходных данных. Причем преобразование
Фурье понимается в смысле обобщенных функций, а преобразование Фурье в обычном смысле может не существовать. В настоящем параграфе приводятся выражения для используемого преобразования Фурье, позволяющие при построении численных алгоритмов использовать метод, изложенный в предыдущем параграфе. Устанавливаются также некоторые соотношения между результатами
Пусть заданы функция f(x) = f(x1, x2, x3) , точка S = (s1, s2, s3) и вектор a = (a 1, a 2, a 3). Лучевым преобразованием функции f(x) будем называть функцию
[pic], являющуюся интегралом от f(x) вдоль луча, исходящего из точки S в направлении вектора a .
Наряду с функцией [pic]в некоторых ситуациях рассматривается функция
[pic], являющаяся интегралом по всей прямой или, что тоже самое, суммой интегралов вдоль лучей из точки z в направлениях a и -a .
Множество точек S, для которых известно лучевое преобразование обычно является множеством точек, принадлежащих некоторой кривой, являющейся траекторией движения источника излучения.
Пусть задана кривая, по которой движется источник, Ф(l ) = (Ф1(l ), Ф2(l
), Ф3(l )) параметр l пробегает некоторый интервал Щ действительной прямой.
Для любого a = (a 1, a 2, a 3) и l О Щ определим функцию
[pic].
Функция g(a ,l ) есть интеграл от функции f(x) вдоль проходящего через точку Ф(l ) в направлении вектора a . Отметим, что при любом фиксированном l функция [pic]является l однородной функцией a степени -1:
[pic]. (2.1.1)
Для функций, имеющих финитный носитель, в [101] получена формула:
[pic]. (2.1.2)
При фиксированном l функция G+(b ,l ) есть преобразование Фурье от функции [pic]по переменной a , b = (cosq cosf , sinq cosf , sinf ). В формуле (2.1.2) l зависит от x и b и выбирается из условий: скалярное произведение (b , x) равно (b ,f (l )), но (b ,Ф(l )) не равно нулю.
Значение функции f(x) может быть восстановлено в точке x, если такое l существует для любого b . Геометрически это означает, что любая плоскость, пересекающая точку x носителя функции, пересекает кривую Ф(l ) так, что знаменатель в (2.1.2) не обращается в нуль. Примером кривой, удовлетворяющей условиям Кириллова-Туя, является совокупность двух единичных окружностей, лежащих во взаимно перпендикулярных плоскостях, если носитель лежит в единичном шаре. Для цилиндрических объектов можно использовать винтовую линию.
В формулу (2.1.2) входит G+(b ,l ) - преобразование Фурье от функции
[pic], однако преобразование Фурье, понимаемое в обычном смысле:
[pic], в данном случае не существует, так как [pic]является однородной и имеет на бесконечности порядок 1/к a к . Преоразование Фурье здесь понимается в смысле обобщенных функций. Поскольку [pic]однородная функция, то при любом фиксированном l исходные данные, полностью определяются своими значениями на поверхности к a к =1. Переход к функции, заданной во всем пространстве
R3 при использовании преобразования Фурье приводит к обобщенным функциям.
Преобразование Фурье в смысле обобщенных функций является линейным функционалом над соответствующим пространством. Подробнее об этом будет сказано в следующих параграфах. Здесь нам важно отметить, что не любой функционал задается с помощью регулярной функции. Для того, чтобы использовать формулы типа (2) для построения алгоритмов, необходимо показать, что [pic]задается с помощью регулярной функции и иметь для нее выражения через функцию [pic]. В работе [101] дается выражение, связывающее
[pic], при x отличном от нуля с помощью регулярных операций с искомой функций f(x), то есть фактически показано, что функционал [pic]задается с помощью регулярной функции. Однако для построения алгоритмов томографической реконструкции нужно [pic]выразить не через искомую функцию f(x), а через исходные данные [pic].
Итак, перейдем к нахождению [pic]. Мы будем использовать то, что
[pic]является однородной функцией по a фиксированном l . В [95] доказано следующее
Утверждение: Пусть [pic]есть преобразование Фурье в смысле обобщенных функций от однородной функции [pic], тогда
[pic]. (2.1.3)
Строгое доказательство требует существенного использования аппарата обобщенных функций, понимаемых как линейные функционалы над соответствующим пространством. Здесь мы ограничимся изложением основных моментов доказательства. В частности, замену переменных в расходящихся интегралах мы будем делать по тем же правилам, что и в обычных.
Представим [pic]в виде
[pic],
(поскольку параметр l фиксирован, его на данном этапе можно опустить).
Как уже отмечалось выше, интеграл является расходящимся, тем не менее, переходя к сферическим координатам по обычным правилам, получаем:
[pic], где b = b (j ,q ) = (cosq cosj , sinq cosj , sinj ), j О [-p /2, p /2], q
О [0, p ].
Учитывая, что [pic], а также то, что интегрирование по углам j и q соответствует интегрированию по единичной сфере, приходим к выражению
[pic].
Интеграл по r есть преобразование Фурье от r ++. Используя таблицы для преобразования Фурье обобщенных функций [19], приходим к выражению (2.1.3).
Для действительных функций f(x) в формуле (2) нужна мнимая часть [pic]:
[pic].
Используя обобщенные функции, сосредоточенные на поверхности [19], получаем следующее следствие:
[pic].
Здесь S(x ) = {g О S2Ѕ (x , g ) = 0), [pic]v производная по направлению x
. Подставляя в (2.1.2) функции [pic]и [pic], зависящие от параметра l , получаем формулу обращения, пригодную для построения численных алгоритмов:
[pic](2.1.4)
Здесь S(x ) v окружность, являющаяся пересечением единичной сферы и плоскости P(b ). Плоскость P(b ) проходит через начало координат ортогональна вектору b . Символ W (x ) означает интегрирование по окружности. Оператор L(b , D) означает дифференцирование функции [pic]в направлении вектора b :
[pic], при этом l , зависящее от b и x, остается фиксированным.
Как и выше, b = b (q , j ) = (cosq cosj , cosq sinj , sinq ), l = l (q , j ) = l (x, b ) такое, что скалярное произведение (x, b ) равно (b , g (l
)) и (b , g /(l )).
В формуле (4) используются регулярные функции, и она пригодна для построения численных алгоритмов.
Замечание. А.С. Денисюком независимо и другим методом, без явного использования преобразования Фурье обобщенных функций, получены формулы обращения функции g+ в Rn . При n = 3 формулы А.С. Денисюка и формулы, получаемые изложенным способом из формулы Туя, совпадают.
Выше были получены формулы, позволяющие строить численные алгоритмы восстановления функции f(x) = f(x1, x2, x3) по ее лучевому преобразованию
[pic]
Далее мы будем опускать символ f и использовать обозначение [pic].
При фиксированном S функция [pic]является функцией в трехмерном пространстве, но в силу ee однородности существуют поверхности, такие что
[pic]полностью определяется своими значениями на них (поверхности расположения приемников излучения).
Исходные данные в виде функции [pic]удобно использовать, если матрица приемников расположена на сфере. Однако в реальных ситуациях матрицу приемников обычно располагают на плоскости или поверхности цилиндра. В этих случаях удобно использовать несколько иной вид исходных данных.
Плоский детектор.
Мы будем предполагать, что для источника, находящегося в точке S = (s1, s2, s3), исходные данные регистрируются в плоскости P, определяемой уравнением xs1 + ys2 + zs3 = -Ѕ SЅ . Плоскость P, определяется следующими условиями: плоскость P перпендикулярна лучу, соединяющему источник с началом координат; плоскость P проходит через точку S= (s1, s2, s3.)
Расстояние D между плоскостью регистрации и источником равно удвоенному расстоянию от источника до начала координат. В плоскости регистрации будем использовать прямоугольную систему координат (p1, p2), начало которой находится в точке пересечения с лучем, соединяющим источник с точкой (0, 0,
0). Таким образом, если источник находится в точке S = (s1, s2, s3), то начало системы координат (p1, p2) в плоскости наблюдения находится в точке с трехмерными координатами -s1, -s2, -s3 =- S.
При реконструкции в конусе лучей наиболее распространенными примерами траекторий источника являются винтовая линия и совокупность двух окружностей лежащих в пересекающихся плоскостях.
Траектория в виде двух окружностей.
Рассмотрим окружность, лежащую в плоскости z =0.
Направление оси p2 в плоскости регистрации будет совпадать с направлением оси z.
Ось p1 системы координат возьмем на линии пересечения плоскости регистрации с плоскостью, содержащей окружность, по которой движется источник. Для окончательного определения системы координат необходимо выбрать одно из двух возможных направлений оси p1. Если s3 = 0, s1 = r cosl
, s2 = r sinl (источник движется в плоскости z =0), то положительный единичный вектор на оси p1 выберем так, чтобы он совпадал с вектором (cos(l
+p /2), sin(l +p /2), 0) = (-sinl , cosl , 0) = (-s2/Ѕ SЅ , s1/Ѕ SЅ , 0).
Точка, имеющая в плоскости регистрации координаты (p1, p2), имеет следующие пространственные координаты: x = -p1 sinl - r cosl = -p1 s2 /Ѕ SЅ - s1 , y = p1 cos l - r sinl = p1 s1 /Ѕ SЅ - s2 , z = p2.
В случае плоского детектора, исходными данными являются интегралы по лучам, соединяющим точки (p1, p2) в плоскости регистрации с источником S.
Регистрируемая функция gr(p1, p2, l ) есть интеграл от искомой функции f(x) = f(x1, x2, x3) вдоль луча исходящего из точки S = (s1, s2, s3) =
(rcosl , r sinl , 0) в направлении точки
P = (-p1 sin l - rcosl , p1 cosl - r sinl , p2 ) = (-p1 s2/Ѕ SЅ v s1, p1 s1/Ѕ SЅ v s2, p2).
Интегральная форма регистрируемой функции имеет вид:
[pic]
При t = 0 луч проходит через точку S = (rcosl , rsinl , 0), при t = 1 v через точку P = (p1, p2) = (-p1 sin l - rcosl , p1 cosl - r sinl , p2).
Итак, мы имеем соотношение между функциями gr(p1, p2, l ) и [pic]:
[pic],
[pic].
Наряду с обозначением gr(p1, p2, l ), мы будем использовать обозначения gr(p1, p2, S(l )), gr(p1, p2, S) и gr(P, S) , здесь S(l ) точка на траектории источника, соответствующая параметру l , P = (p1, p2). Мы выразили функцию gr(p1, p2, l ) через функцию [pic]= g+ (x , l ).
В формуле обращения лучевого преобразования используется функция g+ (x , l ) =[pic] для того, чтобы использовать gr(p1, p2, l ), регистрируемую в случае плоского детектора, нужно выразить g+ (x , l ) используя gr(p1, p2, l ).
Для дальнейшего нам потребуются координаты (p1, p2) (в системе координат плоскости регистрации) точки пересечения плоскости регистрации данных с лучем (S +tx ) = (s1 + tx 1, s2 + tx 2, s3 + tx 3). Эти координаты имеют вид:
[pic]
[pic].
[pic].
Теперь мы можем выразить [pic]используя gr(p1, p2, l ):
[pic]= g+ (x , l ) = gr(2 Ѕ S(l )Ѕ (s2(l )x 1 v s1(l )x 2) /[pic], -2Ѕ
S(l )Ѕ 2x 3 /[pic],l ), если [pic]< 0, [pic]= 0, если [pic]і 0.
Итак, мы имеем следующее соотношение между функциями: g+ (P, l ) и [pic]= g+ (x , l ); P = (p1, p2), x = (x 1, x 2, x 3,);
[pic]= g+ (x , l ) =
= gr(2 Ѕ S(l )Ѕ (s2(l )x 1 v s1(l )x 2) /[pic], - 2Ѕ S(l )Ѕ 2x 3 /[pic],l
), если [pic]< 0,
[pic]= 0, если [pic]і 0.
При переходе от функции g+ (x , l ) = [pic]к функции gr (P, S) интегрирование по окружности S(l ) в трехмерном пространстве заменяется на интегрирование по прямым линиям в плоскости регистрации. Отметим, что формулы обращения лучевого преобразования, использующие интегрирование вдоль прямых в плоскости регистрации.

Страницы: 1, 2, 3, 4, 5, 6