Введение

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

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

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

Основные операции бинарной морфологии

Для того, чтобы провести операцию бинарной морфологии, нам потребуются

Исходное бинарное изображение (B)
Структурирующий элемент (S)

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

Ниже представлены некоторые примеры структурирующих элеметнов в виде прямоугольника, диска и кольца

$Box(3, 5) = \begin{matrix}*&*&*&*&*\\*&*&*&*&*\\*&*&*&*&*\end{matrix}$

$Disk(5, 5) = \begin{matrix}0&*&*&*&0\\*&*&*&*&*\\*&*&*&*&*\\*&*&*&*&*\\0&*&*&*&0\end{matrix}$

$Ring(5, 5) = \begin{matrix}0&*&*&*&0\\*&0&0&0&*\\*&0&0&0&*\\*&0&0&0&*\\0&*&*&*&0\end{matrix}$

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

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

Рассмотрим основные операции бинарной морфологии.

Наращивание бинарного изображения B структурирующим элементом S обозначается B ⊕ S и задается следующим образом

$B \oplus S = \underset{b\in B}{\cup} S_b$

Здесь оператор объединения применяется ко всем окрестным пикселям. Структурирующий элемент S применяется ко всем пикселям бинарного изображения. Каждый раз, когда начало координат структурирующего элемента совмещается с единичным бинарным пикселом, ко всему структурирующему элементу применяется перенос и последующее логическое сложение ИЛИ с соответствующими пикселями бинарного изображения.

Рассмотрим пример такого преобразования.

Figure 1: Исходное бинарное изображение B

$\begin{matrix} 0&0&0&0&0&0&0&0\\ *&*&*&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&*&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&*&*&0&0&0&0\\ 0&0&0&0&0&0&0&0 \end{matrix}$

Figure 2: Структурирующий элемент S

$\begin{matrix} *&*&*\\ *&*&*\\ *&*&* \end{matrix}$

Figure 3: Наращивание B ⊕ S

$\begin{matrix} *&*&*&*&*&*&*&*\\ *&*&*&*&*&*&*&*\\ *&*&*&*&*&*&*&*\\ 0&*&*&*&*&*&*&*\\ 0&*&*&*&*&*&*&*\\ 0&*&*&*&*&*&*&*\\ 0&*&*&*&*&*&*&*\\ 0&*&*&*&*&0&0&0\\ \end{matrix}$

Разберем код F#, реализующий алгоритм наращивания.

let (>.) value checkValue = if value > checkValue then checkValue else value
let (<.) value checkValue = if value < checkValue then checkValue else value

let subArrayOfMaks array2d mask (centerX, centerY) (maskCenterX, maskCenterY) = 
        let (rows, cols) = Matrix.sizes {values = array2d}
        let (maskRows, maskCols) = Matrix.sizes {values = mask}

        let x1 = centerX - maskCenterX
        let y1 = centerY - maskCenterY
        let x2 = x1 + maskCols - 1
        let y2 = y1 + maskRows - 1

        (x1 <. 0, x2 >. (cols - 1), y1 <. 0, y2 >. (rows - 1))

let upbuilding array2d mask (maskCenterX, maskCenterY) =

        let sub source mask (x, y) (maskCenterX, maskCenterY) =
            let (x1, x2, y1, y2) = subArrayOfMaks source mask (x, y) (maskCenterX, maskCenterY)
            source.[y1..y2, x1..x2] <- mask

        let copy = (Matrix.cloneO {values = array2d}).values
        array2d |> Array2D.iteri (fun y x v -> if v = 1 then sub copy mask (x, y) (maskCenterX, maskCenterY))

        copy

Наращивание реализовано в функции upbuilding

let upbuilding array2d mask (maskCenterX, maskCenterY) =

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

Разберем новые операторы >. и <.

let (>.) value checkValue = if value > checkValue then checkValue else value
let (<.) value checkValue = if value < checkValue then checkValue else value

Они сравнивают два числа
value — рабочее число
checkValue — пограничное число (с которым сравнивается рабочее число)

