Рандомизированный алгоритм при обработке данных ультразвуковых исследований

Сенин Иван Игоревич

Рандомизированный алгоритм при обработке данных ультразвуковых исследований

Ультразвуковая томография нашла широкое применение в медицинской практике. По мере развития технологий стало возможным использование большего количества датчиков для получения более качественного изображения, что привело к большому объему обрабатываемой информации. Однако, характер получаемых при томографии изображений имеет разреженный профиль, из чего следует, что необходимое для восстановления количество информации очень мало. Значительная избыточность данных оставляет возможность для разработки более оптимальной технологии их сбора и анализа в процессе томографического исследования. В этой работе рассмотрена возможность применения техники ``Compressive Sensing'' для улучшения алгоритма реконструкции. Метод показал свою эффективность, позволив снизить объем используемых данных без потери в качестве.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36495f1be77c40d58b85

UUID: 9ccfe61a-e943-46ca-8332-187321b3afbd

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 6

Сенин Иван Игоревич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 353484 bytes

Поделиться работой

Санкт-Петербургский государственный университет Математическое обеспечение и администрирование информационных систем Системное программирование Сенин Иван Игоревич Рандомизированный алгоритм при обработке данных ультразвуковых исследований Бакалаврская работа Научный руководитель: д. ф.-м. н., профессор Граничин О. Н. Рецензент: разработчик ПО ООО “Рэйдикс” Тюшев К. И. Санкт-Петербург 2016

SAINT-PETERSBURG STATE UNIVERSITY Software and Administration of Information Systems Software Engineering Ivan Senin Randomised algorithm in ultrasonic investigation precessing Bachelor’s Thesis Scientific supervisor: professor Oleg Granichin Reviewer: Software Developer at RAIDIX Kirill Tyushev Saint-Petersburg 2016

Оглавление Аннотация 4 1. Введение 1.1. Актуальность . . . . . . . . . . . . . . . . . . . . . . . . . 1.2. Обзор литературы и ключевых работ в предметной области 1.2.1. Travel-time томография . . . . . . . . . . . . . . . . 1.2.2. Compressive Sensing . . . . . . . . . . . . . . . . . . 1.3. Предварительный пример . . . . . . . . . . . . . . . . . . 1.4. Краткий обзор результатов работы . . . . . . . . . . . . . 1.5. Структура работы . . . . . . . . . . . . . . . . . . . . . . . 5 5 6 6 9 13 14 15 2. Постановка задачи 16 3. Основной результат 3.1. Исследование разреженности изображения . 3.2. Ожидания по масштабированию задачи . . . 3.3. Проектирование универсальной матрицы измерений . . . . . . . . . . . . . . . . . . . . 3.4. Алгоритм реконструкции при использовании Compressive Sensing . . . . . . . . . . . . . . . 3.5. Архитектура программного решения . . . . . . . . . . . . . . . . . . . 17 17 19 . . . . . . . 22 . . . . . . . . . . . . . . 24 24 4. Моделирование 4.1. Результаты эксперимента на маломасштабной модели . . 26 27 5. Заключение 30 Список литературы 31 3

Аннотация Ультразвуковая томография нашла широкое применение в медицинской практике. По мере развития технологий стало возможным использование большего количества датчиков для получения более качественного изображения, что привело к большому объему обрабатываемой информации. Однако, характер получаемых при томографии изображений имеет разреженный профиль, из чего следует, что необходимое для восстановления количество информации очень мало. Значительная избыточность данных оставляет возможность для разработки более оптимальной технологии их сбора и анализа в процессе томографического исследования. В этой работе рассмотрена возможность применения техники “Compressive Sensing” для улучшения алгоритма реконструкции. Метод показал свою эффективность, позволив снизить объем используемых данных без потери в качестве. 4

1. Введение 1.1. Актуальность Ультразвуковая томография по качеству и разрешающей способности достигла сопоставимого с МРТ уровня и активно применяется для исследования мягких тканей [7]. Преимуществами УЗИ являются относительно низкая стоимость оборудования и обслуживания, безопасность для организма, неинвазивность техники исследования. Для повышения качества и разрешающей способности получаемого в ходе исследования изображения требуется увеличение количества используемых датчиков и частоты дискретизации сигнала. Всё это ведёт к значительному увеличению объема передаваемых данных, что усложняет как производственный процесс, так и проведение диагностики. Кроме того особенности проведения исследования допускают наличие помех и искажений [28]. Опухоли преимущественно имеют более высокую скорость прохождения ультразвука, чем окружающие ткани. Это делает возможным реконструкцию плотностей тканей в исследуемой зоне с помощью уравнений с участием предполагаемых путей распространения сигнала и временем его прибытия на датчики кругового массива (travel time tomography) [25]. Современные системы основаны именно на таком подходе и получили применение в диагностике рака молочных желез [26][30]. Для метода томографии travel-time типична квадратичная от числа сенсоров зависимость получаемых “сырых” данных для анализа: последовательно с каждого датчика пускается ультразвуковой импульс, который принимают остальные k −1 сенсоров. Это приводит к серьезному повышению требований к вычислительной части устройства томографа, а также к увеличению времени обработки: современные методы реконструкции работают итеративно с временной сложностью итера5

