§ 33. Реконструкция в частотной плоскости, спекл-изображения и метод обратной проекции

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

По теореме отсчетов (см. § 10) функция полностью определяется отсчетами на квадратной сетке точек где диапазоны изменения целых чисел достаточно большие для того, чтобы включить все заметно отличающиеся от нуля значения Если то формула (10.6) сводится к равенству

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

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

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

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

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

Кроме того, очень небольшое число угловых координат точек квадратной сетки являются рациональными комбинациями числа, и это показывает, что вычислительная эффективность алгоритма БПФ не может быть использована при расчете суммы в первом соотношении (9.7). Эту трудность можно обойти путем реконструкции изображения с помощью второй зависимости (9.7) с последующим использованием соотношений (9.8), (9.13) и (9.10), которые могут быть реализованы на эквидистантных точках растра по осям направлении о (предполагая, что заданные проекции эквидистантны по углу о и регулярным образом дискретизованы по оси Здесь плохо то, что по-видимому. не существует какого-либо алгоритма, эффективность которого была бы сравнима с эффективностью алгоритма БПФ при вычислении правой части равенства (9.13). Хотя такая процедура интерполяции более предпочтительна, чем рассмотренный ранее подход, основанный на использовании соотношения (10.4), она все же очень сложна в вычислительном отношении.

Конечно, можно использовать и грубые интерполяционные пропс луры. Самая простая из них — интерполяция по методу ближайшего соседа, при которой оценка значения приравнивается величине которая является ближайшей в том смысле, что при определенных значениях величина минимальна для Такая интерполяция может быть осушествлсна очень эффективно. Несколько более медленно реализуется линейная интерполяция. при которой величина усредненное значение двух ближайших соседних значений взятых с весовыми коэффициентами, обратно пропорциональными соответствующим интервалам между точками радиальной сетки При использовании таких интерполяционных процедур удобно написать

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

где истинный отсчет, преобразованная ошибка интерполяции.

Каждая величина занимает один элемент изображения в частотной плоскости, так что она имеет характер двумерной дельта-функции. Из формул (6.13) и (6.20) или, что эквивалентно, формул (12.2) и (12.4) тогда явствует, что в каждое значение вносит свой вклад каждая величина Из указанного выше псевдослучайного характера последней величины следует, что величина фактически подчиняется распределению Релея. (Это еще одно из многих замечаний, разбросанных по всей книге, которые, хотя и не являются строго обоснованными, но тем не менее имеют существенное практическое значение.) Таким образом, ясно, что значения должны перекрывать реконструированное изображение “пятнами”, которые весьма “родственны" спеклам, о свойствах и причинах образования которых подробно говорится в § 34. Если величины в среднем достаточно большие, то “уровень спеклов” может быть исключительно высоким. В обычной медицинской практике ВТ в настоящее время вообще избегают грубой интерполяции в частотной плоскости, поскольку получающиеся в этом случае артефакты обычно неприемлемы.

Если бы интерполяцию можно было перенести с частотной плоскости на плоскость изображения, то грубая процедура интерполяции могла бы быть удовлетворительной, поскольку максимальное искажение соответствовало бы тогда истинным значениям элементов изображения, сдвинутым на половину интервала дисктеризации. Поэтому мы рассмотрим ниже одну из процедур интерполяции в плоскости изображения. Это метод модифицированной обратной проекции, который в настоящее время является стандартным в практике медицинской рентгеновской ВТ. Его практическая реализация иллюстрируется в примере 5. Теоретическое обоснование этого метода удобно начинать с формул (6.11) и (6.13), которые дают

так как и по определению

Поскольку в силу определения координат (см. § 9) мы имеем из формул (33.5) и (9.5) следует, что

где модифицированная проекция под углом которая определяется одномерной сверткой:

причем результат применения одномерного преобразования Фурье к величине т. е.

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

Спектр пространственных частот реконструированного изображения всегда конечен, что явным образом указывалось ограничением (29.1), налагаемым на вариант спектра, который может быть рассчитан на основе известных данных. Это означает, что практические функции фильтра получаются в результате применения процедуры “аподизаиии” (см. § 15) к величине в интеграле (33.9). На детальную разработку соответствующих процедур аподизации были затрачены большие усилия (см. литературу, цитируемую в вводных замечаниях к данной главе). Каждая модифицированная проекция вычисляется по соответствующей заданной проекции либо заданием соответствующей программы [первое соотношение (9.6) с последующим расчетом результата применения одномерного преобразования Фурье к выражению либо, что эквивалентно, путем выполнения операции свертки в выражении (33.8) на специальном процессоре.

Мы отсылаем читателя к третьему — шестому абзацам примера 5, где дается наглядная иллюстрация к операции обратной проекции и поясняются различия между обратным проецированием обычных и модифицированных проекций. Здесь же мы дадим прямую интерпретацию формулы (33.7), которая представляет собой математическое выражение связи, существующей между истинным изображением и его модифицированными проекциями. Последний интеграл в формуле (33.7) можно рассматривать как некое предписание для последовательного вывода (на дисплей) изображения. Заметим, что подынтегральное выражение определено на всей плоскости изображения, поскольку переменная изменяется, вообше говоря, от до Таким образом, при любом конкретном значении амплитуда регистрируется в каждой точке на прямой, идущей параллельно оси на расстоянии от нее. Это означает, что функция распространяется подобно “разрыву” по всей плоскости изображения с “волновым фронтом", всегда параллельным оси . Эта операция называется операцией обратной проекции. Полное реконструированное изображение формируется в результате обратного проецирования всех заданных проекций после их модификации. Поскольку изображение определено только при как это следует из условий (25.1), не имеет смысла пытаться выполнить операцию обратной проекции вне окружности с уравнением

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

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

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

Хотя метод модифицированной обратной проекции был предложен инженерами-электриками, увлекшимися астрономией, в настоящее время он редко применяется в астрономических исследованиях, поскольку в них измеренные данные обычно лежат на эллиптических, а не на прямых линиях в частотной плоскости. Интерферометрическая задача (см. § 25) в обшем случае решается с использованием алгоритма БПФ после выполнения предельно детальной дискретизации в частотной плоскости. Теперь почти стандартной практикой стало устранение возникающих при этом артефактов путем использования алгоритма “очистки” Хегбома (см. § 17) или его соответствующей модификации (см. § 17, последний абзац).