Как стать автором
Поиск
Написать публикацию
Обновить

Бинарная сегментация изображений методом фиксации уровня (Level set method)

Время на прочтение10 мин
Количество просмотров13K
Алгоритмы*Обработка изображений*Математика*Matlab*
Из песочницы
Сегментация изображений является задачей разбиения цифрового изображения на одну или несколько областей, представляющих интерес. Это фундаментальная проблема в области компьютерного зрения, которая решается многими различными способами, каждый из которых обладает своими преимуществами и недостатками.

В этой статье я кратко рассмотрю понятие метода фиксации уровня и неявно заданных динамических поверхностей (level set method). Также рассмотрю роль этого метода в бинарной сегментации с введением и определением математических конструкций, таких как SDT (Signed Distance Transforms), маркированной карты расстояний.

Слева — исходное изображение, справа — сегментированное

Зачем все это?


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

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

Со второй проблемой обычно борются применением различного рода фильтров. Эта технология отлично описана здесь.

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

Метод фиксации уровня


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

Метод фиксации уровня неявно преобразует контур (в двух измерениях) или поверхность (в трех измерениях), используя функцию более высокой размерности, называемую функцией фиксации уровня $φ(x,t)$. Преобразованный контур или поверхность можно получить как кривую $\Gamma(x,t)=\{φ(x,t)=0\}$.
Преимущество использования этого метода состоит в том, что топологические изменения, такие как слияние и разделение контура или поверхности, решаются неявно, как показано ниже на рисунке.

Можно видеть взаимосвязь между функцией фиксации уровня (снизу) и контуром (сверху). Видно, что изменение поверхности разбивает контур.

Уравнение модификации уровня


Задание контура или поверхности определяется уравнением модификации уровня. Численное решение этого дифференциального уравнения в частных производных мы получаем путем обновления $φ $ для каждого временного слоя. Общий вид уравнения модификации уровня имеет вид:

$ \frac{∂φ}{∂t} = -F|∇φ| \qquad(1)$


где:
$| · |$ — евклидова норма,
$φ$ — функция фиксации уровня,
$t$ — время,
$F $ — коэффициент скорости.

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

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

Для приложений в сегментировании будем задавать $F $ так:

$F = αD(I)+(1-α)∇\cdot\frac{∇φ}{|∇φ|} \qquad(2)$


где:

$D(I)$ — функция данных, склоняющая решения к желаемому сегментированию,
$I$ — значение интенсивности пиксела,
$∇\cdot\frac{∇φ}{|∇φ|}$ — средний член кривизны функции $φ$, который будет сохранять эту функцию гладкой,
$\alpha$ — весовой параметр.

Функция данных $D(I) $ действует как главная «сила», которая управляет сегментацией. Делая $D $ положительным в желаемых областях или отрицательным в не желаемых, модель будет стремиться к желаемому сегментированию. Вот простая функция данных:

$ D(I)= ε-|I-T| \qquad(3) $


Ее работа приведена ниже на рисунке:


То есть если значение пиксела попадает в интервал $(T - ε, T+ ε)$, то диапазон модели будет расширяться, иначе — сужаться.

Поэтому три пользовательских параметра, которые необходимо указать для сегментации, — это $T,$ ε и $\alpha$.
Например, если нам нужно выделить на изображении объекты с низким уровнем яркости (черные области), то $T$ выбирается равным нулю.


Слева выделяем черные области, справа — желаемая область имеет более большую яркость, поэтому выбираем $T = 45$

Также требуется задание начальной маски $m$. От того, какой мы ее выберем, будет зависеть начальное приближение $φ$. Для этого введем понятие маркированной карты расстояний

Маркированная карта расстояний (Signed Distance Transforms)


Карта расстояний – это матрица, которая имеет такую же размерность, что и наше изображение. Строится она таким образом: для каждого пиксела нашего изображения высчитывается значение, равное минимальному расстоянию от этого пиксела до ближайшего пиксела на границе объекта (ов). Математическое определение функции расстояния $d: R^3 \to R $ для множества S:

$d(r,S)=min|r-S| \qquad для всех \qquad r∈R^3$


Маркированная карта расстояний — это матрица, имеющая положительный знак перед значениями пикселов, которые находятся вне объекта, и отрицательный для тех, которые внутри него. Это соглашение о знаке, которое будет выполняться при реализации, однако можно было бы использовать и противоположную конвенцию знака. Следует отметить, что значения ячеек зависят от выбранной метрики: некоторые общие метрики расстояния — это евклидово расстояние, расстояние Чебышёва ( Chebyshev distance, chessboard distance) и расстояние городских кварталов (Taxicab geometry, city block distance, Manhattan distance).


Слева выбранная нами начальная маска. Это объект любой формы, находящийся внутри интересующего нас объекта. Чаще всего это прямоугольник, потому что его довольно просто задать. После определенного количества итераций функция справа (функция модификации уровня) меняется таким образом, чтобы ее пересечение с нулевой плоскостью давало интересующую нас кривую


Дискретизация


Уравнение (1), уравнение уровня, должно быть дискретизировано для последовательного или параллельного вычисления. Это делается с использованием разностной схемы с разностями против потока (upwind).

Метод первого порядка точности для дискретизации по времени уравнения (1) задается прямым методом Эйлера:

$\frac{( φ^{t+∆t}-φ^t)}{∆t}+F^t\cdot∇φ^t=0 \qquad (4) $


где:

$φ^t$ — текущие значения $φ $ в момент времени $t $
$F^t$ — вектор скорости в момент времени $t $
$∇φ^t$ — значения градиента $φ$ в момент времени $t $

При вычислении градиента необходимо проявлять большую осторожность в отношении пространственных производных $φ$. Это лучше всего видно, если рассмотреть расширенный вид уравнения (4):

$\frac{(φ^{t+∆t}-φ^t)}{∆t}+u^t φ_x^t+v^t φ_y^t+w^t φ_z^t=0 \qquad (5)$


где $u, v, w$ — компоненты $x, y, z$ из $F$. Для простоты рассмотрим одномерную форму уравнения (5) в конкретной точке сетки $x_i$:

$\frac{(φ^{t+∆t}-φ^t)}{∆t}+u_i^t \cdot(φ_x )_i^t=0 \qquad (6)$


где $(φ_x )_i^t $ — производная $φ$ по пространству в точке $x_i$. Для того, чтобы понять, какую разностную схему использовать, прямую или обратную, будем основываться на знаке $u_i$ в точке $x_i$. Если $u_i>0$, значения $φ$ перемещаются слева направо, и поэтому следует использовать обратные разностные схемы ( $D_x^-$ в уравнениях 7). Напротив, если $u_i<0$, то для приближения $φ_x$ следует использовать прямые разностные схемы ($D_x^+$ в уравнениях 7). Именно этот процесс выбора аппроксимации для производной $φ$ по пространству, основанный на знаке $u_i$, известен как расчет по потоку.

В двумерном случае для вычисления $φ$ нам требуются производные, расписанные ниже. Стоит упомянуть, что мы делаем предположение изотропности изображения.

$\begin{matrix} D_x=\frac{φ_{i+1,j}-φ_{i-1,j} }2\qquad D_x^+=\frac{φ_{i+1,j}-φ_{i,j}}2\qquad D_x^-=\frac{φ_{i,j}-φ_{i-1,j}}2\\ D_y=\frac{φ_{i,j+1}-φ_{i,j-1} }2\qquad D_y^+=\frac{φ_{i,j+1}-φ_{i,j}}2\qquad D_y^-=\frac{φ_{i,j}-φ_{i,j-1}}2 \end{matrix} \qquad (7)$


$∇φ$ аппроксимируется с использованием схемы upwind:

$∇φ_{max}= \begin{bmatrix} \sqrt{max(D_x^+,0)^2+max(-D_x^-,0)^2} \\ \sqrt{max(D_y^+,0)^2+max(-D_y^-,0)^2} \end{bmatrix} \qquad (8)$


$∇φ_{min}= \begin{bmatrix} \sqrt{min(D_x^+,0)^2+min(-D_x^-,0)^2} \\ \sqrt{min(D_y^+,0)^2+min(-D_y^-,0)^2} \end{bmatrix} \qquad (9)$


Наконец, в зависимости от того, является ли $F_{i,j,k}>0$ или $F_{i,j,k}<0$, $∇φ$ представляется, как:

$∇φ = \begin{cases} ‖∇φ_{max} ‖_2 & \quad \text{если } F_{i,j,k}>0 \\ ‖∇φ_{min} ‖_2 & \quad \text{если } F_{i,j,k}<0 \\ \end{cases} \qquad (10)$


$φ(t+∆t)=φ(t)+ ∆tF|∇φ|\qquad (11)$


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

Кривизна


Будем вычислять ее на основе значений текущего уровня, установленных с использованием приведенных ниже производных. В двух измерениях требуется только $D_x^{+y},D_x^{-y},D_y^{+x},D_y^{-x}$, наряду с определенными ранее производными:

$\begin{matrix} D_x^{+y}=\frac{(φ_{i+1,j+1}-φ_{i-1,j+1})}2 \qquad D_x^{-y}=\frac{(φ_{i+1,j-1}-φ_{i-1,j-1})}2 \\ D_y^{+x}=\frac{(φ_{i+1,j+1}-φ_{i+1,j-1})}2 \qquad D_y^{-x}=\frac{(φ_{i-1,j+1}-φ_{i-1,j-1})}2 \end{matrix} \qquad (12)$



Используя разность нормалей, кривизна вычисляется с использованием указанных выше производных с двумя нормалями $n^+$ и $n^-$.

$n^+= \begin{bmatrix} \frac{D_x^+}{\sqrt{(D_x^+)^2+(\frac{D_y^{+x}+D_y}{2})^2}} \\ \frac{D_y^+}{\sqrt{(D_y^+)^2+(\frac{D_x^{+y}+D_x}{2})^2}} \\ \end{bmatrix} \qquad (13)$


$n^-= \begin{bmatrix} \frac{D_x^-}{\sqrt{(D_x^-)^2+(\frac{D_y^{-x}+D_y}{2})^2}} \\ \frac{D_y^-}{\sqrt{(D_y^-)^2+(\frac{D_x^{-y}+D_x}{2})^2}} \\ \end{bmatrix} \qquad (14)$


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

$H=\frac{1}{2}∇\cdot\frac{∇φ}{|∇φ|} =\frac{1}{2}((n_x^+-n_x^- )+(n_y^+-n_y^- ))\qquad (15)$


Устойчивость


Устойчивость вытекает из условия Куранта-Фридрихса-Леви (CFL), которое утверждает, что скорость движения информации в вычислительной схеме должна быть такой, что за время $∆t$ вы не выходите за рамки одной ячейки, то есть $\frac{∆x}{∆t}>|u|$. Имеем:

$∆t<\frac{∆x}{max{|u|}} \qquad (16)$


Последовательная реализация


Реализуем данный алгоритм бинарной сегментации в системе Matlab. Сделаем это, в первую очередь, для наглядности, так как Matlab прекрасно справляется с работой с таблицами. Также там есть очень простые функции, позволяющие считывать изображения и выводить работу алгоритма на экран. Благодаря этому код будет сжатым и понятным. Естественно, реализация, например на C будет куда выгоднее с точки зрения скорости. Но в данной статье мы не гонимся за скоростью, а пытаемся разобраться в методе.

