Содержание к диссертации
Введение
1. Обзор основных мрт-систем и принципов генерирования изображений 10
1.1. Краткий исторический ракурс теории параллельной реконструкции МРТ-изображений 10
1.2. Обзор и классификация МРТ-систем 11
1.3. Основные технические характеристики МРТ-систем 13
1.4. Импульсные последовательности и их влияние на качество изображения 18
1.5. Принципы построения МРТ-изображений 21
1.6. Обзор тестовых объектов для проведения МРТ 26
1.7. Постановка задачи 27
1.8 Выводы по главе 28
2. Параллельная реконструкция изображений в МРТ 29
2.1. Исследование программных библиотек для реконструкции и обработки МРТ-
изображений 34
2.2. Генерирование многоканальных тестовых изображений 35
2.3. Исследование основных методов параллельной реконструкции МРТ-изображений
46
2.5. Разработка математических моделей и метода параллельной реконструкции МРТ-изображений 72
2.6. Разработка алгоритма параллельной реконструкции МРТ-изображений 79
2.7. Оценка погрешности метода параллельной реконструкции МРТ-изображений 80
2.8. Основные результаты и выводы 88
3. Разработка приемной системы сбора данных МРТ 90
3.1. Разработка приемного тракта MP-сигналов 98
3.2. Ввод данных в ЭВМ 107
3.3. Основные результаты и выводы 109
4. Реконструкция мрт-изображений на многопроцессорной системе
4.1. Исследование пакета программ Gadgetron 112
4.2. Реализация алгоритма параллельной реконструкции на основе графического процессора 115
4.3. Тестирование производительности реконструкции с использованием графического процессора 115
4.4. Основные результаты и выводы 122
Заключение 123
Библиографический список
- Обзор и классификация МРТ-систем
- Генерирование многоканальных тестовых изображений
- Ввод данных в ЭВМ
- Реализация алгоритма параллельной реконструкции на основе графического процессора
Обзор и классификация МРТ-систем
В истории МРТ идея многоканальности и параллельного сбора МР-сигналов впервые была предложена в 1988 году М. Хатчинсон (США) и У. Ра. После этого теория получила дальнейшее развитие в работах Д. Квит и С. Эйнав (Израиль) в 1991 году. Главный принцип многоканальной МРТ заключается в том, чтобы область чувствительности каждой катушки была мала и уникальна. В 1989 году Келтон (США) предложил идею сокращения расстояния между выборками К-пространства в направлении фазового кодирования с помощью использования многоканальных систем.
Впоследствии появились разнообразные методы параллельной реконструкции изображений:
Предложен в 2002 Марком Грисволдом (США) [24]. Главной проблемой для реконструкции МРТ-изображений подобных техник является работа при высоких уровнях шума и при высоких коэффициентах акселерации. Под термином высокий коэффициент акселерации далее в работе будет пониматься случай, когда число недостающих шагов фазового кодирования более 50% от значения по теореме Котельникова-Найквиста-Шенона.
Таким образом, перед исследователями стоит следующая задача: сокращение уровня геометрических искажений и повышение устойчивости методов параллельной реконструкции изображений МРТ. Для достижения этой задачи в работах [41, 47, 48, 49] предлагается метод SPIRiT. Он является разновидностью автокалибровочного метода реконструкции. Во многом он похож на GRAPPA- реконструкцию, но обладает лучшими характеристиками по уровню искажений изображений. Особенно сильно заметна разница при реконструкции изображений с высоким коэффициентом акселерации[110].
Из проведенного обзора можно заключить, что исследование методов параллельной реконструкции при высоких степенях акселерации и разработка устойчивых методов является актуальной задачей в настоящее время.
Обзор и классификация МРТ-систем Как известно, все МРТ-системы медицинского назначения можно разделить на системы закрытого и открытого типа. Принципиальное отличие состоит в конструкции магнита, который создает базовое поле ВО для работы МРТ. Магнит может быть как сверхпроводящим, так и постоянным. Следует отметить, что в настоящее время сверхпроводящие магниты (то есть МРТ закрытого типа) имеют напряженность поля от 1 до 7 Тесла и выше, в то время как МРТ открытого типа имеют напряженность, как правило, не более 0,5 Тесла.
Низкая напряженность магнитного поля ВО приводит к низкой частоте ЯМР. Известно, что увеличение резонансной частоты приводит к улучшению качества МРТ-изображений. Поэтому, как правило, в клинической практике применяются МРТ с напряженностью поля более 1 Тесла. Хотя для исследования головного мозга и позвоночника в большинстве случаев достаточно МРТ открытого типа, чего нельзя сказать о качественной визуализации брюшной полости.
Из общего числа МРТ в России порядка 30 % оборудования изготовлено фирмой Siemens [82]. Среди моделей МРТ Siemens наиболее часто используются Siemens Impact, Harmony, Symphony, Espree и т. д.
Помимо Siemens парк томографов в России представлен следующими фирмами производителями: Philips[71], GE, Philips, Toshiba. Некоторые из них представлены на фотографии 1.2.
Генерирование многоканальных тестовых изображений
Целью исследований, результаты которых приведены в таблице 2.1, является исследование коэффициента сигнал/шум при различных результатах реконструкции. Исследуя полученные изображения (оценка которых представлена в таблице 2.1), можно сделать выводы: - увеличение числа линий автокалибровки при использовании метода параллельной реконструкции изображений GRAPPA ведет к уменьшению фактора геометрических искажений, но несколько снижает коэффициент сигнал-шум; - использование четырех каналов совместно с методом GRAPPA приводит к видимым неравномерным геометрическим искажениям на изображении (в виде неравномерного шума); - в подавляющем числе экспериментов реконструкция изображений методом дробного преобразования Фурье показывает менее качественные результаты, чем параллельная реконструкция методом GRAPPA.
Базовое аналитическое выражение для перехода из пространства изображений в К-пространство для протонной плотности Р(х У имеет следующий вид [121]:
Ввиду того что в (2.5) S(t) представляет собой волновой процесс, то вводится множитель e r(g, gyy+g2z)t учитывающий фазовые различия интегрируемых точек.
При выполнении реконструкции в параллельных МРТ-системах важным фактором является чувствительность катушки, сканирующей участок изображения. В этом случае уровень сигнала в К-пространстве, получаемого с катушки с чувствительностью Ь(х,у) имеет вид [68, 96]:
Выражение (2.6) является основополагающим для ряда методов реконструкции в МРТ-томографии. Дальнейшим развитием данного выражения является использование недекартовой (картезианской) системы координат. В частности, в работе [6] говорится о преимуществе реконструкции МРТ-изображения в формате 3D при использовании ядра Кайзера — Бесселя. Альтернативой методу Кайзера-Бесселя в работе [23] предлагается минимаксная интерполяция для случая оптимизации преобразования Фурье от неравномерных данных (NUFFT).
Важным вопросом работы многоканальных систем МРТ является оценка профиля чувствительности приемных систем сбора данных. Все методы параллельной реконструкции можно разделить на автоматически калибруемые методы и методы с предварительной калибровкой. Методы автокалибровки по произвольным траекториям К-пространства представлены в работе [56].
Зачастую реконструкция изображения МРТ выполняется от подвижных и динамических объектов, таких как сердце или кровеносные сосуды. Обычная Фурье-реконструкция в таких случаях сильно искажена (ввиду того что сбор данных занимает длительное время). Существуют специальные техники для синхронизации с процессами в объекте исследования. Теоретическое обоснование влияния движения объектов во время МРТ-сканирования представлено в работах [14, 120]. В работе [52] предлагается метод специальной модификации К-пространства перед выполнением Фурье-реконструкции.
Проведенный обзор литературных источников и экспериментальные исследования показали, что в основе большинства методов реконструкции лежит двумерное преобразование Фурье. Именно оно легло в основу более современного класса методов параллельной реконструкций изображений. Более подробно математическое обоснование параллельной реконструкции МРТ изображений будет дано в следующих главах.
Исследование программных библиотек для реконструкции и обработки МРТ-изображений
В настоящее время известно более десятка фундаментальных алгоритмов параллельной обработки МРТ-сигналов. В работе [85] сравниваются методы параллельной реконструкции на базе сопоставления теоретических моделей и практических экспериментов.
Группа исследователей [55] разработала специальную библиотеку на базе вейвлет-преобразований для декомпозиции и пост-обработки изображений. В частности, ими была показана высокая эффективность при обработке МРТ-снимков для функциональной диагностики головного мозга (f-MRI).
Наглядно изучить свойства, преимущества и недостатки данных методов позволяет пакет библиотек PULSAR Toolbox [37] для Matlab. Принцип работы математической библиотеки представлен на рисунке ниже:
Формирование данных Модуль реконструкции Реальные данные Оценка профиля чувствительности Оценка искажений SMASH ч г Оценка производительности GRAPPA Рис. 2.2. Структура библиотеки PULSAR Библиотека может работать как с экспериментальными данными, так и с фантомными изображениями. Более подробно работа с данной библиотекой будет проводиться для оценки погрешностей в главе 4.
Генерирование многоканальных тестовых изображений Одним из важнейших тестовых изображений, применяющихся в практической деятельности, является фантом Шеппа — Догана. Общий вид фантома Шеппа — Догана представлен на рис. 2.3 [10].
На рис. 2.3 видно, что фантом можно описать совокупностью эллипсов, ориентированных по отношению к выбранной системе отсчета под разными углами. При этом плотность вещества в пределах каждого из эллипсов можно считать постоянной величиной. В таком случае плотность изображения для фантома, состоящего из R эллипсов, будет иметь следующий вид: где г — радиус-вектор в системе отсчета, связанной с фантомом, описывающий положение точки на плоскости; Р , — область, ограниченная эллипсом, и значение плотности в ней соответственно.
Как видно из (2.5) и (2.7), основной особенностью генерации фантома Шеппа — Догана будет выполнение преобразования Фурье кусочно-постоянной функции плотности с учетом геометрии ограничений, накладываемых ею. Аналитическое представление результатов интегрирования может существенно повысить точность реконструкции и дать качественные предпосылки в дальнейшем исследовании.
Перспективным подходом к выполнению интегрирования может быть переход от интегрирования по эллиптической области к интегрированию по области с более простой геометрией. В работах [3, 98] предложен переход от эллипса к единичной окружности с центром в начале координат:
Ввод данных в ЭВМ
Задачу можно решать различными способами. Наиболее распространенными из них являются прямые методы, позволяющие переформулировать (2.69) как задачу метода наименьших квадратов. Это, кроме гибкой содержательной постановки, делает применение различных методов решения систем линейных уравнений. Например, метод сингулярного разложения (SVD-метод) и метод регуляризации Тихонова, заключающиеся в приведении системы линейных уравнений к некоторому специальному виду, что дает возможность решить их с высокой точностью при, как правило, низкой обусловленности главной матрицы системы.
Оператор А является матрицей, связывающей вектор измерений s(x) и неизвестное решение а. В этом случае CS-постановка примет следующий вид:
Требование редкости для а может носить и содержательный характер, обусловленный тем, что при сканировании может иметь место децимация данных. В этом случае требование редкости с учетом потерь информации помогает найти корректное решение задачи реконструкции.
Отметим, что кроме (2.71) возможна альтернативная постановка, заключающаяся в изменении характера целевой функции и выбранных ограничений. Она заключается в том, чтобы, как и в методе наименьших ІІФх - vll квадратов, минимизировать норму отклонения " Л2, но при условии, что вектор х должен быть достаточно редким. В этом случае, с учетом того, что редкость вектора мы будем оценивать как норму И 1% задачи (2.70), (4.71) могут быть переформулированы следующим образом: Фх-у\\2 -» min, XJ T (2.73) где т - показатель, учитывающий задаваемую в задаче редкость вектора. Задавая различную величину г , можно контролировать степень точности решения задачи СS-реконструкции. При этом постановку (2.65) иногда называют лассо-постановкой, что обусловлено аналогией с задачей о методе наименьших квадратов.
Для решения выражений (2.74)-(2.76) необходимо применить специальные методы поиска. Довольно подробно данный вопрос изучен в работе [64]. В работе [21] подробно рассмотрено несколько алгоритмов оптимизаций методом Брегмана. Учебное пособие [104] содержит описания большого числа методов оптимизации, некоторые из которых использовались при исследовании методов решения уравнений реконструкции МРТ-изображений. В диссертации использовался метод базиса преследования. 2.6. Разработка алгоритма параллельной реконструкции МРТ-изображений
В ходе дальнейшей работы было написано несколько тестовых программ для реализации акселерации метода SPACE-RIP. Первая тестовая программа основана на операторах линейной алгебры и прежде всего псевдоинверсии Мура — Пенроуза. Вторая тестовая программа основана на использовании CS постановке метода SPACE-RIP. При выполнении реконструкции с применением CS-постановки был выбран фреймворк SPGL1[91], осуществляющий CS-реконструкцию в случае задач большой размерности с высокой степенью разреженности данных. На основе данной математической модели был разработан алгоритм параллельной реконструкции МРТ-изображений, который показан на рис. 2.19
Алгоритм реконструкции начинает работу с инициализации параметров: в переменные загружаются такие параметры реконструкции, как число каналов, размер матрицы изображения для реконструкции, данные К-пространства реконструированного изображения и т. д.
После этого производится расчет профилей чувствительности приемных систем каждого канала. Данный расчет выполняется на базе тестовых изображений, выполненных в процессе предварительного (калибровочного)сканирования. НАЧАЛО
Сама реконструкция выполняется итерационно по каждому каналу последовательно. Каждая итерация заканчивается добавлением нового фрагмента изображений к результирующему.
Алгоритм реконструкции завершается процедурой сохранения реконструированного изображения в DCM-файл. Исходный текст алгоритма написан на языке MATLAB. Исходный текст программы представлен в приложении 1. где R(x,y)- эталонное изображение, Т(х,у)- тестовое изображение размерностью М х N. Помимо формул оценки (2.77-2.79) существуют и другие[114]. Например, в докторской диссертации [33] автор утверждает, что для некоторых случаев метод наименьших квадратов некорректен. В замен метода наименьших квадратов автор работы [33] предлагает метод CASE-PDM. Однако данный метод недостаточно документирован и его использование для оценки оказалось весьма затруднительно. Поэтому в данной диссертационной работе использовались более простые классические показатели качества изображений МРТ.
Сравнивая различные методы параллельной реконструкции, было замечено, что в отличие от обычного метода SENSE, отношение сигнал/шум в методе SMASH или AUTO-SMASH ниже, но оно менее зависит (ухудшается) при росте фактора ускорения R.
Реализация алгоритма параллельной реконструкции на основе графического процессора
Gadgetron [25] — набор инструментов для обработки данных, используемый для реконструкции медицинских изображений МРТ. Этот пакет предназначен для исследовательских целей. Набор разрабатывается с оглядкой на то, чтобы сделать прототипирование, отладку и добавление новых алгоритмов восстановления изображений наиболее простым.
Набор содержит несколько приложений для восстановления изображений, которые сразу готовы к использованию. Кроме того, набор содержит большой выбор общих структур данных и алгоритмов, которые позволяют создавать новые алгоритмы на основании уже существующих. Эти компоненты можно использовать в своих разработках как подключаемые библиотеки.
При составлении листинга программ в диссертации сборка Gadgetron инструментов изначально была начата на Windows 7 х64, т. к. операционная система семейства Windows наиболее привычна для использования. 64-битная версия была выбрана для наиболее полного использования возможностей современных компьютеров. В руководстве Gadgetron явно указывается, что сборка под ОС Windows весьма не тривиальная и требует много шагов и внимательности к деталям. Такая сложность вытекает из того, что Gagetron зависит от многих сторонних компонентов, каждый из которых должен быть установлен корректно на систему и доступен для других приложений.
Кроме сборки фреймворка, преследовалась цель использовать легальное программное обеспечение (т. е. использовать только программы и приложения, которые свободно доступны для использования). По этой причине использовался компилятор из свободной среды разработки Visual Studio 2010 Express. Эта среда разработки позволяет редактировать исходные тексты, компилировать и отлаживать относительно простые приложения под Windows. Существенным для этого проекта ограничением версии Express является отсутствие множества системных библиотек для х64 исполняемых файлов (набор доступных к использованию функций очень сильно ограничен). Этих ограничений нет в полной (платной) версии среды разработки Visual Studio 2010. Для сборки приложений, исполняемых на архитектуре х64, был дополнительно установлен набор разработчика Microsoft Windows 7 SDK. В этот набор входят компилятор, линкер, заголовочные файлы и множество библиотек, что позволяет создавать 64-битные исполняемые файлы и пользовательские библиотеки. В этом наборе нет среды разработки и средств отладки приложений. Вместе с этим, набор совместим с Visual Studio Express, что позволяет использовать среду разработки Visual Studio, а компилировать при этом 64-битный код.
После установки Visual Studio и Windows SDK появилась возможность собирать приложение, и дальнейшие шаги производились в соответствии с руководством Gadgetron. Поэтому были установлены следующие приложения:
К сожалению, достичь успешной сборки в данной конфигурации на W7-64 не удалось. Компилятор nvcc (собирающий код для устройств CUDA) проверяет версию компилятора, с которым его сопрягают, и в случае использования компилятора из набора SDK 7.1 завершается с ошибкой (несоответствия версий). Для того чтобы обойти эту ошибку, необходимо или переконфигурировать проект и задать другой компилятор (VC 100) или в уже существующем проекте изменить набор инструментов с SDK 7.1 на V 100. И в том, и в другом случае успешно завершить сборку Gadgetron не удаётся. Ошибка nvcc исчезает, т. е. зависимость компилятора, с которым сопрягается nvcc, удовлетворена, но линковка в дальнейшем не возможна, т. к. используется набор инструментов, неспособный полноценно собирать 64-битное приложение, а часть библиотек была предварительно собрана для 64-битной среды.
Для выхода из этой ситуации в диссертации сборка и дальнейшая работа с фреймворком Gadgetron была произведена в ОС Linux. Было подготовлено несколько тестовых программ. Экспериментальная оценка их приводится в последующих главах.
Реализация алгоритма параллельной реконструкции на основе графического процессора В диссертационной работе для реализации параллельных алгоритмов реконструкции была выбрана видеокарта NVIDEA и среда MATLAB. В данной среде, начиная с 2010b версии, возможны следующие варианты использования ресурсов видеокарты: - вызова пользователями CUDA-ядер (kernels); - вызов специальных процедур, встроенных в MathLab. Во втором способе MATLAB берет на себя все заботы по инициализации GPU, передаче данных из обычной оперативной памяти в память графической карты. Это удобно прикладным специалистам в своей области, которые не являются профессиональными программистами.
Для исследования возможности акселерации в диссертационной работе было написано несколько программ и проведено экспериментальное исследование.
Тестирование производительности реконструкции с использованием графического процессора Оценка производительности алгоритмов реконструкции МРТ-изображений производилась на нескольких компьютерах с разными конфигурациями.
Сборка фреймворка Gadgetron была произведена на компьютере с процессором Intel Dual Core El400 (2 ядра на 2 ГГц) и графической картой Geforce GTX 460 SE (288 CUDA ядер, частота ядер 1300 МГц, частота синхронизации графики 650 МГц, частота памяти 1700 МГц, тип памяти 1 Гб GDDR5). Данная видеокарта использовалась для всех последующих экспериментов. Оценка времени выполнения GRAPPA реконструкции на 116 каналах исходных данных матрицы размером 256x256 показала, что при использовании графического адаптера время сократилось с 1,5 секунд (без графического адаптера, т.е. при вычислении на центральном процессоре) до 0,2 секунд (при вычислении на центральном процессоре с использованием вызовов функций графического адаптера).
Перейдем к рассмотрению возможности акселерации в среде Matlab. Результат экспериментальной оценки реконструкции метода SENSE в среде Matlab представлен в таблице 4.1. Откуда видно, что время выполнения на компьютере с встроенным видео несколько выше, чем при использовании отдельного графического адаптера.
