Разработка программы выделения в реальном времени гладких объектов в трехмерных медицинских изображениях

Размещено на http://www.allbest.ru/

Введение

В работе рассматривается выделение объектов в трёхмерных медицинских изображениях.

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

После получения изображения нужной части тела пациента, одной из задач является выделение замкнутых объектов. Такими объектами могут являться органы пациента, местоположение и объём которых может изменяться при некоторых заболеваниях, и поэтому должен быть определен при исследовании. Также такими объектами могут быть различные образования и опухоли. Основываясь на изображениях таких образований, может быть принято решение о проведении дополнительных процедур для определения злокачественности образования.

Полученные после проведения медицинских процедур данные представляются в стандартизированном формате DICOM, содержащем серии изображений пациента. Поскольку объекты трёхмерны, одного среза недостаточно для правильного выделения. Поэтому серии изображений рассматриваются как трёхмерное изображение, которое может быть представлено трёхмерным массивом численных значений. В работе рассматриваются способы анализа и обработки именно таких сущностей.

Процесс выделения сегментов, содержащих замкнутые объекты, или сегментация, похож на трёхмерную вариацию алгоритма заливки двумерной области цветом. Существуют и другие подходы к сегментации, но в данной работе за основу был взят именно этот по причинам, рассмотренным далее. Также, как и двумерная заливка, эта операция сталкивается с проблемой «протекания». Это означает, что алгоритм выходит за пределы выделяемого объекта, и заполняет области, не относящиеся к нему, но имеющие сходную яркость. Есть несколько причин такого протекания. Во-первых, неизбежно присутствие дефектов, вызванных шумом изображения, полученного с медицинских аппаратов. Вторая причина - сложность и разнородность биологических систем. На медицинских изображениях яркость точки является репрезентацией плотности ткани тела. Человеческий организм состоит из множества различных тканей, и нет гарантии, что на изображении не окажется тканей со схожей плотностью, но не имеющих отношения к выделяемому объекту.

В связи с этим в работе [6] был предложен критерий для предотвращения протекания алгоритма сегментации. Исходя из того, что выделяемые объекты имеют органическую природу, было сделано предположение о том, что в большинстве случаев они имеют гладкую поверхность и высокий радиус кривизны поверхности. Поэтому критерий для предотвращения протекания был основан на вычислении значения кривизны поверхности в точке, которая является кандидатом на добавление к выделяемому объекту. Были проведены эксперименты с двумерными и трёхмерными искусственными объектами, рассмотрены несколько методов сегментации с контролем кривизны. Методы были оценены по критериям применимости, настраиваемости и эффективности сегментации.

Изученные методы применимы в разной степени. Асимптотическая вычислительная сложность для каждого из них линейно зависит от числа рассматриваемых точек, поскольку каждый из методов использует для подсчёта кривизны фиксированное количество точек. В абсолютном измерении самым вычислительно простым способом является способ подсчёта соседей, проверяющий значения 26 точек, а наиболее сложными являются способы подсчёта через производные, поскольку они рассматривают окрестность 5х5х5 выбранной точки, то есть 124 соседних точки при подсчёте производных второго порядка. Однако, данные подсчёты производятся при помощи операции свёртки, эффективная реализация которой позволяет оптимизировать скорость вычислений.

С точки зрения качества получаемых результатов наилучшим образом показывают себя два последних рассматриваемых метода: усреднение по кубу и метод гессиана. При использовании этих методов на различных двумерных и трёхмерных формах удаётся успешно выделить поверхность гладко, без выбросов. Также эти методы позволяют наиболее тонко настраивать жёсткость критериев отсечения как пользователю, так и в автоматическом режиме при помощи динамического подсчёта порога кривизны.

Из рассмотренных и разработанных методов в результате сравнения был выбран наиболее подходящий - метод подсчёта кривизны через производные при помощи модифицированного оператора Собеля, обозначенный как метод усреднения по кубу. Этот метод использован в работе для построения на его основе полноценного инструмента для сегментации трёхмерных изображений.

В работе [6] внимание уделялось в основном применимости алгоритма и качеству сегментации. В меньшей мере рассматривалась сложность алгоритмов и производительность программы. В настоящей работе эти вопросы являются приоритетными, и значительно влияют на список поставленных задач. Важной является оптимизация времени работы программы и требуемой памяти.

Для повышения производительности программы существует несколько подходов. Первый - разработка эффективной с точки зрения памяти и скорости работы операций доступа структуры данных для хранения точек, которые рассматриваются как кандидаты на добавление к объекту на очередной итерации. Второй - введение быстро работающего критерия, определяющего наиболее очевидно подходящие или не подходящие точки, без запуска более тяжеловесной полноценной проверки кривизны. Наконец, имеет смысл исследовать возможность работы алгоритма на изображениях с пониженным разрешением, и влияние такого снижения на качество работы.

В рамках данной выпускной квалификационной работы была поставлена цель - разработать программу выделения в реальном времени гладких объектов в трехмерных медицинских изображениях. Для достижения этой цели необходимо решить следующе задачи:

1. Изучить возможности эффективной работы с трёхмерными матрицами

2. Изучить формат хранения медицинских снимков DICOM

3. Разработать структуру данных для ресурсно-эффективного хранения точек, претендующих на добавление к объекту

4. Реализовать алгоритм сегментации трёхмерных изображений с контролем кривизны