ции O(N log N ) [11]. На настоящий момент использование ультразвуковой томографии осложнено высокими требованиями к вычислительному комплексу и временем на обработку [24]. Уменьшение объема данных позволит эффективно использовать FPGA вычислители для задач определения времени прибытия сигнала на датчик и последующей реконструкции изображения. Compressive Sensing также позволит сократить время на сбор данных путем уменьшения количества производимых проекций, что частично сократит зашумление данных из-за колебаний пациента в процессе обследования. 1.2. Обзор литературы и ключевых работ в предметной области 1.2.1. Travel-time томография Техника томографии по времени прибытия сигнала уже хорошо изучена и освещена в работах [20], [25], [7]. В том числе рассмотрены алгоритмы, устойчивые к присутствию шума (помех) [26]. Далеко не последнюю роль в реконструкции изображения играет точное определение времени прибытия сигнала [4]. Общий процесс обработки данных томографии по времени прибытия показан на рис.1. Основная цель акустической томографии – восстановить параметры неизвестной среды, изучая характеристики распространения звука в ней. Во-первых, для этого требуется точная модель, хорошо описывающая лежащую в её основе физическую систему, и, во-вторых, высокоточные измерения. Тогда решением обратной задачи составляется оценка неизвестной модели. Корректность физической модели, точность измерений и выбор метода решения обратной задачи имеют прямое влияние на качество реконструкции. Распространение энергии акустического сигнала в неоднородной сре6

Рис. 1: Общий вид процесса ультразвуковой томографии времени прибытия. де хорошо описывается дифференциальным уравнением в частных производных второго порядка ∇2 p(r, t) − 1 ∂ 2 p(r, t) = s(r, t), F 2 (r) ∂t2 где p(r, t)  давление в r в момент времени t, F (r)  неизвестная модель распространения скорости и s(r, t)  первоначальный сигнал. Зная исходный сигнал s(r, t), решая обратную задачу, находим модель F (r), которая лучше всего описывает измерения p(r, t)|Ω , записанные на известных границах Ω. Моделирование прямой задачи по уравнению распространения волны вычислительно весьма трудоемко, вследствие чего в travel-time томографии используются принципы геометрической акустики. При пред7

положении, что частота звукового сигнала достаточно высока, можно найти его путь распространения используя принципы Ферма [27]. Тогда обратная задача состоит в реконструкции распределения скорости звука F (r) исходя из времени сигнала в пути, получаемое с показаний датчиков. Следует отметить, что в отличии от задач томографии с помощью рентгена, где сигнал распространяется по прямой, распространение ультразвукового сигнала в неоднородной среде происходит иначе и зависит от распределения скорости. В то время как акустическая томография по времени прибытия сигнала берет свои корни из сейсмологии [16], на сегодняшний день уже множество исследований показало, что ультразвуковая томография имеет большой потенциал в области диагностирования рака груди [15]. Постановка обратной задачи по реконструкции изображения Время сигнала в пути от передатчика к приемнику вдоль траектории его распространения можно представить как ∫ 1 Y = ds, (1) Γ F (r) где Y  время прибытия сигнала, Γ  путь его распространения, F (r)  скорость звука в точке r. Следует заметить, что проходимый сигналом путь Γ зависит от распределения скорости в среде F (r). Существует нелинейная зависимость между временем сигнала в пути и скоростью звука в среде. Можно представить уравнение (1) в дискретной форме, наложив сетку на интересующую для восстановления область с константными значениями скорости звука в клетках общим количеством N ячеек: Y = A(F ) · F, (2) где F  вектор [N × 1], представляющий интересующее распределение скоростей, а A(F )  матрица [M ×N ], представляющая путь, по которо8

му проходит сигнал и t  вектор [M × 1], содержащий время прибытия сигнала, полученное анализом показаний с датчиков. Так, имея в установке k датчиков, получаем M = k · (k − 1) всевозможных траекторий кратчайшего прохождения сигнала. Целью обратной задачи является нахождения такой оценки распределения скоростей m̂ , которая лучше всего описывает время прохождения пути сигналом в уравнении (2). Одним из традиционных алгоритмов решения является метод LASSO с нормой ℓ1 и вариацией для регуляризации [26]. Тогда задача представляется как задача минимизации min ∥A(F ) · F − Y ∥22 + λ∥ΨT F ∥1 + λT V T V (F ), (3) F где Ψ  некоторый базис, переводящий F в разреженное представление и λ, λT V  положительные весовые коэффициенты. Этот метод будем в дальнейшем называть исходным. 1.2.2. Compressive Sensing С начала третьего тысячелетия объемы получаемой и обрабатываемой информации существенно возросли. Главным образом это связано с увеличением размерностей передаваемых сигналов (от одномерных “1-D” к дву- и трех- мерным). При этом рост данных происходит по экспоненциальному закону относительно размерности. Также в современных системах происходит постоянный рост требований к частоте измерений (теорема Шеннона-Найквиста-Котельникова). Все это влечет необходимость сжатия данных для хранения и пересылки, стоимость суммарной обработки оказывается очень дорогостоящей [2]. Интересующую информацию можно представить как x ∈ X, которую предстоит получить с помощью набора наблюдений y ∈ Y. Между x и y существует закономерность явления Φ : X → Y. Если Φ линейный обратимый оператор, то известно, что x = Φ−1 y. Для систем реального мира более типичным является случай, когда результаты наблюдений подвержены 9

