9.2. ФУРЬЕ-АЛГОРИТМ РЕКОНСТРУКЦИИ

Из выражений (9.3) и (9.6) видно, что для функции двух полярных координат имеет место равенство вида

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

Основная трудность в использовании фурье-алгоритма возникает при его применении к реальным исходным данным. Рассмотрим этот процесс поэтапно.

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

Здесь мы используем тот факт, что функция при .

Интеграл в правой части (9.9) необходимо вычислить при некоторых выбранных значениях Замена интеграла римановой интегральной суммой, взятой в точках отсчета, дает

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

Использование в этой формуле такого же обозначения что и в соотношении (8.41) и последующем анализе разд. 8.6, не является случайным. Нетрудно показать, что функция, определенная с помощью равенства (9.11), тождественна введенной в гл. 8 функции Как это уже обсуждалось в разд. 8.6, функцию можно аппроксимировать фурье-образом при но даже в этом интервале пространственных частот она может быть неточной, особенно вблизи границ интервала. На пространственных частотах вероятно, лучше полагать равным нулю, чем пользоваться формулой (9.10).

Мы также должны сделать предположение относительно точек в интервале в котором мы собираемся вычислять значения . С точки зрения выполнения вычислений желательно (по причинам, которые мы здесь упомянем) отобрать величины которые кратны значениям Заметим, что последнее не является реальным ограничением на плотность расположения точек отсчета, поскольку вычисления по формуле (9.10) не изменяются, если мы пожелаем воспользоваться ей при больших [Напомним, что при ].

С учетом приведенных выводов можно утверждать, что первым этапом фурье-алгоритма является вычисление величины

по формуле

при .

Сумма в правой части (9.12) с точностью до представляет собой так называемый дискретный фурье-образ по первой переменной от дискретизированной функции Введем следующее обозначение:

Расчет значений точках при и гораздо менее трудоемок, чем это может показаться на первый взгляд. Заметим прежде всего, что как и в сверточном алгоритме, здесь мы имеем дело с каждым ракурсом отдельно. различных значений при фиксированном и переменной и вычисляются по значениям функции при тех же значениях тип. Так как величины можно вычислить заранее и занести в пямять ЭВМ (так как они не зависят от исходных данных), то вычисления по формуле (9.13) при фиксированном и изменяющемся требуют выполнения операций умножения. (Множитель 2 появляется из-за того, что каждый член суммы вещественных чисел умножается на комплексное число, и это требует двух операций умножения.) Однако имеется гораздо более эффективный способ вычислений по формуле (9.13) причем метод, рассчитывающий данное выражение раздельно для каждого значения Упомянутый алгоритм носит название быстрого преобразования Фурье Его рассмотрение выходит за рамки данной книги, однако здесь необходимо проанализировать характеристики данного алгоритма, поскольку они весьма существенны при оценке эффективности фурье-алгоритма реконструкции изображений.

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

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

рассмотренные в нашем примере величины, экономия вычислений может быть еще более существенной.

Таким образом, с использованием алгоритма БПФ первый этап вычислений по формуле (9.13) при и можно реализовать с относительно малыми затратами машинного времени.

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

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

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

Очевидно, что необходимо найти другой метод расчета величин по значениям в точках с координатами

Один из путей решения этой проблемы состоит в разбиении вычислений по формуле (9.14) на два этапа, а именно:

а) для каждого лежащего в интервале вычисляем при с помощью выражения

б) функцию вычисляем по формуле

Последняя операция включает в себя интерполяцию значений по величинам .

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

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

Трудности, связанные с этим способом вычислений, состоят в том, что точки, в которых рассчитываются значения и точки, в которых необходимо взять значения должны лежать на эквидистантной прямоугольной сетке [например, в центрах ячеек -элементной сетки (разд. 4.1)]. Для наших приложений это не является ограничением, если дело касается выходного изображения, поскольку нам необходимо вычислить величину для находящихся в центрах ячеек. Иначе обстоит дело с входными данными. Первая стадия фурье-алгоритма дает нам значения в точках с координатами лежащих на концентрических окружностях, а не на прямоугольной сетке (рис. 9.2).

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

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

Основная трудность обусловлена тем, что алгоритм БПФ позволяет вычислить лишь двумерный обратный фурье-образ. Соотношение между выходными данными после БПФ и истинными значениями обратного фурье-образа аналогичны соотношению, связывающему функцию определенную соотношением (9.11), и фурье-образ . В частности, если расстояние между точками, в которых производится интерполяция значений равно 6 (рис. 9.2), то выходные данные после БПФ могут быть достаточно точными в прямоугольной области, ограниченной соотношениями Более того, вблизи границ области возможно ухудшение данных вследствие возникновения ложных пространственных частот (разд. 8.6). Таким образом, для получения реконструируемых изображений приемлемого качества нам необходимо выбрать величину 6 достаточно малой, так, чтобы обратная ей величина Последнее требование приводит к двум нежелательным для вычислений последствиям.

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

Рис. 9.3. Билинейная интерполяция в точке .

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

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

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

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

Обозначим через РЁ (соответственно через длину дуги, соединяющую точки (соответственно ). Тогда

Характеристики фурье-алгоритма реконструкции при использовании стандартных проекционных данных, полученные в параллельном пучке, иллюстрируются рис. 9.4 и 9.5, а также табл. 9.1. Были произведены две реконструкции, в одной из которых двумерный обратный фурье-образ вычислялся с помощью алгоритма БПФ в точках, на первой стадии было принято а число лучей в одном ракурсе составляло При второй реконструкции была предпринята попытка уменьшить эффект возникновения ложных частот путем выбора массива из точек для алгоритма обратного двумерного БПФ. Размер увеличенной области реконструкции стал, таким образом, равным 51,888 см (напомним, что размер элиза составляет , поэтому Мы стремились подобрать значение таким, чтобы для получения требуемой точности интерполяции на второй стадии. Мы использовали обоих реконструкциях перед интерполяцией производилось поточечное перемножение вычисленной величины со значением функции обобщенного «окна» Хэмминга при Это позволяет сравнивать данные результаты с результатами, полученными

Таблица 9.1 (см. скан) Мера расстояния между изображениями для реконструированных изображений рис. 9.4

(см. скан)

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