Семинар "Суперкомпьютерные технологии в науке, образовании и промышленности", 29 октября 2013 г .

Суперкомпьютерное моделирование в обратных задачах ультразвуковой томографии

Гончарский А.В.

Романов С.Ю.

(докладчик)

Московский государственный университет, НИВЦ

Постановка задачи волновой томографии

проходящая волна регистрируемый сигнал проходящая волна приемник-передатчик

•

Постановка задачи во временной области

c

(

u c

( (

r r

) •

Волновое уравнение

r

)

u tt

(

r

,

t

)  

u

(

r

,

t

)    

u



t

0 )  ,

t c

0  0 )

const

(

r

,

t

, при

r



R

-

(

r

  0

q

) 

f

(

t

)

q



X

положение источников известна вне неоднородности

u

(

r

,

t

) 

U

(

r

,

t

),

r



Y

,

t

 0

известна на приемниках

•

Требуется найти с(r) в R. Обратная задача нелинейная.

Приложение к ультразвуковой томографии в медицине

•

Более 40 тысяч женщин России ежегодно заболевают раком молочной железы. Доля лиц с поздними стадиями заболевания среди первичных больных превышает 40%.

•

Ультразвуковые томографы высокого разрешения позволят осуществлять раннюю диагностику рака .

•

Каждое 4 е онкологическое заболевание – рак груди.

•

Почему не X-Ray томография и не МРТ?

Невозможность проведения безопасных регулярных рентгеновских обследований из-за дозы получаемого излучения.

МРТ – высокая стоимость.

•

Почему не использовать стандартные ультразвуковые сканеры?

Все УЗИ не являются томографическими. В обычных УЗИ, информация получается с одного или нескольких ракурсов. В томографах используется обзор с 360 градусов. За счет этого высокая разрешающая способность.

В США и Европе ведутся интенсивные работы по разработке программного обеспечения и макетов УЗ томографов

Постановка задачи волновой томографии.

Градиентные методы решения

•

Прямая задача

c u

(

r

(

r

)

u tt

,

t

 (

r

, 0 )

t

)  

u t



u

(

r

,  

n u

|

ST



p

(

r

,

t

)

t

)  0 )  (

r

 0 

r

0 ) 

f

(

t

) •

Сопряженная задача

c

(

r

)

w w



n

(

r

,

t w

|

ST tt

 (

r T

,

t

) ) 

u

 

w

(

r

,

t

| 

ST w



t

(

r U

,

t

) 

T

 ) 0  0 •

Обратная задача - минимизация функционала невязки

 (

u

(

c

)) 

u

|

ST



U

2 •

Градиент



C

(

u

(

c

) ,

функционала имеет вид:

T c

)  

w t

(

r

, )

u t

(

r

,

t

)

dt

0

здесь

u

(

r

,

t

)

w

(

r

,

t

)

решения основной и сопряженной задач

•

При проведении расчетов прямой задачи брали граничное условие неотражения на границе области расчетов

c

(

r

) 

n u

|

ST

 

t u

|

ST

•

Т – выбирали достаточным для прохождения основных прямых и отраженных волн.

Численный алгоритм

• •

Аппроксимации частных производных 2-го порядка точности Явная разностная схема для расчета по временным слоям распространения звуковой волны (расчет «в прямом и обратном времени»)

u ijlk

 1 

u ijlk

( 2  8  2

с ijl h

2 )  (

u i

 1

jlk



u i

 1

jlk

)  2

с ijl h

2  (

u ij

 1

lk



u ij

 1

lk с ijl h

2 )  2  (

u ijl

 1

k



u ijl

 1

k с ijl h

2 )  2 

u ijlk

 1

Условие неотражения на границе



n u

|

ST

 

c

0

.

5 

t u

|

ST

Условие устойчивости Куранта

c

 0 .

5  

h

3

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

grad ijl



k m

  0

u ijlk

 1  

u ijlk w ijlk

 1  

w ijlk



Итерационная последовательность: 1. В качестве начального приближения используется С 0 =const. 2. Для С 0 с помощью явной разностной схемы решается прямая задача, находим волну u(r,t) на каждом из детекторов. 3. Для функции u(r,t), решается сопряженная задача по явной схеме, получаем w(r,t). 4. Используя полученные значения u(r,t) и w(r,t), вычисляется градиент .

5. Зная градиент в точке С 0 , находится минимум функционала по параметру k в области k>0. 6. Точка минимума функционала принимается за С 1

Φ

(

c (

0

)



k Φ



C

(

c (

0

)

))

. Процесс возвращается к пункту 1.