различным помехам y = Φx + ξ. При значительных помехах задача решается статистически, используя m >> N при x ∈ RN . Известно, что такая задача о восстановлении x может быть решена даже в случае нецентрированных коррелированных помех за счет случайного выбора матрицы Φ [19]. Установлено, что рандомизация позволяет не только устранить эффект смещения, но и уменьшить количество итераций алгоритма оценивания x [1]. На практике интересно рассмотреть возможности восстановления x ∈ RN при m << N , что, очевидно, невозможно в общем случае. Однако, оказалось, что задача может быть решена с достаточной точностью, при выполнении некоторых условий [17]. Подобная парадигма получения и восстановления данных называется Compressive Sensing или Опознание со сжатием. Постановка задачи о получении и восстановлении данных в общем виде выглядит следующим образом. Пусть x  интересующая информация. Она проявляется себя через сигнал f = Ψx. С помощью каких-либо приборов регистрации получают наблюдения y = Af . По ним исследователю необходимо восстановить исходную информацию x, формируя оценки x̂. Проектирование матрицы измерений Получение m << N наблюдений можно представить как скалярное произведение с множеством некоторых векторов {ai }m i=1 , или умножением на матрицу размерности m × N : y = Af . Тогда из постановки задачи следует y = Af = AΨx = Φx, где Φ  матрица m × N . Необходимыми и достаточными условиями, предъявляемыми теорией опознания со сжатием, являются: 10

• сигнал f является s−разреженным: ∑ существует такой базис Ψ, что f = N j=1 x[j]ψj , где только s коэффициентов x[j] ̸= 0, • выполняется “свойство ограниченной изометрии” (Restricted Isometry Property, RIP): RIP(δ, m) = √ 1−δ ≤ ∥Φz∥2 √ ≤ 1 + δ, ∥z∥2 где δ ∈ (0, 1) и z  произвольный ненулевой m-разреженный вектор, • свойство “некогерентности” (малой взаимной зависимости) A и Ψ: µ(A, Ψ) = √ N max i,j |⟨ai , ψj ⟩| . ∥ai ∥2 Установлено, что случайные матрицы A с высокой вероятностью некогерентны с любым фиксированным базисом Ψ. Из [9] следует достаточное условие для высокой вероятности точного восстановления s−разреженных векторов по m наблюдениям: m ≥ cµ(A, Ψ)2 s log N, (4) где c  некоторая положительная постоянная. В замечании из [8] обращают внимание, что на практике часто достаточно m ≈ 4s. Также известно, что условие RIP(δ, 2s) является недостаточным для восстановления в присутствии помех или если в сжимаемом сигнале N − s компоненты малы, но не равно нулю. Для робастности же достаточным условием будет RIP(δ, 3s). Явное построение матрицы измерений A: Φ = AΨ требует проверки условия RIP для каждой комбинации положений s ненулевых компонент вектора z длиной N (т. е. CNs ). Однако, оказывается, что случайная матрица m × N с элементами a[i, j] ∼ N (0, m1 ) имеет следующие полезные свойства [10]: 11

• если m ≥ c1 s log Ns при 0 < δ < 1, тогда P{A ∈ RIP(δ, m)} ≥ 1 − 2e−c2 m , где P  вероятность, c1 , c2 > 0  малые постоянные, зависящие от δ. И тогда, следовательно, s−разреженные сигналы длины N могут быть восстановлены по m << N измерениям с высокой вероятностью. • если матрица A, удовлетворяет указанным условиям, то и Φ = AΨ также будет являться матрицей с независимо нормально-распределенными элементами, что влечет Φ ∈ RIP(δ, m) с той же вероятностью и для любого ортонормированного базиса Ψ. Для использования опознания со сжатием помимо проектирования матрицы измерений A требуется разработать алгоритм реконструкции сигнала, задачей которого восстановить по y ∈ Rm , матрице измерений A и базису Ψ сигнал f ∈ RN , или эквивалентное ему спектральное разреженное представление x. При m < N для s-разреженных сигналов имеется бесконечное число x′ : Φx′ = y, т. к. если Φx = y =⇒ Φ(X + r) = y для любого r ∈ Ker Φ. Таким образом, алгоритм реконструкции должен найти вектор разреженного представления сигнала в (N − m)-размерном подпространстве H = Ker Φ + x. Методы решения обратной задачи Традиционным подходом в решении подобных обратных задач является поиск в H вектора с минимальной энергией (нормой) ℓ2 : x̂ = argmin ∥x′ ∥2 : Φx′ = y, 12

