3.7. Итерационные методы реконструкции
Идея итерационных методов реконструкции состоит в последовательной коррекции приближенных значений рентгеновской плотности отдельных ячеек с целью их согласования со значениями измеренных лучевых сумм [2].
Сначала задается произвольный начальный набор значений
, например постоянная по изображению плотность. Затем по этим начальным значениям рассчитываются проекции. Если полученная расчетом лучевая сумма меньше ее измеренного значения, то плотность каждой ячейки, дающей вклад в данную лучевую сумму, увеличивается на определенную величину. После проведения такой процедуры для всех ячеек и лучей первая итерация считается выполненной. Количество последовательных итераций определяется степенью необходимой точности и скоростью убывания ошибки. Обычно весовые коэффициенты при итерационной реконструкции не вычисляют по строгим формулам, а принимают равными 1 или 0 в зависимости от того, лежит ли центр ячейки в полосе Dl, занимаемой лучом, или нет. Процесс последовательных коррекций аналогичен процедуре обратного проецирования и может быть представлен как
где 
— плотность i-й ячейки до очередной итерации; mk — плотность i-й ячейки после итерации; 
— коррекция, вносимая в i-ю ячейку от j-го луча.
Обычно используются два варианта расчета коррекций : аддитивный или мультипликативный. В первом случае коррекция пропорциональна величине весового коэффициента aij:
(3.31)
где Dpj k = рj – рj k — разность между измеренным и рассчитанным после очередной итерации значением лучевой суммы.
В случае мультипликативной коррекции величина 
пропорциональна плотности ячейки 
, определенной при предыдущей итерации:
. (3.32)
В зависимости от организации процедуры итераций обычно выделяют три наиболее употребительных метода: одновременной коррекции ILST (iterative Least — Squares Technique); поточечной коррекции SIRT (Simultaneous Iterative Reconstruction Technique) и получевой коррекции ART (Algebraic Reconstruction Technique).
В методе ILST все коррекции рассчитываются одновременно в начале итерации и затем корректируются сразу все ячейки; число необходимых умножений на одну итерацию составляет при прежних условиях pn3. При методе SIRT каждая ячейка корректируется одновременно для всех лучей, проходящих через нее, причем при коррекции последующих точек учитываются результаты предыдущих коррекций. Число необходимых умножений на одну итерацию даже больше, чем в методе обращения матрицы и составляет p2n4.
При реконструкции по методу ART коррекция осуществляется одновременно для всех точек, дающих вклад в отдельный луч. Затем процедура повторяется для следующего луча и т.д. При этом учитываются результаты предыдущих коррекций. По затратам машинного времени (число умножений порядка pn3) этот метод наиболее быстрый.
Рассмотрим элементарный пример реконструкции двумерного объекта, когда привлечения ЭВМ вообще не требуется.
Возьмем в качестве объекта квадрат ABCD (рис. 3.10), разделенный на 9 равных клеток. Числа от 1 до 9, разбросанные по клеткам, соответствуют плотности или какой-нибудь другой характеристике, находимой томографически. Пусть известны 4 проекции, определяемые направлениями сторон квадрата АВ и AD и его диагоналей АС и BD. Если в каждой проекции взять по 3 лучевые суммы, то в первых двух случаях вклад внесут все клетки, а в двух других – лишь по 7 клеток из 9.
Таким образом, мы исходим из 12 значений лучевых сумм и ищем 9 структурных элементов объекта, т.е. решаем переопределенную задачу.
Рис. 3.10. Тест-объект из девяти элементов Стрелками указаны 12 используемых лучевых сумм в четырех проекциях
Выберем для восстановления алгоритм ART. Начнем с проекции, образованной лучами, параллельными стороне АВ, т.е. P(0°)=(15,15,15). Каждое значение лучевой суммы разделим на число пересекаемых клеток и припишем этим клеткам найденную величину (в данном случае 5). Будем считать полученный результат первой итерацией (рис. 3.11,а).
Как видим, в нашем случае для взятой проекции объект представляется совершенно однородным. Если погрешность восстановления оценивать по формуле
(3.33)
где индекс k нумерует клетки, i обозначает номер итерации, 
— исходные значения на рис. 3.10, то s1 = 45,8 %.
Перейдем к следующей проекции (лучи идут вдоль AD). Теперь для каждого луча следует скорректировать сумму чисел, получаемых после первой итерации, на известную лучевую сумму данной проекции. Так, сумму в первом столбце (равную 15) следует, очевидно, уменьшить на 2, т.е. вычесть из каждого числа по 2/3; в третьем столбце нужно, наоборот, добавить в каждую клетку по 2/3; во втором столбце изменений нет. Видно, что в нашем примере вторая итерация фактически оказывается неинформативной: предыдущее однородное распределение лишь слегка деформируется, создавая небольшой градиент вдоль АВ и совершенно не выявляя структуры объекта (рис. 3.11,б). Погрешность восстановления, оцениваемая по (3.33), даже несколько возрастает по сравнению с первой итерацией: s2 — 48,7 %.
а) б)
в) г)
д) е)
Рис.
3.11. Различные стадии (а — е) восстановления тест - объекта, изображенного на рис. 3.10, итерационным методомТретья итерация (рис. 3.11,в, лучи идут параллельно диагонали АС) заметно меняет дело, поскольку в значениях лучевых сумм неоднородность объекта проявляется отчетливо. Принцип выравнивания остается прежним, только соответствующие разности следует равномерно распределять по двум или трем клеткам в зависимости от того, какая лучевая сумма принята в расчет; при этом получается s3 = 28,7 %. Такая же ситуация возникает и при четвертой итерации (рис. 3.11,г), завершающей первый цикл процедуры; теперь объект уже напоминает исходный, s4 = 14,6 %.
Далее можно вновь привлечь первую проекцию и начать второй цикл выравнивания. На рис. 3.11,д,е показаны результаты, получаемые после шестой (s6 = 9,94 %) и восьмой (s8 = 4,25 %) итераций. Процесс можно продолжать и дальше (его сходимость гарантирована), но уже и из проделанных вычислений ясно, что 7—8 итераций позволяют получить неплохой результат восстановления.
В контрольных вопросах к этой главе предложена задача (7), которую рекомендуем решить при помощи калькулятора по вышеизложенной методике.
Несмотря на приближенный характер решений, при достаточном числе итераций современные алгебраические методы по точности реконструкции не уступают аналитическим методам, но требуют больших затрат машинного времени.
Общий вывод, который можно сделать из всего сказанного, сводится к тому, что каждый из упомянутых алгоритмов имеет свои положительные качества и свои недостатки. Поэтому было бы неверным отдать предпочтение какому-то одному из них, объявив остальные не заслуживающими дальнейшего развития.