Содержание к диссертации
Введение
Глава 1. Реконструкция томографических изображений 15
1.1. Присутствия шума в результатах измерений сигналов 15
1.2. Алгоритмы реконструкции распределения плотности 19
1.3. Структура данных о распределении протонной плотности 24
1.4. Качество томографических изображений 28
1.5. Зависимость соотношения сигнал/шум от выбранных параметров сканирования 32
1.6. Классификация искажений томографических изображений 36
Глава 2. Математические методы обработки изображений 41
2.1. Преобразование Фурье 41
2.2. Преобразование Габора и оконное преобразование Фурье 43
2.3. Основные положения теории вейвлет-анализа 46
2.4. Вейвлет-преобразования двухмерных сигналов 50
2.5. Сравнение вейвлет-преобразования и преобразования Фурье при обработке медицинских изображений 55
Глава 3. Коррекция изображений 58
3.1. Методы обработки изображений, улучшающие их визуальное восприятие 58
3.2. Фильтрация изображений 64
3.2.1. Ранговая фильтрация 64
3.2.2. Медианная фильтрация 65
3.2.3. Адаптивная фильтрация 66
3.3. Вейвлет-анализ изображений 67
3.3.1. ВейвлетХаара 68
3.3.2. Вейвлет Добеши 69
% 3.3.3. Симлеты 70
3.3.4. Койфлеты 71
3.3.5. Другие виды вейвлетов 72
3.4. Исследование влияния различных вейвлет-фильтров на точность передачи геометрических размеров 74
3.5. Предлагаемый способ вейвлет-коррекции изображений с НЧ фильтрацией пространства 77
3.6. Постановка задачи коррекции профиля чувствительности поверхностной катушки 79
Глава 4. Решение задач шумоподавления 82
4.1. Обоснование выбора исследуемых объектов 82
4.2. Описание экспериментальной части 84
4.3. Предлагаемый алгоритм шумоподавления на основе вейвлетов 87
4.4. Анализ результатов коррекции 89
4.5. Алгоритм коррекции профиля чувствительности катушки 94
Глава 5. Сравнительный анализ оценки качества коррекции изображений 98
5.1. Оценка визуального качества 98
5.2. Количественные критерии качества изображений 100
5.3. Результаты расчетов количественных оценок 102
5.4. Практическое использование полученных результатов 106
5.5. Разработка рекомендаций по эксплуатации томографов 112
Выводы по главе 5 113
Заключение 114
Литература 116
Приложение
- Алгоритмы реконструкции распределения плотности
- Основные положения теории вейвлет-анализа
- Фильтрация изображений
- Предлагаемый алгоритм шумоподавления на основе вейвлетов
Введение к работе
Актуальность темы. В последние годы стремительными темпами развиваются технические устройства, позволяющие проводить исследования внутренней структуры объектов. Зачастую одни и те же изменения можно проводить с помощью устройств, основанные на различных принципах действия, при этом достоверность полученных данных будет сопоставима. В подобных условиях на первое место выходит информационная составляющая исследований.
На данном этапе наиболее информативным методом исследования внутренней структуры объектов является томография, на полноту и достоверность результатов которой влияет целый ряд факторов, зависящих, в первую очередь, от принципов реализации метода [74, 83].
Наиболее широко распространена компьютерная томография (КТ), позволяющая по ряду снимков путём математической обработки реконструировать плотность исследуемого вещества в ряде сечений. Популярность КТ связана с относительно низкими эксплуатационными затратами и короткой продолжительностью исследования. Это наиболее универсальный метод, позволяющий исследовать внутреннюю структуру как биологических, так и промышленных объектов, содержащих элементы из пластика, керамики, стекла и пр. Современные компьютерные томографы позволяют проводить измерения плотности вещества в диапазоне от -1024 до +3071 HU (плотность большинства биологических тканей от -150 до +100 HU). Возможности современных КТ-сканеров определяются не столько аппаратной частью, сколько математическим и программным обеспечением, обеспечивающим реконструкцию изображения.
Качество изображения зависит от многих факторов: системы коллимации (параллельности лучей и толщины среза), алгоритма реконструкции, физических процессов, участвующих в сборе данных и т.д. При этом интенсивность пиксела реконструированного изображения не всегда соответствует истинным значениям коэффициентов ослабления рентгеновских
лучей исследуемым объектом. В этом случае на изображениях появляются различные артефакты и искажения изображений. Другим фактором, влияющим на достоверность получаемых данных, является температура окружающей среды, перепады которой приводят к повышению уровня шума на изображения и снижению точности измерений.
Другим существенным ограничением информативности КТ-изображений является получение данных только в поперечной плоскости, хотя современное программное обеспечение позволяет реконструировать изображения в любой другой плоскости, но при этом снижается разрешение изображений и точность проводимых измерений.
Другим широко распространенным методом в силу своей высокой информативности является магнитно-резонансная томография (МРТ), используемая преимущественно для исследований биологических объектов и медицинской диагностики. По сравнению с КТ, МРТ имеет более высокую разрешающую способность и контрастность изображений, даёт более чёткое представление об объёме и неравномерности распространения тканей [19]. Её несомненным преимуществом по сравнению с другими видами томографии является возможность непосредственного получения изображений в любой произвольной плоскости. Наиболее эффективны в эксплуатации высокопольные сверхпроводящие МР-томографы [75], позволяющие получать изображения с высоким пространственным разрешением и проводить функциональные исследования.
Соотношение сигнал/шум магнитно-резонансного (MP) изображения зависит от целого ряда параметров сканирования (последовательности импульсов, матрицы изображения, типа РЧ катушки) и свойств вещества (протонной плотности, процессов релаксации, диффузии). Различные последовательности импульсов позволяют регулировать вклад того или иного параметра в интенсивность регистрируемого сигнала для получения на изображении оптимального контраста между нормальными и измененными тканями [80]. К недостатками МРТ можно отнести высокую продолжи-
тельность исследования. Препятствием для проведения МРТ является наличие в отображаемом объекте металла, который вносит существенные искажения в однородность магнитного поля и вызывает появление артефактов изображений. Точность измерений и качество изображений в МРТ зависит от целого ряда факторов, связанных с настройками самого сканера, алгоритмами обработки сигнала или свойствами отображаемого объекта.
Исследования промышленных объектов с помощью МРТ затруднено, т.к. любой отображаемый объект должен содержать протоны водорода и не должен содержать металлических частиц. С другой стороны, МРТ является незаменимым методом неразрушающего контроля качества продуктов питания и проведения спектрального анализа различных веществ и тканей.
Довольно часто при анализе томографических изображений приходится сталкиваться со снимками, характеризующимися повышенным уровнем шума или содержащими помехи или артефакты, которые могут привести к постановке неверного диагноза. Цифровая форма представления таких данных позволяет проводить их дополнительную компьютерную обработку, позволяющую повысить соотношение сигнал/шум и качество измерений.
Шумоподавление является важной задачей томографии, поскольку шум определяет погрешность реконструкции распределений плотности вещества и точность определения геометрических размеров.
Шумом называют беспорядочные случайные колебания различной физической природы, отличающиеся сложностью временной и спектральной структуры [94]. Уменьшить случайный шум на изображениях можно увеличивая мощность рентгеновского излучения в КТ, что приводит к росту эксплуатационных затрат, или увеличивая число измерений (что увеличивает продолжительность исследования) в МРТ. В то же время, сокращение времени сканирования позволяет устранить ошибки измерений, вызванные изменением положения исследуемого объекта в пространстве, повысить качество исследований пациентов с нарушениями сознания и увеличить пропускную способность аппарата.
В настоящее время практически ценной и актуальной является задача разработки методов и средств устранения шума и артефактов изображений в томографии, новые пути решения которой стали возможны благодаря цифровой форме представления результатов измерений и использованию протокола DICOM 3.0, выполняющему функцию стандартизации медицинской графической информации [13, 65]. DICOM-файл хранит не только сведения об условиях проведения исследования, положении пациента в момент сбора данных и другие сведения, требуемые при анализе изображений объекта, полученных в разное время и в разных условиях, но и матрицу истинных значений интенсивностей регистрируемых сигналов. Последнее свойство DICOM является его отличительной особенностью, позволяющей проводить эффективную компьютерную коррекцию изображений, поскольку прочие форматы представления данных (например, jpeg) поддерживают лишь 256-уровневую шкалу интенсивностей, что приводит к потере части информации. Данная диссертационная работа посвящена разработке методов и алгоритмов шумоподавления на основе вейвлет-анализа изображений, представленных в формате DICOM.
Цель диссертационной работы. Целью данной работы является анализ и совершенствование существующих, а также разработка новых методов и средств шумоподавления.
Достижение данной цели позволяет решить важную научно-техническую проблему повышения качества изображений внутренней структуры объектов при сокращении времени исследования и лучевой нагрузки в различных видах томографии.
Задачи исследования. Для достижения поставленной цели были решены следующие задачи:
-анализ и сравнение существующих методов повышения качества изображений;
- исследование влияния вейвлет-фильтров на точность передачи геометрических размеров и на соотношение сигнал/шум;
-разработка методов и средств обработки томографических изображений, позволяющих повысить соотношение сигнал/шум и устранить наиболее характерные артефакты при сокращении времени сканирования;
- разработка и тестирование алгоритмов и программ шумоподавления, работающих с данными в исходном формате томографа, позволяющих снизить временные затраты на обработку данных и объем требуемой машинной памяти;
-разработка методики коррекции профиля чувствительности регистрирующей сигнал радиочастотной катушки;
-нахождение количественных оценок результатов шумоподавления, их статистический анализ и сравнение полученных методов шумоподавления с существующими.
Основные положения, защищаемые в диссертации:
использование ранговой, медианной и адаптивной фильтрации хотя и снижает уровень шума, приводит к размыванию изображений и потере мелких структур объектов;
разработанный метод вейвлет-анализа позволяет повысить качество изображений промышленных и биологических объектов при сокращении времени исследования, хотя, как известно, вейвлет-фильтры не снижают точность передачи геометрических размеров, но при этом позволяют увеличить соотношение сигнал/шум;
созданные программные средства, использующие данные в исходном формате томографа, позволяют повысить эффективность шумоподавления;
предложенная методика коррекции профиля чувствительности радиочастотной катушки позволяет повысить соотношение сигнал/шум и устранить ряд артефактов изображений;
проведенный статистический анализ качества коррекции изображений предложенным методом вейвлет-анализа позволяет сделать вывод об эффективности разработанных в диссертации методов и средств шумоподавления.
Методы исследования. Основные результаты работы получены с применением вейвлет-преобразовований данных и преобразования Фурье. В работе использованы методы ранговой и медианой фильтрации, фильтрации Винера, методы коррекции исходных данных к-пространства, вейв-лет-анализ, методы статистической обработки экспериментальных данных.
Научная новизна работы состоит в следующем:
исследовано влияние применения различных вейвлет-фильтров на точность передачи геометрических размеров в томографии и соотношения сигнал/шум изображений;
на основе проведенного системного анализа признаков артефактов томографических изображений предложены алгоритмы улучшения их качества;
разработан алгоритм коррекции искажений, вызванных профилем чувствительности катушки;
разработано программное обеспечение, реализующее подавление шума изображений и работающее с данными в исходном формате прибора;
исследована зависимость качества откорректированных изображений от используемого алгоритма фильтрации и его параметров.
Совокупность представленных в работе результатов может рассматриваться как фундаментальные основы инженерных наук, касающиеся теории и эффективности функционирования томографических измерительных систем. Внедрение технических решений, принятых на основании исследований, проведенных в диссертационной работе, вносит значительный вклад в совершенствование томографических исследований.
Достоверность научных результатов, полученных в работе, обеспечивается строгостью постановки задач, применяемых математических методов, статистической обработкой полученных результатов и их сравнением, где это возможно, с известными данными. Обработка экспериментальных данных проводилась на базе кафедры измерительных технологий и компьютерной томографии СПбГУ ИТМО и кафедры рентгенологии Санкт-
Петербургской медицинской академии последипломного образования. Томографические изображения были получены в Санкт-Петербургской медицинской академии последипломного образования на МР-томографе "Signa Infinity" ("General Electric") с магнитным полем В = 1,5 Тл и на компьютерном томографе LightSpeed Plus ("General Electric"), в СПб ГУЗ «Городская Покровская больница» на МР-томографе "Signa Infinity" ("General Electric") с магнитным полем В = 1,0 Тл и компьютерном томографе "Somatom AR.SP" ("Siemens"), в Санкт-Петербургском научно-исследовательском психоневрологическом институте им. В.Н. Бехтерева на МР-томографе "Universal-Max" с магнитным полем В = 0,15 Тл.
Практическая ценность результатов работы заключается в создании средства шумоподавления, реализуемого с помощью программного продукта, эффективное применение которого возможно для широкого класса томографических исследований, требующих высокого качества изображений и точности измерений. Выполненные исследования и статистический анализ полученных результатов позволяют выработать рекомендации для пользователей медицинских томографов в областях, требующих быстрого получения изображений с высоким соотношением сигнал/шум. Метод улучшения качества и сжатия томографических изображений с помощью вейвлет-фильтров может применяться в клинической практике в диагностических центрах, в том числе при проведении удаленных консультаций (телемедицина) [58].
Реализация работы. Разработанные средства шумоподавления были использованы в ФГУ «Российский научно-исследовательский нейрохирургический институт им. проф. А.Л. Поленова» в рамках федеральной НИР «Некоторые фундаментальные вопросы современной нейрохирургии». Результаты работы были использованы при обработке результатов томографических исследований мелких биологических структур головного мозга, полученных с высоким пространственным разрешение.
Предложенные методики повышения качества томографических изо-
бражений позволило улучшить их качество и повысить информационную составляющую результатов исследований, проводимых на магнитно-резонансном томографе Universal-Max, установленном в Научно-исследовательском психоневрологическом институте им. В.Н. Бехтерева.
Разработанные алгоритмы шумоподавления в томографии, позволяющие повысить соотношение сигнал/шум и устранить характерные для быстрого сбора данных в MP-томографии артефактов, прошли апробацию и были внедрены на магнитно-резонансном и компьютерном томографах СПб ГУЗ «Городская Покровская больница».
Результаты диссертационной работы были внедрены в ГОУ ВПО «Санкт-Петербургская медицинская академия последипломного образования» при коррекции артефактов изображений, вызванных взаимодействием биологического объекта с поверхностью радиочастотной катушки. Предложена методика снижения артефактов на основе коррекции профиля чувствительности катушки и вейвлет-анализа сигнала.
Опубликованные результаты работы внедрены в учебный процесс на факультете точной механики и технологий СПбГУ ИТМО. Материалы были использованы при чтении лекций и проведении лабораторных работ по дисциплинам ЕН.Ф.06 «Физические основы получения информации», СД.Ф.02 «Теория измерений», СДМ.В.02 «Современные виды томографии», СД.Р.05 «Конструирование медицинских томографов», СД.Р.08 «Лучевая диагностика в клинической медицине», а также в учебном пособии [81], получившем гриф «Рекомендовано УМО по образованию в области приборостроения и оптотехники в качестве учебного издания для студентов высших учебных заведений, обучающихся по направлению подготовки 200100 "Приборостроение"».
Практическое использование результатов диссертационной работы подтверждено соответствующими документами.
Работа получила развитие и поддержку Российского фонда фундаментальных исследований по проекту 05-08-65468а как работа в области фун-
даментальных основ инженерных наук.
Апробация работы. Результаты работы докладывались, обсуждались и получили положительную оценку на кафедре измерительных технологий и компьютерной томографии и кафедре мехатроники; на I-III конференции молодых ученых (2003 г., 2004 г., 2005 г., СПб); на IX международной конференции «Региональная информатика» (2004 г., СПб); 4th European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS (Finland, 2004), конференции «Оптика и образование» (2004 г., 2006 г., СПб); XXXIV научной и учебно-методической конференции СПбГУ ИТМО (2005 г., СПб); Политехническом симпозиуме «Молодые ученые - промышленности Северо-Западного региона» (2005 г., СПб); конференции «Информатика и управление в медицинских системах» (2006 г., СПб).
Публикации. По материалам диссертации опубликовано 13 работ, в том числе одно учебное пособие с грифом УМО.
Структура и объем диссертации. Диссертация состоит из введения, 5 глав, заключения, 3 приложений, библиографического списка из 95 наименований. Объем диссертации 124 страницы, 42 рисунка и 10 таблиц.
В первой главе рассмотрены принципы получения информации о распределении плотности в сечениях объекта с помощью рентгеновского излучения и явления магнитного резонанса ядер. Рассматриваются виды задач реконструкции плотности, структура регистрируемых данных, причины появления шума изображений. Проведена классификация возникающих артефактов изображений, проанализированы причины их возникновения и проявление на изображениях, приведены существующие пути устранения.
Во второй главе рассмотрены математические методы и алгоритмы, используемые при реконструкции изображений в томографии, приведены особенности их реализации, выявлены их достоинства и недостатки применительно к конкретным измерительным задачам. Рассмотрены линейные преобразования в анализе сигналов, проведен анализ некоторых хорошо известных математических методов обработки сигналов. Приведены ос-
новные положения теории вейвлет-анализа, проанализированы возможные способы применения вейвлет-анализа к обработке изображений.
В третьей главе рассмотрены различные методы шумоподавления и повышения качества изображений. Рассмотрены вопросы поиска вейвлет-базисов, наилучшим образом подходящих для аппроксимации изображений, рассмотрен ряд детализирующих вейвлет-функций. Проведен анализ влияния применения различных методов шумоподавления на точность передачи геометрических размеров. Показано, что исследуемые вейвлет-фильтры не оказывают существенного влияния на точность измерений.
В четвертой главе описываются разработанные методы и средства шумоподавления применительно к конкретным измерительным задачам. Приводится описание ряда модельных экспериментов, обоснован выбор исследуемых объектов, приведены результаты обработки их изображений. Показано, что теория вейвлет-анализа является адекватным математическим аппаратом для выявления и устранения шума изображений. Предложен способ обработки данных на основе сочетания коррекции профиля чувствительности поверхностной катушки и применения вейвлет-фильтрации сигнала.
В пятой главе проанализированы возможности применения различных вейвлет-фильтров относительно существующих задач шумоподавления изображений промышленных и биологических объектов. Рассмотрены критерии оценки качества полутоновых монохромных изображений. Проведен статистический анализ результатов коррекции изображений и количественная оценка их качества.
Таким образом, в данной работе решена научная проблема, заключающаяся в разработке новых методов и средств шумоподавления, повышающих качество изображений в томографии. Как следует из результатов работы, такие средства могут быть разработаны на основе вейвлет-анализа изображений. В силу универсальности, предложенные в работе методы шумоподавления охватывают широкий спектр современных прикладных измерительных задач и эксплуатации томографических систем.
Алгоритмы реконструкции распределения плотности
Компьютерная томография является одним из методов исследования внутренней структуры материальных объектов. Она получила широкое применение в различных областях науки: в медицине для исследования внутренних органов, в промышленности для контроля качества изделий, в геофизике для исследования мантии Земли, в физике для диагностике плазмы и т.д. Промышленные рентгеновские компьютерные томографы позволяют обнаружить дефекты размером порядка 1 мкм, благодаря чему удается повысить качество и надежность выпускаемой продукции. Срок окупаемости затрат на промышленные рентгеновские томографы во многих случаях в 5-Ю раз меньше срока окупаемости технологического оборудования. К достоинствам вычислительной КТ относятся высокая скорость и невысокая стоимость исследования.
В рентгеновской КТ внутренняя структура объекта характеризуется распределением плотности. В данном случае под плотностью понимается плотность электронов, т.е. количество электронов в единице объема. Тонкий пучок рентгеновских лучей сканирует сечение объекта, испытывает частичное поглощение веществом, и изменение интенсивности рентгеновского излучения, являющееся интегральной характеристикой плотности исследуемого сечения, фиксируется детектором излучения. Затем с помощью компьютера по полученным данным производится реконструкция исследуемого распределения плотности в сечении объекта. Качество и скорость реконструкции распределений плотности зависят от конструктивных особенностей томографа и от используемого математического аппарата.
Решение математических задач компьютерной томографии сводится к решению операторных уравнений 1-го рода. При нахождении их приближенных решений необходимо использовать методы регуляризации, позволяющие учитывать дополнительную информацию о решаемой задаче [90, 92]. Одна из главных проблем, возникающих при решении математических задач томографии, - выбор оптимального алгоритма, критерием отбора ко # торого может служить, например, качество изображения [54, 59].
Рассмотрим основные математические соотношения, на которых базируются современные методы вычислительной томографии.
Пусть на плоскости (х, у) в прямоугольной системе координат заданадвухмерная функция Дх,у), интегрируемая по всем возможным прямым,лежащим в данной плоскости (рис. 1). Всякая прямая может быть описана уравнениемгде s- расстояние от начала координат до рассматриваемой прямой; д угол, образованный с осью х перпендикуляром, опущенным на прямую изt начала координат.s и ф .Поэтому результат R интегрирования функции Дх,у) по некоторой прямой будет зависеть от этих же параметров (R = R(s, ф)):где 8 - дельта-функция Дирака.
Подобное интегрирование можно рассматривать как некоторое преобразование, которое для функции Дх,у) на плоскости {х,у} ставит в соответствие R(s, ф) на множестве всех прямых. Это преобразование называется преобразованием Радона, а функцию R(s, p) называют образом функции Дх,у) в пространстве Радона. В томографии ставится математическая задача поиска неизвестной функции Дх,у), если известна функцияR(s, ф), являющаяся образом функции f(x,y) в пространстве Радона.
Различные алгоритмы восстановления различаются между собой способом учёта технических особенностей; степенью детальности учёта флук-туационных явлений; объёмом априорных сведений. Проекция изображения формируется объединением набора измерений, проведенных вдоль параллельных линий. Существует соотношение, определяющее связь между преобразованием Фурье этих функций - теорема о центральном сечении.
Пусть R{a ,(p) - одномерное преобразование Фурье (спектр Фурье) функции R{s,q ) по переменной s, a F(u,v) - двумерное преобразование Фурье (пространственный спектр) функции f(x,y) по переменным х и у:
Теорема о центральном сечении говорит, что если функция f(x,y) и ее радоновский образ R(s, (р) имеют преобразования Фурье, то одномерное преобразование Фурье радоновского образа R(s, ср) по переменной s равно функции, описывающей центральное сечение двумерного преобразования Фурье, соответствующее тому значению (р, при котором вычисляетсяпреобразование Фурье функции R{s,qj). Математическая формулировкатеоремы о центральном сечении имеет вид:Задача восстановления изображения базируется на теореме о центральном сечении. Функцию f(x,y) можно найти по двумерному преобразованию Фурье F(u, v):Перейдём в плоскости (M,V) К полярным координатам а,(р:и = со cos ср, v = со sin ср. Тогда уравнение (1) будет иметь вид:
Равенство (2) является искомой формулой обращения, позволяющей найти функцию f(x,y).Алгоритм фоновой проекции относительно прост для параллельной схемы сканирования, но реконструкция занимает много времени. Веерное сканирование намного быстрее, но алгоритм для него более сложен [95].
Регистрируемые детектором данные это результат взаимодействия рентгеновского излучения и исследуемого вещества [49, 69]. Коэффициент поглощения фотонов при похождении через материал зависит от свойств материала:интенсивность испускаемого и регистрируемого излучения, соответственно, ц. - коэффициент линейного ослабления материала.
Интенсивность излучения, прошедшего через исследуемый объект: I = I0 рассеивания излучения вокселами. Рис. 2. Схема получения данных при компьютерной томографии Найти коэффициенты поглощения для каждого воксела можно с помощью метода обратного проецирования, предполагающего получение информации о характере поглощения излучения во многих ракурсах.Для слоя, состоящего из четырех вокселов (рис. 2), коэффициенты ослабления излучения Рис. 3. Проекционные данные (а) и реконструированное КТ-изображение (б)
В современных томографах матрица изображения чаще всего имеет размер 512x512 или 1024x1024 пикселов и восстанавливать коэффициенты рассеивания приходится, решая соответствующее количество уравнений.
Выходные данные КТ сканера измеряются в единицах Хаунсфилда (HU), и для современных томографов лежат в диапазоне от -1024 до +3071 HU. Соотношение между коэффициентом линейного ослабления иссле-дуемого вещества (д,) и соответствующей единицей Хаунсфилда (//): Н = М,-М А. ,1Q00_
Значения плотности некоторых веществ известны: плотность воздуха составляет -1000 HU, воды 0 HU, пластика +100 HU. Плотность металла выходит за предел нормального диапазона значений, который томограф может отобразить, поэтому компьютер присваивает им наибольшее возможное значение. Значения плотности большинства биологических тканей лежат в достаточно узком диапазоне (табл. 1). 4
Основные положения теории вейвлет-анализа
В отличие от Фурье-анализа, при котором размер частотно-временного окна остается постоянным, вейвлет-анализ позволяет анализировать сигнал в изменяющемся частотно-временном окне. Это дает возможность проводить анализ непериоических и значительно изменяющих свои частотные характеристики во времени сигналов.
Термин «вейвлет» (дословный перевод - маленькая волна, всплекск) был введен Гросманом и Морле в середине 80-х применительно к анализу сейсмических и акустических сигналов. В настоящее время семейство анализаторов, названных веивлетами, применяется в задачах распознавания образов, при обработке и синтезе различных сигналов [17] (например, речевых), при анализе изображений различной природы [6] (например, фотографии, томограммы, спутниковые изображения и др.). Также вейвлеты используют для сжатия больших объемов информации.
Вейвлет преобразование одномерного сигнала состоит в его разложении по базису, сконструированному из обладающей определенными свойствами функции (вейвлета) путем применения к ней преобразования переноса и изменения масштаба. Каждая из функций этого базиса характеризует как определенную пространственную (временную) частоту, так и ее локализацию в пространстве (времени).
Таким образом, в отличие от преобразования Фурье, вейвлет преобразование обеспечивает двумерную развертку исследуемого одномерного сигнала, причем частота и координата рассматриваются как независимые переменные. В результате появляется возможность анализировать свойства сигнала одновременно в физическом (время, координата) и в частотном пространствах [60, 88].
В данной работе излагается суть применения метода вейвлет-анализа к фильтрации изображений (шумоподавлению) на примере томограмм.
В теории вейвлетов сигнал s(t) представляют в виде взвешеннойсуммы простых составляющих - базисных функций ц/к it) (вейвлетов), ум ф ноженных на коэффициенты Ск:
Функцию у/ є L (4R) называют базисным вейвлетом, если она удовлетворяет следующему условию:где у/ - преобразование Фурье функции у/.
Т.к. базисные функции у/к (t) предполагаются заданными как функция вполне определенного вида, то только коэффициенты Ск содержат ин f формацию о конкретном сигнале. Таким образом, можно говорить о воз можности представления сигналов на основе рядов с различными базисными функциями. Например, ряд Фурье использует в качестве базисных функций синус и косинус.
Вейвлеты характеризуются своим временным и частотным образами. Временной образ определяется некоторой функцией времени y/(t), а частотный образ определяется ее Фурье-образом у/(со) = F(co), который задает огибающую спектра вейвлета. Фурье-образ определяется выражением:
Одна из основополагающих идей вейвлет-представления сигналов заключается в разбивке приближения к сигналу на две составляющих - грубую (аппроксимирующую) и уточненную (детализирующую) - с последующим их уточнением итерационным методом. Каждый шаг такого уточнения соответствует определенному уровню декомпозиции и реставрации сигнала. Это возможно как во временной, так и в частотной областях.
В основе непрерывного вейвлет-преобразования лежит использхова-ние двух непрерывных и интегрируемых по всей оси t (или х) функций: - вейвлет-функция y/(t) с нулевым значением интеграла ( \y/(t)dt = 0),определяющая детали сигнала и порождающая детализирующие коэффициенты; - масштабирующая функция cp(t) с единичным значением интеграла ( fa(t)dt = 1), определяющая грубое приближение (аппроксимацию) сигна -00ла и порождающая коэффициенты аппроксимации.Аппроксимирующие функции q (t) присущи не всем вейвлетам, а
только тем, которые относятся к ортогональным.
Детализирующая функция y/(t) создается на основе той или иной базисной функции if/(f(t), которая, как и y/{t), определяет тип вейвлета. Базисная функция должна удовлетворять тем же требованиям, что и функция \f/(t) и должна обеспечивать выполнение двух основных операций:- смещение по оси времени t-i//0(t-b) при b є $R;- масштабирование -а 12у/0\— 1 при а О и а є 9?+ - {0}.
Параметр а задает ширину этого пакета, а Ъ - его положение. Слеа Таким образом, для заданных значений а и Ъ функция y/(t) и естьвейвлет. Вейвлеты, обозначаемые как y/(t), иногда называют «материнскими вейвлетами», поскольку они порождают целый ряд вейвлетов определенного типа.
Дискретное вейвлет-преобразование сигнала определяется соотношением:Выражение (7) является основой для вейвлет-анализ сигналов. На нем базируются методы фильтрации сигналов, кодирования сигналов, кратно-масштабного анализа [50, 17].
Идея кратномасштабного анализа была разработана Мейером. В процессе проведения кратномасштабного анализа, исследуемая функция подвергается разложению на составляющие с различным разрешением [61].
Пусть функция состоит из медленно меняющихся и быстро меняющихся сегментов. Если мы хотим представить эту функцию на каком-тоодном уровне аппроксимации, то мы должны использовать шаг дискретизации h, определяемый из быстро меняющегося сегмента функции. Этоприведет к большему количеству лишней информации, которую нам придется хранить и обрабатывать. Представив функцию с использованием нескольких разных шагов дискретизации, мы можем значительно уменьшитьколичество информации, необходимое для точного восстановления функции. Например, можно представить функцию на четырех уровнях дискретизации, каждый раз удваивая шаг. С каждым удвоением шага, мы получаем аппроксимацию функции и ее детализацию (рис. 8).
Фильтрация изображений
При поиске оптимальных алгоритмов обработки сигнала неизбежно приходится опираться на некоторые статистические модели сигналов и шумов. Чаще всего при формировании этих моделей используются концепции линейности, стационарности и нормальности (гауссовости). Наиболее часто для обработки полутоновых изображений используют алгоритмы ранговой, медианной и адаптивной фильтрации.
Преимуществом ранговых алгоритмов в отличие от методов линейной фильтрации, является отсутствие пространственной инерционности, которая заключается в том, что при использовании линейных фильтров влияние отдельных деталей изображения проявляется на результирующем изображении на расстоянии порядка размеров апертуры фильтра. В этом случае пиксели исходного изображения, соответствующие ненулевым элементам маски фильтра, переупорядочиваются в вариационный ряд в порядке возрастания. Пикселю исходного изображения, соответствующему центральному элементу маски, присваивается значение с заданным порядковым номером (рангом). Изменяя положение центрального пикселя маски, свойства ранговой фильтрации можно менять в широких пределах [62, 66].
Ранговая фильтрация может подчеркивать дефекты (или нужные осо і,. бенности) и способна эффективно отфильтровывать мелкие детали, на пример, шумы. В последнем случае, метод ранговой фильтрации практически равноценен медианной фильтрации. Алгоритм ранговой фильтрации позволяет создавать эффект расфокусировки (эрозии) и фокусировки (уточнения, наращивания фона и т.д.) обрабатываемого изображения. В первом случае центральный пиксель маски надо выбирать как пиксель минимальный, а во втором - как пиксель максимальный (по порядку в маске).Рис. 12. Ранговая фильтрация с окном 3 3: (а) исходное изображение, (б) ранг 1; (в) ранг 9
Кроме применений для сглаживания, усиления детальности, выделения деталей изображений и границ деталей, ранговые алгоритмы можно употреблять для решения многих других задач обработки изображений, в частности, для диагностики статистических характеристик искажений видеосигнала, стандартизации изображений, определения статистических характеристик видеосигнала и измерения текстурных признаков. 3.2.2. Медианная фильтрация
Медианная фильтрация это метод низкочастотной фильтрации изображений, позволяющий убрать резкие выбросы, и являющийся частным случаем ранговой фильтрации [41, 62]. Медианная фильтрация подавляет «импульсный» шум, для которого является более эффективным средством, чем фильтрация на основе операций свертки.
Медианный фильтр реализуется как процедура локальной обработки скользящим окном mxn (маской), которое включает нечетное число отсчетов. Для каждого положения окна все значения пикселов сортируются в порядке возрастания, и за новое значение фильтруемой точки принимается % значение центральное элемента маски (медиана).
В результате применения медианного фильтра наклонные участки и резкие перепады значений яркости на изображениях не изменяются. Это является достоинством метода, поскольку на изображениях контуры несут основную информацию. В то же время, любой пиксел, значение которого намного отличается от его соседей, исключается, т.к. если в точке был выброс, то она оказывается на краю окна и не попадает в отфильтрованное изображение. Медианная фильтрация намного более чувствительна к резко отличающимся значениям пикселов, чем усреднение, поэтому лучше удаляет подобные выбросы не снижая резкости (точности) изображения.
Медианный фильтр является нелинейным. Недостатком метода является относительно слабая эффективность при фильтрации флуктуационно-го шума. Кроме того, при увеличении размера маски происходит размытие контуров изображения, и, как следствие, снижение четкости изображения. Фильтр эффективен при подавлении импульсных помех при подборе размеров маски, когда известны минимальные размеры объектов и максимальные размеры искаженных помехой локальных областей. Поэтому на практике для обработки томографических изображений данный вид Метод адаптивной фильтрации Винера обычно применяется для подавления на изображениях гауссовского белого аддитивного шума [63, 89]. При осуществлении данной фильтрации учитываются статистические особенности изображения в пределах выбранного окна, в частности среднее значение яркости пикселов и ее среднеквадратическое отклонение.
Если в выбранном окне дисперсия яркостей большая, то фильтр Винера дает небольшое сглаживание; если дисперсия мала, то фильтр дает большее сглаживание изображения.
Данный метод часто -дает лучший результат, чем линейная фильтрация. Адаптивный фильтр более избирательный, чем сопоставимый линей щ Применение новых эффективных способов обработки изображений с помощью вейвлет-преобразований позволяет легко выявлять все локальные особенности функций, сигналов и изображений с привязкой их ко времени [1] или координатам пространства. Вейвлет-преобразование сигнала состоит в его разложении по базису, сконструированному из обладающей определенными свойствами функции (вейвлета) посредством масштабных изменений и переносов [52]. Каждая из функций этого базиса характеризует как определенную пространственную (временную) частоту, так и ее локализацию в физическом пространстве (времени). В результате появляется возможность анализировать свойства сигнала одновременно в физическом (время, координата) и частотном пространствах.
С каждым вейвлетом можно связать несколько фильтров. При разложении сигнала на низкочастотную и высокочастотную составляющую используются два фильтра разложения, которые определенным образом свя заны с масштабирующим фильтром W. При реконструкции сигнала используются еще два фильтра, которые также вычисляются по масштабирующему фильтру W. К масштабирующим фильтрам относятся вейвлет-фильтры Добеши, симлеты, койфлеты, биортоганальные и обратные био-ртогональные сплайновые фильтры [91]. Выбор конкретного вейвлета во многом определяется исследуемым сигналом и задачами анализа. Учитывая характерные особенности различных вейвлетов во временном и в частотном пространстве, можно выявлять в анализируемых сигналах те или иные свойства, которые незаметны на их графиках, особенно в присутствии сильных шумов [32, 34, 40]. Рассмотрим особенности некоторых вейвлетов.
При анализе любого сигнала надо прежде всего выбрать соответствующий базис, т.е. систему функций, которые будут играть роль «функциональных координат». В большинстве случаев приходится иметь дело с сигналами, которые представлены квадратично-итегрируемыми функциями, определенными на вещественной оси (или квадратично-суммируемыми последовательностями комплексных чисел [3]). Выбор анализирующего вейвлета не определен заранее, его следует выбирать в соответствии с решаемой проблемой.
Вейвлет Хаара - простейший ортогональный дискретный вейвлет, порождающий ортонормированный базис [56] и широко используется в функциональном анализе. Для декомпозиции сигнала используются две функции. Одна это аппроксимирующая функция, которая у вейвлета Хаара равна +1 на всем компактном носителе, и детализирующая функция, определяемая соотношением
Предлагаемый алгоритм шумоподавления на основе вейвлетов
На основе анализа существующих методов шумоподавления и результатов обработки с их помощью изображений, предлагаемый алгоритм шумоподавления в общем случае выглядит следующим образом:- выбор исходного изображения;- считывание из DICOM-файла метаданных (модальность оборудования, параметры исследования, тип исследуемой области);- считывание из DICOM-файла матрицы изображения (например, абсолютные значения интенсивности MP-сигнала, значения плотности в единицах Хаунсфилда) и преобразование данных;- выбор базисной вейвлет-функции или семейства вейвлет-функций;- выбор уровня N вейвлет-разложения данных;- пороговая обработка детализирующих коэффициентов;- реконструкция изображения;- оценка качества коррекции;- выбор коэффициента предварительной гамма-коррекции (при необходимости);- определения размера и типа маски повышения резкости изображения (при необходимости).
Выбор исходного изображения осуществляется оператором, исходя из визуальной оценки качества полученных данных. При этом в качестве формата представления обрабатываемых данных используется DICOM, хранящий информацию о результатах и условиях проведения измерений.
Считывание метаданных из DICOM-файла позволяет получить информацию о модальности диагностического устройства и выбрать оптимальные параметры обработки, данные об исследуемой анатомической области или фантоме также содержат информацию, используемую в коррекции. В случае компьютерной томографии необходимо учитывать радиус сканирования сканера, а для МРТ - тип импульсной последовательности.
Матрица изображения содержат значений интенсивностей пикселов, реконструированные сканером из исходных проекционных данных.
Выбор базисной вейвлет-функции или семейства базисных вейвлет-функций определяется типом исходного изображения и видом измерений.
Выбор уровня N вейвлет-разложения данных определяет порядок разложения сигнала и количество рассчитываемых детализирующих и аппроксимирующих коэффициентов.
В случае выбора семейства вейвлет-функций (Добеши, симлетов, койфлетов и т.д.) необходимо выбрать порядок вейвлета, задаваемый его номером, и определяющий вид аппроксимирующей и детализирующей функций.
Поскольку шумовая компонента представляет собой сигнал меньший по модулю, чем основной, поэтому при удалении шума выбирают нулевые значения детализирующих коэффициентов, меньших некоторого порогового значения т. Выбор порога полезного сигнала определяется степенью зашумленности исходного изображения. При жесткой обработке сохраняются все коэффициенты, большие или равные по абсолютной величине порога т, а меньшие коэффициенты обращаются в нуль. При мягкой пороговой обработке наряду с обращением в нуль коэффициентов, меньших по модулю чем т, происходит уменьшение по модулю остальных коэффициентов на величину т. От выбора порогового уровня фона (оценки дисперсии шума), зависит качество шумоподавления сигнала, оцениваемое в виде соотношения сигнал/шум. Установка малых значений порога сохраняет фон в коэффициентах детализации и поэтому приводит лишь к незначительному увеличению отношения сигнал/шум. При больших значениях порога можно потерять коэффициенты, которые несут существенную информацию. Поиск оптимального значения т означает отыскание такого порога, который при наименьшем смещении восстановленного сигнала обеспечивает наибольшее значение отношения сигнал/шум.
Вейвлет-реконструкция изображения производится на основании первоначальных аппроксимирующих коэффициентов уровня N и модифициро ванных детализирующих коэффициентах уровней от 1 до N.
Качество проведенной коррекции позволяют оценить количественные критерии, рассчитываемые в подпрограмме для исходного и откорректированного изображений.
Для изменения яркости пикселов при отображении (без изменения абсолютных значений их интенсивности), используется гамма-коррекция изображения, позволяющая изменить контраст между самым ярким и самым темным участками изображения.
В случае если после подавления шума изображения требуется повы-сить его резкость, возможно использование повышающих резкость фильтров.
Коррекция томографических изображений осуществлялась с помощью программного обеспечения, написанного для работы в среде MATLAB (приложение 2) и состоящего из нескольких подпрограмм, позволяющих осуществлять обработку данных с различных модальностей и использующих входные параметры из DICOM-файла. Среда MATLAB имеет удобный интерфейс и используется для программирования и обра-ботки результатов научных исследований [18]. Каждое из полученных изображений обрабатывалось с помощью 98 фильтров (адаптивного, рангового и медианного с различным окном, вейвлет-фильтров Добеши (со 2 по 20 порядок), симлетов (с 1 по 20 порядок), койфлетов (с 1 по 50 порядок), биортогональных и обратных биортогональных вейвлетов различных порядков, вейвлетов Хаара, Мейера, мексиканская шляпа, Морле). Получаемые данные записывались в виде файла-изображения в отдельную директорию.
При коррекции изображений однородных фантомов в компьютерной и магнитно-резонансной томографии результат был приблизительно одинаков для различных типов вейвлетов. При визуальной оценке полученных изображений оставалась заметна зернистая структура, но она становилась
