ГЛАВА 7. ВОССТАНОВЛЕНИЕ ИЗОБРАЖЕНИЙ НА ЭВМ

7.1. БЫСТРЫЕ АЛГОРИТМЫ ЦИФРОВЫХ ПРЕОБРАЗОВАНИЙ

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

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

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

где матрица коэффициентов искомый вектор принятое изображение При для решения системы (7.1) достаточно найти обратную матрицу:

Алгебраические методы можно разделить на прямые и итерационные. Важнейшим из прямых методов является метод исключения Гаусса. В целом, основную идею гауссовых схем исключения можно представить как разложение матрицы на произведение двух треугольных матриц . В таком представлении система записывается в виде и сводится к двум системам, которые легко решаются. Отметим, что если вычислены матрицы то для решения системы требуется порядка операций умножения и деления; для нахождения решения также требуется операций. Однако предшествующие вычисления и требуют соответственно операций умножения и деления, поэтому часто оказывается выгодней использовать вычисление чем [76].

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

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

Дальнейшее сокращение объема вычислений достигается для уравнения свертки, когда значения весовой функции зависят только от разности аргументов. В этом случае ленточная матрица приобретает так называемую теплицеву структуру: все элементы на главной диагонали, а также на любой диагонали, параллельной главной, одинаковы. Для решения системы уравнений с теплицевыми матрицами существуют эффективные численные схемы решения [62]. Так, при треугольном разложении вида для элементов нижней треугольной матрицы справедливо равенство если Элементы этой матрицы на ненулевых диагоналях вычисляют по алгоритму Холецкого:

где Далее решают основное матричное уравнение в два этапа с помощью формул

Будем называть приведенной операцией среднее время выполнения операции сложения действительных чисел.

Машинная операция умножения эквивалентна по времени в среднем полутора операциям сложения. В дальнейшем будем оценивать объем вычислений по приведенным операциям. Лучшие современные мини-ЭВМ выполняют одну приведенную операцию за времена от 200 не до что соответствует приведенных операций в секунду. Количество приведенных операций, требуемых для получения оценки с использованием алгоритма Холецкого равно

Для метода регуляризации А. Тихонова по формуле (2.12) требуется приблизительно приведенных операций. Применительно к алгоритму Холецкого это число составляет приведенных операций. Так как обычно то выигрыш в объеме вычислений оказывается значительным. Для решения систем уравнений высокого порядка вместо гауссовских вычислительных схем, дающих «точное» решение, часто удобнее использовать приближенные итерационные методы, которые в целом несколько экономичнее прямых методов. Их достоинством является простота вычислительных процедур и возможность проведения вычислений с помощью корректирующихся процессов. Одна ошибка в вычислениях у таких процессов не отражается на окончательном результате. Кроме того, преимуществом итерационных методов является то, что одновременно решается только одно уравнение. Поэтому системы большой размерности с помощью итерационных методов могут быть решены на малых ЭВМ, обладающих достаточным объемом внешней памяти.

Из итерационных методов широкое применение нашли метод простой итерации, метод Гаусса — Зейделя и их разновидности [112], (см. §§ 4.1, 4.6). Детальное исследование итерационных методов решения систем линейных уравнений приведено в [30].

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

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

Остановимся теперь на одном из центральных вопросов цифровой обработки сигналов — вопросе вычисления

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

где X — интервал дискретизации функции номер дискретного отсчета функции интервал дискретизации спектра — номер дискретного отсчета спектра. ДПФ переводит последовательность из значений в последовательность также содержащую значений. Формула (7.5) определяет периодическую последовательность чисел с периодом Запишем (7.5) в следующем виде:

где .