или, что тоже самое, в форме метода наименьших квадратов x̂ = ΦT (ΦΦt )−1 y. Поскольку ℓ2 норма отражает энергию сигнала, а не разреженность, результат будет содержать много ненулевых компонент. С задачей поиска разреженного решения хорошо справляется “ℓ0 −норма”. Однако, задача ℓ0 оптимизации невыпукла и относится к комбинаторному типу, процедуры ее решения, в общем случае, N P -сложные [23]. Существует два подхода для приближенного решения таких задач: “жадные” стратегии оптимизации [22] или переход к ℓ1 -норме. Хорошо известно, что минимизация ℓ1 нормы способствует разреженности [18] и позволяет хорошо приближать разреженные сигналы используя только m ≥ c1 s log Ns случайных измерений [17][10]. Это задача выпуклой оптимизации, сводимая к известной задаче линейного программирования “выбор базиса” (Basis Pursuit Denoising) [12], эффективные методы решения которого имеют временную сложность O(N log N /s) [5]. Техника обработки сигналов Compressive Sensing может быть использована для решения проблемы передачи и обработки большого количества информации с массива датчиков. Этот метод уже успешно применяется в магнитно-резонансной томографии [21, 13], а также в ультразвуковых исследованиях [14]. Однако, данные, используемые для получения изображения в ультразвуковой томографии имеют иную природу, чем в классическом (B-mode) ультразвуковом исследовании, так как после замера сигнала из него извлекается время прибытия звуковой волны, с которой и происходит вся дальнейшая работа по реконструкции снимка. Также этот факт отличает ультразвуковую томографию от МРТ, где все необходимые для реконструкции данные получают уже в разреженном представлении. 13

1.3. Предварительный пример В целях получения наилучшего изображения требуется большое количество датчиков и высокая частота дискретизации сигнала. Всё это ведёт к необходимости обработки серьезного объема данных. Оценить количество данных, поступающих непосредственно с датчиков можно: V = k 2 zYmax ν, где V  общее число измерений сигнала, k  количество датчиков, z  количество исследуемых срезов, Ymax  время прохождения сигнала с учетом затухания эха, ν  частота дискретизации. Примером современного коммерческого ультразвукового томографа является The SoftVue [6], который имеет следующие технические характеристики системы: • Мастер сервер: 2 процессора quad-core Intel Xeon E5620, 192 ГБ ОЗУ • Сервер реконструкции: 2 процессора quad-core Intel Xeon E5620, 96 ГБ ОЗУ, 2 ГП Nvidia Tesla M2070 Частота дискретизации (МГц) 10 12 14 Число датчиков 256 512 1024 0.21 0.86 3.44 0.26 1.03 4.13 0.30 1.20 4.81 Таблица 1: Объем исходных данных ГБ за один срез в зависимости от числа датчиков в массиве и частоты дискретизации сигнала [6]. Так в таблице 1 производители описали количество данных, получаемых с датчиков за один срез. На обработку такого среза при использовании только одного сервера для реконструкции требуется несколько минут, а при использовании всего вычислительного комплекса – 20 секунд. Всего производится 70 срезов [6]. Таким образом время, за14

трачиваемое только лишь на вычисление томографических снимков, составляет приблизительно 23 минуты. 1.4. Краткий обзор результатов работы В ходе выполнения работы был разработан метод, позволяющий при разработке ультразвукового томографа внедрить в применяемые алгоритмы по сбору и обработке данных парадигмы опознания со сжатием. Эксперименты, проведенные на компьютерных моделях, показали состоятельность метода. В частности, даже в маломасштабной модели оказалось возможным сокращение используемых данных до 68% от общего их объема. 1.5. Структура работы Сначала в разделе 3 проверяется и демонстрируется факт разреженности восстанавливаемого сигнала, для чего было найдено подходящее трансформирующее кодирование. После чего рассматриваются различные варианты построения рандомизированной матрицы измерений, а также процесс реконструкции и архитектура программного решения. Далее приведены оценки необходимого для реконструкции объема данных при дальнейшем увеличении числа датчиков и разрешения снимка. Раздел 4 описывает проведенные эксперименты, где также анализируются полученные результаты. Наконец, подведение итогов проделанной работы происходит в разделе 5. 15

2. Постановка задачи Цели работы: • Исследовать возможность применения техники Compressive Sensing в задаче ультразвуковой томографии. • Разработать рандомизированный алгоритм для сбора и обработки данных ультразвуковой томографии. • Провести эксперименты на программной симуляции. Для достижения поставленных целей необходимо: 1. Определить характеристики разреженности сигнала s, а также их изменения при дальнейшем увеличении разрешения получаемого изображения. 2. Для сбора измерений и последующего восстановления изображения следует спроектировать рандомизированную матрицу измерений. Согласно теории Compresive Sensing допустимы различные варианты ее построения. 3. Разработать модификацию алгоритма реконструкции томографического снимка. 4. На основании результатов программной симуляции выяснить эффективность полученного алгоритма. 16

Рис. 2: Карта скоростей, восстанов- Рис. 3: Представление карты скоро- ленная из сигнала. стей после вейвлет-преобразования. 3. Основной результат 3.1. Исследование разреженности изображения Для сравнения результатов алгоритмов реконструкции с исходной моделью в качестве метрики будем использовать среднеквадратичную ошибку (RMSE). Эта метрика уже была использована в работах по Compressive Sensing. Анализ характера разреженности сигнала Известно, что томографические снимки имеют разреженный характер. К категории эффективных трансформирующих преобразований, позволяющих привести восстановленную карту скоростей F к разреженному представлению, относятся вейвлет-преобразования. Исследования Compressive Sensing в сфере МРТ показали эффективность вейвлета Добеши db6 [21]. Пример использования вейвлета db6, примененного к результату реконструкции искусственной маломасштабной модели и результат преобразования X = Ψ−1 F показаны на рис. 2 и рис. 3 соответственно. 17