5. Провести эксперименты по выделению объектов на реальных медицинских изображениях

6. Разработать программу для сегментации трёхмерных медицинских изображений

7. Разработать техническую документацию

Работа структурирована следующим образом: в первой главе рассмотрены более ранние работы, посвященные обзору методов сегментации трёхмерных изображений, в том числе применительно к медицинским целям, а также использование кривизны поверхности для анализа изображений, структуры данных для хранения разреженного набора точек и прочие методики, которые могут быть полезны при разработке. Во второй главе подробно описываются алгоритмы и структуры, выбранные для реализации в данной работе. В третьей главе описываются используемые технологии и особенности разработки программы. В приложении А приведено техническое задание на разработку программы выделения в реальном времени гладких объектов в трехмерных медицинских изображениях.

Глава 1. Кривизна и анализ дискретных поверхностей

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

1.1 Способы сегментации изображений

Сегментация - это разделение изображение на несколько непересекающихся сегментов. В материале [7] проведен достаточно подробный сравнительный анализ различных алгоритмов, которые осуществляют сегментацию, и применяются в анализе медицинских изображений. Описаны достоинства и недостатки каждого из способов, заданы границы области применимости алгоритмов. Также в этой работе предоставляется общая последовательность действий при сегментации изображения, которая может быть использован в дальнейшем.

Для начала зададим определение сегментации изображения так, как оно использовано в [7]

Способы сегментации делятся на 3 типа по числу учитываемых параметров:

· Один параметр - интенсивность. Алгоритм использует только яркость точки для определения принадлежности к сегменту

· Два параметра - интенсивность и положение. Алгоритм использует яркость точки и информацию о её окружении

· Три параметра - интенсивность, положение и форма сегмента. Алгоритм использует яркость, окружение точки и заранее известную информацию о предполагаемой форме выделяемого объекта, например, математическое описание этой формы

В общем случае можно сказать, что чем больше информации использует алгоритм, тем он может быть точнее. Однако в области выделения объектов на медицинских изображениях, предсказать форму можно лишь для определенного набора хорошо известных и изученных объектов, тогда как для использования того же алгоритма для новой задачи необходимо строить математическое описание формы нового объекта. Такие подходы используются во многих аналогах. Например, в работах [8] и [9] приводится описание сегментации с целью выделения лёгочных узелков. Однако для построения программного решения, позволяющего выполнять более широкий спектр задач, от такого подхода следует отказаться.

Также алгоритмы сегментации можно классифицировать по скорости:

· Мгновенные - простые и неточные методы сегментации;

· Интерактивные - выполняемые в пределах сессии работы пользователя с программой;

· Отложенные - выполняемые на распределенных машинах, состоящие из множества фаз, сочетающие несколько алгоритмов.

По степени автономности:

· Ручные - постоянное вмешательство пользователя;

· Полуавтоматические - пользователь меняет параметры алгоритма, задает первое приближение сегмента;

· Автоматические - полностью выполняемые без участия пользователя.

Мгновенные алгоритмы не подходят для выполнения качественной сегментации медицинских изображений, поскольку низкая точность не позволяет выполнять поставленные задачи диагностики. Отложенные алгоритмы теоретически могут быть применимы, но нежелательны, поскольку если алгоритм требует много времени и больших вычислительных мощностей, то быстрее и менее затратно провести сегментацию вручную. Таким образом, для выполнения поставленных задач необходим интерактивный алгоритм. Это обычно означает, что алгоритм должен состоять только из одной основной фазы обработки, а дополнительные методы и преобразования должны относиться к мгновенным алгоритмам и быть отключаемыми.

Говоря о степени автономности алгоритма, нужно учитывать, что, для того чтобы использовать полностью автоматически алгоритм, нужно быть уверенным, что он сработает с достаточной точностью на любых входных параметрах. Это возможно опять же при наличии узкой задачи выделения одного определенного типа объектов, но затруднительно при большем разнообразии входных данных. Ручные же алгоритмы не дают большого преимущества перед выполнением диагностики полностью без средств компьютерной обработки изображений. Поэтому имеет смысл рассматривать полуавтоматические алгоритмы. Здесь стоит оговориться, что даже полуавтоматические алгоритмы имеют свои ограничения. Например, задача выделения тазобедренного сустава выполнима лучше всего только инструментами ручной сегментации, о чём также сказано в [7].

В работе [10], как и в [7], в качестве двух основных подходов к сегментации определяются анализ цветового пространства изображения и наращение региона. Одним из самых простых примеров первого типа является пороговая фильтрация изображения. Пример результата такой фильтрации приведен на рис. 1.

В отношении монотонных изображений это означает, что точки с яркостью ниже определенного порога считаются частью фона, а точки с яркостью выше этого порога являются частью объекта. Этот метод является очень простым как для реализации, так и для выполнения. Однако, данный метод не учитывает никаким образом пространственную информацию о точках. За счёт скорости данного алгоритма, существует возможность его использования в качестве вспомогательного, без потери производительности. Автор [10] в своей работе предлагает особый метод пороговой фильтрации изображения, основанный на рассмотрении сегментов изображения как нечётких множеств второго типа, и вычисления порогового значения путём максимизации нечеткости. В несколько упрощенном виде, этот алгоритм может быт использован для бинаризации изображения или повышения его контрастности.

а)

б)