Если value больше (или меньше) checkValue, то оператор возвращает checkValue, иначе — value. Это напоминает оператор объединения по nil языка Swift (??), только он позволяет задать конкретное пограничное число, а также можно задать направлене сравнения (больше или меньше).
Например, следующее выражение

2 <. 3 // Вернет 3

Или такое выражение

3 <. 2 // Вернет 3

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

Следующая вспомогательная функция subArrayOfMaks

let subArrayOfMaks array2d mask (centerX, centerY) (maskCenterX, maskCenterY) =

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

array2d — Исходная двумерная матрица (бинарное изображение)
mask — Структурирующий элемент (маска) в виде двумерной матрицы
(centerX, centerY) — Координата пикселя в бинарном изображении, который нуждается в текущей обработки
(maskCenterX, maskCenterY) — Начало координат в структурирующем элементе. Обычно в качестве начала координат задают центр маски, но вы можете задать любую другую корректную координату. Следует понимать, что выбор центра координат маски может влиять на конечный результат.

Функция наращивания upbuilding на вход принимает следующие параметры

array2d — Бинарное изображение
mask — Структурирующий элемент
(maskCenterX, maskCenterY) — Начало координат структурирующего элемента

let upbuilding array2d mask (maskCenterX, maskCenterY) =

        // вспомогательная внутренняя функция для получения субматрицы, относительно обрабатываемого пикселя бинарного изображения и заданной маски
        let sub source mask (x, y) (maskCenterX, maskCenterY) =
            // получаем граничные координаты субматрицы
            let (x1, x2, y1, y2) = subArrayOfMaks source mask (x, y) (maskCenterX, maskCenterY)
            // наполняем субматрицу заданной маской
            source.[y1..y2, x1..x2] <- mask

        // создаем копию бинарного изображение, заполненную нулями
        let copy = (Matrix.cloneO {values = array2d}).values
        // заполняем копию алгоритмом наращивания оригинального изображение
        array2d |> Array2D.iteri (fun y x v -> if v = 1 then sub copy mask (x, y) (maskCenterX, maskCenterY))

        // возвращаем копию как результат
        copy

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

Эрозия бинарного изображения B структурирующим элементом S обозначается B ⊖ S и задается выражением

$B \ominus S = { b | b + s \in B ∀s \in S }$

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

Пример выполнение операции эрозии показан ниже (относительно исходного изображения, представленого в Figure 1).

Figure 4: Эрозия B ⊖ S

$\begin{matrix} 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&*&*&0&0\\ 0&0&0&0&*&*&0&0\\ 0&0&0&0&*&*&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ \end{matrix}$

Реализация F#. Здесь используется материал, рассмотренный работы с матрица на F#
Функциональное программирование на примере работы с матрицами из теории линейной алгебры

// вспомогательная функция для проверки, содаржит ли оригинальная матрица array2d структурирующий элемент mask?
let contains array2d mask =
        if Matrix.isEquallySized (Matrix.ofArray2D array2d) (Matrix.ofArray2D mask) then
            not (array2d
                 |> Array2D.mapi (fun x y v -> if mask.[x, y ] = 0 then 1
                                               elif mask.[x, y] = v then 1
                                               else 0)
                 |> Seq.cast<int>
                 |> Seq.contains 0)
        else false

// функция, реализующая алгоритм эрозии
let erosion array2d mask (maskCenterX, maskCenterY) =

        let sub source (dest: int [,]) mask (x, y) (maskCenterX, maskCenterY) =
            let (x1, x2, y1, y2) = subArrayOfMaks source mask (x, y) (maskCenterX, maskCenterY)
            let subArray2d = source.[y1..y2, x1..x2]

            if contains subArray2d mask then
                dest.[y, x] <- 1

        let copy = (Matrix.cloneO {values = array2d}).values
        array2d |> Array2D.iteri (fun y x v -> if v = 1 then sub array2d copy mask (x, y) (maskCenterX, maskCenterY))

        copy

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

Следующие две операции являются синтезом уже рассмотренных.

Замыкание бинарного изображения B структурирующим элементом S обозначается B • S и задается выражением

$B • S = (B \oplus S) \ominus S$

Реализация F#

let closure array2d mask (centerX, centerY) = 
        let step1 = upbuilding array2d mask (centerX, centerY)
        let step2 = erosion step1 mask (centerX, centerY)

        step2

Обработанное изображение выглядит так
Figure 5: Замыкание B • S

$\begin{matrix} 0&0&0&0&0&0&0&0\\ 0&*&*&*&*&*&*&0\\ 0&0&*&*&*&*&*&0\\ 0&0&*&*&*&*&*&0\\ 0&0&*&*&*&*&*&0\\ 0&0&*&*&*&*&*&0\\ 0&0&*&*&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ \end{matrix}$

Размыкание бинарного изображения B структурирующим элементом S обозначается B ◦ S и задается выражением

$B ◦ S= (B \ominus S) \oplus S$

Реализация F#

let opening array2d mask (centerX, centerY) = 
        let step1 = erosion array2d mask (centerX, centerY)
        let step2 = upbuilding step1 mask (centerX, centerY)

        step2

Figure 6: Размыкание B ◦ S

$\begin{matrix} 0&0&0&0&0&0&0&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&*&*&*&*&0\\ 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ \end{matrix}$

Пример использования эрозии для поиска связанных структур в бинарном изображении

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

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

Figure 7: Оригинальное бинарное изображение для операции условного наращивания

$\begin{matrix} 0&0&0&0&0&0&0&0&0&0\\ 0&0&*&0&0&0&0&0&0&0\\ 0&*&*&*&0&0&0&*&0&0\\ 0&0&0&0&0&0&0&0&*&0\\ 0&0&0&0&*&0&0&0&0&0\\ 0&0&0&0&*&*&0&0&*&0\\ 0&0&0&*&*&*&*&0&*&0\\ 0&0&0&*&0&0&0&0&*&0\\ 0&0&0&0&0&0&0&*&*&0\\ 0&0&0&0&0&0&0&0&0&0\\ \end{matrix}$

Figure 8: Структурируюший элемент

$\begin{matrix} *\\ *\\ *\\ \end{matrix}$

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

Figure 9: Результирующее изображение

$\begin{matrix} 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&*&0&0&0&0&0\\ 0&0&0&0&*&*&0&0&*&0\\ 0&0&0&*&*&*&*&0&*&0\\ 0&0&0&*&0&0&0&0&*&0\\ 0&0&0&0&0&0&0&*&*&0\\ 0&0&0&0&0&0&0&0&0&0\\ \end{matrix}$

Реализация такого поиска на F# приведена ниже

// задаем оригинальное бинарное изображене
let origin = array2D [[0;0;0;0;0;0;0;0;0;0]
                      [0;0;1;0;0;0;0;0;0;0]
                      [0;1;1;1;0;0;0;1;0;0]
                      [0;0;0;0;0;0;0;0;1;0]
                      [0;0;0;0;1;0;0;0;0;0]
                      [0;0;0;0;1;1;0;0;1;0]
                      [0;0;0;1;1;1;1;0;1;0]
                      [0;0;0;1;0;0;0;0;1;0]
                      [0;0;0;0;0;0;0;1;1;0]
                      [0;0;0;0;0;0;0;0;0;0]]

// задаем матрицу
let mask = array2D [[1]
                    [1]
                    [1]]

// ищем связные компоненты и маркируем их
let marked = Algorithms.markingOfConnectedComponents origin
// используем эрозию
let erosion = Algorithms.erosion origin mask (0, 1)

printfn "origin = \n %A" origin
printfn "marked = \n %A" marked
printfn "erosion =\n %A" erosion

// создаем множество для хранения найденных меток связных компонентов. используется для уменьшения итераций обработки изображения (мемоизация)
let set = erosion |> Seq.cast<int> |> Seq.filter ((<>)0) |> Set.ofSeq

// создаем пустую копию исходного изображения, в которую будем складывать результат
let copy = (Matrix.cloneO {values = marked}).values

// обрабатываем эрозийное изображение с целью восстановления цельных связанных компонентов
erosion
|> Array2D.iteri (fun row col v -> if v = 1 then
                            // получаем метку в найденном участке
                            let label = marked.[row, col]
                            if not (set.Contains label) then
                                // перебираем промаркированное изображение и если метки совпадают, то по полученным координатам заполняем копию единичками
                                marked |> Array2D.iteri (fun row1 col1 v1 -> if v1 = label then copy.[row1, col1] <- 1)
                         )

printfn "copy =\n %A" copy

Результаты работы алгоритма ниже

origin = 
 [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 1; 0; 0; 0; 0; 0; 0; 0]
 [0; 1; 1; 1; 0; 0; 0; 1; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 1; 0]
 [0; 0; 0; 0; 1; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 1; 1; 0; 0; 1; 0]
 [0; 0; 0; 1; 1; 1; 1; 0; 1; 0]
 [0; 0; 0; 1; 0; 0; 0; 0; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 1; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]
marked = 
 [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 1; 0; 0; 0; 0; 0; 0; 0]
 [0; 1; 1; 1; 0; 0; 0; 3; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 4; 0]
 [0; 0; 0; 0; 5; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 5; 5; 0; 0; 6; 0]
 [0; 0; 0; 5; 5; 5; 5; 0; 6; 0]
 [0; 0; 0; 5; 0; 0; 0; 0; 6; 0]
 [0; 0; 0; 0; 0; 0; 0; 6; 6; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]
erosion =
 [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 1; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]
copy =
 [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 1; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 1; 1; 0; 0; 1; 0]
 [0; 0; 0; 1; 1; 1; 1; 0; 1; 0]
 [0; 0; 0; 1; 0; 0; 0; 0; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 1; 1; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]

Объединение, пересечение, дополнение, разность

Рассмотрим вышеупомянутые операции и их реализацию на F#

Figure 10: Исходные матрицы

$A= \begin{matrix} 0&*&0\\ 0&*&*\\ 0&*&0\\ \end{matrix} , B= \begin{matrix} 0&0&0\\ *&*&*\\ 0&0&0\\ \end{matrix}$

Объединение двух матриц A и B

$A \cup B$

Figure 11: Объединение

$\begin{matrix} 0&*&0\\ *&*&*\\ 0&*&0\\ \end{matrix}$

Пересечение двух матриц A и B

$A \cap B$

Figure 12: Пересечение

$\begin{matrix} 0&0&0\\ 0&*&*\\ 0&0&0\\ \end{matrix}$

Дополнение матрицы A

$display$

Figure 13: Дополнение

$\begin{matrix} *&0&*\\ *&0&0\\ *&0&*\\ \end{matrix}$

Разность двух матриц A и B

$display$

Figure 13: Разность

$\begin{matrix} 0&*&0\\ 0&0&0\\ 0&*&0\\ \end{matrix}$

Реализация на F#

// объединение
let union array2d1 array2d2 = 
        if Matrix.isEquallyArraySized array2d1 array2d2 then
            array2d1 |> Array2D.mapi (fun x y v -> v ||| array2d2.[x, y])
        else failwith "array2d1 is not equal to array2d2"

// пересечение
let intersection array2d1 array2d2 = 
        if Matrix.isEquallyArraySized array2d1 array2d2 then
            array2d1 |> Array2D.mapi (fun x y v -> v &&& array2d2.[x, y])
        else failwith "array2d1 is not equal to array2d2"

// дополнение (по своей сути есть инверсия)
let inverse array2d =
        array2d |> Array2D.mapi (fun _ _ v -> v ^^^ 1)

let complement array2d = 
        inverse array2d

// разность
let difference array2d1 array2d2 = 
        if Matrix.isEquallyArraySized array2d1 array2d2 then
            array2d1 |> Array2D.mapi (fun x y v -> if v <> array2d2.[x, y] then v else 0)
        else failwith "array2d1 is not equal to array2d2"

Выделение периметра связанного компонента

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

Figure 14: Исходное бинарное изображения для выделения периметра

$\begin{matrix} 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&*&*&*&*&0&0\\ 0&*&*&*&*&*&0&0\\ 0&*&*&*&*&*&*&1\\ 0&*&*&*&*&*&*&*\\ 0&0&*&*&*&*&*&*\\ 0&0&0&0&0&0&0&0 \end{matrix}$

Figure 15: Структурирующий элемент

$\begin{matrix} 0&*&0\\ *&*&*\\ 0&*&0\\ \end{matrix}$

Реализация на F#

let borderAllocation array2d mask (maskCenterX, maskCenterY) = 
        let e = erosion array2d mask (maskCenterX, maskCenterY)
        let border = difference array2d e
        border

В результате обработки получаем следующее изображение

$\begin{matrix} 0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0\\ 0&0&*&*&*&*&0&0\\ 0&*&0&0&0&*&0&0\\ 0&*&0&0&0&0&*&1\\ 0&*&0&0&0&0&0&*\\ 0&0&*&*&*&*&*&*\\ 0&0&0&0&0&0&0&0 \end{matrix}$

Нахождение самого круглого объекта в бинарном изображении

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

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

$\begin{matrix} 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&*&*&*&*&0&0\\ *&*&*&*&0&0&0&0&0&*&*&*&*&*&*&0\\ *&*&*&*&0&0&0&0&*&*&*&*&*&*&*&*\\ *&*&*&*&0&0&0&0&*&*&*&*&*&*&*&*\\ *&*&*&*&0&0&0&0&*&*&*&*&*&*&*&*\\ *&*&*&*&0&0&0&0&0&*&*&*&*&*&*&0\\ *&*&*&*&0&0&0&0&0&0&*&*&*&*&*&0\\ *&*&*&*&0&0&0&0&0&0&0&0&0&0&0&0\\ *&*&*&*&0&0&0&0&0&0&0&0&0&0&0&0\\ *&*&*&*&0&0&*&*&*&0&0&0&0&0&0&0\\ *&*&*&*&0&0&*&*&*&0&0&0&0&0&0&0\\ *&*&*&*&0&0&*&*&*&0&0&0&0&0&0&0\\ *&*&*&*&0&0&0&0&0&0&0&0&0&0&0&0 \end{matrix}$

F#

// исходная матрица с несколькими изображениями
    let multyBinImage = array2D [[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]
                                 [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]
                                 [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]
                                 [0;0;0;0;0;0;0;0;0;0;1;1;1;1;0;0]
                                 [1;1;1;1;0;0;0;0;0;1;1;1;1;1;1;0]
                                 [1;1;1;1;0;0;0;0;1;1;1;1;1;1;1;1]
                                 [1;1;1;1;0;0;0;0;1;1;1;1;1;1;1;1]
                                 [1;1;1;1;0;0;0;0;1;1;1;1;1;1;1;1]
                                 [1;1;1;1;0;0;0;0;0;1;1;1;1;1;1;0]
                                 [1;1;1;1;0;0;0;0;0;0;1;1;1;1;0;0]
                                 [1;1;1;1;0;0;0;0;0;0;0;0;0;0;0;0]
                                 [1;1;1;1;0;0;0;0;0;0;0;0;0;0;0;0]
                                 [1;1;1;1;0;0;1;1;1;0;0;0;0;0;0;0]
                                 [1;1;1;1;0;0;1;1;1;0;0;0;0;0;0;0]
                                 [1;1;1;1;0;0;1;1;1;0;0;0;0;0;0;0]
                                 [1;1;1;1;0;0;0;0;0;0;0;0;0;0;0;0]]

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

// маркируем связанные пиксели
let markedMultyImages = Algorithms.markingOfConnectedComponents mulyBinImage
printfn "marked = \n %A" markedMultyImages

Результат

marked = 
 [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 1; 1; 1; 1; 0; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1]
 [2; 2; 2; 2; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1]
 [2; 2; 2; 2; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 0; 1; 1; 1; 1; 0; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
 [2; 2; 2; 2; 0; 0; 5; 5; 5; 0; 0; 0; 0; 0; 0; 0]
 [2; 2; 2; 2; 0; 0; 5; 5; 5; 0; 0; 0; 0; 0; 0; 0]
 [2; 2; 2; 2; 0; 0; 5; 5; 5; 0; 0; 0; 0; 0; 0; 0]
 [2; 2; 2; 2; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]

Алгоритм нашел три связанные группы и промаркировал их тремя метками (1, 2, 5)

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

// группируем по меткам

// groups key - метка, value - массив координат всех пикселей, помеченных этой меткой
let groups = Dictionary<int, (int*int) list>()
markedMultyImages
|> Array2D.iteri (fun row col v -> if v <> 0 then
                                       if not (groups.ContainsKey(v)) then
                                            groups.Add(v, [row, col])
                                       else
                                            groups.[v] <- (row, col)::groups.[v])

printfn "groups = \n %A" groups

Результат

groups = 
 seq
  [[1, [(9, 13); (9, 12); (9, 11); ... ]];
   [2, [(15, 3); (15, 2); (15, 1); ... ]];
   [5, [(14, 8); (14, 7); (14, 6); ... ]]]

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

Площадь фигуры (это сумма всех пикселей)
Длину периметра (сумма пикселей по периметру фигуры)
Большая формула (будет позже)

let mapAreas = groups
               |> Seq.map (|KeyValue|)
               |> Map.ofSeq

// вычисляем площади фигур
let areas = mapAreas
            |> Map.map (fun k v -> v |> List.length)
printfn "areas =\n %A" areas

areas — словарь, у которого ключи — метки фигур, а значения — вычисленные площади
Результат

areas =
 map [(1, 44); (2, 48); (5, 9)]

Видно, что фигура, помеченная меткой 1 имеет площать 44 пикселя, фигура с меткой 2 — 48 пикселей, фигура с меткой 5 — 9 пикселей.

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

// вырезаем отдельные фигуры помеченные метками
    let subarrays = mapAreas
                    |> Map.map (fun k v ->
                                    let minRow = v |> List.map (fst) |> List.min
                                    let maxRow = v |> List.map (fst) |> List.max
                                    let minCol = v |> List.map (snd) |> List.min
                                    let maxCol = v |> List.map (snd) |> List.max

                                    let dimRow = maxRow - minRow + 2 + 1 // add two zeros on left and right sides
                                    let dimCol = maxCol - minCol + 2 + 1

                                    let array2d = Array2D.create dimRow dimCol 0
                                    v |> List.iter (fun elem -> (array2d.[(fst elem) - minRow + 1, (snd elem) - minCol + 1] <- multyBinImage.[fst elem, snd elem]))

                                    array2d
                               )

    printfn "subarrays =\n %A" subarrays

Результат

subarrays =
 map
  [(1, [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
        [0; 0; 0; 1; 1; 1; 1; 0; 0; 0]
        [0; 0; 1; 1; 1; 1; 1; 1; 0; 0]
        [0; 1; 1; 1; 1; 1; 1; 1; 1; 0]
        [0; 1; 1; 1; 1; 1; 1; 1; 1; 0]
        [0; 1; 1; 1; 1; 1; 1; 1; 1; 0]
        [0; 0; 1; 1; 1; 1; 1; 1; 0; 0]
        [0; 0; 0; 1; 1; 1; 1; 0; 0; 0]
        [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]);

   (2, [[0; 0; 0; 0; 0; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 1; 1; 1; 1; 0]
         [0; 0; 0; 0; 0; 0]]);

   (5, [[0; 0; 0; 0; 0]
        [0; 1; 1; 1; 0]
        [0; 1; 1; 1; 0]
        [0; 1; 1; 1; 0]
        [0; 0; 0; 0; 0]])]

Теперь мы готовы выделить периметры каждой фигуры и найти их длины. Длина периметра — это сумма ненулевых пикселей полученной фигуры.

// выделяем периметры фигур и вычисляем их длины
    let erosionMaks = array2D [[0;1;0]
                               [1;1;1]
                               [0;1;0]]

    let borders = subarrays
                  |> Map.map (fun label array2d -> Algorithms.borderAllocation array2d erosionMaks (1, 1))
    let bordersLength = borders |> Map.map (fun label elem -> elem |> Seq.cast<int> |> Seq.filter ((<>)0) |> Seq.length)

    printfn "borders =\n %A" borders
    printfn "borders length =\n %A" bordersLength

Результат

borders =
 map
  [(1, [[0; 0; 0; 0; 0; 0; 0; 0; 0; 0]
        [0; 0; 0; 1; 1; 1; 1; 0; 0; 0]
        [0; 0; 1; 0; 0; 0; 0; 1; 0; 0]
        [0; 1; 0; 0; 0; 0; 0; 0; 1; 0]
        [0; 1; 0; 0; 0; 0; 0; 0; 1; 0]
        [0; 1; 0; 0; 0; 0; 0; 0; 1; 0]
        [0; 0; 1; 0; 0; 0; 0; 1; 0; 0]
        [0; 0; 0; 1; 1; 1; 1; 0; 0; 0]
        [0; 0; 0; 0; 0; 0; 0; 0; 0; 0]]);

(2, [[0; 0; 0; 0; 0; 0]
      [0; 1; 1; 1; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 0; 0; 1; 0]
      [0; 1; 1; 1; 1; 0]
      [0; 0; 0; 0; 0; 0]]);

   (5, [[0; 0; 0; 0; 0]
        [0; 1; 1; 1; 0]
        [0; 1; 0; 1; 0]
        [0; 1; 1; 1; 0]
        [0; 0; 0; 0; 0]])]

borders length =
 map [(1, 18); (2, 28); (5, 8)]

Длина периметра фигуры 1 равна 18 пикселей, фигура 2 — 28 пикселей и фигуры 5 — 8 пикселей.

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

// ищем центры каждой фигуры
    let centers = borders
                   |> Map.map (fun k v ->
                                 let centerRow = v |> Array2D.mapi (fun r _ _ -> float r) |> Seq.cast<float> |> Seq.average // list.average requires float
                                 let centerCol = v |> Array2D.mapi (fun _ c _ -> float c) |> Seq.cast<float> |> Seq.average
                                 (int centerCol, int centerRow)
                              )
    printfn "centers =\n %A" centers

Результат

centers =
 map [(1, (4, 4)); (2, (2, 6)); (5, (2, 2))]

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

$C=\frac{μR}{σR}$

μR — Среднее радиальное расстояние
σR — Среднеквадратичное отклонение радиального расстояния

Среднее радиальное расстояние вычисляется по следующей формуле

$μR = \frac{1}{K}\sum\limits_{k=1}^{K-1}||(r_k, c_k) - (\overline{r}, \overline{c})||$

(r_k, c_k) — координата центра фигуры
(r, c) — текущаяя координата

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

// вычисляем расстояние между точками
    let diff point1 point2 = 
        let square x = x * x
        let s1 = square ((fst point1) - (fst point2))
        let s2 = square ((snd point1) - (snd point2))
        sqrt(double (s1 + s2))

Вычисляем коэффициенты μR каждой фигуры.

// получаем словарь, где key - метка, value - вычисленный коэффициент (double)
let mr = borders
         |> Map.map (fun label v ->
                        let aver = v
                                    |> Array2D.mapi (fun r c v -> diff (r, c) centers.[label])
                                    |> Seq.cast<double>
                                    |> Seq.sum
                        aver / (double v.Length)
                 )

Коэффициент σR вычисляется по следующей формуле

$σR = (\frac{1}{K}\sum\limits_{k=1}^{K-1}[||(r_k, c_k) - (\overline{r}, \overline{c})|| - μR]^2)^{1/2}$

let qr = borders
     |> Map.map (fun label v ->
                    let aver = v
                            |> Array2D.mapi (fun r c v ->
                                                let t1 = diff (r, c) centers.[label]
                                                let t2 = (t1 - mr.[label]) * (t1 - mr.[label])
                                                t2
                                            )
                            |> Seq.cast<double>
                            |> Seq.sum
                            |> sqrt
                aver / (double v.Length)
         )

// вычисляем округлость. чем выше коэффициент, тем круглее фигура
    let roundness = mr |> Map.map (fun label v -> v / qr.[label])
    printfn "roundness =\n %A" roundness

Результат

roundness =
 map [(1, 25.06054796); (2, 20.37134197); (5, 13.43281947)]

Как видно, фигура с меткой 1 имеет самый высокий коэффициент округлости. Этой фигурой является овал и это правильный ответ.

Заключение

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

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

F#, морфология бинарных изображений