Обратные задачи УЗ томографии.

Модельные расчеты

Модельный объект.

Восстановленная структура томографического слоя

• • •

Начальное приближение = const = C 0 Количество источников – 8, количество приемников – 320.

Время - 4 ч на 512 ядрах, 700 итераций, невязка в 100 раз.

• • •

Сетка по XY – 1000x1000.

Разрешение 1-3 мм.

Длина импульса 5 мм.

Реконструкция по независимым данным

•

Проблема: Провести реконструкцию по независимым данным.

•

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

•

Хорошее согласие с численными расчетами. График. Вид преломленной волны на границе. Аналитические и численные расчеты.

Модель Результаты реконструкции

Томографические схемы на прохождение и на отражение

Почему важно использовать полный набор данных?

полные данные на прохождение на отражение стандартный УЗИ сканер При полном наборе данных:

•

высокое разрешение,

•

абсолютные значения скорости

Томографические схемы.

Послойная реконструкция

Послойная реконструкция:

•

Используется в рентгеновской томографии.

• • •

Исходная трехмерная задача разбивается на N (N – число слоев) независимых двумерных обратных задач.

•

Сетка по XY – 1000x1000.

Решать двумерные задачи проще.

Огромный объем вычислений

• • •

Сечений – 40.

Количество источников – 8.

Количество приемников – 320.

Структура распараллеливания

Источник Приемники Область независимых вычислений

•

Область обменов Задача в целом Слой 1 Слой 2 … Слой N Источник 1 Источник 2 … Источник K Часть 1 Часть 2 … Часть M

•

Максимальное число ядер = число слоев х число источников х число частей

• • • •

Подход позволяет легко практически неограниченно распараллеливать вычисления.

Эффективность распараллеливания 60% (ускорение в ~ 12000 раз) на 20480 ядрах. Высокая масштабируемость для десятков тысяч процессов.

Подход позволяет свободно маневрировать, при наличии ограничений на число процессоров.

•

Масштабируемость программы

• • •

Выравнивание загрузки ядер -> разбиваем матрицу на блоки равного размера.

Минимизации межпроцессорных обменов -> минимизируем периметр блоков по отношению к их площади. Варианты разбиения сетки Время, сек 1*64 692 2*32 386 4*16 254 8*8 220

•

Таблица - Сравнение вариантов разбиения сетки на блоки, размер сетки – 1002х1002, число процессов (блоков) – 64.

Диаграмма - Повышение масштабируемости на одной матрице по мере модернизации программы

250 200 150 100 basic opt_v1 opt_v2 opt_v3 opt_v4 50 0 144 196 256 324

Число процессов Таблица - Результаты тестирования на суперкомпьютере «Ломоносов» Число процесс ов N

proc

Время 15 итерац ий сек 1 1893 2*2 3*3 4*4 903 370 291 5*5 6*6 162 121 7*7 67 8*8 64 9*9 44 10*10 40 11*11 41 12*12 42

400

13*13 44

484

14*14 45

Эффективность программы

•

E

.

E



T

0

T



K N proc



T

• •

N proc

• • • • •

E -

эффективность;

N proc T 0

число процессов; время расчетов на N

proc

производительностью; с пиковой T = T(N

proc

) экспериментальное время расчета на N

proc

; К = const - время расчетов с одним процессом с пиковой производительностью.

•

Зависимость эффективности E от N одной матрице: теоретическая (сплошная ) и экспериментальная (пунктир).

proc

на Уменьшение Е из-за увеличения времени обменов по отношению к вычислениям.

Локальный максимум Е из-за увеличением доли вычислений в кэше.

При построении теоретического графика использована простейшая параметрическая модель и типичные значения параметров для процессоров общего назначения.

Производительность (ускорение) программы

Пр

Пр = const / T;

• •

N proc

число процессов;

T = T(N proc

) время расчетов N

proc N proc

Зависимость производительности системы от числа процессоров на одну матрицу (пунктир).

•

Производительность растет с линейной скоростью вплоть до 100 процессов на матрицу.

•

Знание зависимости эффективности и производительности от N на число доступных ядер.

proc

позволяет выбирать стратегии распределения ядер при наличии и отсутствии ограничений

Решение задачи на 20480 ядрах суперкомпьютера «Ломоносов»

•

Впервые проведен масштабный вычислительный эксперимент на 20 480 ядрах суперкомпьютера по послойной реконструкции модельного 3D объекта в УЗ томографии .

Верхний ряд модельные сечения 3D объекта.

Нижний ряд восстановленные сечения.

• • •