Ясно, что для вычисления одного из коэффициентов Фурье в (7.6) требуется операций умножения и столько же операций сложения комплексных чисел, следовательно, для вычисления ДПФ последовательности из отсчетов требуется операций. До сих пор мы рассматривали только одномерные преобразования. Если же учесть, что двумерное изображение упорядочено в кадр квадратного формата со стороной то объект содержит отсчетов. Пусть формат разложения изображения составляет 256X256 отсчетов, что составляет 65 536 точек. Если одно комплексное умножение производится за а это соответствует производительности быстрой мини-ЭВМ, то получим, что время вычисления преобразования Фурье по строкам и столбцам составит более десяти часов. Даже для выполнения простейшей операции линейной фильтрации в этом случае необходимы десятки часов машинного времени.

Разработанные алгоритмы быстрого преобразования Фурье позволяют на несколько порядков уменьшить время вычисления . В настоящее время существует несколько вариантов вычисления быстрого преобразования Фурье — «стандартное» БПФ, алгоритм простых множителей и алгоритм Винограда. Рассмотрим кратко

основные принципы построения этих алгоритмов [34]. Обратимся к (7.6) и представим себе, что общая длина последовательности может быть разбита на произведение двух сомножителей:

Произведем замену индексов:

В новых значениях для выражение (7.6) запишется следующим образом:

где

Теперь процедуру вычисления ДПФ можно свести к следующим этапам:

— вычислению функции

— умножению на

— вычислению

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

Получаем четную и нечетную последовательности а выражение для принимает вид

Если длина последовательности равна 2, то нам необходимо сделать ровно шагов, причем каждый шаг

преобразует ДПФ длины которых у нас длины Смысл о писанного способа вычисления ДПФ заключается в разбиении общей длины последовательности на подпоследовательности до тех пор, пока в подпоследовательностях не останется по одному числу. Эти числа являются первичным массивом для подстановки в систему рекуррентных соотношений. На каждой итерации выполняется операций комплексного умножения и столько же сложений. Учитывая, что всего шагов, получаем что общее число операций умножения и сложения составляет приблизительно Если при прямом вычислении ДПФ общее число операций было равно то выигрыш составит и при больших будет очень велик. Вернемся к задаче вычисления БПФ кадра 256X256 элементов. Если мы попытаемся выполнить два БПФ описанным методом (по строкам и столбцам изображений), то общее число операций равно Если время комплексного умножения равно то время вычисления БПФ кадра составит около 10 секунд. Таким образом, получается огромный выигрыш в объеме вычислений.

Сравнительно недавно появилась новая модификация алгоритма БПФ, называемая также алгоритмом простых множителей [33, 103]. Для ее реализации необходимо, чтобы общая длина последовательности была равна произведению взаимно простых чисел.

В алгоритме на основе метода простых множителей число умножений по сравнению со стандартным БПФ уменьшается еще примерно в 1,2-1,5 раза [34]. Наконец, алгоритм Винограда пользуется представлением БПФ из точек в качестве многомерного ДПФ из точек, где взаимно простые [33]. Численная эффективность алгоритма Винограда примерно в 1,5 раза выше, чем у алгоритма простых множителей, однако лрограммная реализация его оказывается сложнее, чем в стандартном алгоритме

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

последовательности нулевых значений. После этой операции находим преобразование Фурье:

Найдя преобразования Фурье весовой функции, получим;

При обратном вычислении преобразования функции от произведения спектров находим свертку

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

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

Важны также требования к памяти ЭВМ. При обработке изображения элементов и 8 бит информации, на элемент (256 уровней яркости) необходимы, по крайней мере, байт бит) памяти для изображения, байт для весовой функции системы и байт для запоминания действительных и мнимых частей комплексных функций. Сравнительный анализ различных методов обработки для одномерного сигнала длины приведен в табл. 7.1. Максимальные оценки соответствуют применению алгоритмов обычного прямой свертки с вычислением весовой функции через передаточную функцию системы и решению системы алгебраических уравнений методом исключения Гаусса, минимальные оценки — применению стандартных алгоритмов БПФ, прямой свертки с формированием инверсного ядра и решению системы алгебраических уравнений с учетом ленточной структуры матриц. Анализ данных табл. 7.1 показывает, что если ограничиваться линейными

Таблица 7.1. (см. скан)

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

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