12.5. Фурье-алгоритм восстановления томограмм

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

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

По двумерному фурье-образу  можно найти саму функцию , используя обратное двумерное преобразование Фурье в полярных координатах (12.22). Подставив (12.23) в (12.22), получим

.           (12.25)

Таким образом, в фурье-алгоритме восстановления томограмм вначале вычисляются одномерные фурье-образы  по проекциям  (см. 12.17). Получаем двумерный спектр томограммы в полярной системе координат . Затем выполняем обратное двумерное преобразование Фурье также в полярной системе координат.

Основная трудность в использовании фурье-алгоритма возникает при его применении к реальным данным. Напомним, что проекционные данные имеют дискретный характер, так как вычисляются для дискретных значений ,  и представляют собой массив значений проекций , ;. При реализации фурье-алгоритма по дискретным данным используются процедуры одномерного и двумерного ДПФ или БПФ. Однако эти преобразования осуществляются на прямоугольной сетке отсчетов, а фурье-образы проекций мы получим на полярной сетке (рис. 12.10).

Рис. 12.10. Прямоугольная и полярная сетки отсчетов в пространстве частот

Поэтому для оценки значений  в узлах прямоугольной сетки (где ; ) по известным значениям  в узлах полярной сетки необходима та или иная интерполяция. Здесь можно использовать как простейшие методы интерполяции (например, интерполяция по ближайшему значению), так и метод, реализуемый с помощью БПФ [12.6], который применяется для интерполяции в радиальном направлении.

При использовании БПФ каждую проекцию в пространственной области предварительно дополняют нулями. Затем находят фурье-образ таких проекций. Чем больше нулей добавлено в исходную — точечную проекцию, тем меньше шаг дискретизации по частоте  (почему?) и, следовательно, тем больше число интерполированных значений . Интерполяцию с помощью БПФ иллюстрирует рис. 12.11.

Рис. 12.11. Интерполяция за счет дополнения нулями: а — исходная последовательность, содержащая шестнадцать отсчетов, одиннадцать из которых равны единице, остальные - нулю; б — модуль БПФ исходной последовательности; в — двукратно увеличенная исходная последовательность за счет дополнения нулями; г — модуль ее БПФ

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

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

2) интерполяция значений  на прямоугольной сетке по значениям  на полярной сетке;

3) дискретное или быстрое обратное преобразование Фурье на прямоугольной сетке, которое и дает оценку изображения  также на прямоугольной сетке.

На рис. 12.12 показаны результаты восстановления изображения «Фантом» с помощью фурье-алгоритма. Видно, что в результате применения интерполяции устранено характерное для алгоритмов восстановления томограмм появление артефактов (ложных объектов и линий), которые особенно заметны на краях изображения. Использование 180 проекций и интерполяции с помощью БПФ (число отсчетов в исходной проекции было увеличено в четыре раза за счет дополнения нулями) позволило практически идеально восстановить изображение «Фантом».

Обсудим выбор параметров . Очевидно, что величина шага дискретизации по пространственным переменным ,  и  должна выбираться по теореме Котельникова исходя из величины верхней частоты  фурье-образов. Из теоремы о центральном сечении следует, что одномерный фурье-образ проекции определяет сечение («спицу», рис. 12.10) двумерного фурье-образа восстанавливаемой томограммы. Поэтому верхние частоты у одномерного фурье-образа проекции и двумерного фурье-образа томограммы равны.

Следовательно, шаги дискретизации  и  по пространственным переменным  и  могут быть одинаковыми, при этом число отсчетов . Если используется БПФ при реализации фурье-алгоритма, то дополнительным условием при выборе числа отсчетов  является его кратность степени 2. Увеличение шага дискретизации  в пространственной области приводит к наложению спектров, а уменьшение  (при этом число отсчетов  увеличивается) к возрастанию дозы облучения при получении дополнительных данных и объема вычислений, связанных с их обработкой.

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

Рис. 12.12. Результат восстановления изображения фантом фурье-алгоритмом: а…в - восстановление с интерполяцией по ближайшему значению; г…е - восстановление с интерполяцией за счет дополнения нулями (исходная последовательность была увеличена в 4 раза)

При эквидистантном расположении проекций расстояние между «спицами» (см. рис. 12.10), на которых проекции задают  фурье-образ томограммы, равно , где частота  определяет радиус окружности, на которой расположены отсчеты. Следовательно, интервал между отсчетами пропорционален  и достигает максимума при . При большом количестве проекций . Здесь использовано то, что  при малых углах . Теорема Котельникова требует, чтобы , поэтому . Таким образом, число лучей  зависит от верхней частоты  и размера  восстанавливаемой области, а число проекций  должно быть больше  примерно в полтора раза. На самом деле обеспечить такое большое количество проекций достаточно трудно, так как увеличение  приводит к увеличению времени сбора проекционных данных. Поэтому на практике вместо одномерной интерполяции часто применяют более сложные двумерные интерполяционные процедуры, использующие не только ближайшие соседние отсчеты, расположенные вдоль «спицы», но и отсчеты, расположенные на соседних «спицах». Из рис. 12.10 видно, что низкочастотная область заполнена отсчетами значительно плотнее, чем область высоких частот. Поэтому качество интерполяции значений  на низких частотах будет лучше, чем на высоких. В непрерывном случае этот недостаток компенсируется умножением  на  (см. (12.25)). В фурье-алгоритме с использованием ДПФ или БПФ процедура умножения фурье-образов проекций на , по сути, заменяется процедурой интерполяции. Это приводит к снижению четкости восстанавливаемого изображения . Следовательно, альтернативой фурье-алгоритму может быть алгоритм для восстановления томограмм, использующий дискретный аналог соотношения (12.25).