Содержание к диссертации
Введение
Глава 1. Томографические методы восстановления изображений с параллельной системой регистрации проекционных данных 16
1.1 Основные определения 17
1.2 Аналитические алгоритмы на основе интегральных преобразованиях в задачах томографии 21
1.2.1 Метод Фурье - синтеза 21
1.2.2 Метод фильтрации и обратного проецирования 23
1.3 Итерационные томографические методы 27
1.3.1 Итерационный алгоритм алгебраической реконструкции 28
1.3.2 Методы максимума энтропии и мультипликативного ART 30
1.3.3 Итерационный алгоритм Гершберга - Папулиса 34
1.3.4 Алгоритм реконструкции на основе разложения в ряд Неймана (NDSL) 44
1.4 Некорректность задач томографии и методы их регуляри зации 46
1.5 Выводы 47
Глава 2. Численная реализация двух итерационных алгоритмов компьютерной томографии с параллельной системой сканирования 49
2.1 Параллельный итерационный алгоритм Гершберга-Папулиса 50
2.1.1 Библиотека элементарных фантомов для моделирования томографического эксперимента 52
2.1.2 Выбор метода интерполяции в алгоритме Г-П 57
2.1.3 Влияние ширины полосы влияния в полосовой интерполяции на сходимость алгоритма Г-П 64
2.1.4 Выбор начальной ширины полосы влияния в полосовой интерполяции и комбинированная интерполяция 68
2.1.5 Влияние весовых коэффициентов на сходимость алгоритма Г-П 70
2.1.6 Двумерная фильтрация проекционных данных в Фурье пространстве 73
2.1.7 Критерии останова итерационного процесса по минимуму невязки 77
2.1.8 Результаты алгоритма Г-П с использованием четвертого критерия останова итерационного процесса по минимуму невязки 81
2.1.9 Одномерная фильтрация проекционных данных в Фурье пространстве 86
2.2 Параллельный итерационный алгоритм Шеппа - Логана и численные результаты 92
2.2.1 Сравнение итерационного алгоритма Г-П и классического алгоритма Шеппа - Логана 92
2.2.2 Итерационный алгоритм NDSL 95
2.2.3 Выбор параметра релаксации для NDSL 96
2.2.4 Реконструкции томограмм, полученные итерационными алгоритмами Г-П и NDSL 101
2.3 Решение задачи стеганографии с помощью итерационного отделения возмущения на синограмме 105
2.3.1 Алгоритм отделения возмущения от проекционных данных 106
2.3.2 Эксперимент по отделению возмущения от фонового изображения 108
2.3.3 Решение задачи стеганографии с помощью итерационного отделения возмущения от синограммы 115
2.4 Выводы 124
Глава 3. Веерный итерационный алгоритм Гершберга - Папулиса 125
3.1 Двумерная задача томографии в веерной постановке 127
3.1.1 Примеры деформации 131
3.1.2 Две модификации итерационного веерного алгоритма Гершберга - Папулиса 133
3.1.3 Численная реализация итерационного веерного алгоритма Гершберга - Папулиса 136
3.1.4 Сравнение результатов итерационного алгоритма Г- П для параллельной и веерной систем сбора данных 156
3.1.5 Численные результаты веерного Г-П при увлечении разрешения томограмм 162
3.2 Выводы 167
Заключение 169
Приложение I 180
Приложение II 184
- Метод фильтрации и обратного проецирования
- Влияние ширины полосы влияния в полосовой интерполяции на сходимость алгоритма Г-П
- Эксперимент по отделению возмущения от фонового изображения
- Две модификации итерационного веерного алгоритма Гершберга - Папулиса
Введение к работе
Актуальность работы. Во многих областях науки, таких как медицина, геофизика, астрофизика, промышленная дефектоскопия, диагностика плазмы и других, возникает проблема определения внутренней структуры объектов. Для решения данной задачи во многих случаях являются неприемлемыми прямые методы исследования, связанные с разрушением объекта, поэтому создаются специальные системы получения данных. Физический принцип этих систем состоит в использовании воздействия, представляющего собой физический процесс произвольной природы (излучение, волновое поле и.т.д.), и последующей регистрации отклика этого процесса на объект.
В некоторых важных случаях восстановление изображений является вообще практически единственным средством получения информации. Это относится, например, к томографии - получению изображения по набору проекций. Наиболее распространена томография в медицинской диагностике. Однако различные методы вычислительной томографии, позволяющие исследовать внутреннюю структуру объекта не разрушая его, применяются сейчас во многих областях, таких как электронная микроскопия, биохимия, физика Земли, радиоастрономия, исследования океана и космоса.
Вычислительная томография применима в тех случаях, когда внутренняя структура объекта может быть исследована с помощью некоторого вида излучения, которое распространяется с интенсивностью, убывающей по формуле:
Л До = ехр{- / g(x)dx},
где /о - начальная интенсивность излучения луча L, до его прохождения через исследуемый объект, а Д - интенсивность после прохождения. Здесь д(х) - коэффициент ослабления. Решая задачу восстановления функции д(х) можно восстановить требуемые внутренние характеристики объекта. Например, при рентгеновском просвечивании коэффициент поглощения связан с плотностью биотканей, а при диагностике плазмы - с пространственным распределением температур, концентраций электронов и ионов.
Так как не существует точных формул восстановления по конечному набору проекций, решение задачи восстановления на практике ищется в виде приближенного решения. Для задач, когда число ракурсов мало, необходимым условием для нахождения такого решения является применение априорной информации о восстанавливаемом объекте. Использование методов регуляризации при наличии шумовых компонент в проекционных данных способствует нахождению устойчивого решения задачи. Так как в каждом конкретном случае имеем свою постановку задачи, то основной проблемой становится выбор наилучшего алгоритма для ее решения. Основными критериями качества алгоритма рекон струкции являются качество восстановленного изображения и время его получения.
Целью работы является разработка и исследование новых итерационных методов в задачах малоракурсной вычислительной томографии, а также их применение в цифровой обработке изображений; создание и реализация интерактивного диалогового программного комплекса на универсальной ЭВМ.
В соответствии с поставленной целью в диссертационной работе решаются следующие задачи:
- модернизация и исследование итерационного алгоритма Гершберга -Папулиса для параллельной геометрии сбора данных. Введение в алгоритм новых, вспомогательных параметров, способных увеличить скорость сходимости итерационного процесса и уменьшить итоговую погрешность реконструкции. Подробное исследование вопросов регуляризации;
- для сравнения с алгоритмом Гершберга - Папулиса создание второго итерационного алгоритма для параллельной геометрии, направленного на применение в задачах малоракурсной томографии (разложение обратного оператора Радона в ряд Неймана). На основе полученного алгоритма решение задачи стеганографии (скрытие одного изображения в другом);
вывод теоремы о центральном сечении для веерной геометрии сканирования с помощью деформирующего преобразования и ее численная проверка на тестовых моделях. На основе полученной теоремы: создание, разработка и исследование итерационного алгоритма Гершберга Папулиса для веерной системы сбора данных. Сравнение между собой всех исследуемых в диссертации итерационных алгоритмов; - создание удобного в применении пользователем ЭВМ интерактивного диалогового приложения, в котором существует возможность проводить вычислительный эксперимент для любой из поставленных в диссертации томографических задач.
Методы исследования. Основные результаты работы получены с использованием интегральной геометрии PI вычислительных методов, а также методов цифровой обработки изображений и математического моделирования. Для численного моделирования и программной реализации разработанных алгоритмов использовались методы прикладного программирования.
Научная новизна работы состоит в следующем:
1. Разработана модифицированная версия алгоритма Гершберга - Папулиса для параллельной системы сбора данных. Посредством применения нового типа интерполяции (комбинированной) и применения весовых множителей удается ускорить сходимость и получать более точные результаты реконструкции. В итерационный процесс введены новые критерия его останова.
2. Представлена новая теорема о центральном сечении применительно к веерной геометрии сбора данных, на основе которой становится возможным перенос многих алгоритмов с параллельной геометрии на веерную, в частности впервые создан итерационный алгоритм Гершберга - Папулиса для решения томографической задачи с веерными проекционными данными.
3. Создан итерационный алгоритм отделения помехи от синограммы на основе разложения обратного оператора Радона в ряд Неймана, с помощью которого исследована задача стеганографии применительно к томографической постановке. Впервые осуществлено скрытие одного изображения в другом в пространстве Радона, а затем его итерационное восстановление.
Основные положения, выносимые на защиту:
1. Предложены модифицированные томографические методы для восстановления изображений с параллельной системой регистрации проекционных данных, направленные на задачи с малым числом ракурсов.
2. Разработан новый итерационный алгоритм Гершберга - Папулиса для веерной геометрии сбора данных.
3. Решена задача стеганографии с помощью математического аппарата томографии.
Практическая ценность работы. Предложенный в диссертации итерационный алгоритм Гершберга - Папулиса может быть применен для обработки экспериментальных данных в физической томографии с малым числом параллельных проекций. Впервые предложен новый алгоритм Гершберга - Папулиса для обработки веерных экспериментальных данных.
Достоверность результатов полученных результатов и выводов подтверждается анализом разработанных численных алгоритмов и проведением численных экспериментов.
Апробация работы. Основные результаты диссертации докладывались и обсуждались на международной конференции "Математические методы в геофизике - ММГ" (Новосибирск, 2003 г.), международном симпозиуме "International Symposium on Computed Tomography and Image Processing for Industrial Radiology" (Берлин, Германия, 2003 г.), международной конференции "Perspectives in Inverse Problems" (Хельсинки, Финляндия, 2004 г.), международной конференции "Applied Inverse Problems" (Сайренсэстэр, Англия, 2005 г.), 4-ом международном конгрессе "4h World Congress on Industrial Tomography" (Эйзу, Япония, 2005 г.), международной конференции "Review of Progress in Quantitative NDE" (Брунсвик, США, 2005 г.), международной конференции "Обратные и некорректные задачи" (Новосибирск, 2007 г.), конференции молодых ученых ИВМиМГ СО РАН (Новосибирск, 2008 г.), международной конференции IEEE Region 8 Intl. Conf. SIBIRCON 2008 (Новосибирск, 2008 г.), международном конгрессе Nuclear Science Symposium, Medical Imaging Conference IEEE (Dresden, 2008г.), а также на семинарах лаборатории обработки изображений ИВМиМГ СО РАН, ИВТ СО РАН, ИМ им.С.Л. Соболева СО РАН и ИТПМ им. С.А. Христиановича СО РАН. Результаты исследований, проводимых в рамках работы над диссертацией, отмечены первой премией конференции молодых ученых ИВМиМГ СО РАН (2008г.).
Личный вклад автора. Основные научные и практические результаты диссертации получены автором лично. Вывод теоремы о центральном сечении для веерной геометрии принадлежит д.ф.-м.н. В.В. Пикалову и д.ф.-м.н. В.П. Голубятникову, автору принадлежит ее численная проверка. Из печатных работ, опубликованных диссертантом в соавторстве, в диссертацию вошли только те результаты, в получении которых он при нял непосредственное творческое участие.
Публикации. Результаты исследований по теме диссертационной работы опубликованы в 11 печатных работах.
Структура и объем работы.
Диссертация состоит из введения, трех глав, заключения, списка литературы из 69 наименований и двух приложений. Содержание основного текста диссертации изложено на 170 страницах, содержит 74 иллюстрации, 2 таблицы и 2 блок-схемы.
Основное содержание работы
Во введении обоснована актуальность проблем, рассматриваемых в диссертации, определены основные цели и задачи исследования, показана его научная новизна и практическая ценность, сформулированы выносимые на защиту положения и представлен краткий обзор содержания работы.
Первая глава посвящена постановке задачи вычислительной томографии, а также обзору и анализу методов, используемых для ее решения. В разделе 1.1 даются основные определения, такие как понятие томографии, систем сбора данных, методов реконструкции, преобразований Фурье и Радона. В разделе 1.2 дается определение теоремы о центральном сечении в параллельной системе сбора данных и описаны основные аналитические алгоритмы томографии: метод Фурье - синтеза и метод фильтрации и обратного проецирования. В разделе 1.3 внимание заострено на малоракурсных томографических постановках и методах решения данных задач. Подробно рассмотрены следующие итерационные томографические методы: итерационный алгоритм алгебраической реконструкции (ART), метод максимума энтропии (MENT) и мультепли-кативного ART (MART), итерационный алгоритм Гершберга- Папулиса (Г-П) и алгоритм реконструкции на основе разложения обратного оператора Радона в ряд Неймана (NDSL, Neumann Decomposition with Shepp - Logan filter). Последние два итерационных метода рассмотрены более подробно, так как их изучению посвящена данная диссертация. Приводится доказательство сходимости итерационного процесса для алгоритма Г-П. В разделе 1.4 рассматривается проблема некорректности постановок в томографических задачах и необходимость применения методов регуляризации.
В главах 2-3 содержится материал, представляемый к защите.
Вторая глава посвящена исследованию и численной реализации итерационных алгоритмов Гершберга - Папулиса PI алгоритма NDSL для параллельной геометрии сбора данных в задачах с ограниченным числом ракурсов. В разделе 2.1 приведены результаты исследования алгоритма Гершберга - Папулиса для параллельной геометрии, и исследуются разнообразные его модификации. Рассматривается проблема выбора процедуры интерполяции в методе и предложен новый способ - комбинированная интерполяция. Исследуется сходимость итерационного процесса Г-П при наличии шумовых компонент в проекционных данных. Особое внимание уделено проблемам регуляризации и фильтрации. В разделе 2.2 приводятся результаты по второму итерационному алгоритму для малоракурсной томографии - методу реконструкции на основе разложения обратного оператора Радона в в ряд Неймана (NDSL). Демонстрируются результаты сравнения алгоритмов Г-П и NDSL. Ана логично экспериментам с алгоритмом Гершберга - Папулиса исследовано влияние фильтрации на зашумлснные проекционные данные и выбор релаксационного параметра для управления сходимостью итерационного процесса.
В разделе 2.3 предметом рассмотрения является задача стеганографии (скрытие одного изображения в другом) применительно к томографической постановке. Описывается и реализуется численно алгоритм отделения возмущения от фона, где налагаемым возмущением является скрываемый текст или изображение. Алгоритм отделения возмущения от фона является итерационным и базируется на методе разложения в ряд Неймана, описанным в предыдущем разделе. Приведенные в разделе 2.3 численные результаты демонстрируют успешное применение томографического аппарата к задаче стеганографии.
Третья глава посвящена новому итерационному алгоритму Гершберга - Папулиса, для веерной системы сбора данных. В разделе 3.1 главы вводится теорема о центральном сечении для веерной геометрии. В ее основе лежит хорошо известное в проективной геометрии нелинейное координатное преобразование, переводящее задачу веерной томографии в задачу параллельной томографии, но для деформированного пространства, причем для каждого ракурса наблюдения такая деформация - своя. На основе данной теоремы реализуется численный алгоритм Гершберга - Папулиса для веерной геометрии для малоракурсной томографии. Алгоритм представлен в двух модификациях, между которыми проведен сравнительный анализ. Использовались интерполяции различных порядков точности для улучшения конечного результата реконструкции, такие как бикубическая и интерполяция эрмитовыми сплайнами. Исследовано поведение алгоритма на зашумленных проекциях, а также применение к ним сглаживания. Сравнены между собой алгоритмы Гершберга - Па-пулиса и NDSL для параллельной геометрии и алгоритм Гершберга -Папулиса для веерной. Представлены рекомендации по ускорению вычислительного времени работы алгоритма.
В приложении I описаны созданные прикладные программы для моделирования томографического эксперимента. Весь пакет программ объединен в одно интерактивное диалоговое приложение (FBTA, Fan-Beam Tomography Algorithms), созданное на языке Matlab.
В приложении II представлены две блок-схемы: пакета FBTA и алгоритма Г-П для параллельной геометрии сбора данных. В блок - схемах наглядно объясняется система работы пакета FBTA и томографических алгоритмов.
Автор выражает благодарность своим научным руководителям: доктору физико-математических наук Валерию Владимировичу Пикалову и доктору технических наук профессору Валерию Павловичу Пяткину, за постоянное внимание к работе, практические рекомендации, организационную помощь и моральную поддержку, а также всему коллективу лаборатории обработки изображений ИВМиМГ СО РАН за полезные обсуждения результатов работы на лабораторных семинарах и отчетных сессиях.
Метод фильтрации и обратного проецирования
Рассмотрим метод фильтрации и обратного проецирования (ФОГТ) [48], так как он является наиболее применяемым алгоритмом в медицинской томографии с параллельной системой сбора данных. В данной диссертации ФОП использовался в итерационном алгоритме разложения в ряд Немана (NDSL), а также как один из способов получения начального приближения искомой функции в итерационном процессе. Пусть двумерное обратное Фурье - преобразование для нахождения неизвестной функции д(х,у) представлено в виде (1.6). В частотной области декартовы координаты (их, иу) выражаются через полярные координаты (ир. в) следующим образом.
Уравнение (1-Ю) представляет собой операцию фильтрации проекционных данных, где \vp\ импульсная характеристика фильтра (частотный отклик). Таким образом, fl{vp) является фильтрованной проекцией. Уравнение (1.9) отвечает за операцию обратного проецирования. Результат фильтрации проекций, а затем их обратное проецирование, дает искомую функцию д(х, у).
Поскольку ф попадает в класс обобщенных функций, то для ее численной замены производят замену точной сингулярной функции на некоторую приближенную, в итоге появляется проблема выбора фильтра ф(р).
Особую известность при обработке алгоритмов медицинской и физической реконструктивной томографии приобрел фильтр Шеппа - Логана [64], который в дискретном виде можно записать так: де/1р) = (j -1) к = 0 ±1 ±2 - (1Л4) где hp - шаг по прицельному параметру р.
По рассмотренным аналитическим алгоритмам можно сделать следующие выводы. При выводе алгоритма Фурье - синтеза, Фурье - преобразование было записано в декартовой системе координат, а при выводе ФОП - в полярной. Несмотря на общую основу, практическая реализация этих алгоритмов различна. Алгоритм ФОП оказался наиболее применимым методом в компьютерных томографах. Его преимуществом является то, что операции свертки и обратного проецирования для каждого ракурса могут выполняться независимо, а результирующее изображение представляет собой сумму изображений, полученных для каждого из ракурсов.
Большинство аналитических методов реконструкции предполагает наличие достаточно большого числа проекционных данных. В случае недостаточного их количества, эта группа методов является неспособной к получению удовлетворительных результатов. Поэтому для задач восстановления при ограниченном (малоракурсная томография) или неполном наборе проекций используются различные итерационные алгоритмы.
Итерационные методы имеют ряд преимуществ перед аналитическими алгоритмами, в частности: 1) возможность создания нелинейных алгоритмов с применением априор но известных свойств искомого решения - априорной информации. Априорная информация может включать в себя положительность решения, его пространственную и частотную ограниченность, гладкость и.т.д. 2) при реализации этих методов возможна работа в интерактивном режиме, позволяющая сделать компромиссный выбор между качеством восстановления и временем обработки.
Итерационные методы характеризуются медленной сходимостью, однако разнообразные их модификации позволяют в той или иной степени решать эту проблему. Перейдем к рассмотрению итерационных алгоритмов. На методах, которым посвящена данная диссертация, остановимся более детально.
Решение системы (1.17) с прямым обращением матрицы А может быть записано как: д = A f, однако такое решение требует огромных вычислительных затрат и в большинстве случаев для томографических задач не может быть использовано. Проблемой является и то, что число неизвестных обычно больше, чем имеющееся число уравнений. Многообразие способов решения СЛАУ (1.17) дано, например, в [16].
Впервые применение алгебраических методов реконструкции было описано в работе Гордона и др. [45] в 1970 г. Впоследствии было установлено, что предложенный способ совпадает с алгоритмом решения СЛАУ (1.17), предложенным Качмажем [30] в 1937 г. Гордон и др. предложили итерационную процедуру решения (1.17), в которой используется некоторое начальное приближение вектора элементов изображения, д Є RJ. Полученная итерационная процедура была названа методом ART (Algebraic Reconstruction Technique) и суть его заключалась в следующем.
Рассмотрим (1.17) в виде: /j = ]Cj=iau5r Здесь пиксели изображения gj индексируются по столбцам, а элементы синограммы /; индексируются последовательно по лучам в каждой проекции идущей одна за другой. На /с-й итерации текущее приближение дк уточняется (корректируется), что дает новое приближение дк+1. При этом рассматривается только один луч, например г-й, а изменению подвергаются только те компоненты вектора дк, которые соответствуют пикселям, пересекаемым данным лучом. Величина невязки между измеренным значением /г-и приближенной величиной проекции J2j=i aij9j полученной из текущего приближения дк, перераспределяется между пикселями, расположенными вдоль г -го луча, пропорционально их весам а - в луче.
Влияние ширины полосы влияния в полосовой интерполяции на сходимость алгоритма Г-П
Начальная полуширина Sw выбрана равной 1.8 для модели No. 200 (см. Рис. 2.11) и 1.6 для модели No. 216 (см. Рис. 2.12). Множитель Sa = 0.8 одинаков для обоих экспериментов. Исходя из полученных результатов видно, что "замораживание" полос на малых номерах итераций, ведет к расходимости итерационного процесса, а следовательно, и к накоплению погрешности реконструкции. Поэтому, чтобы не допустить расходимости алгоритма и получить наименьшую итоговую ошибку реконструкции, на каждой итерации необходимо уменьшение ширины полосы. В процессе сужения ширины полосы происходит уточнение процедуры полосовой интерполяции. То есть с каждой итерацией в полосу, где производится перенос значений с радиальной линии, попадает меньшее количество точек. Зависимость количества точек, попадающих в полосу, от номера итерации показана на Рис. 2.136, а сужающаяся в итерациях ширина полосы приведена на Рис.2.13а. На Рис. 2.14 продемонстрированы результаты реконструкции модели No. 200 (см. Рис 2.14а), полученные алгоритмом Г-П с использованием полосовой интерполяции. На Рис. 2.146 представлена реконструкция на 120-ой итерации с постоянным уменьшением Sw в итерациях по правилу (2.49), Ai = 6.5%. Реконструкция с "замораживанием" ширины полосы на 3-ей итерации приведена на Рис. 2.14в, Ді = 27.4%.
Это обусловлено тем, что при узкой ширине полосы 2Sw, в процедуру переноса значений не попадают с первых итераций многие близкие к радиальной линии точки. И наоборот, при завышенном значении ширины полосы (Sw — 2.5), наблюдается медленная сходимость итерационного процесса на первых десятках итераций. Причиной этого служит захват удаленных точек на первых итерациях алгоритма, когда ширина полосы является излишне широкой и для этих точек интерполяция ближайшего соседа является слишком грубой. Исходя из полученных результатов, начальная полуширина полосы была выбрана равной 1.8. Оперируя значениями Sw = 1-6 — 1.8 удалось добиться наилучшей сходимости итера ционного процесса и минимальной погрешности на последних итерациях.
Однако эксперимент с комбинированной интерполяцией дает более точный результат, чем полосовая интерполяция с оптимально выбранным начальным значением Sw С эмпирически подобранными параметрами ІІ и Sw Для билинейной и полосовой интерполяций было обнаружено, что при использовании комбинированной интерполяции есть возможность уменьшить ошибку реконструкции на первых итерациях на 3 -10%. На последних итерациях кривая ошибок реконструкции для комбинированной рштерполяции почти сливается с кривой для полосовой, но здесь комбинированная интерполяция проявляет себя лучше полосовой на 1-2%.
После серии проведенных экспериментов по реконструкции разных элементарных моделей был сделан вывод о целесообразности использования процедуры комбинированной интерполяции в алгоритме Г-П, так как данный вид интерполяции проявил себя лучше полосовой.
Следующий эксперимент проведен с целью уменьшить итоговую погрешность реконструкции и увеличить скорость сходимости итерационного процесса. Это осуществлялось посредством введения в алгоритм Г-П весовых коэффициентов. С помощью данных параметров осуществлялся учет значений, полученных в Фурье пространстве на предыдущей итерации. Весовые множители Wprev и Wexp участвуют в алгоритме следующим образом: Wprev учитывает значения полученные в процедуре ин тєрполяции на предыдущем шаге, a Wexp учитывает экспериментальные данные на текущей итерации.
Зависимость ошибок реконструкции Ді(%) в зависимости от числа итераций п для модели No. 216, проверка влияния весовых параметров на сходимость алгоритма с параметрами: Nx = Ny = 128, Np = 128, К = 7; кривая No.l: (W , W v) = (0,1); No.2: (0.2,0.8); No.3: (0.5,0.5); No.4: (0.8,0.2). алгоритм Г-П весовых параметров и показана зависимость ошибок реконструкции от числа итераций на моделях No. 201 и 216. Было рассмотрено четыре крайних случая, когда веса были равны (И .еи, Wexp) = (0,1), (0.2,0.8), (0.5,0.5) и (0.8,0.2). Проведенные эксперименты показали, что применение весов (Wprev, Wexp) = (0.5, 0.5) дает результат в большинстве случаев лучший, чем при исходном значении весов (И »б, WeXp) (0,1) (когда значения интерполяции полученные на предыдущем шаге не рассматриваются).
Как уже упоминалось в Главе 1, задача томографического восстановления изображения является некорректно поставленной, так как содержит в проекционных данных случайные шумовые компоненты и количество данных ограничено. Поэтому данная задача требует применение некоторых приемов регуляризации, например, частотной фильтрации.
Для исследования устойчивости алгоритмов к ошибкам во входных данных были использованы две модели случайных шумов в проекционных данных. Первая - нормально распределенный шум с нулевым средним и переменной дисперсией, имеющей значение к, в процентах от вс-73 личины проекции в рассматриваемой точке (тип 1). Вторая - тот же шум, но с постоянной дисперсией, равной к, в процентах от величины максимума проекции (тип 2). Далее в экспериментах в основном используется 1-ый тип шума.
Перед тем как перейти к демонстрации результатов, полученных с использованием процедуры регуляризации, введем несколько критериев останова итерационного процесса по минимуму нормы невязки. Это является необходимым условием, так как в реальном эксперименте точное решение неизвестно, следовательно, нет возможности сравнения истинной томограммы и ее реконструкции, поэтому нет возможности отслеживать поведение кривой реконструкции и задавать число итераций, опираясь на результаты по данной кривой. В описанных ранее экспериментах использовалось 120 итераций, так как примерно к 120-ой итерации кривая ошибки реконструкции выходила на насыщение, таким образом можно было прервать итерационный процесс в реальном времени. В реальном эксперименте значения погрешности реконструкций отсутствуют, но известны значения нормы невязок. Используя зависимость норм невязок от номера итерации, можно принимать решения по останову итерационного процесса, так как не всегда процесс обладает сходимостью и требуется останов задачи реконструкции до того, как началась расходимость.
Эксперимент по отделению возмущения от фонового изображения
В данном разделе в численном эксперименте исследуется итерационный алгоритм (2.61) по разделению суммы двух изображений. Результат наложения возмущения на фоновое изображение показан на рис. 2.42г (возмущение помещено в левый верхний угол). На этом рисунке размер возмущения был небольшой, поскольку использованное изображение Lena имело размер 64 х 64. Модельная томограмма вычислялась внутри единичного круга, ее размерность - 256 х 256. Проекции вычислялись в диапазоне углов от 0 до 7Г с изменением прицельного параметра в диапазоне [—1,1]. Результаты отделения возмущения от фонового изображения за 20 итераций показаны на рис. 2.43. Отделение возмущения от фона произошло неплохо, поэтому оставшиеся от разделения артефакты на фоновом изображении обладают достаточно малой амплитудой. Зависимость погрешностей разделения Afw и Aw (в %) от номера итерации изображена на рис. 2.43в. В этом эксперименте параметры итерационного процесса (2.61) следующие: а = 1,/? = 1. Координаты вставки возмущения (номера соответствующих индексов в матрице фонового изображения, в процессе реконструкции они полагались известными): гго = 30, т/о = 15 амплитуда возмущения - hdist = 0.1. В алгоритм на каждой итерации вносится априорная информация в виде положительности возмущения и фонового изображения.
По качеству отделенного в итерациях изображения Lena (см. рис. 2.43а) видно, что некоторые участки изображения не были отделены от фона (верхняя часть изображения). На восстановленном фоновом изображении (см. рис. 2.436) слабо видны те части возмущения, которые были не полностью отделены от фона в итерационном процессе. Зависимости ошибок Afw и Aw на рис. 2.43в демонстрируют медленную сходимость и высокую итоговую погрешность реконструкции возмущения, Aw = 61%. Несмотря на столь значительную погрешность изображение восстановленного возмущения, представленное отдельно на рис. 2.43д, демонстрирует то, что его основная часть восстановилась неплохо, а погрешности накопились по периферии изображения, на областях
Реконструкция возмущения W; б) реконструкция фоновой сииограммы fw\ в) зависимость погрешностей реконструкции (в %) от числа итераций для фона (кривая No.l) и возмущения (2); г) реконструкция томограммы по восстановленной синограмме; д) восстановленное возмущение размерностью Nx — Ny = 64. с малой амплитудой, слабо отличающихся от окружающей возмущение синограммы. При этом качество восстановления синограммы достаточно высокое (кривая 1, рис. 2.43в), что показывает и низкий уровень артефактов на восстановленной томограмме (рис. 2.43г).
Чтобы исследовать влияние шумовых компонент на процесс разделения фонового изображения и возмущения, вводился 1, 3 и 5% шум в итоговые изображения, результаты которого приведены на Рис. 2.43. Набрасыванию шума подвергается суммарное изображение (см. Рис. 2.46 (а-в)). Результаты выделения возмущения продемонстрированы на Рис. 2.46 (ж-и). Внешний вид полученых возмущений и высоких величин А0 (см. Рис. 2.46 (г-е)) указывает на отрицательное воздействие шума в эксперименте, даже при его весьма низком уровне (к — 1%).
Расположим скрываемое возмущение в центре фонового изображения (XQ = 90, уо = 60) и увеличим его размерность до Np — 128. К = 180. Фоновое изображение при этом сохранено с размерностями Np = 256, К = 360. Амплитуда добавляемого возмущения равна hdist = 0.1. Теперь возмущение находится в центре фонового изображения. Разделение двух изображений дало улучшенный результат, восстановленное возмущение в виде изображения Lena (см. рис. 2.49е) является более качественным (Aw — 28%, рис. 2.50) и обладает лучшей четкостью из-за увеличения размерности. Восстановленная томограмма на рис. 2.49д также свидетельствует о качественном разделении, с низким уровнем амплитуды артефактов.
Для исследования влияния шумовых компонент в задаче стеганографии (с последним набором параметров), вводился 1, 3 и 5% шум в суммарное изображение (см. Рис. 2.52 (а-в)). Результаты выделения возмущения продемонстрированы на Рис. 2.52 (ж-и). Судя по внешнему виду полученых при отделении возмущений, видно что уже 3% уровень шума и выше вносит в эксперимент ошибки, которые сказываются на качестве восстановленного возмущения.
Изображение с фоном в виде синограммы от модели и наложенным на него возмущением (синограммы от изображения Lena) Np = 128, К = 180, hdist = 0.4; б) восстановленное возмущение; в) восстановленный фон; г) восстановленная томограмма; д) зависимость погрешностей реконструкции от числа итераций для фона (1) и возмущения (2); е) результат обратного преобразования Радона от восстановленной синограммы, Nx = Ny = 128.
Суммарные изображения с добавлением 1, 3 и 5% шума (тип 1), где фоном является модельная синограмма размерностью Np = 256, К = 360 и возмущением на ней - синограмма от изображения Lena размерностью Np = 128, А" = 180, hdist = 0.4 (а-в); зависимость погрешности реконструкции (в %) от числа итераций для фона (1) и возмущения (2) (г-е); востановленные возмущения (ж-и).
По полученным результатам можно сделать вывод, что разделение двух изображений в пространстве Радона может быть выполнено наилучшим образом с помощью алгоритма (2.61), если выполнено несколько условий: 1) возмущение в виде синограммы для скрываемого изображения должно иметь меньшую размерностью, чем фоновое изображение, - необходимо, чтобы не все фоновые проекции были испорчены налагаемым возмущением (систематической ошибкой); 2) возмущение располагается в центре фонового изображения; 3) амплитуда возмущения составляет не менее 10-20% от максимума фонового изображени; 4) уровень случайного шума в экспериментальных данных составляет не более 1%.
Соблюдая данные условия, можно решать разнообразные задачи стеганографии предложенным итерационным методом. В данной работе было рассмотрено скрытие изображения в виде преобразования Радона от него, однако существуют большое количество методов, каким образом можно "спрятать"изображение. Например, одним из способов может служить скрытие изображения на нескольких масштабах вейвлет - преобразования фоновой синограммы. При этом номера масштабов служат дополнительными управляющими параметрами процессов кодирования и декодирования.
Две модификации итерационного веерного алгоритма Гершберга - Папулиса
Численная реализация веерного алгоритма Г - П была осуществлена двумя различными способами (первая и вторая модификации) [22]. Опишем подробно каждую из модификаций, для этого будем использовать операторное представление координатных преобразований (3.65), (3.66).
С помощью последовательного применения указанных операторов можно перейти к координатам (и, v), затем двумерное прямое преобразование Фурье дает возможность перейти к частотным координатам (уи, i/v), в пространстве которых может быть использована теорема о центральном сечении. После чего, как и в алгоритме Г-П для параллельной геометрии, коэффициенты Фурье заменяются на известные полученные применением одномерного Фурье-преобразования к исходным проекционным данным. Отметим, что подобная замена проводится в повернутой на угол Р системе координат, в которой текущий угол веерного наблюдения всегда равен нулю. Как следствие, замена спектра в Фурье-пространстве необходима лишь на оси абсцисс ии, и требует лишь одномерной интерполяции. Обратное двумерное Фурье-преобразование дает функцию в координатах (u,v), и, чтобы вернуться к исходным координатам (ж, у), необходимо осуществить обратную деформацию и поворот на угол /5.
Как видно из таблицы 3.1, погрешность интерполяций меняется в зависимости от угла поворота объекта, а также от расстояния до фокальной точки (D). На каждом ракурсе такое количество интерполяций дает высокую итоговую погрешность реконструкции. Происходит это в следствии накопления ошибок от интерполяции в итерациях. Чтобы избавиться от нескольких интерполяций, а следовательно уменьшить погрешность реконструкции, была разработана еще одна модификация веерного алгоритма Г-П (модификация 2). Целью служило избавиться от лишних интерполяций с помощью уменьшения числа поворотов объекта.
Рассмотрим более подробно виды интерполяций применяемые в модификациях 1 и 2. Одномерная интерполяция, которая входила только в модификацию 1 (блок прямой и обратной деформации) осуществлялась двумя типами интерполяций: линейной и эрмитовыми сплайнами. В диссертации эрмитовы сплайны использовались в модификации 1 в одномерной интерполяции, а бикубическая интерполяция использовалась во второй модификации в двумерном случае. Для интерполяции сплайнами использовались встроенные в пакет Matlab функции [8].
Эксперименты, описываемые далее, основаны на результатах, полученных в ходе исследования итерационного алгоритма Г-П на параллельной геометрии в предыдущей главе. Большая часть модификаций при реализации алгоритма Г-П для параллельной геометрии в главе 2 была реализована для дальнейшего их применения в алгоритме Г-П для веерной геометрии. К основным модификациям, улучшающим скорость сходимости и качество реконструкции, относятся: выбор ширины пропускания Sw в полосовой интерполяции, множителя Ss и периода то, введение весовых множителей Wi, Wjj, выбор начального приближения томограммы до использование априорной информации, включающую в себя положительность, гладкость и пространственную ограниченность томограммы, а также ограничение по частоте Найквиста.
В эксперименте включена априорная информация о положительности томограммы и зануление спектра вне частоты Найквиста. Здесь результаты применения одномерной линейной интерполяции в процедуре деформации томограммы показаны на кривых No. 1-3, а интерполяция эрмитовыми сплайнами продемонстрирована на кривых No. 4-6. Поворот изображсния осуществлялся с использованием процедуры билинейной интерполяции. Интерпретируя полученные результаты, можно сделать одно Рис. 3.57: Реконструкции модели No. 224 на 20-ой итерации, априорная информация включала в себя положительность и ограничение по частоте Найквиста: (а-б) визуализация реконструкции для кривой No. 2 на Рис. 3.56, Ді = 15.3%, линейная интерполяция; (в-г) визуализация реконструкции для кривой No. 6 на Рис. 3.56, Ді = 10.3%, интерполяция эрмитовыми сплайнами.значный вывод о том, что от точности интерполяции напрямую зависит итоговая ошибка реконструкции Д]. Варьирование ширин полосы, периода и множителя, не сказывается на ошибке реконструкции на последней итерации, в отличие от алгоритма Г-П в параллельной геометрии. И только на кривых No. 1, 2, 4 и 5 можно увидеть незначительное кратковременное улучшение Ai, что обусловлено тем, что на данной итерации ширина полосы становится меньше одного пикселя. Кривые No. 4-6, которые отвечают за реконструкцию с использованием интерполяции эрмитовыми и кубическими сплайнами, демонстрируют 5-ти процентное улучшение в ошибке реконструкции на 20-ой итерации, но итерационный процесс в этом случае также не сходится. Такое поведение итерационного процесса (медленная сходимость) и высокие итоговые погрешности вызваны большим количеством интерполяций на каждой итерации (четыре интерполяции в каждом ракурсе, здесь для К = 13 имеем 48 интерполяций (тринадцатый ракурс совпадает с первым) на каждой итерации). Такое количество интерполяций на каждой итерации вносит искажения в итоговую реконструкцию и ведет к расходимости итерационного процесса.
На Рис. 3.57 приведены примеры реконструкции модели No. 224 на 20-ой итерации. На двумерных реконструкциях между (а-б) и (в-г) нет существенного различия по внешнему виду восстановленного объекта и количеству артефактов. Одномерная интерполяция эрмитовыми сплайнами в первой модификации (см. Рис. 3.57(в-г)) дала улучшение результата на 5% по сравнению с линейной, однако полученная погрешность Ai=10.3%, все равно является высокой для данной элементарной модели (гауссианы). Этот вывод следует из экспериментов по параллельному Г-П, где ошибка реконструкции составляет Ai = 5% в отсутствии шума в проекциях.