Рисунок 1 Пороговая фильтрация. а): Исходное изображение б): результат пороговой фильтрации

Алгоритмы наращивания региона - это группа алгоритмов, к которым и относится метод сегментации с контролем кривизны, который берётся за основу в данной квалификационной работе. Эти алгоритмы итеративны, на каждой итерации к выделяемому объекту добавляются новые точки из числа точек-кандидатов. Демонстрация принципа работы этих алгоритмов приведена на рис. 2. Алгоритмы из этой группы работают по следующему принципу: На начальном этапе определяется некоторый стартовый регион сегментации, который считается принадлежащим к искомому сегменту. Чаще всего такой регион задается пользователем. На каждой итерации алгоритма точки, находящиеся рядом с точками стартового региона, назначаются точками, претендующими на добавление к региону. Основная часть алгоритма состоит в том, чтобы определить, следует ли добавлять точку-претендент к региону. После того, как точки-претенденты обработаны, некоторые из них добавлены к региону, другие отвергнуты, назначается новый набор точек-претендентов уже на основе расширенного региона и алгоритм начинает новую итерацию.

Рисунок 2 Принцип работы алгоритмов наращивания регионов [7]

Автор [10] приводит в качестве недостатка этой группы алгоритмов чрезмерную сегментацию, то есть разделение изображения на слишком большое число сегментов. Поскольку в данной работе используется выделение только одного сегмента, с инициированием пользовательской базовой областью, эта проблема полностью исключается.

Наконец, вводятся критерии оценки алгоритма сегментации и рекомендации по созданию качественного алгоритма. Качественный алгоритм должен быть

· Надежным

· Быстрым

· Точным

И для построения такого алгоритма следует рассматривать два подхода

· Комбинация нескольких методов и алгоритмов

· Максимальное использование заранее имеющейся информации об объектах

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

1.2 Кривизна в сегментации изображений

Метрики кривизны поверхности нередко используется в задачах сегментации и обработки трёхмерных изображений, в том числе медицинских. Различные проекты используют разные способы подсчёта кривизны, в зависимости от формата представления данных и решаемой задачи. Методы подсчёта необходимых метрик, выбранные для ВКР, описаны в соответствующем разделе. В работе [8] на основе кривизны построен алгоритм снижения шума изображения. В работе [9] кривизна лёгочных узелков используется для дифференциации доброкачественных и злокачественных образований. В работе [11] метод, основанный на измерении кривизны, используется для детектирования изображений костей на компьютерных томограммах.

1.3 Структуры хранения точек в трёхмерном изображении

В различных задачах, связанных с трёхмерными объектами, их визуализацией и изменением, их чаще всего представляют в виде триангулированных объектов. [12] Это упрощает перемещение, аффинные преобразования объектов и прочие операции. Также это уменьшает размер информации, которую необходимо хранить о каждом объекте.

Однако такое представление объектов не слишком применимо для задач сегментации. Важно сохранять информацию если не о каждой точке изображения, то о каждой рассматриваемой на конкретной итерации алгоритма. Таким образом, необходимо иметь структуру данных, позволяющую эффективно работать с трехмерным массивом точек, часто разреженным, при этом по возможности не держа в памяти большое число копий матрицы изображения целиком, и имея возможность быстрого доступа к соседям точки. Поэтому были рассмотрены различные примеры структур данных, хранящих информацию о смежных точках в трёхмерном пространстве. Это часто древовидные структуры или построенные по принципу многомерного связного списка. Также применяется подход с использованием kd-деревьев [13] и графоподобная структура DCEL [14]. Структуры данных различаются по количеству дополнительно хранимой информации, скорости доступа и возможностям пространственной ориентации, то есть получения информации о взаимном расположении точек. DCEL предполагает сильно разреженные точки и модель, сходную с полигональной, с дополнительными связями между вершинами. Это удобно для различных целей, но не подходит для работы с изображениями, заданными трехмерным массивом дискретных точек, так как хранимая дополнительная информация становится избыточной. Структуры на основе многосвязных списков также несколько избыточны за счет того, что каждый элемент хранит данные о нескольких других, однако это позволяет хранить данные об искривленной поверхности в трёхмерном пространстве. Однако, для реализации алгоритма сегментации с контролем кривизны необходима возможность получить соседей точки-кандидата, не принадлежащей к выделенной поверхности. Поэтому, пространственные данные многосвязного списка оказываются мало востребованными. Древовидные структуры и, в особенности, kd-деревья, позволяют структурировать набор точек по различным принципам, сохраняя быстрый доступ по координатам и возможность использования множества существующих алгоритмов работы с деревьями.

1.4 Общий порядок операций при сегментации

Для сегментации медицинских изображений имеет смысл взять за основу общий порядок действий, описанный в [7]:

· Повышение качества изображения

· Уменьшение размера данных

· Основная сегментация

· Уточнение сегментации

· Повышение качества сегментации.

· Удаление изолированных точек

· Заполнение дыр

Такая последовательность действий позволяет решить многие проблемы, возникающие при сегментации. Повышение качества изображения позволяет снизить уровень шума, чтобы чётче отделить объект от окружающего фона. Разумеется, алгоритм сегментации в любом случае должен быть устойчив к шуму, поскольку полностью избавиться от него невозможно. Уменьшение размера данных - один из способов повышения производительности алгоритма. Балансируя между качеством и скоростью работы, можно найти наилучший уровень сжатия изображения. Основная сегментация - в данной работе, процесс выделения объекта с контролем кривизны. Также, в данном проекте планируется не использовать алгоритмов уточнения сегментации, но позволить пользователю добавлять или отменять итерации алгоритма, для его уточнения. Наконец, повышение качества сегментации - это удаление образовавшихся дефектов и сглаживание выделенного объекта. В случае необходимости, с этой задачей можно справиться при помощи операций математической морфологии.