Ошибка (RMSE), м/с 30 20 10 0 0 100 200 300 400 500 600 Число использованных коэффициентов результата вейвлет-трансформации Рис. 4: Ошибка восстановления карты скоростей искусственной модели при использовании неполного набора коэффициентов вейвлетпреобразования. Амплитуда 400 300 200 100 0 0 10 20 30 40 50 60 70 Модуль коэффициентов вейвлет преобразования в отсортированном порядке Рис. 5: Набор коэффициентов вейвлет-преобразования. Анализ разреженности искусственной маломасштабной модели Согласно маломасштабной модели, общее количество коэффициентов составляет N = 4096. Легко убедиться в разреженности сигнала: рис. 4 демонстрирует, что при использовании только лишь 100-200 коэффициентов происходит восстановление близкое к полному. В контексте опознания со сжатием это означает, что характеристику сигнала s для маломасштабной модели можно выбрать на интервале 18

от 100 до 200 без существенной потери качества. Опираясь на отсортированный по амплитуде набор коэффициентов вейвлет преобразования (рис. 5) в дальнейшем для искусственной модели будем полагать s = 112 = 121, что составляет около 3% от всех коэффициентов N . Анализ разреженности томографического снимка Для того, чтобы оценить влияние сложности деталей, содержащихся на изображении, рассмотрим свойства снимка, полученного при МРТ. Результаты анализа его вейвлет-преобразования представлены на рис. 6 и 7. Таким образом для томографического снимка можем взять характеристику разреженности smri = 252 = 625, что составляет около 3.8% от размера сигнала N = 128 × 128 = 16384. Результат восстановления при использовании smri коэффициентов показан на рис. 8. Ошибка (RMSE), м/с 100 80 60 40 20 0 400 800 1,200 1,600 Число использованных коэффициентов результата вейвлет-трансформации Рис. 6: Ошибка восстановления томографического снимка при использовании неполного набора коэффициентов вейвлет-преобразования. 3.2. Ожидания по масштабированию задачи Характер разреженного сигнала, выражающийся параметром s, по большей мере зависит от сложности изображения, чем от его разрешения. Так при неизменном параметре s было проведено масштабирова19

3,000 Амплитуда 2,000 1,000 0 −1,000 −2,000 0 20 40 60 80 100 Порядковый номер коэффициента 120 140 Рис. 7: Набор коэффициентов вейвлет-преобразования. Рис. 8: Результат восстановления МРТ снимка при использовании smri = 252 коэффициентов вейвлет-преобразования, что составляет 3.8% от их общего объема. ние искусственной модели и вычислены значения среднеквадратичного отклонения (см. табл 2). Результаты показывают, что при использовании только 225 из 20482 коэффициентов вейвлет-преобразованного изображения среднеквадратичная ошибка составляет 0.16% от амплитуды скоростей в изучаемой модели. Это дает повод полагать, что при увеличении требуемого разрешения изображения томографической реконструкции характеристика разреженности s должна увеличится незначительно, и этим фактом можно пренебречь в дальнейших теоретических оценках. 20

N 642 1282 2562 5122 10242 20482 RMSE, м/с 0.0458 0.3131 0.8941 1.4651 1.6269 1.7754 Таблица 2: Значение среднеквадратичного отклонения от длины сигнала N (числа пикселей) восстановленного изображения на искусственной модели (описанной в разделе 4) при характеристике разреженности s = 152 . В качестве трансформирующего кодирования был использован вейвлет db6. Можем полагать, что разрешение изображения N ≈ k 2 , где k  количество датчиков в устройстве. Тогда число уравнений, получаемых в результате максимально возможного числа проекций томографического обследования составляет M = k · (k − 1). Важным результатом теории опознания со сжатием является тот факт, что для восстановления разреженного изображения достаточным является использование только m измерений: m ≈ 4s log N . s (5) Так по рис. 9 и табл. 3 видно преимущество, получаемое при использовании парадигмы опознания со сжатием: количество уравнений, необходимых для реконструкции требуемого разрешения, растет логарифмически по сравнению с традиционным способом использования полного набора M уравнений. M 64 × 64 128 × 128 256 × 256 512 × 512 1024 × 1024 m 1705 2376 3047 3718 4389 Таблица 3: Количество измерений достаточных для реконструкции (m) и полное количество произведенных измерений M . 21

m/M × 100% 40 30 20 10 0 0 500 1,000 1,500 2,000 Количество датчиков в устройстве Рис. 9: Процент уравнений, достаточных для реконструкции изображения, от общего объема получаемых данных. 3.3. Проектирование универсальной матрицы измерений Согласно [3] хороший результат дает матрица со следующим распределением элементов:  √   с вероятностью 16 + m3   Ai,j = 0 с вероятностью 23  √   − 3 с вероятностью 1 m 6 Такая матрица с высокой вероятностью будет удовлетворять требованиям, описанным в разделе 1.2.2. Вариант 1: базовый вариант матрицы Эта же матрица будет использована для тестирования базового метода опознания со сжатием. Это позволит сократить объем данных для последующих вычислений, однако требуется полный объем сведений о времени прибытий сигналов. Таким образом время на предваритель22

