Диссертация (1145332), страница 9
Текст из файла (страница 9)
является то, что из двух наборов коэффициентов {µ1 , k1 , J1 }, {µ2 , k2 , J2 } утверждается совпадение только µ1 и µ2 , а обостальных ничего не говорится. Из формулы (1.34) видно, что неизвестныефункции k, J не влияют на процедуру определения функции µ. В физическихтерминах это означает, что влияние рассеяния и наличие внутренних источников подавляется за счет выбора соответствующих внешних источников.2.3.Алгоритм восстановления коэффициента ослабления наоснове многократного облученияРассмотрим обратную задачу, являющуюся частным случаем задачи томографии, рассмотренной в §2.2. Требуется определить функцию µ из соотношений1ω·∇r f (r, ω, α)+µ(r)f (r, ω, α) =4πZk(r, ω·ω 0 )f (r, ω 0 , α)dω 0 +J(r, ω), (1.41)Ωf (ξ, ω, α) = h(ξ, ω, α),(ξ, ω, α) ∈ Γ−0 × [0, 1],(1.42)f (η, ω, α) = H(η, ω, α),(η, ω, α) ∈ Γ+0 × [0, 1],(1.43)в которых известны только функции h и H.
Эту задачу для краткости будемназывать обратной задачей (1.41)-(1.43).Пусть функция h удовлетворяет условиям (h1 )-(h4 ). Нетрудно видеть,что все результаты, полученные для рассмотренной в §2.2 обратной задачи,будут справедливы и для задачи (1.41)-(1.43), которую в целях экономии ресурсов ЭВМ мы и будем рассматривать для тестирования алгоритма.Обратимся к формуле (1.34), которая для задачи (1.41)-(1.43) выглядит55следующим образом:0d(η,−ωZ )µ(η − tω 0 )dt = ln0∂0∂α h(ξ, ω , α)limα→0 ∂ H(η, ω 0 , α)∂α!≡ ge(η, ω 0 ),(η, ω 0 ) ∈ Γ+0,(1.44)где ξ = η − d(η, −ω 0 )ω 0 . Фиксируем сечение η3 = const, для которого справедлива формула (1.44), и в этом сечении введем полярные координатыr1 = ρ cos γ,r2 = ρ sin γ,ρ ∈ [0, ∞), γ ∈ [0, 2π).00Пусть ω 0 = (− sin γ, cos γ, 0), ω⊥= (cos γ, sin γ, 0), ω⊥· ω 0 = 0, тогда, учиты-вая, что µ(r) = 0 вне G, равенство (1.44) запишется в видеZ+∞0µ(ρω⊥(Rµ)(ρ, γ) ≡(γ) − tω 0 (γ))dt = g(ρ, γ),(1.45)−∞0где g(ρ, γ) ≡ ge(η(ρ, γ), ω 0 (γ)), η(ρ, γ) = r + d(r, ω 0 )ω 0 , r = ρω⊥(γ).
Эта фор-мула справедлива почти в каждом горизонтальном сечении. Таким образом,в каждом таком сечении нам известны интегралы от следа функции µ по по0чти всем прямым, проходящим через точки r = ρω⊥(γ) в направлении ω 0 (γ),ρ ∈ [0, ∞), γ ∈ [0, 2π). То есть, задано двумерное преобразование Радона(Rµ)(ρ, γ) от функции µ(r). Получается, что задача определения функции µсвелась к задаче обращения преобразования Радона в горизонтальных сечениях области G.Выпишем основные этапы тестирования алгоритма решения обратной задачи (1.41)-(1.43).1. Задание области G, разбиения G0 и функций µ, k, J.
Задание функции h(r, ω, α), удовлетворяющей условиям (h1 )-(h4 ), и сечения, в которомвосстанавливается функция µ.2. Численное решение прямой задачи (1.41), (1.42) и вычисление функции∂0∂α H(η, ω , α)при достаточно малом положительном числе α.563. Вычисление функцииge(η(ρ, γ), ω 0 (γ)) = ln−1 !∂∂h(ξ, ω 0 , α)H(η, ω 0 , α)∂α∂αпри достаточно малом положительном числе α.4. Решение уравнения(Rµ)(ρ, γ) = g(ρ, γ)относительно неизвестной функции µ.5.
Графическое представление распределения функции µ(r) в некоторомсечении как полутоновое изображение оттенков серого цвета.Определенные проблемы возникают на этапе 2. Ввиду многомерности задачи (1.41), (1.42) наиболее приемлемым методом ее решения является методМонте-Карло, для которого характерны медленная сходимость и, как следствие, продолжительное время вычислений на ЭВМ. Поэтому, чтобы обеспечить удовлетворительный временной интервал вычислений, были выбраныне самые "тяжелые" тесты.Для проверки алгоритма решения обратной задачи (1.41)-(1.43) нам требуется вычислить решение прямой задачи в наборе точек (η, ω), причем сдостаточно хорошей точностью.
Для этого была выбрана одна из версий метода Монте-Карло, называемая методом сопряженных блужданий [62,67]. Сутьэтого метода заключается в том, что решение прямой задачи ищется в видеусеченного ряда Неймана (1.13)f (r, ω) =mX((AS)n f0 )(r, ω),n=0который представляет собой сумму интегралов((AS)n f0 )(r, ω),n = 0, 1, ..., m,(1.46)57гдеd(r,−ω)Z1exp(−τ (r, ω, t))4π(ASf0 )(r, ω) =0Zk(r − tω, ω · ω 0 )f0 (r − tω, ω 0 )dω 0 dt,Ωf0 (r, ω) = h(r − d(r, −ω)ω, ω) exp(−τ (r, ω)) + (AJ)(r, ω).Для вычисления этих многомерных интегралов применяются вероятностныеметоды [30, 78], причем вычисления в методе сопряженных блужданий организованы так, что все слагаемые усеченного ряда Неймана (1.46) вычисляются одновременно.При выполнении четвертого этапа применялся алгоритм свертки и обратной проекции [89,148]. Этот алгоритм в ряду всех методов, используемых дляобращения преобразования Радона, является одним из самых распространенных [148].Метод сопряженных блужданий и алгоритм свертки и обратной проекциибыли программно реализованы и протестированы.Для описания вычислительных экспериментов, на которых тестировалсяалгоритм решения обратной задачи (1.41)-(1.43), нам понадобятся следующиеобозначения:m − количество членов усеченного ряда Неймана;n − число траекторий (итераций), используемых для вычисления выходящего излучения H в точке (η, ω);p, q − числа, характеризующие сетку при вычислении функции H(η(ρi , γj ),ω 0 (γj )), i = 1, ..., p, j = 1, ..., q;scat(r) − величина, определяющая в процентах уровень рассеяния, онаравна альбедо однократного рассеяния, умноженному на 100;A − величина, определяющая во сколько раз максимально ослабляетсянерассеянное излучение при прохождении через среду G.
Она равна exp(τ ),где τ - оптический диаметр области G.Для тестирования алгоритма решения обратной задачи (1.41)-(1.43) былореализовано два численных эксперимента, в которых область G являлась58единичным шаром с центром в начале координат. Восстановление функцииµ осуществлялось в плоскости r3 = 0.Tecт 1. Рассматриваемый тестовый пример называется фантом Кормака[112]. Множество G0 является объединением четырех подобластей Gi :G1 = {r = (r1 , r2 , r3 ) : 0.81 < r12 + r22 + r32 < 1},G2 = {r = (r1 , r2 , r3 ) : (r1 − 0.5)2 + r22 + r32 < 0.04},G3 = {r = (r1 , r2 , r3 ) : r12 + r22 + r32 < 0.01},G4 = G\3[Gi .i=1Функция J(r, ω) ≡ 0, а коэффициенты µ, k постоянны в Gi и определяютсяследующим образом:(µ(r) =(k(r) =3, r ∈ G1SG2SSG2G3 ,1, r ∈ G4 ,0.9µ(r), r ∈ G1µ(r),SG3 ,r ∈ G4 .Интенсивность входящего излучения определяется равенствомh(r, ω, α) = (|r3 | + |ω3 | + α)0.5 .(1.47)Интенсивность выходящего излучения H аппроксимировалась суммой первых 10 членов ряда Неймана.
Для вычисления функции H(r, ω) в одной точке(r, ω) было реализовано 5000 итераций метода Монте-Карло. Функция∂∂α Hрассчитывалась в 10000 точках.На рисунке 1.1 представлены результаты численного эксперимента приследующих значениях параметров:n = 5000, m = 10, p = 30, q = 200, 90% ≤ scat(r) ≤ 100%, A ' 50.59Tecт 2. Множество G0 состоит из четырех областей Gi :G1 = {r = (r1 , r2 , r3 ) :r1 − 0.60.042+ r 220.3 r 23+< 1},0.04G2 = {r = (r1 , r2 , r3 ) : (r1 + 0.5)2 + (r2 − 0.2)2 + r32 < 0.0025},G3 = {r = (r1 , r2 , r3 ) : (r1 − 0.1)2 + (r2 + 0.6)2 + r32 < 0.0016},G4 = G\3[Gi .i=1Функция J(r, ω) ≡ 0, а коэффициенты µ, k постоянны в Gi :(µ(r) =0, r ∈ G1SG2SG3 ,5, r ∈ G4 .Индикатриса рассеяния k(r) = 0.95µ(r).
Интенсивность входящего излучения определяется равенством (1.47). Интенсивность выходящего излученияH аппроксимировалась суммой первых 3 членов ряда Неймана. Для вычисления функции H(r, ω) в одной точке (r, ω) использовалось 2000 итерацииметода Монте-Карло. Функция∂∂α Hвычислялась в 12000 точках.На рисунке 1.2 представлены результаты второго численного эксперимента при следующих значениях параметров:n = 2000, m = 3, p = 20, q = 200, scat = 95%, A ' 450000.Представленные результаты численных экспериментов демонстрируют эффективную работу разработанного алгоритма.60Рисунок 1.1 — Графическое представление множества значений коэффициента ослабления µ в плоскости r3 = 0.
В верхней части рисунка представленоригинал функции, а в нижней – ее реконструкция, полученная методом многократного облучения (тест 1 – фантом Кормака).61Рисунок 1.2 — Графическое представление множества значений коэффициента ослабления µ в плоскости r3 = 0. В верхней части рисунка представленоригинал функции, а в нижней – ее реконструкция, полученная методом многократного облучения (тест 2).622.4.Алгоритмы параллельных вычислений реконструкцииструктуры трехмерного объектаВосстановление структуры трехмерного объекта по данным радиационного просвечивания является важной задачей компьютерной томографии. Одним из подходов в решении указанной задачи является сканирование исследуемого объекта по отдельным сечениям.
При этом восстановление изображения в отдельном сечении в итоге сводится к решению задачи обращенияпреобразования Радона. Достаточно популярный для решения указанной задачи алгоритм свертки и обратной проекции характеризуется хорошим быстродействием. К примеру, реконструкция изображения размером 400 × 400пикселей, полученная по данным просвечиваний по 10201 прямым, требуетоколо 2 секунд машинного времени при использовании процессора Intel(R)Xeon(R) CPU E5345 @ 2.33GHz.Однако, восстановление трехмерного изображения объекта по набору сечений в более высоком разрешении и с большим количеством просвечиванийзаймет существенное время.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