1.5 Прочие полезные методики

Стоит отметить работу [15], цели и задачи которой достаточно схожи с теми, что поставлены в данной работе. Эта работа посвящена сегментации томограмм мозга с использованием алгоритма наращивания регионов. В ней можно найти несколько полезных идей, которые могут быть применены. Одна из ключевых идей для упрощения работы пользователя с программой и повышения эффективности алгоритма - это использование пользовательских стартовых областей в качестве ядра для наращивания. Это означает, что пользователь, например, мышью выделяет некую область изображения, которая автоматически становится частью будущего сегмента. На первой итерации алгоритма наращивание начинается именно с этой области. Кроме того, алгоритм из [15] сохраняет статистическую информацию об уже выделенных точках объекта. Это позволяет анализировать распределение яркости в объекте без хранения всех его точек в памяти.

В первой главе был проведен обзор материалов, посвященных сегментации трёхмерных изображений, выделены критерии оценки качества алгоритмов сегментации. Также была рассмотрена базовая последовательность действий при сегментации. Помимо этого, были изучены структуры данных, предназначенные для хранения точек в трёхмерном пространстве, и прочие методики, которые могут помочь в разработке.

Глава 2. Описание выбранных методов

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

2.2 Метод выделения объектов с контролем кривизны

В предыдущей работе [6] рассмотрены различные методы сегментации изображений путём наращивания региона с контролем кривизны поверхности.

· Метод количества соседей

· Метод математической морфологии

· Метод итеративного подсчёта производных

· Метод подсчёта производных по Собелю

· Метод усреднения по кубу

· Метод гессиана

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

Главным преимуществом данного способа является его простота как в реализации, так и вычислительная - для каждого из соседей проверяется его значение один раз и сумма этих значений сравнивается с заданным порогом. Также, результаты работы данного метода легко представимы в наглядном виде, а численные значения метрики возможно анализировать с точки зрения здравого смысла (число соседей является понятной на бытовом уровне метрикой).

Основной недостаток такого способа - высокая дискретизация получаемых значений. Вариантов количества соседей всего 8 в двухмерном случае и 26 в трёхмерном, что противоречит потребности в плавном и гладком регулировании поверхности выделяемого объекта. Эксперименты показали низкую точность данного способа, и от его использования было принято решение отказаться.

Метод математической морфологии основан на том, чтобы измерять радиус кривизны поверхности при помощи структурных сферических элементов различного радиуса. Для этого необходимо производить для каждой точки несколько раз операцию морфологии, последовательно увеличивая радиус. Это оказалось слишком затратно с точки зрения вычислительной эффективности. Однако, операции математической морфологии показали свою работоспособность для постобработки выделенных сегментов с целью устранения дефектов и сглаживания. Поэтому, метод имеет своё применение, и подробнее описан в соответствующем разделе ВКР.

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

Следующий метод, решающий эти проблемы, это метод подсчёта производных по Собелю.

Метод также основан на определении кривизны через производные. Общая формула подсчёта кривизны выглядит таким образом:

Для трёхмерного случая, подсчитываются минимальная и максимальная кривизна по всем направлениям, эти значения называются главными кривизнами, и на их основе вычисляются метрики средней и полной кривизны (кривизны Гаусса)

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

Для подсчёта производных в точках дискретной поверхности в данном методе используется модифицированный оператор Собеля. Обычно данный оператор используется в обработке двумерных изображений для поиска границ объектов, поскольку высокое значение производной является признаком резкого перепада яркости точек. [16]

Принцип работы оператора Собеля заключается в том, что область изображения подвергается свёртке с матрицей оператора и результирующая матрица содержит значения производных в каждой из точек по направлению, заданному матрицей оператора. Оригинальная матрица Собеля для производной по X:

На рис. 3 приведен принцип работы свёртки на примере функций одной переменной. Считаем, что функции дискретны, и, следовательно, являются векторами значений. Рисунок показывает последовательное применение свертки раз функции y с функциями h. В первой строке изображена исходная кусочно-линейная функция , ядро свёртки и результат свёртки . На второй строке полученная функция подвергается свертке с ядром и так далее. Прочие обозначения на изображении условны, и служат для иллюстрации того, каким образом изменяются значения функции при свёртке.

Рисунок 3 Принцип работы операции свёртки на примере векторов значений

Подсчёт производной операцией свёртки эффективнее итеративного подхода, можно экспериментировать с различными значениями коэффициентов в матрице, а получаемые значения более разнообразны, что позволяет тоньше настраивать порог кривизны.

Одним из методов, показавшим наиболее высокую точность результатов и возможность наиболее тонкой настройки является метод, модифицирующий метод производной по Собелю, названный «Усреднением по всему кубу».

Для целей сегментации трёхмерных изображений, этот оператор был видоизменен с целью адаптации к трёхмерному пространству.

В методе усреднения по кубу используется трёхмерное ядро свёртки со следующими коэффициентами:



В данном случае значения в каждой точке ядра считаются как , где М(c.v) - манхэттенское расстояние до центральной точки, в котором подсчитывается производная. Для каждой точки изображения I с координатами i, j, k свёртка с такой маской подсчитывается как

В данном примере индексация элементов ядра свёртки идётся с единицы для удобства.

Если для какой-либо точки не существует соседей: необходимых для подсчёта (например, точки с координатами i-1,j,k), то есть точка находится на краю изображения хотя бы по одному из измерений, то существует несколько способов разрешения такой проблемы. Отзеркаливание, то есть принятие несуществующих точке равными существующим с другой стороны от рассматриваемой точки. Закругление, то есть принятие несуществующих точек равными точкам с другой стороны изображения по данной оси. Игнорирование, то есть не вычисление свёртки для крайних точек. Первые два используются при обработке изображений, но не могут являться источником достоверной информации о кривизне поверхности за пределами изображения, то есть в области, не попавшей в сканирование. Поэтому, принято решение выбрать способ игнорирования, и не учета граничных точек в алгоритме. Поскольку для подсчета кривизны требуется вычислить вторую производную, описываемый оператор свертки сначала применяется ко всем точка в 3x3x3-окружении точки, для которой идёт подсчёт, а затем применяется повторно к самой рассматриваемой точке.

Общая схема работы алгоритма выглядит следующим образом:

Создается сферическая базовая область с центром в точке, указанной пользователем. Вычисляется начальный набор точек-кандидатов на добавление к объекту. Для каждой из точек-кандидатов вычисляется метрика кривизны поверхности, и если она оказывается выше порогового значения, точка добавляется к объекту. Список точек кандидатов пересчитывается для объекта с добавленными точками. Алгоритм продолжается до тех пор, пока на очередной итерации не оказывается, что ни одна из точек-кандидатов не подходит для добавления.

Эксперименты показали, что данный метод успешно справляется с выделение объектов различных форм, без заполнения дефектных областей. Именно этот метод был выбран для реализации проекта ВКР.

Также было разработано два подхода к определению порога кривизны для алгоритма.

Поскольку метрика кривизны не является нормализованной, возможные значения кривизны варьируются в зависимости от коэффициентов матрицы свёртки от нуля до некоторого максимума. Для выбранных коэффициентов этот максимум равен 12500.

Максимальное значение метрика принимает в точках, со всех сторон, окруженных уже выделенными точками, а нулевое значение принимается при полном их отсутствии. Отсутствие нормализации позволяет производить все вычисления в целых числах. Первый подход к выбору порога кривизны подразумевает ввод порогового значения пользователем, что увеличивает настраиваемость алгоритма под размер значимых дефектов. При использовании второго подхода порог кривизны вычисляется автоматически на каждой итерации как 40-процентный квантиль значений кривизны в точках-кандидатах и точках, добавленных на предыдущей итерации. Такой подход более универсален и позволяет избавить пользователя от заботы о подборе оптимального порогового значения. Этот подход был назван динамическим.

2.2 Метод пороговой фильтрации

В качестве метода предобработки изображения с целью уменьшения шума и повышения контрастности был выбран метод пороговой фильтрации, основанный на нечетких множествах 2 типа. Это метод подробно описан в работе [10]. Сегмент изображения, то есть выделяемый объект, рассматривает в данном методе как нечеткое множество. Выбрав функцию принадлежности к данному множеству, можно рассчитать функцию от значений принадлежности всех точек изображения, называемую функцией неопределенности. Для этого строится гистограмма яркости изображения, и функция неопределенности высчитывается на основе значений этой гистограммы. Затем функция принадлежности сдвигается вдоль уровня яркости точек, и находится сдвиг, при котором функция неопределенности изображения принимает максимальное (или минимальное) значение. Функция принадлежности должна быть симметричной, и, следовательно, иметь центр симметрии по оси яркости. В положении функции принадлежности, в котором достигается экстремальное значение, запоминаем положение центра симметрии функции. Это значение уровня яркости и становится пороговым для проведения фильтрации изображения. Поскольку используются нечёткие множества второго типа, вместо одного значения функции принадлежности, используются верхняя и нижняя границы «отпечатка неопределенности» (footprint of uncertainty). Они подсчитываются как

где µ_U - верхняя граница,

µ_L - нижняя граница,

а µ(x) - значение функции принадлежности в точке. Итоговая функция неопределенности изображения выглядит как

Где А - трёхмерное изображение,

М, N и K - размеры изображения по трём измерениям,

L - число градаций яркости (например, 256 для восьмибитных изображений),

h(g) - значение гистограммы яркости изображения для уровня g.

2.3 Метод постобработки изображения

