Содержание к диссертации
Введение
1 Физика описываемых процессов 13
1.1 Томография 13
1.1.1 Классификация 13
1.1.2 Преобразование Радона 15
1.2 Синхротроны и лабораторные томографы 22
1.2.1 Монохроматор 25
1.2.2 Сцинтиллятор 27
1.2.3 ПЗС-камера 29
1.2.4 Оптическая система 31
1.3 Артефакты 32
2 Алгоритмы и численные методы 34
2.1 Вид кольцевых артефактов 34
2.2 Априорная информация 45
2.3 Одномерные проекции 51
2.4 Двумерные проекции 62
2.5 Нерегулярные кольцевые артефакты 67
2.6 Волновые артефакты 79
3 Программный комплекс 88
3.1 Современные вычислительные архитектуры 88
3.1.1 Классические процессоры 88
3.1.2 Графические процессоры 91
3.2 Кольцевые артефакты 92
3.2.1 Программы для обработки изображений 92
3.2.2 Программы для подавления кольцевых артефактов 96
Заключение 107
Список литературы 109
- Синхротроны и лабораторные томографы
- Априорная информация
- Нерегулярные кольцевые артефакты
- Программы для обработки изображений
Введение к работе
Диссертационная работа посвящена задаче уменьшения кольцевых артефактов, возникающих в рентгеновской компьютерной томографии. Для ее решения изучаются физические процессы, вызывающие данные артефакты. В зависимости от типа основного физического процесса предлагаются численные алгоритмы решения, запрограммированные для использования на обычных и графических процессорах. Техника распараллеливания и специализированные библиотеки программ позволяют производить коррекцию артефактов за время, малое по сравнению со временем реконструкции трехмерного объекта или его сечения.
Актуальность темы
Во многих отраслях науки, медицины и промышленности требуется определить внутреннюю структуру объекта, не нарушая его целостности. Компьютерная томография является одним из самых распространенных методов, используемых для решения таких задач. В данной работе в качестве источника излучения используются рентгеновские лучи, создаваемые лабораторными томографами и синхротронами.
Для восстановления структуры объекта в томографии за последние 30 лет были предложены различные методы реконструкции. Но в ходе эксперимента всегда имеются ошибки, вызванные как разнообразными дефектами компонентов установки, так и сложностью физических явлений, возникающих в процессе прохождения рентгеновских лучей через изучаемый объект, их преобразования в видимое излучение, а затем при регистрации изображений. Все это приводит к случайным и систематическим ошибкам во входных данных, в результате чего возникают артефакты при реконструкции объекта. Поэтому, даже если некоторый метод и позволяет реконструировать объект точно при отсутствии ошибок, то при использовании реальных данных получается сильно искаженная картина. В связи с этим задачи устранения различных артефактов в компьютерной томографии представляют собой очень большой интерес.
Среди разнообразия артефактов выделяются кольцевые артефакты, как наиболее распространенные и сильно затрудняющие анализ объекта. Они представляют собой тонкие концентрические кольца или дуги, наложенные на обычное реконструируемое изображение — сечение объекта. В литературе предложено множество способов уменьшения таких артефактов, но они в основном базируются на сглаживании входных данных или реконструированного изображения. Такое сглаживание приводит к потере некоторых (обычно малых по размеру)
деталей. Кроме того, многие методы не полностью убирают указанные артефакты и даже могут привести к появлению новых.
Цель работы
На основе имеющейся априорной информации о физических процессах, происходящих во время сканирования объекта, разработать методы, более эффективно подавляющие кольцевые артефакты.
Положения, выносимые на защиту
Получена формула, задающая общий вид кольцевого артефакта при параллельной геометрии рентгеновских лучей. Использование данной формулы совместно со знанием коэффициента поглощения в некоторых областях двумерного сечения объекта позволяет быстро удалить артефакты на всем сечении с точностью, приемлемой для практических приложений.
Сформулирована задача подавления кольцевых артефактов в виде задачи минимизации функционала Тихонова. Для случая одномерных проекций была предложена аналитическая формула для нахождения точки минимума данного функционала, в случае же двумерных проекций была предложена формула для быстрого нахождения градиента функционала.
Разработан метод подавления нерегулярных кольцевых артефактов, интенсивность которых меняется в зависимости расположения относительно оси вращения.
Предложен метод коррекции кольцевых и волновых артефактов, возникающих в результате колебаний рентгеновского пучка (дрожания ПЗС-камеры, оптической системы, сцинтиллятора), основанный на использовании информации о характере колебаний и дополнительных измерениях интенсивности излучения.
Создан программный комплекс для подавления кольцевых артефактов на обычных и графических процессорах с применением многопоточных библиотек компаний Intel и nVidia.
Научная новизна
В работе был использован метод регуляризации Тихонова для подавления кольцевых артефактов. Полученная аналитическая формула для корректировки регулярных артефактов позволяет существенно увеличить скорость и точность корректировки кольцевых артефактов по
сравнению со стандартными численными методами. Предложенные алгоритмы позволяют использовать параллельные вычисления и были реализованы как для обычных многоядерных компьютерных процессоров, так и для графических карт.
Практическая ценность
Результаты, полученные в диссертации, могут быть применены как в лабораторных, в том числе медицинских, томографах, так и на синхротронах для улучшения качества реконструируемых объектов.
Личный вклад автора
Основные результаты диссертации, в том числе положения, выносимые на защиту, были получены автором лично. Проведение экспериментов по сканированию трехмерных объектов на синхротронах и лабораторных установках, а также интегрирование комплекса программ, предназначенного для подавления кольцевых артефактов, в общие комплексы программ, используемых для реконструкции объектов, было произведено соавторами статей.
Апробация работы
Основные результаты диссертационной работы были представлены: на семинарах "Сотрудничество для успеха" по проблемам компьютерной томографии, организованных при поддержке Научного совета по инженерным и физическим наукам (Великобритания) в Уилмслоу (графство Чешир), 9-11 января 2008 г. и 15-17 марта 2009 г.; в лаборатории 112 синхротрона "Даймонд" (Diamond Light Source, Дидкот, графство Оксфордшир, Великобритания), 31 июля 2008 г.; на Британском коллоквиуме по прикладной математике в Университете Ноттингема (Великобритания), 7-9 апреля 2009 г. [6]; на ежегодной конференции "Женщины в математике", проводимой Лондонским математическим сообществом, 24 апреля 2009 г.; в Университете Манчестера (Великобритания) на семинарах школ математики и наук о материалах, а также на конференции студентов и аспирантов по математике, 18 мая 2009 г. [7]; в Аргоннской национальной лаборатории (США) в марте 2010 г.; в медицинском центре Университета Чикаго (США) на семинаре отделения радиологии, 18 марта 2010 г.; в Университете Шеффилда (Великобритания), в школе математики и статистики, 9 июня 2010 г.; на конференции "SPIE Optics & Photonics 2010" в Сан Диего (Калифорния, США), 1-5 августа 2010 г. [8]; на научном семинаре "Обратные задачи математической физики" в НИВЦ МГУ под научным руководством д.ф.-м.н.
проф. А. Г. Яголы, А. Б. Бакушинского, А. В. Тихонравова, 29 сентября 2010 г.; на научном семинаре кафедры математики физического факультета МГУ 29 сентября 2010 г.
Публикации
По теме диссертации было опубликовано 8 работ, из которых 2 тезиса конференций [6, 7], 1 статья в трудах конференции [8] и 5 статей в рецензируемых отечественных и зарубежных научных журналах [1-5].
Структура работы
Диссертация написана на 121 странице, состоит из титульного листа, оглавления, введения, трех глав, заключения и списка литературы (101 наименование).
Синхротроны и лабораторные томографы
Пусть имеется задача: определить внутреннее строение объекта, не нарушая его целостности. Решениями таких задач занимается томография (от греческого слова "тодоо-" — сечение). Задачи томографии имеют широкую область применения в разных областях науки. Наиболее широко распространенными задачами являются изучение человеческого тела в медицине [15,68], растений и животных в биологии [48], горных пород в геологии [66], океана [9,71], Земли [75], ионосферы [14], плазмы в физике [19,20], структуры Солнца в задачах астрофизики [53]. Задачи томографии также встречаются в атомной энергетике [5], сей- смике [10,12], при обнаружении подводных объектов [23]. Для решения этих задач было разработано множество математических методов, см. например [1,2,8,11,18,28-31,34,35,62,64,73]. Приведенные монографии не претендуют на полноту обзора задач томографии в соответствующих областях науки.
В томографии выделяют несколько основных типов зондирующего излучения: ультразвуковые волны [81]; электромагнитное излучение (гамма-кванты в позитронно-эмис- сионной томографии [98], рентгеновские лучи [59], лучи оптического диапазона [15,48]); электромагнитные поля (магнитно-резонансная [77] и злектро- импедансная томография [56]); элементарные частицы (нейтроны [36], протоны [61], нейтрино [50]). В данной работе рассматривается только рентгеновская томография, что конечно не ограничивает возможности применения полученных результатов в других областях, например в оптической томографии. По типу регистрируемого излучения можно выделить два основных типа томографии: трансмиссионную, когда происходит измерение внешнего излучения после прохождения и ослабления внутри объекта, и эмиссионную, когда регистрируется излучение, созданное внутренними источниками под воздействием внешнего излучения. Возможны следующие основные типы геометрии сканирования, см. рис. 1: параллельная, когда все лучи, проходящие через объект являются параллельными (или почти параллельными), например, на синхротронах расстояние от источника излучения обычно около 50 метров, а размер объекта не превышает 5 миллиметров, поэтому угловое расхождение не превышает 0,01, такой же порядок наблюдается для лабораторных установок, используемых в нанотомографии, когда расстояние между источником и объектом порядка одного метра, а размер самого объекта — несколько сотен микрометров; коническая, когда лучи, испускаемые рентгеновской трубкой, выходят из одной точки, данная геометрия характерна для лабораторных установок, применяемых в микротомографии, при этом проекции точек пересечения лучей и объекта на плоскость, перпендикулярную центральному лучу, могут находится на расстоянии в несколько процентов (иногда до 20%) от размера объекта; спиральная, когда источники излучения и детекторы вращаются вокруг оси, вдоль которой движется объект (в этом случае в стандартной системе отсчета, связанной с объектом, источники и детекторы движутся по спирали), данная геометрия популярна в медицинской томографии. Для каждой из этих геометрий разработаны свои методы реконструкции объекта. Наиболее популярными методами являются метод обратного проецирования фильтрованных проекций (ОПФП) для параллельной геометрии, см. например [17,73], и метод ФДК (сокращение от имен авторов: Фелдкамп, Дэвис, Кресс [51]) для конической геометрии. Отметим, что не всякая геометрия сканирования позволяет полностью восстановить объект, см. более подробно в [64]. К примеру в случае конической геометрии точная реконструкция невозможна, а метод ФДК является приближенным методом. В диссертации рассматриваются только задачи трансмиссионной томографии для параллельной и конической геометрий. Пусть Ь — луч, вдоль которого распространяется монохроматическое рентгеновское излучение. Тогда если /о — интенсивность излучения до объекта, то интенсивность излучения сразу после прохождения объекта равна от » где р = у ц{х)д.х — оптический путь, л{х) — линейный коэффициент поглощения, а интегрирование производится вдоль луча Ь. Отметим, что формула записана без учета волновой природы света и в общем случае является только первым приближением. Для более точной аппроксимации необходимо учесть дифракцию света на объекте, а в для получения точного значения решить соответствующее волновое уравнение. В классической компьютерной томографии стараются измерять интенсивность света сразу за объектом. Учет же поправок, связанных с волновой природой света, привел к новому направлению в томографии, фазово-контрастной томографии, которая в основном используется при нахождении границ сред с близкими коэффициентами поглощения [38].
Рассмотрим классический вариант сканирования в рентгеновской томографии с параллельной геометрией лучей, см. рис. 2. Пусть ось вращения объекта вертикальна, т.е. параллельна оси Ох системы координат, все рентгеновские лучи параллельны оси Оу, а плоскость детектора параллельна плоскости хОх. Тогда рассматриваемый трехмерный объект делится на множество двумерных срезов, параллельных плоскости хОу и независимых друг от друга. Таким образом, для каждого и фиксированного значения % детекторы измеряют синограмму — двумерное изображение одномерных проекций среза объекта как функцию проекционного угла (проекционный угол в изображается по оси ординат, а линейная проекционная координата х изображается по оси абсцисс). Для простоты рассмотрения считаем, что г = 0 и точка О принадлежит пересечению проекции оси вращения на плоскость детектора.
Введем также систему координат, связанную с вращающимся объектом, с центром в точке О пересечения вертикальной оси и рассматриваемого среза, при этом при 9 = 0 координатная ось 0 сонаправлена с осью Ох, а ось О г) — с осью Оу.
Априорная информация
В многослойных монохроматорах используется принцип интерференционного отражения. Имеется система параллельных плоскостей с малыми коэффициентами отражения, почти весь свет, падающий на каждую плоскость, проходит внутрь без потерь, в результате чего образуется система параллельных лучей. Разность хода между двумя отраженными лучами равна Поэтому при выполнении условия
Брэгга-Вульфа все лучи будут в одной фазе. При фиксированном п условие Брэгга-Вульфа может выполняться только для волн определенной длины А. В результате при освещении системы белым светом интерференционное отражение наблюдается только для некоторых сравнительно узких спектральных линий [22]. Отметим, что монохроматический пучок, полученный с помощью многослойного монохроматора, обычно в 100 раз интенсивнее пучка, полученного с помощью кристаллического монохроматора [46]. Поэтому в связи с улучшением технологий производства многослойных рентгеновских покрытий наблюдается переход к многослойным монохроматорам. К сожалению, у многослойных моно- хроматоров есть и недостаток: ширина спектра энергий в 100 раз шире, чем в случае кристаллического монохроматора, например на станции 2-ВМ относительная ширина спектра для многослойного монохроматора АЕ/Е равна Ю-2, а для кристаллического — Ю-4.
Отметим, что монохроматор часто состоит из двух частей, наклон которых относительно падающего пучка и расстояние между которыми регулируется с помощью высокоточных моторов. Для уменьшения вибраций вся система охлаждается до сверхнизких температур, часто с помощью жидкого азота [44]. К сожалению, небольшие угловые колебания частей монохроматора возможны, что с учетом больших расстояний до образца (порядка 30 м) и малых размеров сенсора камеры (до 10 мкм) приводит к изменению профиля интенсивности с течением времени (обычно в вертикальном направлении). Окно из бериллия тоже испытывает некоторые колебания в связи с нагреванием, вызываемым рентгеновских излучением, и охлаждением, осуществляемым жидким азотом.
Сцинтилляторами называют вещества, способные излучать свет под действием ионизирующего излучения. В томографии рентгеновское излучение выбивает электрон с внутренней орбиты атома, что вызывает переход другого электрона на вакантное место с испусканием фотона меньшей энергии (часто видимого излучения). Сцинтилляторы могут быть как неорганическими, так и органическими. Для синхротронов и лабораторных томографов используются неорганические кристаллы. Некоторые характеристики, учитывающиеся при выборе сцинтиллятора: световыход — количество высвечиваемых фотонов на единицу поглощенной энергии, спектр высвечивания, время высвечивания — время, в течение которого пологощенная энергия рентгеновского излучения преобразуется в видимый свет (для задач сверхбыстрой томографии, когда регистрируются более 1000 изображений в секунду, очень важно, чтобы каждое новое изображение не зависело от ранее зарегистрированных), радиационная прочность, так как под действием падающего излучения происходит нагрев сцинтиллятора, формирование дефектов в кристаллической структуре, толщина. Отметим, что декартову систему координат выбирают таким образом, что при параллельной геометрии рентгеновских лучей координатная ось Оу параллельна этим лучам, а ось Ог — оси вращения объекта. Для увеличения разрешения экспериментальной установки необходимо уменьшить путь, который проходят рентгеновские лучи через сцинтиллятор. Поэтому сцинтиллятор выравнивается перпендикулярно этим лучам, т.е. параллельно плоскости хОх.
Пусть Ксщтт (х, г) — функция распределения интенсивности видимого излучения от точечного источника рентгеновского излучения единичной интенсивности. В большинстве случаев эта функция является радиально-симметричной. Для грубого сравнения различных оптических элементов, к которым, конечно, относится и сцинтиллятор, используется понятие полуширины функции Ксцижт_(х, г) (ширина функции на половине максимальной высоты) [37].
При прохождении фотона рентгеновского излучения через сцинтиллятор может происходить цепное высвечивание фотонов меньшей энергии, при чем не всегда вдоль направления луча. Поэтому при увеличении толщины сцинтиллятора увеличивается и полуширина функции КсЦинт.(я,2)» что приводит к уменьшению разрешающей способности системы, но при этом увеличивается световыход. Поэтому очень важно подобрать толщину сцинтиллятора для выбранной задачи. Улучшение технологии выращивания кристаллов в некоторых случаях позволяет увеличить толщину сцинтиллятора без уширения функции /Ссцинт.(ж, г). Например в [69] описано использование кристаллических игл иодида цезия (СБ1): диаметр каждой иглы примерно 5 мкм, и за счет явления полного внутреннего отражения луч видимого света не выходит через боковую поверхность, таким образом полуширина функции Ксщшт{х, г) будет равна диаметру иглы. Для некоторых задач такое разрешение системы бывает вполне удовлетворительным, для других же, например задач нанотомографии, такая полуширина в десятки, а иногда и сотни раз превышает допустимую.
Разработка теоретических моделей, позволяющих найти функцию Ксципт.(х, г), — очень сложная задача [101], а проведение эксперимента по измерению этой функции сопровождается техническими трудностями. Поэтому зачастую делаются некоторые предположения, например функция не зависит от положения рентгеновского луча и имеет радиальную симметрию, т.е.. Во время эксперимента используют свинцовую пластину с очень ровным краем [45] и затем, решая интегральное уравнение, находят функцию А;сцинт.(г) или ее образ Фурье.
Нерегулярные кольцевые артефакты
Существующие методы коррекции Кольцевые артефакты являются очень распространенными, поэтому не удивительно, что существует довольно много различных методов их коррекции. Большинство таких методов опубликованы в различных монографиях, статьях или представлены в качестве докладов или тезисов на конференциях. Другая же часть методов была разработана в различных компаниях по производству томографического (в том числе медицинского) оборудования и представляет собой коммерческую тайну, так как позволяет получить существенные преимущества при продаже оборудования вместе с программным обеспечением. К сожалению, даже для опубликованных методов невозможно найти соответствующие программы, чтобы сравнить преимущества и недостатки того или иного метода, — такие программы обычно являются частью программных комплексов, разработанных для больших лабораторий, и дают конкурентные преимущества при привлечении пользователей.
Отметим, что задача коррекции кольцевых артефактов является некорректной, а поэтому без знания дополнительной информации о решении, а также уровня погрешностей входных данных, невозможно найти приближенное решение, сходящееся к точному при стремлении погрешностей к нулю. Даже если для реальных задач уровень погрешностей известен, то он имеет какое-то фиксированное значение и не стремится к нулю, поэтому невозможно оценить близость найденного приближенного решения к точному [99].
В некоторых случаях очень сложно понять, является ли та или иная особенность на изображении артефактом или реальным объектом. Это особенно сложно для объектов с цилиндрической симметрией, например, при сканировании биологических объектов они помещаются в цилиндрические пробирки или трубки, а ось вращения пробирки зачастую совпадает или близка к оси вращения всего объекта. В этом случае проекция реального объекта (пробирки) выглядит также как и в случае кольцевого артефакта. В других случаях могут наблюдаться сильные изменения коэффициента поглощения вблизи оси вращения объекта. Например, при изучении пористой структуры стебля в питательную среду добавляют частицы металлов (часто золота), которые потом оседают внутри клеток стебля. Такая частица, когда она находится близко к центру вращения объекта, дает большой скачок на синограмме. В той же мере это относится и к пустотам внутри стебля. Поэтому без наличия какой-то дополнительной информации о структуре объекта или проведения дополнительных измерений отличить артефакт невозможно. К сожалению, проведение дополнительных измерений в некоторых случаях невозможно, так как восстановление структуры объекта часто происходит после проведения эксперимента.
Существующие методы коррекции кольцевых артефактов могут быть условно поделены на три группы. 1) Подавление артефактов происходит до реконструкции изображения. Часто предполагается некоторая гладкость синограммы. В [40] для каждого п-ого столбца синограммы рп{г) находятся средние значения у (г), после чего применяется специальный движущийся фильтр, который заменяет значение у (г) значением 2/5(г), равным усредненному значению у (г) для 2N + 1 соседних значений. В результате исходная синограмма заменяется нормализованной по формуле рп(г) = рп(г) 2/5(г)/?/(г). В [97] интенсивность, измеренная в каждом пикселе, заменяется усредненным значением в восьми соседних пикселях. В [70] применяется специальный фильтр, основанный на использовании вейвлетов и быстрого преобразования Фурье. В [80] описывается применение итеративного взвешенного медианного фильтра. 2) Подавление артефактов осуществляется после реконструкции объекта [84,100]. В [84] восстановленное изображение преобразуется в полярную систему координат, где происходит обработка на основе медианного фильтра. 3) Изменение процедуры сканирования или измерительных приборов. В [47] детектор постоянно сдвигается, таким образом артефакт как бы "размазывается" по различным областям ПЗС-матрицы. В [74] предложен специальный гибридный сенсор. Изучение статистики генерируемых электронов приводит к алгоритму нахождения дополнительных корректирующих множителей для регистрируемых проекций. В [60] был разработан специальный твердотельный детектор, а также предложен метод коррекции послесвечения детектора, характерного для многих сцинтилляторов. Этот метод основан на интеративной фильтрации регистрируемых изображений и также используется для подавления возникающих кольцевых артефактов. Конечно, некоторые методы используют идеи из нескольких групп. Отметим также, что кольцевые артефакты могут быть вызваны разными физическими процессами и поэтому требуют различных методов коррекции. Например, кусок грязи или пыли на сцинтилляторе уменьшает интенсивность проходящего через него луча света стандартным способом (как экспоненциальная функция от толщины куска) и не зависит от общей интенсивности рентгеновского излучения, падающего на сцинтиллятор. С другой стороны при наличии трещины внутри сцин- тиллятора интенсивность излучаемого видимого света часто пропорциональна всей интенсивности падающего излучения. Это приводит к тому, что при сканировании объекта, вытянутого вдоль какой-нибудь горизонтальной оси, интенсивность наблюдаемого кольцевого артефакта зависит от угла поворота объекта. Поэтому сложно, если вообще возможно, разработать общий метод коррекции кольцевых артефактов, который мог бы быть применен как для различных по структуре объектов, так и для различных экспериментальных установок.
Программы для обработки изображений
Видно, что использование библиотеки MKL дает наибольшее увеличение скорости работы программы. Можно достичь еще большего увеличения скорости, если воспользоваться библиотекой ТВВ (Threading Building Blocks) компании Intel, которая позволяет автоматически создавать многопоточные программы. Так, например, команда выполняет некоторую операцию с помощью нескольких потоков команд. Для заданного индекса в "цикле" используется один поток, разные индексы могут соответствовать разным потокам, последовательность перебора индексов не гарантирована, т.е. может быть 1,2,3,4.
У пользователя есть возможность как самому задать число используемых потоков, так и предоставить это системе. Примеры зависимостей, найденных для разного числа потоков, показаны на рис. 24. Отметим, что увеличение пх и числа потоков приводит к увеличе нию размера "скачков" на кривых. Конечно, имеются некоторые флуктуации при вычислении времен, так как ЦП работает на своей максимальной мощности (его загруженность достигает 99%) и ему иногда необходимо отвлекать ресурсы на работу как самой операционной системы, так и других программ. Но многократный повтор вычислений и использование большого числа (1000) корректируемых сечений позволяет убрать эти флуктуации. Таким образом, положение и глубина "провалов" является фиксированной для заданной системы и определяются как самой библиотекой МКЬ, так и используемыми процессором, памятью и материнской платой компьютера. При обработке нескольких сечений объекта каждый поток работает с одним из сечений не зависимо от других потоков. Поэтому следует допускать, что в какой-то момент времени два разных потока обратятся к одним и тем же элементам матрицы 5ут для последующей работы на (возможно) разных ядрах ЦП. В этом случае возможен некий конфликт и один из потоков подождет, пока другой поток закончит использование данных элементов. Тогда воизбежание естественного простоя одного из потоков наиболее логичным является увеличение копий одной и той же матрицы так, что каждый поток работал с одной из копий. Но для использованного четырехядерного процессора такой подход наоборот привел к замедлению работы программы, см. рис. 25. Это может быть объяснено тем, что вместо последовательного доступа к памяти, когда система в основном обращается к близким элементам памяти, происходит увеличение "случайного" доступа к элементам, находящимся на удалении друг от друга. На рис. 26 показано увеличение скорости работы программы при использовании библиотек 1РР и МКЬ. Видно, что применение библиотеки МКЬ вместе с автоматическим распараллеливанием с помощью ТВВ дает максимальное увеличение скорости вычислений. Так, при больших значениях пх использование IPP дает увеличение в 1,13 раз, MKL (од- нопотоковая программа) — 1,97, a MKL + ТВВ — примерно в 3 раза. Для малых значений пх (пх 1000) увеличение скорости более заметно и может быть объяснено правильной организацией памяти при использовании библиотек компании Intel. Теперь обратимся к графическим процессорам. Компания nVidia предоставляет библиотеку cuBLAS, содержащую некоторые функции из набора подпрограмм BLAS, которые были реализованы для использования на ГП. В качестве графической карты использовалась видео-карта GeForce GTX 280, содержащая 240 ядер, объединенных в 30 мультипроцессоров. На рис. 27 показаны времена работы программ на центральном и графическом процессорах. Видно, что для графических процессоров тоже характерно появление "пиков" и "провалов" во временах работы. Здесь справедливы те же аргументы, что были приведены для централь-Увеличение скорости работы программы при использовании библиотек 1РР и МКЬ (в случае, когда используется один поток или число потоков выбирается системой автоматически) ного процессора. Можно также отметить, что у графических процессоров также имеются дополнительные особенности такие, как ограниченное число регистров (8192) на один мультипроцессор, которые и производят непосредственно арифметические вычисления; объединение потоков в группы по 32 потока и ограничение числа активных потоков и групп для одного мультипроцессора; большое время (от 400 до 600 тактов процессора), необходимое для извлечения данных из глобальной или локальной памяти; необходимость копирования данных из центральной памяти компьютера в память графической карты. В любом случае времена работы программ на обеих архитектурах сопоставимы. При замене графической карты на более мощную (например GeForce GTX 295 или 480) время работы программы на ГП должно стать меньше, чем на ЦП. Обратимся теперь к методу, изложенному в разделе 2.4. Будем считать, что регистрируемое изображение — квадрат, т.е. nz = пх. Дискретный функционал Тихонова минимизировался с помощью градиентных методов: метода наискорейшего спуска (мне) и метода сопряженных градиентов (мсг). На рис. 28 показано, как зависит время, затрачиваемое на одну итерацию, от размера изображения. Видно, что для при пх 4900 программа выполняется быстрее на ЦП, чем на ГП. При увеличении пх графическая карта обеспечивает большую скорость, чем центральный процессор.