40 сечений, шаг между сечениями 5 мм. Источников – 8, приемников – 320, импульс - 5 мм.

• •

Сетка по XY – 1000x1000.

Разрешение 1-3 мм.

Время расчета 4 ч.

Послойная реконструкция vs. непосредственно 3D реконструкция

Послойная томография 3D томография

•

Послойная реконструкция: исходная трехмерная задача разбивается на N (N – число слоев) независимых двумерных обратных задач.

•

3D реконструкция: решается трехмерное уравнение.

1 • • • •

Решать двумерные задачи проще.

Послойная реконструкция хорошо подходит для X-ray, т.к. коэффициент рефракции ~ 1.

Для ультразвука существенны эффекты дифракции, рефракции, переотражения, излучение не сохраняется в слое.

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

Послойная реконструкция vs. непосредственно 3D реконструкция

•

Проблема: Насколько качественна послойная реконструкция 3D объекта в УЗ томографии?

•

Модельный эксперимент для сферы.

•

Экспериментальные данные получены из точных аналитических формул в 3D.

•

В случае < 0.5R: качество реконструкции хорошее.

• •

В случае > 0.5R: искажение геометрической формы, появление артефактов, размывание вдоль оси z.

Для получения высокого разрешения необходимо решать 3D.

1

Верхний ряд модельные сечения 3D шара радиуса 6см.

Z = 0 см Z = 3 см Z = 4 см Z = 5 Нижний ряд восстановленные сечения см Z = 6 см

Структура распараллеливания в 3D

Задача в целом Источник 1 Источник 2 … Источник K параллелепипед 1 параллелепипед 2 … параллелепипед M

• •

Число источников = числу карт GPU.

Графические карты NVIDIA Tesla X2070 суперкомпьютер «Ломоносов», С/С++/OpenCL.

•

Обмен данными между узлами по MPI занимает менее 2% от времени расчета каждой карты Рис. Схема параллельных расчётов для нескольких источников одного шага градиентного спуска

Структура распараллеливания в 3D

X Y Z Рис. Процедура вычислений по явной по времени трёхмерной разностной схеме Рис. Схема расчета одной подзадачи (для одного источника).

• • • •

Явная разностная схема – высокая степень параллелизма. Около 10 8 вычислимых операций параллельно Область размера 400х400 в плоскости ХУ разбивается на прямоугольные блоки размера 32х16 точек Каждый мультипроцессор (MP) GPU обрабатывает отдельный блок 32*16*Zn точек сетки последовательно по слоям, где Zn – размер области вдоль оси Z. 32х16=512 потоков При последовательных послойных расчетах автоматически активно используются регистры и кэш-память GPU.

• •

В быстрой памяти GPU хранятся данные 2 шагов по времени -> высокая производительность.

Во внешней памяти - данные для границ. Ввиду малого объема, время обмена по медленному каналу с внешней памятью составляет около 2% от времени расчетов.

Полный диапазон данных в 3D

•

Впервые разработаны алгоритмы и создан комплекс масштабируемых программ для суперкомпьютеров на GPU в 3D задачах УЗ томографии в медицине.

•

Проведены модельные расчеты на сетках до 500х500х500 для прямой и обратной задачи с полным и неполным диапазоном углов данных, разрешение порядка 2-3 мм.

Фантом скорости Схема эксперимента Восстановленная скорость

Неполный диапазон данных в 3D

Решение задачи на суперкомпьютере «Ломоносов» 17 графических карт NVIDIA Tesla X2070.

    

Источников -17 Время ~ 2 часов Длина импульса - 5 мм. Сетка 420х420х420. Область 10х10х10 см Схема эксперимента 1 Сечение скорости Восстановленная скорость Схема эксперимента 2

Модельные расчеты на GPU

5 итераций; 15 итераций; 120 итераций Результаты модельных расчетов

10  (u(c (n) )) 10 • • •

Невязка за 120 итераций уменьшилась в 30 раз.

Невязки 2.3*10 -4 соответствует уровню ошибки ~ 0.23%.

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

Число итераций Функционал невязки 6 источников Время расчета (мин) 5 6.3*10 -3 5 15 1.6*10 -3 15

10 10 0 20 40 60 80 100 120 n

График невязки от числа итераций 120 2.3*10 -4 120

•

Задача волновой томографии с учетом поглощения

В задачах ультразвуковой томографии на частотах > 0.5 Мгц поглощение более 50%. Частота 20 Мгц – только для подкожного слоя.

Разработан новый подход, учитывающий поглощение в среде.

•

c u



Основная задача Модель1 (поглощение не зависит от частоты)