На случай, если при сегментации изображения не удалось избежать мелких дефектов, и объект получился с отверстиями, граница поверхности недостаточно гладкая, имеет смысл разработать метод постобработки изображения, повышающий качество сегментации. Также он может помочь избавиться от лишних шумовых пикселей вне объекта, но при использовании алгоритма сегментации через наращивание региона такие пиксели не могут образовываться. Для этих целей хорошо подходят операции математической морфологии. Бинарная математическая морфология определяет операции над двоичными изображениями, преобразования при помощи структурных элементов. Структурный элемент также является двоичным изображением, меньшего размера, чем преобразуемое изображение. Результатом операций математической морфологии является изображение, совпадающее по размеру с исходным. Для структурного элемента вводится особая точка, которая называется начальной, и имеет ключевое значение в обработке. Часто используются структурные элементы в виде прямоугольника или круга. В таких элементах начальной точкой является геометрический центр фигуры. В процессе обработки, называемой «зондированием», структурный элемент совмещается начальной точкой поочередно с каждой точкой изображения. Для каждого такого наложения проверяется некоторое условие, определяемое типом операции, и если условие выполняется, то в результирующее изображение по координатам этой точки ставится бинарная единица (белый цвет), в противном случае ставится ноль (черный цвет). Существуют две основные операции математической морфологии - наращивание и эрозия, а также производные операции, являющиеся комбинацией последовательно применяемых основных операций. В описании операций под «белой» точкой подразумевается точка бинарного изображения, в которой значение равно единице, а «чёрной» точкой называется точка с значением ноль.

Наращивание (dilation) - операция, используемая для расширения исходного объекта с сохранением его контуров. Для каждой точки, при наложении на неё структурного элемента, выполняется следующая операция:

Если начальная точка структурного элемента наложена на белую точку

То закрасить все точки, которые покрыты структурным элементом

Эрозия (erosion) - операция, используемая для сокращения исходного объекта. Она является противоположной к наращиванию. Для каждой точки, при наложении на неё структурного элемента, выполняется следующая операция:

Если все точки, покрываемые структурным элементом, белые

То закрасить точку, на которую наложена начальная точка структурного элемента

Операции открытия (open) и закрытия (close) являются производными, и определяются как комбинации наращивания и эрозии. Для изображения I и структурного элемента S определения выглядят как



Морфологические операции со структурным элементом в виде круга с начальной точкой в центре позволяют производить сглаживание изображения и удаление дефектов после сегментации.

2.4 Используемые структуры данных

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

Первое предположение, которое было сделано, состоит в том, что в единственном экземпляре изображение всё же необходимо хранить полностью. Потоковая обработка не представляется возможной из-за постоянной необходимости получения соседних точек, а обеспечение случайного доступа к точкам без хранения в памяти приводит к огромному числу операций обращения к жёсткому диску, что замедляет работу алгоритма. Таким образом, исходим из того, что существует объект с исходным изображение, к которому можно обратиться по координатам за O(1) или получить сегмент изображения размером MxNxK за O(MNK).

Наиболее простой способ хранения списка точек-кандидатов - это запись максимальных и минимальных координат точек-кандидатов по всем трём измерениям. Таким образом, хранение всего шести чисел позволяет обозначить верхнюю границу числа точек-кандидатов. Существенным недостатком в данном случае является то, что фактически рассматриваются не только точки-кандидаты, но все точки параллелепипеда, описанного вокруг массива этих точек. Таким образом, проверке подвергаются не только точки-кандидаты, но и точки, уже добавленные к объекту, а также, при сложной форме объекта, и множество лишних точек. Это не сказывается на качестве алгоритма. Но увеличивает число операций в нём. Однако такой способ позволяет быстро вычислить границы области точек-кандидатов и требует минимального размера памяти.

Другой способ хранения точек, претендующих на добавление, это хранение всех их координат в некой структуре данных. Таким образом, в памяти хранятся все значения точек-кандидатов. Простой пример такой структуры данных - стек, основанный на множестве без повторов. Для заполнения такого стека достаточно добавлять туда все точки, граничащие с недавно добавленными, а стек автоматически обработает дублирования. Это позволяет хранить только нужные точки, и уменьшает число проверок. Однако за счёт проверки на уникальность, сложность вставки в такой стек O(logN). Кроме того, такое хранение ничего не говорит о том, как именно сгруппированы точки-кандидаты. Эта информация была бы полезна для определения приоритетных направлений роста алгоритма. Приоритетными направлениями считаются области изображения, в которых на предыдущих итерациях алгоритма было выделено большее число точек. Высока вероятность того, что в таких областях не находится граница объекта, и нет необходимости проводить полную трудоемкую процедуру проверки кривизны. Это позволяет ускорить заполнение основной части объекта, проводя контроль кривизны только ближе к поверхности, для предотвращения протекания. Поэтому имеет смысл рассмотреть третий способ хранения точек - в виде kd-дерева. Kd-дерево - это пространственная структура данных для организации точек в многомерном пространстве. По своей сути, это бинарное дерево, в котором каждый лист является точкой пространства, а каждый не-листовой узел является разделяющей гиперплоскостью в данном пространстве. На каждом уровне дерева происходит разделение пространства по одной из осей. Таким образом, поиск точки по её координатам, и добавление новой точки, имеют логарифмическую сложность. При этом можно оценить плотность расположения точек в той или иной части пространства, что и помогает выделить приоритетные направления роста сегмента. Структура такого дерева приведена на рис. 6. Это позволит быстрее определять плотность пространственной группировки точек, обрабатывать приоритетные направления, а также иметь информацию о наличии среди точек-соседей других кандидатов.

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

В данной главе были описаны основные алгоритмы, необходимые для осуществления поставленных задач: выделения объекта с контролем кривизны, пред- и постобработки изображения. Также рассмотрены вариации структур данных для хранения оперативной информации о состоянии алгоритма.

Глава 3. Разработка программы

В этой главе описываются используемые технологии и особенности разработки программы.

3.1 Язык и платформа разработки