ную обработку (извлечение из измерений датчиков времени прибытия звуковой волны) не будет уменьшено. Вариант 2: равновероятный отсев данных В отличие от варианта 1, (M − m) случайных столбцов матрицы будут заполнены нулями. Это далает возможным сократить также время на предварительную подготовку данных (извлечение времени прибытия звуковой волны), так как часть датчиков может не производить измерения. Вариант 3: отсев проекций Теперь, в отличие от варианта 2, зануляться будут группы столбцов  датчики, принявшие сигнал от какого-либо датчика, пославшего импульс. Уменьшение числа совершаемых проекций позволяет не только уменьшить объем данных для обработки, но и сократить время на непосредственно исследование пациента, что ведёт к уменьшению зашумления данных от движений тела пациента. Вариант 4: точечный выбор данных В каждой строке матрицы находится только один ненулевой коэффициент. Все ненулевые коэффициенты равны. Заметим: в то время как вариант 1 опирается на полный объем данных, варианты 2,4 используют только m из M измерений. Вариант 3 округляет m вверх до ближайшего значения, позволяющего использовать полное число (kRx ) принимающих датчиков. 23

3.4. Алгоритм реконструкции при использовании Compressive Sensing Реконструкция будет производиться методом сопряженных градиентов, минимизируя: min ∥ÂA(F )F − ÂY ∥22 + λ∥ΨT F ∥1 . F (6) Здесь Â является одним из вариантов рандомизированных матриц измерений, в качестве Ψ используется вейвлет db6, остальные переменные уже описаны в разделе 1.2.1. В процессе реконструкции перемножение Â · Y выполняется единожды: в начале вычислений. Также произведение Â · A(F ) можно произвести неявно: производя вычисления матрицы путей A(F ) только для заведомо используемых данных. Исключением будет являться вариант-1 построения матрицы измерений Â, так как он использует полный объем данных. 3.5. Архитектура программного решения В ходе работы было реализовано программное решение, позволяющее адаптировать алгоритм реконструкции к ультразвуковым томографам с различными характеристиками. Общая структура программы показана на диаграмме 10. Рис. 10: Архитектура программы, оптимизирующую реконструкцию томографического снимка. 24

Исходными данными для программы оптимизации алгоритма реконструкция являются сущности: TomographProperties, SignalData и TissueModel. TomographProperties содержит такие характеристики томографа, как количество датчиков, диаметр кольца, на котором расположен массив датчиков. TissueModel представляет собой модель опухоли, описывающей параметры скорости распросранения звука в среде и затухания сигнала. По модели проводится программная симуляция, результаты которой в виде записанных показаний датчиков располагаются в объекте SignalData. Предварительная обработка данных заключается в извлечении из показаний датчиков времени прибытия звуковой волны. Эта нетривиальная задача решается методом, предложенным в [4] сущностью TraveltimePicker. По исходной модели также необходимо определить характеристику разреженности s. Это происходит в SparsityEstimation посредством анализа результата вейвлет-трансформации модели с использованием модификации алгоритма First Break Picking (“STA/LTA” метод). ImageReconstruction содержит реализацию алгоритма получения снимка, описанном в разделе 3.4. Для определения наилучших параметров алгоритма реконструкции используется метод стохастической оптимизации “имитация отжига”, сущность SimulatedAnnealing. В качестве функции стоимости используется среднеквадратичное отклонения от заданной модели. За вычисление функции стоимости отвечает QualityAssessment. Наконец, ReconstructionSetting содержит параметры, необходимые для алгоритма реконструкции, такие как весовые коэффициенты, величину достаточной сходимости, максимальное число итераций, разрешение томографического снимка. 25

4. Моделирование Описание маломасштабной модели Для исследования по внедрению парадигмы опознания со сжатием в сфере ультразвуковой томографии была взята следующая модель: • Исследуемая область заполнена водой. В центр помещена опухоль диаметром 5 мм. • Скорость прохождения звука в воде fw = 1500 м/с Скорость прохождения звука через опухоль: ft = 2600 м/с. • Диаметр кольца массива датчиков d = 4 см. • Полное число датчиков массива, равномерно распределенных по периметру кольца: k = 100. • Число датчиков, используемых для получения сигнала, kRx = 41 ·k. • Сетка реконструкции N = 64 × 64. Модель представлена на рис. 11. Рис. 11: Схема маломасштабной модели. Для симуляции модели был использован специальный инструмент для ультразвуковых исследований в среде Matlab “K-Wave”. K-Wave неоднократно использовался в научных работах1 [29]. 1 http://www.k-wave.org/publications.php 26