(

n

(

r r u

)

u tt

,

t

|

ST

(

r

  ,

t

0 ) )

p

(  

a u

(

t r

( )

u t r

, (

t r

, 

r

,

t

)

t

)  0 ) 

u

(  0

r

,

t

)   •

Обратная задача – минимизация функционала невязки

 (

u

(

c

) ,

c

) 

u

|

ST



U

2 (

r



q

) 

f

(

t

) •

Сопряженная задача

c

(

r

)

w

 (

n w r w tt

,

t

|

ST

(

r

, 

t

) 

T u

) 

a(r)w

 |

ST w t



U

(

t r

(

r

, ,

t t

)  

T



w

( ) 

r

, 0

t

)  0 •

Градиент



C

(

u

(

c

) ,

c

функционала имеет вид:

T

)  

w t

)

u t

(

r

)

dt

где

u

(

r

,

t

)

w

(

r

, 0

t

) 

a

(

u

(

c

,

a

))  

T w t

(

r

,

t

)

u

(

r

,

t

)

dt

Основная задача Модель2 (поглощение квадратично зависит от частоты стоксовское поглощение в вязких средах)

c

(

r

)

u tt

(

r

,

t

) 

a

(

r

)( 

u

)

t

 

u

(

r

,

t

)   (

r



q

) 

f

(

t

)

Обратные задачи УЗ томографии с учетом поглощения. Модельные расчеты

• •

Удалось в модельных расчетах реконструировать одновременно 2 функции (скорости и поглощения). Скорость восстанавливается лучше.

Фантом скорости 500 итераций 3500 итераций Фантом поглощения 64 источника 64 источника Точное и восстановленное сечения Восстановленное поглощение

•

Обратные задачи УЗ томографии с учетом поглощения. Модельные расчеты

• •

Восстановить скорость можно без учета поглощения. Восстановить поглощение при неправильной модели – проблематично.

Фантом скорости Фантом поглощения Решение на основе Модели 1 по данным из Модели 2 4 источника без учета поглощения

Суперкомпьютерное моделирование

•

Проблема: Насколько влияет плотность на реконструкцию в задачах УЗ томографии в медицине.

•

Из уравнения движения, непрерывности и состояния можно получить



u

(

r

,

q

,  ) 

c

2  2 (

r

)

u

=

f

(

r

,

q

,  )     (

r

(

r

) ) 

u

(

r

,

q

,  ).

•

Экспериментальные данные получены из точных аналитических формул в 2D для разных плотностей

•

В мягких тканях человека плотность варьирует ~ 20%.

•

Импедансы равны – отражения нет, хуже разрешение

•

Если акустический импеданс противоположен, то появляется двоение границы.

•

Результаты моделирования обнадеживают.

 1   2  1 

График. Вид отраженной и преломленной волн для разных плотностей

 2  1   2

Модель Результаты реконструкции

Суперкомпьютерное моделирование

•

Влияние шага между детекторами .

λ/2 1.5λ 3λ Результаты реконструкции 4 источника

•

Влияние количества источников .

Результаты реконструкции.

Шаг между детекторами 3λ

•

Вывод: Уменьшение числа приемников можно компенсировать увеличением числа источников 8 источников 32 источника

Суперкомпьютерное моделирование

•

Влияние длины импульса .

Результаты реконструкции 4 источника λ=5мм λ=10мм

•

Количество уровней оцифровки сигнала 1024 уровня 128 уровней Результаты реконструкции 4 источника 32 уровня

Основные результаты и выводы

•

Предложены эффективные методы и алгоритмы решения 3D задач волновой томографии как коэффициентных обратных задач через прямое вычисление градиента функционала невязки. Методы предполагают как послойное (двумерное) восстановление 3D объекта, так и решение непосредственно трехмерного уравнения.

•

Решение непосредственно 3D уравнения позволяет повысить качество реконструкции и разрешение. Удается восстанавливать не только форму неоднородности, но и абсолютные значения скоростей.

•

3D задача приводит к сверхбольшому объему данных в нелинейной обратной задаче с числом неизвестных порядка ста миллионов.

•

Разработанные алгоритмы хорошо масштабируются на процессорах общего назначения и на GPU. Предложена структура программы, позволяющая выделить большое количество независимых параллельных процессов в общем потоке задачи.

•

В моделях с учетом поглощения удалось в модельных расчетах реконструировать одновременно 2 функции (скорости и поглощения). Скорость восстанавливается лучше.

•

Работа выполнена при поддержке РФФИ проекты № 12-07-00304-А и №. 13-07-00824 А