Ключевым отличием данной работы, от предшествующей ей состоит в том, что высокое значение придается производительности программы и возможности её реального применения. Для этого были введены требования к времени работы отдельных алгоритмов, а также к наличию минимального графического пользовательского интерфейса (см. приложение А, техническое задание). В связи с этим, от разработки прототипов в среде Matlab было принято решение перейти к разработке на языке C++.

Этот язык отличается возможностью низкоуровневой работы с памятью, а его компиляторы обладают масштабными возможностями по оптимизации кода. Также, наличие ориентированных на производительность библиотек работы с матрицами, с похожим на Matlab программным интерфейсом, позволяет минимизировать сложность переноса кода.

В качестве платформы разработки была выбрана система MS Windows и специализированные для неё среду разработки и компилятор - Visual Studio 2015 и MS Visual C++. Такой выбор платформы связан с тем, что она наиболее используема среди потенциальных пользователей. Также, среда разработки VS предоставляет широкие возможности по профилированию программы и мониторингу состояния памяти процесса, что позволяет эффективно отслеживать производительность программы, и находить узкие места, требующие оптимизации.

3.2 Средства работы с матрицами

Для работы изображениями в формате серий матриц было рассмотрено несколько инструментов. Первый - библиотека Armadillo, реализующая широкий спектр функций линейной алгебры и имеющая открытый исходный код. Эта библиотека показывает превосходство в скорости, благодаря оптимизации операций работы с сегментами матриц, а также минимизации числа временных значений и цепному выполнению операций. [18]

Она основана на объединении двух других библиотек линейной алгебры - BLAS и LAPACK. Это позволяет при линковке программы заменять различные реализации этих двух библиотек, и находить наиболее производительные для данной задачи.

Было проведено сравнение скорости работы трёхмерных массивов Armadillo (называемых в библиотеке Cubes, то есть кубы), и обычных трёхмерных массивов С++, построенных на указателях. Сравнение проводилось на операциях инициализации и заполнения матрицы, на операциях взятия сегмента матрицы, сериализации в файл, а также проведения операции свёртки с небольшим ядром. Сравнение показало, что в среднем функции библиотеки Armadillo работают быстрее, чем самописные, однако на больших кубах, порядка 1000х1000х1000 библиотека переставала справляться с нагрузкой и блокировала работу программы и компьютера. Поэтому было принято решение использовать самописные классы и функции для хранения целого объекта изображения, но операции по сегментации проводить на отдельных частях изображения, преобразованных в Armadillo-массивы. Ещё одним преимуществом такого подхода является то, что исходное изображение можно хранить в компактном формате матриц с целочисленными значениями без знака, а вычисления проводить на Armadillo-кубах небольшого размера, но со значениями с плавающей точкой. Таким образом, без увеличения требуемой памяти можно снизить вероятность влияния ошибок округления. Также в библиотеке отсутствует операция трёхмерной свёртки, однако имеются исследования [19], показывающие, что, если ядро свёртки является разложимым, более эффективно произвести несколько операций свёртки с ядром меньшей размерности, и затем объединить результаты. Поскольку используемые в алгоритмах данной работы ядра свёртки являются разложимыми, было принято решение реализовать самостоятельно трёхмерную свёртку через такое разложение.

Для измерения производительности работы отдельных алгоритмов, операций и программы в целом были использованы средства профилирования и отладки Visual Studio, а также библиотека ctime.

Ещё одна библиотека, связанная с обработкой матриц и изображений, OpenCV, может предоставить несколько полезных функций, в том числе для быстрой визуализации серий изображений. Для использования обеих библиотек, а также собственных методов работы с трёхмерными массивами были написаны вспомогательные функции-конвертеры, позволяющие оперировать изображениями с совокупными возможностями трёх технологий. Проведенные замеры подтвердили, что временные затраты на конвертацию не слишком большие, и оправданы количеством предоставляемых возможностей.

3.3 Средства визуализации

Одним из требований к работе является предоставление пользователю простого интерфейса взаимодействия с программой, для загрузки изображений и запуска работы алгоритмов. Также необходима возможность вывода на экран изображения, представленного в виде матрицы OpenCV на экран, совместно с элементами интерфейса. С этой целью была выбрана платформа WinForms. Эта платформа предоставляет весь необходимый инструментарий для создания несложных графических интерфейсов, в том числе визуальный редактор, поставляемый с Visual Studio. Главную сложность при разработке графического интерфейса составляет взаимодействие модулей программы, написанных на C++/CLI и С++, с передачей данных между объектами в управляемой и неуправляемой памяти. Тем не менее, это позволяет пользоваться скоростью низкоуровневой работы в модулях логики программы, а также удобством управляемых классов в модулях визуального взаимодействия. Также стоит отметить наличие средств высокоуровневого взаимодействия с файловой системой Windows.

Для чтения изображений в формате DICOM, была использована библиотека Imebra, поддерживающая большую часть стандарта. Библиотека позволяет читать файлы DICOM, производить конвертацию цветового поля, считывать метаданные и получать прямой доступ к непосредственным данным изображения. С её помощью возможно преобразовывать данные в любой из используемых форматов трёхмерных массивов и работать как с обычными числовыми данными.

3.4 Особенности реализации

В работе реализованы методы автоматического и полуавтоматического выполнения для стадий пороговой фильтрации и основной сегментации. Пороговую фильтрацию возможно проводить по заданному пользователем порогу, а также автоматически по порогу, подсчитанному через максимизацию функции неопределенности для всего изображения или только для стартового региона сегментации.