Псевдокод алгоритма можете видеть ниже:

  • Вход: Оригинальное изображение I; исходная маска m; порог T; ошибка ε; количество итераций N; число, показывающее через какое количество итераций пересчитываем карту расстояний, ITS.
  • Выход: Результат сегментации – изображение seg.
    Инициализируем φ_0 с помощью Signed Euclidean Distance Transform (SEDT) и маски m.
    Высчитываем D(I)= ε-|I-T|.
  • Цикл i от 1 до N
    • Высчитываем производные первого порядка
    • Высчитываем производные второго порядка
    • Высчитываем значения нормалей $n^+$ и $n^-$
    • Вычисляем значение градиента $∇φ$
    • Вычисляем значение функции $F= αD(I)+(1-α)∇\cdot\frac{∇φ}{|∇φ|}$
    • Обновляем функцию уровня $φ(t+∆t)= φ(t)+ ∆t\cdot F\cdot|∇φ|$
    • Если i % ITS == 0 тогда
      Переинициализировать $φ$ с помощью SEDT
      Конец если

    Конец цикла

Начнем с того, что выберем формат изображений, с которым будем работать. Очень удобным оказался формат PGM. Это открытый формат хранения растровых изображений типа bitmap без сжатия в оттенках серого. У него простой ASCII заголовок, а само изображение есть последовательность однобайтных (для оттенков серого от 0 до 255) целых чисел без знака. Подробнее о формате вы можете прочитать непосредственно в самом стандарте, ссылка приведена в разделе источников.

Разобьем наш код на две части: «пусковую установку» и «ядро», чтобы отделить инициализацию и код обновления набора уровней. Первый файл обрабатывает загрузку и уменьшает разрешение изображения с помощью функций “imread” и “imresize”. Пользователь задает параметры для пороговых значений $T$, диапазона $ε$ и взвешивания кривизны $\alpha$, запускает пусковую установку и затем переходит к рисованию замкнутого многоугольника, который будет формировать начальную маску (обеспечивающую некоторую базовую интерактивность).

Итак, считываем изображение и сразу же создаем нулевую матрицу $m$, которая будет нашей маской:

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')';

Теперь считаем коэффициент скорости $F$

 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. Вот ссылки на исходные коды: стартовый файл и функция, непосредственно отвечающая за сегментацию.

Источники


  1. Статья на википедии про метод
  2. Статья по поводу фильтров
Теги:
Хабы:
+18
Комментарии10
5
Карма
0
Рейтинг
Евгений Волков @22dla

Пользователь

Отправить сообщение

Публикации

Показать лучшие за всё время
Разблокировать темную тему

Истории

Ближайшие события

Серия занятий «Тренировки по алгоритмам 5.0» от Яндекса
Дата1 марта – 19 апреля
Время19:00
Место
Онлайн
Подробнее в календаре
Дизайнеры, выбирайте себе компанию по вайбам на Хабр Карьере
Дата15 – 28 апреля
Место
Онлайн
Подробнее в календаре
Weekend Offer в AliExpress
Дата20 – 21 апреля
Время10:00 – 20:00
Место
Онлайн
Подробнее в календаре
Вебинар «Истории Dodo Brands и Offprice: как компании мигрировали 1С в облако»
Дата23 апреля
Время11:00
Место
Онлайн
Подробнее в календаре
Вебинар «Информационная архитектура компании»
Дата23 апреля
Время19:00
Место
Онлайн
Подробнее в календаре
Вебинар «Специфика синхронных и асинхронных инсертов в ClickHouse»
Дата23 апреля
Время20:00
Место
Онлайн
Подробнее в календаре
Конференция «Я.Железо»
Дата18 мая
Время14:00 – 23:59
Место
Москва • Онлайн
Подробнее в календаре
Firebird Conf: конференция для разработчиков и администраторов СУБД Firebird
Дата6 июня
Время09:00 – 20:00
Место
Москва
Подробнее в календаре

Ваш аккаунт

Разделы

Информация

Услуги

Техническая поддержка
© 2006–2024, Habr