Матрица измерений SNR, dB RMSE, м/с Базовая матрица 23.89 ± 0.06 181.89 ± 1.30 Равновероятный отсев 23.84 ± 0.06 183.06 ± 1.34 Отсев проекций 23.77 ± 0.07 184.59 ± 1.48 Точечный выбор данных 23.80 ± 0.04 184.15 ± 0.85 Исходный алгоритм 23.63 183.36 Таблица 4: Качество реконструкции снимка Compressive Sensing алгоритмами при s = 112 и исходным алгоритмом. 4.1. Результаты эксперимента на маломасштабной модели Алгоритм реконструкции с использованием техники опознания со сжатием был протестирован на четырех вариантах построения рандомизированной матрицы измерений на 50 итерациях каждый. Результаты представлены в таблице 4, лучшие показатели выделены жирным шрифтом. Среди протестированных методов построения матрицы измерений лучшего результата достиг базовый вариант. Это можно объяснить тем фактом, что при построении он исходит из полного объема данных, комбинируя которые можно уменьшить влияние помех. Помехи являются неотъемлемой составляющей ультразвуковой томографии хотя бы потому, что пациент во время обследования не может быть абсолютно неподвижен. Все варианты алгоритмов с опознанием со сжатием показали лучшее качество реконструкции по метрике “отношение сигнал-шум” (SNR), однако по абсолютной среднеквадратичной ошибке варианты с исключением целых проекций и точечного выбора данных оказались хуже исходного алгоритма. “Отсев проекций” отсекает целые группы данных, что сказалось на качестве восстановления. Варианты 1–3 комбинируют различные измерения, что приводит к их большей шумоусойчивости, чего не делает вариант 4 “Точечный выбор данных”, продемонстрировавший худшие результаты. В маломасштабной модели для реконструкции можно использовать до M = kRx · k 2 = 2500 уравнений. Повторив анализ, описанный в разделе 3.1, мы находим s = 112 . Таким образом, применяя уравнение (5), по27

Рис. 12: Результаты реконструкции исходного алгоритма (слева) и Compressive Sensing с матрицей измерений i.i.d. (справа). лучаем достаточное для хорошей реконструкции количество уравнений m ≈ 1705, что составляет около 68% от общего их числа M . На графике 13 отражена зависимость качества реконструкции изображения от количества исходных данных, где наилучшее значение показателя абсолютной среднеквадратичной ошибки достигается при использовании 68-75 % данных и в дальнейшем не претерпевает существенных изменений. 220 Ошибка (RMSE), м/с CS-i.i.d. Legacy 210 200 190 180 10 20 30 40 50 60 70 80 Процент использованных данных 90 100 Рис. 13: Зависимость ошибки реконструкции снимка от объема использованных данных при использовании Compressive Sensing с матрицей измерений с равновероятным отсевом (CS-i.i.d.). Для сравнения показан результат исходного алгоритма, использующий полный объем данных (Legacy). Результаты реконструкций исходного алгоритма, используемом в 28

ультразвуковой томографии, и Compressive Sensing представлены на рис. 12. 29

5. Заключение В ходе работы были достигнуты следующие результаты: • Исследована возможность применения техники Compressive Sensing в задаче ультразвуковой томографии. • Разработан и реализован на языке MATLAB рандомизированный алгоритм для сбора и обработки данных ультразвуковой томографии. • Проведены эксперименты на программной симуляции. В этой работе была исследована возможность применения техники Опознания со сжатием, или “Compressive Sensing”, в области ультразвуковой томографии “времени прибытия”: изучены свойства разреженности результатов полученных изображений, а также влияние увеличения разрешения снимка на числовую характеристику разреженности s. Пренебрежимое увеличение этой характеристики позволяет эффективно проводить обработку данных, даже при увеличении их исходного количества на много порядков. Также была разработана и апробирована на искусственных моделях модификация алгоритма реконструкции томографических снимков с применением техники опознания со сжатием. Эксперименты показали, что возможна реконструкция изображения без существенной потери качества при использовании неполного объема данных. Кроме того, некоторые варианты построения матриц измерений продемонстрировали положительное влияние на шумоустойчивость алгоритма. 30

Список литературы [1] Граничин ОН, Поляк БТ. Рандомизированные алгоритмы оценивания и оптимизации при почти произвольных помехах. –– Наука М., 2003. [2] Граничин Олег Николаевич, Павленко Дмитрий Валентинович. Рандомизация получения данных и ℓ1 -оптимизация (опознание со сжатием) // Автоматика и телемеханика. –– 2010. –– no. 11. –– P. 3– 28. [3] A simple proof of the restricted isometry property for random matrices / Richard Baraniuk, Mark Davenport, Ronald DeVore, Michael Wakin // Constructive Approximation. –– 2008. –– Vol. 28, no. 3. –– P. 253–263. [4] An improved automatic time-of-flight picker for medical ultrasound tomography / Cuiping Li, Lianjie Huang, Nebojsa Duric et al. // Ultrasonics. –– 2009. –– Vol. 49, no. 1. –– P. 61–72. [5] Berinde Radu, Indyk Piotr, Ružić M. Practical near-optimal sparse recovery in the l1 norm // Communication, Control, and Computing, 2008 46th Annual Allerton Conference on / IEEE. –– 2008. –– P. 198– 205. [6] Breast imaging using ultrasound tomography: From clinical requirements to system design / Olivier Roy, Signe Schmidt, Cuiping Li et al. // Ultrasonics Symposium (IUS), 2013 IEEE International / IEEE. –– 2013. –– P. 1174–1177. [7] Breast imaging with 3D ultrasound computer tomography: results of a first in-vivo study in comparison to MRI images / Torsten Hopp, Lukas Šroba, Michael Zapf et al. // Breast Imaging. –– Springer, 2014. –– P. 72–79. [8] Candè Emmanuel J, Wakin Michael B. An introduction to compressive 31