Сегментацию возможно проводить со статическим порогом кривизны, задаваемым пользователем, а также с динамическим порогом, подсчитываемым на каждой итерации на основе анализа значения предыдущей итерации:

xes = Utils::convert_d(image->getRegion(minx, miny, maxx, maxy, minz, maxz));

sample = wBright * Utils::convert_d(res_image->getRegion(minx, miny, maxx, maxy, minz, maxz)) + wCurv * xes;

sample = Utils::conv3(Utils::conv3(sample, kernel), kernel);

sample %= -xes + 1;

values = sort(nonzeros(vectorise(sample)));

thres = values(quantile * values.size());

Для реализации операций математической морфологии необходимо разработать критерий применения условий срабатывания операций. Для математического описания условий, таких как «структурный элемент должен полностью быть наложен на объект» можно использовать сложение матриц, и, поскольку известно, что они бинарные, использовать агрегированные значения для определения срабатывания условия. Так, например, для операции наращивания

for (int i = minx; i <= maxx; i++)

{

for (int j = miny; j <= maxy; j++)

{

for (int k = minz; k <= maxz; k++)

{

cube subcube = body.subcube(j, i, k, size(kernel));

subcube += kernel;

result(j + size_y, i + size_x, k + size_z) =

subcube.max() > 1 ? 1: 0;

}

}

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

Для операции эрозии

for (int i = minx; i <= maxx; i++)

{

for (int j = miny; j <= maxy; j++)

{

for (int k = minz; k <= maxz; k++)

{

cube subcube = body.subcube(j, i, k, size(kernel));

subcube -= kernel;

result(j + size_y, i + size_x, k + size_z) =

subcube.min() < 0 ? 0: 1;

}

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

Ещё одним важным требованием к реализации является предоставление возможности корректирования результатов работы алгоритма при помощи отмены последних итераций или же наоборот искусственного продления работы алгоритма на одну или несколько итераций после того, как был достигнут критерий останова. Для отмены итераций используется стек истории ограниченного размера, позволяющий удалить из выделенного объекта точки, добавленные на последних итерациях. Для продления работы алгоритма по команде пользователя снижается порог кривизны и проводится повторный проход по точкам-кандидатам последней итерации.

В данной главе были рассмотрены технологии, выбранные для реализации программы, библиотеки для эффективной работы с матрицами и визуализации результатов.

Глава 4. Экспериментальная проверка алгоритмов

В этой главе приводится описание практических экспериментов, проведенных в процессе написания ВКР, и результатов этих экспериментов

4.1 Данные для экспериментов

Тестирование предлагаемых алгоритмов проводилось в несколько этапов по ходу разработки. Первоначальные эксперименты по проверке работоспособности основного алгоритма сегментации с контролем кривизны были проведены на искусственно созданных бинарных изображениях различной геометрической формы. К объектам были добавлены дефекты, потенциально способствующие протеканию алгоритма. Основной задачей данных экспериментов было проверить, возможно ли при помощи контроля кривизны нейтрализовать влияние дефектов на результат сегментации. Такой подход позволил изучить работоспособность прототипа алгоритма, до реализации программной инфраструктуры. Эксперименты проводились в среде Matlab.

Первое тестовое изображение является моделью ситуации в реальных данных, когда два объекта схожей плотности соединены тонкими перемычками других тканей или шумовых дефектов. Проверяется, сможет ли алгоритм сегментации, выделив один из объектов, избежать распространения по перемычкам на второй. Второе изображение - куб с несколькими дефектами - служит для проверки работы алгоритма в случае, когда форма выделяемого объекта не близка к сферической, но тем не менее результат должен быть таким же - объект должен быть выделен исключая дефектные выбросы. Эксперименты с этими двумя изображениями, а также другими формами входных данных позволили выбрать наиболее подходящий алгоритм контроля кривизны из изученных. Третье изображение - сфера с дефектами различного размера, направления и расположения относительно центра объекта. Такой объект тестирования позволил изучить зависимость порога кривизны, при котором получаются удовлетворительные результаты, от особенностей предотвращаемых дефектов.

Серый график «Частично» показывает, какого размера выбросы были частично усечены, но оставались признаки того, что эти выбросы были. Небольшие выступы и неровности. График «Полностью» соответственно показывает размер выбросов, которые были удалены полностью без видимых остатков. Отдельный график «3D» показывает отсечение неплоских выбросов.

По полученным в данном эксперименте результатам невозможно однозначно определить вид зависимости. Видно, что зависимость положительная, что говорит об общей применимости критерия, а также что значение порога около 1000 позволяет эффективно отсечь большую часть выбросов. График «3D» показывает более плавный и стабильный рост, поскольку неплоские выбросы имеют более чёткую зависимость между кривизной и размером выброса. В отличие от резко возрастающих графиков, описывающих поведение на плоских дефектах, график «3D» показывает видимую тенденцию. Также, он является наиболее важным, поскольку плоские шумы или побочные объекты редко встречаются в реальной медицинской практике.

Изучение зависимости размера отсекаемых дефектов от порога кривизны не позволило вычислить универсальный порог для все изображений, однако помогло в подборе величин порога для последующих экспериментов, а также в разработке динамического метода вычисления порога.

...