13. Неитерационные алгоритмы реконструкции с использованием разложения функции в ряд

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

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

13.1. РАЗЛОЖЕНИЕ ИЗОБРАЖЕНИЯ НА КОЛЬЦЕВЫЕ ГАРМОНИКИ

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

Рассмотрим приведенный в разд. 12.1 критерий квадратичной оптимизации для случая, когда Последний в точности совпадает с методом наименьших квадратов [выражение (6.34)]. Для сохранения общности анализа предположим также, что хотя, как нетрудно проверить, рассматриваемый метод сохраняет все свои наиболее существенные черты и для любой другой диагональной обращаемой матрицы Как было показано в разд. 12.1, искомое решение является решением с минимальной нормой для системы уравнений вида

[соотношения (12.17), (12.18) и (12.20)]. Совместная система уравнений

(13.1) называется системой нормальных уравнений, связанной с задачей минимизации нормы

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

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

Предположим, что имеются такие матрицы размером , а выбор обозначений станет понятен ниже], что можно записать в виде блочно-диагональной матрицы, т.е.

где — нулевая матрица размером При этом обратная матрица (конечно, если она существует) имеет вид

Таким образом, процедура обращения сводится к инверсии каждой из матриц по отдельности, что требует всего операций перемножения. Если то при этом достигается существенная экономия вычислений. Кроме того, при этом для умножения вектора на обратную матрицу необходимы лишь операций умножения, которые можно производить поблочно. В итоге можно сказать, что если представлена в виде матрицы (13.2), то задачу можно свести к частным задачам и тем самым существенно уменьшить объем вычислений на ЭВМ. Для иллюстрации этого утверждения нами был использован метод Гаусса,

однако данный вывод в принципе не зависит от способа решения системы (13.1).

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

и

где

а — те единственные целые числа, которые в интервалах

и удовлетворяют соотношению

[ср. с выражением (6.8)].

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

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

Для существует единственное значение и, лежащее в интервале также единственное значение для которых справедлива формула

Определим при базисное изображение с помощью соотношения

где для и, лежащего в интервале имеем

Кроме того, при

причем значения определяются в соответствии с выражением (13.10).

Почему такая система базисных функций изображений является разумной? Нам необходимо убедиться в том, что любое интересующее нас изображение можно аппроксимировать линейной комбинацией базисных изображений Заметим, что, если выбрать величину достаточно малой для всех лежащих в интервале и и удовлетворяющих условию разумным было бы предполагать функцию приблизительно равной

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

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

функции определенных формулой (13.13). Усечение упомянутого ряда Фурье с сохранением гармоник в интервале дает приближенное значение функции Комбинируя приближенные выражения для [соотношение (13.15)] и для усеченного ряда Фуре, получим аппроксимированное значение функции в виде линейной комбинации элементарных изображений определяемых соотношением (13.11).

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

Для определения элементов матрицы проекций необходимо вычислить значения во всех точках

Для упрощения обозначений при введем функцию определенную как

В данных обозначениях выражение (6.4) примет вид

Пусть для имеем

С использованием соотношений (13.12) и (13.18) выражение (13.17) можно переписать в виде

Дальнейшего упрощения данного выражения можно достичь с помощью метода, весьма эффективного в наших исследованиях. Утверждается, что имеет место соотношение вида

Важность этого соотношения будет» яснее, если при

определить функцию в виде

При этом соотношение (13.20) можно переписать в виде

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

Доказательство того, что соотношение (13.20) является следствием (13.19), необходимо дать отдельно для случаев однако ниже подробно рассмотрен лишь первый случай, а доказательство для второго случая аналогично. Итак, имеем

Рассуждая аналогично и замечая, что получаем

Справедливость утверждения (13.20) доказывается путем сложения левых и правых частей соотношений (13.23) и (13.24).

Получив выражение (13.22) для любого элемента матрицы проекций мы можем вычислить обобщенный элемент матрицы т.е.

Теперь получим следующее свойство ортогональности, элементарное, но громоздкое доказательство которого здесь опушено.

Пусть произвольное положительное целое, произвольные вещественные числа, определяемые для любых целых значений соотношением вида

Пусть также произвольное целое число, а - неотрицательные целые числа, для которых а . Тогда

и

Из приведенных выше соотношений нетрудно видеть, что при имеет место равенство

где если если При выполнении неравенства имеем

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

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

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