sampling // Signal Processing Magazine, IEEE. –– 2008. –– Vol. 25, no. 2. –– P. 21–30. [9] Candes Emmanuel, Romberg Justin. Sparsity and incoherence in compressive sampling // Inverse problems. –– 2007. –– Vol. 23, no. 3. –– P. 969. [10] Candès Emmanuel J, Romberg Justin, Tao Terence. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information // Information Theory, IEEE Transactions on. –– 2006. –– Vol. 52, no. 2. –– P. 489–509. [11] Chen Chen, Huang Junzhou. Compressive sensing MRI with wavelet tree sparsity // Advances in neural information processing systems. –– 2012. –– P. 1115–1123. [12] Chen Scott Shaobing, Donoho David L, Saunders Michael A. Atomic decomposition by basis pursuit // SIAM review. –– 2001. –– Vol. 43, no. 1. –– P. 129–159. [13] Compressed sensing MRI / Michael Lustig, David L Donoho, Juan M Santos, John M Pauly // Signal Processing Magazine, IEEE. –– 2008. –– Vol. 25, no. 2. –– P. 72–82. [14] Compressed sensing of ultrasound images: sampling of spatial and frequency domains / Céline Quinsac, Adrian Basarab, JeanMarc Girault, Denis Kouamé // Signal Processing Systems (SIPS), 2010 IEEE Workshop on / IEEE. –– 2010. –– P. 231–236. [15] Detection of breast cancer with ultrasound tomography: First results with the Computed Ultrasound Risk Evaluation (CURE) prototype / Nebojsa Duric, Peter Littrup, Lou Poulo et al. // Medical physics. –– 2007. –– Vol. 34, no. 2. –– P. 773–785. [16] Dines Kris A, Lytle R Jeffrey. Computerized geophysical tomography // Proceedings of the IEEE. –– 1979. –– Vol. 67, no. 7. –– P. 1065–1073. 32

[17] Donoho David L. Compressed sensing // Information Theory, IEEE Transactions on. –– 2006. –– Vol. 52, no. 4. –– P. 1289–1306. [18] Donoho David L. For most large underdetermined systems of linear equations the minimal L1-norm solution is also the sparsest solution // Communications on pure and applied mathematics. –– 2006. –– Vol. 59, no. 6. –– P. 797–829. [19] Granichin Oleg. Linear regression and filtering under nonstandard assumptions (Arbitrary noise) // Automatic Control, IEEE Transactions on. –– 2004. –– Vol. 49, no. 10. –– P. 1830–1837. [20] Leonid Kunyansky. Fast reconstruction algorithms for the thermoacoustic tomography in certain domains with cylindrical or spherical symmetries // Inverse Problems and Imaging. –– 2012. –– Vol. 6, no. 1. –– P. 111–131. [21] Lustig Michael, Donoho David, Pauly John M. Sparse MRI: The application of compressed sensing for rapid MR imaging // Magnetic resonance in medicine. –– 2007. –– Vol. 58, no. 6. –– P. 1182–1195. [22] Mallat Stéphane G, Zhang Zhifeng. Matching pursuits with timefrequency dictionaries // Signal Processing, IEEE Transactions on. –– 1993. –– Vol. 41, no. 12. –– P. 3397–3415. [23] Natarajan Balas Kausik. Sparse approximate solutions to linear systems // SIAM journal on computing. –– 1995. –– Vol. 24, no. 2. –– P. 227–234. [24] Ozmen N. Ultrasound Imaging Methods for Breast Cancer Detection. –– TU Delft, Delft University of Technology, 2014. [25] Quan Youli, Huang Lianjie. Sound-speed tomography using firstarrival transmission ultrasound for a ring array // Medical Imaging / International Society for Optics and Photonics. –– 2007. –– P. 651306– 651306. 33

[26] Robust ultrasound travel-time tomography using the bent ray model / Ali Hormati, Ivana Jovanovic, Olivier Roy, Martin Vetterli // Proc. SPIE. –– Vol. 7629. –– 2010. –– P. 76290I. [27] Schuster Arthur. An introduction to the theory of optics. –– E. Arnold, 1904. [28] Shannon Claude E. Communication in the presence of noise // Proceedings of the IRE. –– 1949. –– Vol. 37, no. 1. –– P. 10–21. [29] Treeby Bradley E, Cox Benjamin T. k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields // Journal of biomedical optics. –– 2010. –– Vol. 15, no. 2. –– P. 021314–021314. [30] Ultrasound transmission computed tomography of the breast. / JS Schreiman, JJ Gisvold, James F Greenleaf, RC Bahn // Radiology. –– 1984. –– Vol. 150, no. 2. –– P. 523–530. 34

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв