Сегментация изображений является задачей разбиения цифрового изображения на одну или несколько областей, представляющих интерес. Это фундаментальная проблема в области компьютерного зрения, которая решается многими различными способами, каждый из которых обладает своими преимуществами и недостатками.
В этой статье я кратко рассмотрю понятие метода фиксации уровня и неявно заданных динамических поверхностей (level set method). Также рассмотрю роль этого метода в бинарной сегментации с введением и определением математических конструкций, таких как SDT (Signed Distance Transforms), маркированной карты расстояний.
Слева — исходное изображение, справа — сегментированное
В различных отраслях часто возникает потребность получения сегментированных изображений. Например, в медицине, где существует интерес к получению цифровых 3D моделей органов для дальнейшего анализа. Также сегментация применяется в распознавании лиц, машинном зрении и т.д.
Сегментация изображений является особенно трудной задачей по нескольким причинам. Во-первых, в настоящее время не существует универсальных алгоритмов сегментации, поэтому необходимо создавать специальные алгоритмы, параметры которых подбираются пользователем, исходя из конкретных условий. Во-вторых, искажения изображения, такие как шум, неоднородность и т.д., очень трудно учесть в алгоритмах сегментации без участия пользователя.
Со второй проблемой обычно борются применением различного рода фильтров. Эта технология отлично описана здесь.
Я столкнулся с задачей сегментации изображений в нефтегазовой отрасли. Первоначально мы имели плоские изображения керна (как на рисунке выше). Получив сегментированное изображение, мы можем провести интерполяцию границ. Проведя данную операцию для всех изображений, мы можем получить математическую модель нашего образца горной породы, с которой можно работать и на которой, например, можно задавать течение флюидов (если рассматривать нашу модель как модель механики сплошной среды).
Этот метод был введен американскими математиками Стэнли Ошером и Джеймсом Сетьяном в 1980-е годы. Он стал популярным во многих дисциплинах, таких как компьютерная графика, обработка изображений, вычислительная геометрия, оптимизация, вычислительная гидродинамика и вычислительная биофизика.
Метод фиксации уровня неявно преобразует контур (в двух измерениях) или поверхность (в трех измерениях), используя функцию более высокой размерности, называемую функцией фиксации уровня . Преобразованный контур или поверхность можно получить как кривую .
Преимущество использования этого метода состоит в том, что топологические изменения, такие как слияние и разделение контура или поверхности, решаются неявно, как показано ниже на рисунке.
Можно видеть взаимосвязь между функцией фиксации уровня (снизу) и контуром (сверху). Видно, что изменение поверхности разбивает контур.
Задание контура или поверхности определяется уравнением модификации уровня. Численное решение этого дифференциального уравнения в частных производных мы получаем путем обновления для каждого временного слоя. Общий вид уравнения модификации уровня имеет вид:
где:
— евклидова норма,
— функция фиксации уровня,
— время,
— коэффициент скорости.
Это уравнение в частных производных очень похоже на гиперболическое. Но из-за того, что в правой части уравнения помимо производной по пространству также есть производная по времени оно, строго говоря, таким не является.
При определенном задании вида вектора с точностью до параметров, мы можем направлять линию уровня на разные области или формы.
Для приложений в сегментировании будем задавать так:
где:
— функция данных, склоняющая решения к желаемому сегментированию,
— значение интенсивности пиксела,
— средний член кривизны функции , который будет сохранять эту функцию гладкой,
— весовой параметр.
Функция данных действует как главная «сила», которая управляет сегментацией. Делая положительным в желаемых областях или отрицательным в не желаемых, модель будет стремиться к желаемому сегментированию. Вот простая функция данных:
Ее работа приведена ниже на рисунке:
То есть если значение пиксела попадает в интервал , то диапазон модели будет расширяться, иначе — сужаться.
Поэтому три пользовательских параметра, которые необходимо указать для сегментации, — это ? и .
Например, если нам нужно выделить на изображении объекты с низким уровнем яркости (черные области), то выбирается равным нулю.
Слева выделяем черные области, справа — желаемая область имеет более большую яркость, поэтому выбираем
Также требуется задание начальной маски . От того, какой мы ее выберем, будет зависеть начальное приближение . Для этого введем понятие маркированной карты расстояний
Карта расстояний – это матрица, которая имеет такую же размерность, что и наше изображение. Строится она таким образом: для каждого пиксела нашего изображения высчитывается значение, равное минимальному расстоянию от этого пиксела до ближайшего пиксела на границе объекта (ов). Математическое определение функции расстояния для множества S:
Маркированная карта расстояний — это матрица, имеющая положительный знак перед значениями пикселов, которые находятся вне объекта, и отрицательный для тех, которые внутри него. Это соглашение о знаке, которое будет выполняться при реализации, однако можно было бы использовать и противоположную конвенцию знака. Следует отметить, что значения ячеек зависят от выбранной метрики: некоторые общие метрики расстояния — это евклидово расстояние, расстояние Чебышёва ( Chebyshev distance, chessboard distance) и расстояние городских кварталов (Taxicab geometry, city block distance, Manhattan distance).
Слева выбранная нами начальная маска. Это объект любой формы, находящийся внутри интересующего нас объекта. Чаще всего это прямоугольник, потому что его довольно просто задать. После определенного количества итераций функция справа (функция модификации уровня) меняется таким образом, чтобы ее пересечение с нулевой плоскостью давало интересующую нас кривую
Уравнение (1), уравнение уровня, должно быть дискретизировано для последовательного или параллельного вычисления. Это делается с использованием разностной схемы с разностями против потока (upwind).
Метод первого порядка точности для дискретизации по времени уравнения (1) задается прямым методом Эйлера:
где:
— текущие значения в момент времени
— вектор скорости в момент времени
— значения градиента в момент времени
При вычислении градиента необходимо проявлять большую осторожность в отношении пространственных производных . Это лучше всего видно, если рассмотреть расширенный вид уравнения (4):
где — компоненты из . Для простоты рассмотрим одномерную форму уравнения (5) в конкретной точке сетки :
где — производная по пространству в точке . Для того, чтобы понять, какую разностную схему использовать, прямую или обратную, будем основываться на знаке в точке . Если , значения перемещаются слева направо, и поэтому следует использовать обратные разностные схемы ( в уравнениях 7). Напротив, если , то для приближения следует использовать прямые разностные схемы ( в уравнениях 7). Именно этот процесс выбора аппроксимации для производной по пространству, основанный на знаке , известен как расчет по потоку.
В двумерном случае для вычисления нам требуются производные, расписанные ниже. Стоит упомянуть, что мы делаем предположение изотропности изображения.
аппроксимируется с использованием схемы upwind:
Наконец, в зависимости от того, является ли или , представляется, как:
Коэффициент F скорости, как обсуждалось ранее, основан на значениях интенсивности пикселов и значениях кривизны.
Будем вычислять ее на основе значений текущего уровня, установленных с использованием приведенных ниже производных. В двух измерениях требуется только , наряду с определенными ранее производными:
Используя разность нормалей, кривизна вычисляется с использованием указанных выше производных с двумя нормалями и .
Две нормали используются для вычисления дивергенции, с помощью которой мы вычисляем среднюю кривизну, как показано ниже:
Устойчивость вытекает из условия Куранта-Фридрихса-Леви (CFL), которое утверждает, что скорость движения информации в вычислительной схеме должна быть такой, что за время вы не выходите за рамки одной ячейки, то есть . Имеем:
Реализуем данный алгоритм бинарной сегментации в системе Matlab. Сделаем это, в первую очередь, для наглядности, так как Matlab прекрасно справляется с работой с таблицами. Также там есть очень простые функции, позволяющие считывать изображения и выводить работу алгоритма на экран. Благодаря этому код будет сжатым и понятным. Естественно, реализация, например на C будет куда выгоднее с точки зрения скорости. Но в данной статье мы не гонимся за скоростью, а пытаемся разобраться в методе.
Псевдокод алгоритма можете видеть ниже:
Начнем с того, что выберем формат изображений, с которым будем работать. Очень удобным оказался формат PGM. Это открытый формат хранения растровых изображений типа bitmap без сжатия в оттенках серого. У него простой ASCII заголовок, а само изображение есть последовательность однобайтных (для оттенков серого от 0 до 255) целых чисел без знака. Подробнее о формате вы можете прочитать непосредственно в самом стандарте, ссылка приведена в разделе источников.
Разобьем наш код на две части: «пусковую установку» и «ядро», чтобы отделить инициализацию и код обновления набора уровней. Первый файл обрабатывает загрузку и уменьшает разрешение изображения с помощью функций “imread” и “imresize”. Пользователь задает параметры для пороговых значений , диапазона и взвешивания кривизны , запускает пусковую установку и затем переходит к рисованию замкнутого многоугольника, который будет формировать начальную маску (обеспечивающую некоторую базовую интерактивность).
Итак, считываем изображение и сразу же создаем нулевую матрицу , которая будет нашей маской:
Для создания маски удобно выбирать координаты центра и отклонения от него. Получим прямоугольник, заполненный единицами. Его и будем передавать в ядро.
Теперь рассмотрим функцию
Напишем функцию, которая будет высчитывать кривизну:
Производные вычисляются вычитанием сдвинутых матриц функции . Обратите внимание, что производные не вычисляются поэлементно, так как это было бы менее эффективно.
Теперь считаем коэффициент скорости
После переходим к обновлению уровней. Для этого сначала высчитываем градиент для соблюдения устойчивости:
И изменяем функцию набора уровней :
Будем выводить промежуточный результат на экран каждые 20 итераций:
С помощью функции
Я проводил тесты алгоритма на двух типах изображений: медицинских и изображений горных пород (так как решал эту задачу на дипломе).
Для медицинских изображений это работает вот так:
Мы видим, что алгоритм очень хорошо справляется с данной задачей.
Стоит сказать о том, что данный алгоритм можно реализовать для трехмерного случая. В этой статье я не буду расписывать особенности, но покажу как это работает:
В этом случае начальная маска — это куб (в общем случае, конечно, это любое трехмерное тело)
В сегментации изображений горных пород в первую очередь сталкиваешься с проблемой шума. Это связано со многими причинами: начиная от погрешности томографов, заканчивая ошибками в обратном преобразовании Радона (именно так получаются такие изображения). Именно поэтому сначала нужно изображения отфильтровать. Очень хорошо с этим справляется фильтр анизотропной диффузии Перона и Малика. Ссылка на статью об этом фильтре уже была в начале статьи. Оставлю ее и в конце поста.
Вот, что происходит, если сегментировать зашумленное изображение:
А вот, что с отфильтрованным изображением:
Метод фиксации уровня это отличный метод для проведения бинарной сегментации. Основным его недостатком является большая вычислительная сложность. Время, за которое Matlab на моей машине обрабатывал изображение 512х512 (было проведено 5000 итераций) составляет 2737.84 секунд (45 минут). Это, конечно, неимоверно много.
Но не стоит закрывать статью после прочтения этого факта. Самой главной причиной выбора мной этого алгоритма было то, что он очень хорошо параллелится. Я это делал с помощью CUDA, вы можете использовать что-то другое и написать об этом, например, в комментариях. Если кому-то интересно, то в перспективе буду писать следующую статью на тему распараллеливания этого алгоритма.
P.S. Вот ссылки на исходные коды: стартовый файл и функция, непосредственно отвечающая за сегментацию.
В этой статье я кратко рассмотрю понятие метода фиксации уровня и неявно заданных динамических поверхностей (level set method). Также рассмотрю роль этого метода в бинарной сегментации с введением и определением математических конструкций, таких как SDT (Signed Distance Transforms), маркированной карты расстояний.
Слева — исходное изображение, справа — сегментированноеЗачем все это?
В различных отраслях часто возникает потребность получения сегментированных изображений. Например, в медицине, где существует интерес к получению цифровых 3D моделей органов для дальнейшего анализа. Также сегментация применяется в распознавании лиц, машинном зрении и т.д.
Сегментация изображений является особенно трудной задачей по нескольким причинам. Во-первых, в настоящее время не существует универсальных алгоритмов сегментации, поэтому необходимо создавать специальные алгоритмы, параметры которых подбираются пользователем, исходя из конкретных условий. Во-вторых, искажения изображения, такие как шум, неоднородность и т.д., очень трудно учесть в алгоритмах сегментации без участия пользователя.
Со второй проблемой обычно борются применением различного рода фильтров. Эта технология отлично описана здесь.
Я столкнулся с задачей сегментации изображений в нефтегазовой отрасли. Первоначально мы имели плоские изображения керна (как на рисунке выше). Получив сегментированное изображение, мы можем провести интерполяцию границ. Проведя данную операцию для всех изображений, мы можем получить математическую модель нашего образца горной породы, с которой можно работать и на которой, например, можно задавать течение флюидов (если рассматривать нашу модель как модель механики сплошной среды).
Метод фиксации уровня
Этот метод был введен американскими математиками Стэнли Ошером и Джеймсом Сетьяном в 1980-е годы. Он стал популярным во многих дисциплинах, таких как компьютерная графика, обработка изображений, вычислительная геометрия, оптимизация, вычислительная гидродинамика и вычислительная биофизика.
Метод фиксации уровня неявно преобразует контур (в двух измерениях) или поверхность (в трех измерениях), используя функцию более высокой размерности, называемую функцией фиксации уровня . Преобразованный контур или поверхность можно получить как кривую .
Преимущество использования этого метода состоит в том, что топологические изменения, такие как слияние и разделение контура или поверхности, решаются неявно, как показано ниже на рисунке.
Можно видеть взаимосвязь между функцией фиксации уровня (снизу) и контуром (сверху). Видно, что изменение поверхности разбивает контур.
Уравнение модификации уровня
Задание контура или поверхности определяется уравнением модификации уровня. Численное решение этого дифференциального уравнения в частных производных мы получаем путем обновления для каждого временного слоя. Общий вид уравнения модификации уровня имеет вид:
где:
— евклидова норма,
— функция фиксации уровня,
— время,
— коэффициент скорости.
Это уравнение в частных производных очень похоже на гиперболическое. Но из-за того, что в правой части уравнения помимо производной по пространству также есть производная по времени оно, строго говоря, таким не является.
При определенном задании вида вектора с точностью до параметров, мы можем направлять линию уровня на разные области или формы.
Для приложений в сегментировании будем задавать так:
где:
— функция данных, склоняющая решения к желаемому сегментированию,
— значение интенсивности пиксела,
— средний член кривизны функции , который будет сохранять эту функцию гладкой,
— весовой параметр.
Функция данных действует как главная «сила», которая управляет сегментацией. Делая положительным в желаемых областях или отрицательным в не желаемых, модель будет стремиться к желаемому сегментированию. Вот простая функция данных:
Ее работа приведена ниже на рисунке:
То есть если значение пиксела попадает в интервал , то диапазон модели будет расширяться, иначе — сужаться.
Поэтому три пользовательских параметра, которые необходимо указать для сегментации, — это ? и .
Например, если нам нужно выделить на изображении объекты с низким уровнем яркости (черные области), то выбирается равным нулю.
Слева выделяем черные области, справа — желаемая область имеет более большую яркость, поэтому выбираем
Также требуется задание начальной маски . От того, какой мы ее выберем, будет зависеть начальное приближение . Для этого введем понятие маркированной карты расстояний
Маркированная карта расстояний (Signed Distance Transforms)
Карта расстояний – это матрица, которая имеет такую же размерность, что и наше изображение. Строится она таким образом: для каждого пиксела нашего изображения высчитывается значение, равное минимальному расстоянию от этого пиксела до ближайшего пиксела на границе объекта (ов). Математическое определение функции расстояния для множества S:
Маркированная карта расстояний — это матрица, имеющая положительный знак перед значениями пикселов, которые находятся вне объекта, и отрицательный для тех, которые внутри него. Это соглашение о знаке, которое будет выполняться при реализации, однако можно было бы использовать и противоположную конвенцию знака. Следует отметить, что значения ячеек зависят от выбранной метрики: некоторые общие метрики расстояния — это евклидово расстояние, расстояние Чебышёва ( Chebyshev distance, chessboard distance) и расстояние городских кварталов (Taxicab geometry, city block distance, Manhattan distance).
Слева выбранная нами начальная маска. Это объект любой формы, находящийся внутри интересующего нас объекта. Чаще всего это прямоугольник, потому что его довольно просто задать. После определенного количества итераций функция справа (функция модификации уровня) меняется таким образом, чтобы ее пересечение с нулевой плоскостью давало интересующую нас кривую
Дискретизация
Уравнение (1), уравнение уровня, должно быть дискретизировано для последовательного или параллельного вычисления. Это делается с использованием разностной схемы с разностями против потока (upwind).
Метод первого порядка точности для дискретизации по времени уравнения (1) задается прямым методом Эйлера:
где:
— текущие значения в момент времени
— вектор скорости в момент времени
— значения градиента в момент времени
При вычислении градиента необходимо проявлять большую осторожность в отношении пространственных производных . Это лучше всего видно, если рассмотреть расширенный вид уравнения (4):
где — компоненты из . Для простоты рассмотрим одномерную форму уравнения (5) в конкретной точке сетки :
где — производная по пространству в точке . Для того, чтобы понять, какую разностную схему использовать, прямую или обратную, будем основываться на знаке в точке . Если , значения перемещаются слева направо, и поэтому следует использовать обратные разностные схемы ( в уравнениях 7). Напротив, если , то для приближения следует использовать прямые разностные схемы ( в уравнениях 7). Именно этот процесс выбора аппроксимации для производной по пространству, основанный на знаке , известен как расчет по потоку.
В двумерном случае для вычисления нам требуются производные, расписанные ниже. Стоит упомянуть, что мы делаем предположение изотропности изображения.
аппроксимируется с использованием схемы upwind:
Наконец, в зависимости от того, является ли или , представляется, как:
Коэффициент F скорости, как обсуждалось ранее, основан на значениях интенсивности пикселов и значениях кривизны.
Кривизна
Будем вычислять ее на основе значений текущего уровня, установленных с использованием приведенных ниже производных. В двух измерениях требуется только , наряду с определенными ранее производными:
Используя разность нормалей, кривизна вычисляется с использованием указанных выше производных с двумя нормалями и .
Две нормали используются для вычисления дивергенции, с помощью которой мы вычисляем среднюю кривизну, как показано ниже:
Устойчивость
Устойчивость вытекает из условия Куранта-Фридрихса-Леви (CFL), которое утверждает, что скорость движения информации в вычислительной схеме должна быть такой, что за время вы не выходите за рамки одной ячейки, то есть . Имеем:
Последовательная реализация
Реализуем данный алгоритм бинарной сегментации в системе Matlab. Сделаем это, в первую очередь, для наглядности, так как Matlab прекрасно справляется с работой с таблицами. Также там есть очень простые функции, позволяющие считывать изображения и выводить работу алгоритма на экран. Благодаря этому код будет сжатым и понятным. Естественно, реализация, например на C будет куда выгоднее с точки зрения скорости. Но в данной статье мы не гонимся за скоростью, а пытаемся разобраться в методе.
Псевдокод алгоритма можете видеть ниже:
- Вход: Оригинальное изображение I; исходная маска m; порог T; ошибка ?; количество итераций N; число, показывающее через какое количество итераций пересчитываем карту расстояний, ITS.
- Выход: Результат сегментации – изображение seg.
Инициализируем ?_0 с помощью Signed Euclidean Distance Transform (SEDT) и маски m.
Высчитываем D(I)= ?-|I-T|. - Цикл i от 1 до N
- Высчитываем производные первого порядка
- Высчитываем производные второго порядка
- Высчитываем значения нормалей и
- Вычисляем значение градиента
- Вычисляем значение функции
- Обновляем функцию уровня
- Если
i % ITS == 0тогда
Переинициализировать с помощью SEDT
Конец если
Конец цикла
Начнем с того, что выберем формат изображений, с которым будем работать. Очень удобным оказался формат PGM. Это открытый формат хранения растровых изображений типа bitmap без сжатия в оттенках серого. У него простой ASCII заголовок, а само изображение есть последовательность однобайтных (для оттенков серого от 0 до 255) целых чисел без знака. Подробнее о формате вы можете прочитать непосредственно в самом стандарте, ссылка приведена в разделе источников.
Разобьем наш код на две части: «пусковую установку» и «ядро», чтобы отделить инициализацию и код обновления набора уровней. Первый файл обрабатывает загрузку и уменьшает разрешение изображения с помощью функций “imread” и “imresize”. Пользователь задает параметры для пороговых значений , диапазона и взвешивания кривизны , запускает пусковую установку и затем переходит к рисованию замкнутого многоугольника, который будет формировать начальную маску (обеспечивающую некоторую базовую интерактивность).
Итак, считываем изображение и сразу же создаем нулевую матрицу , которая будет нашей маской:
I = imread('test_output.pgm');
I = imresize(I, 1);
init_mask = zeros(size(I,1),size(I,2));
Для создания маски удобно выбирать координаты центра и отклонения от него. Получим прямоугольник, заполненный единицами. Его и будем передавать в ядро.
x = 400;
y = 800; %координаты центра маски
dx = 10;
dy = 10;
init_mask(y - dy : y + dy, x - dx : x + dx) = 1; %заполняем единицами
seg = simpleseg(I, init_mask, max_its, E, T); %Запускаем ядро программы
Теперь рассмотрим функцию
simpleseg, которая и является ядром кода. Вначале нужно получить SDF из маски init_mask:phi = mask2phi(init_mask); %получаем SDF
function phi = mask2phi(init_a) %функция получения SDF
phi=bwdist(init_a)-bwdist(1-init_a);
Напишем функцию, которая будет высчитывать кривизну:
%Вычисление кривизны
function curvature = get_curvature(phi)
dx = (shiftR(phi)-shiftL(phi))/2;
dy = (shiftU(phi)-shiftD(phi))/2;
dxplus = shiftR(phi)-phi;
dyplus = shiftU(phi)-phi;
dxminus = phi-shiftL(phi);
dyminus = phi-shiftD(phi);
dxplusy = (shiftU(shiftR(phi))-shiftU(shiftL(phi)))/2;
dxminusy = (shiftD(shiftR(phi))-shiftD(shiftL(phi)))/2;
dyplusx = (shiftR(shiftU(phi))-shiftR(shiftD(phi)))/2;
dyminusx = (shiftL(shiftU(phi))-shiftL(shiftD(phi)))/2;
nplusx = dxplus./sqrt(eps+(dxplus.^2 )+((dyplusx+dy )/2).^2);
nplusy = dyplus./sqrt(eps+(dyplus.^2 )+((dxplusy+dx )/2).^2);
nminusx= dxminus./sqrt(eps+(dxminus.^2)+((dyminusx+dy)/2).^2);
nminusy= dyminus./sqrt(eps+(dyminus.^2)+((dxminusy+dx)/2).^2);
curvature = ((nplusx-nminusx)+(nplusy-nminusy))/2;
Производные вычисляются вычитанием сдвинутых матриц функции . Обратите внимание, что производные не вычисляются поэлементно, так как это было бы менее эффективно.
%-- Вычисление производных с помощью сдвига матриц
function shift = shiftD(M)
shift = shiftR(M')';
function shift = shiftL(M)
shift = [ M(:,2:size(M,2)) M(:,size(M,2))];
function shift = shiftR(M)
shift = [M(:,1) M(:,1:size(M,2)-1)];
function shift = shiftU(M)
shift = shiftL(M')';
Теперь считаем коэффициент скорости
alpha=0.5;
D = E - abs(I - T);
K = get_curvature(phi);
F = alpha*single(D) + (1-alpha)*K;
После переходим к обновлению уровней. Для этого сначала высчитываем градиент для соблюдения устойчивости:
dxplus = shiftR(phi)-phi;
dyplus = shiftU(phi)-phi;
dxminus = phi-shiftL(phi);
dyminus = phi-shiftD(phi);
gradphimax_x = sqrt(max(dxplus,0).^2+max(-dxminus,0).^2);
gradphimin_x = sqrt(min(dxplus,0).^2+min(-dxminus,0).^2);
gradphimax_y = sqrt(max(dyplus,0).^2+max(-dyminus,0).^2);
gradphimin_y = sqrt(min(dyplus,0).^2+min(-dyminus,0).^2);
gradphimax = sqrt((gradphimax_x.^2)+(gradphimax_y.^2));
gradphimin = sqrt((gradphimin_x.^2)+(gradphimin_y.^2));
gradphi = (F>0).*(gradphimax) + (F<0).*(gradphimin);
dt = .5/max(max(max(abs(F.*gradphi))));
И изменяем функцию набора уровней :
phi = phi + dt.*(F).*gradphi;
Будем выводить промежуточный результат на экран каждые 20 итераций:
if(mod(its,20) == 0)
showcontour(I,phi,its);
subplot(2,2,4); surf(phi); shading flat;
end
С помощью функции
showcontour:function showcontour(I, phi, i)
subplot(2, 2, 3); title('Evolution') ;
imshow (I);
hold on;
contour(phi, [0 0],'g','LineWidth', 2);
hold off ; title([num2str(i) 'Iterations'] ); drawnow;
Результаты
Я проводил тесты алгоритма на двух типах изображений: медицинских и изображений горных пород (так как решал эту задачу на дипломе).
Для медицинских изображений это работает вот так:
Мы видим, что алгоритм очень хорошо справляется с данной задачей.
Стоит сказать о том, что данный алгоритм можно реализовать для трехмерного случая. В этой статье я не буду расписывать особенности, но покажу как это работает:
В этом случае начальная маска — это куб (в общем случае, конечно, это любое трехмерное тело)
В сегментации изображений горных пород в первую очередь сталкиваешься с проблемой шума. Это связано со многими причинами: начиная от погрешности томографов, заканчивая ошибками в обратном преобразовании Радона (именно так получаются такие изображения). Именно поэтому сначала нужно изображения отфильтровать. Очень хорошо с этим справляется фильтр анизотропной диффузии Перона и Малика. Ссылка на статью об этом фильтре уже была в начале статьи. Оставлю ее и в конце поста.
Вот, что происходит, если сегментировать зашумленное изображение:
А вот, что с отфильтрованным изображением:
Выводы и дальнейшая работа
Метод фиксации уровня это отличный метод для проведения бинарной сегментации. Основным его недостатком является большая вычислительная сложность. Время, за которое Matlab на моей машине обрабатывал изображение 512х512 (было проведено 5000 итераций) составляет 2737.84 секунд (45 минут). Это, конечно, неимоверно много.
Но не стоит закрывать статью после прочтения этого факта. Самой главной причиной выбора мной этого алгоритма было то, что он очень хорошо параллелится. Я это делал с помощью CUDA, вы можете использовать что-то другое и написать об этом, например, в комментариях. Если кому-то интересно, то в перспективе буду писать следующую статью на тему распараллеливания этого алгоритма.
P.S. Вот ссылки на исходные коды: стартовый файл и функция, непосредственно отвечающая за сегментацию.
Источники
Поделиться с друзьями
Hedgehogues
Простите, а что Вы имели ввиду, когда говорили, что уравнение (1) похоже на гиперболическое? В гиперболическом уравнении, как минимум, есть вторые производные.
Hedgehogues
Второй аспект, который мне не очень понятен — это модуль. Как вы собираетесь дифференцировать уравнение?
Поясните пожалуйста, чем в реальной задаче может являться функция данных, параметр скорости и весовой коэффициент? У Вас приводится функционал для D(I). Из каких соображений он выбирался? Почему берётся модуль, а не квадрат, который часто намного проще и удобнее дифференцировать?
Поправте пожалуйста разметку формул. Ужасно неудобно читать.
Hedgehogues
Hedgehogues
Что такое начальная маска m?
Про матрицу расстояний. Что Вы подразумеваете под границей объекта? Откуда объекты берутся. Предполагается, что они уже выделены?
Что означает вот эти фразы? Про какие пикселы ведётся речь непонятно совершенно. Вы говорите о чём-то вроде свёрток?
Поправьте пожалуйста разметку кода. Тоже лють какая-то.
В целом, видно, что человек проделал какую-то работу. Что-то вроде как понимает (или должен понимать), но доносит он до народа это ужасно. Набор формул, которые плохо связаны друг с другом. Мне не нравится. Хотел разобраться и закодить. Но понял, что себе дороже: слишком много времени потрачу на разбор формул, в которых, видимо, есть местами ошибки. Плюс ко всему, они плохо оформлены, что тоже накладывает отпечаток.
Поправьте пожалуйста статью. Спасибо (=
22dla
Начальная маска, если смотреть на рисунке, это кривая, которая полностью находится внутри интересующего нас объекта. То есть это начальное условие для нашего уравнения.
Я не понимаю вашего вопроса про объекты. На изображении невооруженным глазом видны объекты, которые нам интересны: это может быть мозг, печень, да и в принципе что угодно. Их мы и стараемся сегментировать, думаю, вы это понимаете.
То есть пользователь, видя картинку, задает маску так, чтобы она находился внутри объекта (который он видит невооруженным взглядом), потом алгоритм меняет эту маску, расширяя ее или сужая.
22dla
Производные аппроксимируются вычитанием матриц, читайте внимательнее.
В этой задаче функция D выбиралась из соображений задачи сегментации: выделить на изображении нужные мне объекты. Предложили этот функционал американские математики, поэтому вопрос про квадрат лучше задать им. Но, в принципе, можно пробовать разные функционалы в зависимости от задачи.
Статью поправил
Gregivy
Обратите внимание на функцию F. Там за счет скалярного, как я понял, произведения оператора градиента на градиент функции «фи», выходят ее вторые производные.
Hedgehogues
Да, видимо. Я потом об этом тоже подумал. Но дальше без бутылки не разобраться. Автор писал статью для себя. Грустно… Человек, сталкнувшийся с этим алгоритмом первый раз потратит уйму времени, прежде чем поймёт, что до него хотели донести. Картинки, такое ощущение вообще натырены.
22dla
Про картинки вообще необоснованная претензия. Если что-то вам непонятно, вопросы задавайте более конкретно. Я скинул рабочий исходный код, с помощью которого вы можете протестить все вышесказанное.
Hedgehogues
Картинки — это в давесок ко всему остальному. Я задал довольно много вопросов. Но Вы не ответили на них.