Подсчет объектов на бинарном изображении. Часть 1

    Аннотация


    imageРаз, два, три, четыре, пять. Будем в прятки мы играть. В статье рассказывается про алгоритм разметки (или подсчета) объектов на бинарном изображении и о том, как без дополнительного прохода вычисляются (в еще неопубликованной части 2) геометрические характеристики этих объектов. Алгоритмы подобного типа часто используются при распознавании образов на бинарном препарате и показывают свою вычислительную эффективность.
    В завершении статьи, читателям предлагается интересная задачка, грамотное решение которой существует и необходимо, при практической реализации алгоритма. Приводится исходный код, но в отличии от предыдущих моих постов, он выполнен не на языке MatLab а в абсолютно свободной, не менее мощной среде SciLab.

    Введение


    В практических приложениях, рассчитанных на использование в реальном времени распознавание образов часто выполняют в виде следующей последовательности:
    Получение изображения (1) -> Адаптивная бинаризация (2) -> Серия морфологических операций (3) -> Разметка объектов (4) -> Заполнение пространства признаков объектов (5).
    image

    Рисунок 1 — Упрощенная схема распознавания образов на изображении. Схема исключает этапы разделения пространства признаков и принятия решения.

    На рисунке 2 изображены фрагменты изображения после бинаризации и серии морфологических операций (рис.2а) и собственно результат разметки (рис.2б).
    image

    Рисунок 2 — а) — бинаризованное изображение, б) — изображение с размеченными объектами.

    В статье мы будем не только размечать объекты, но подсчитывать геометрические характеристики такие как площадь и периметр. Их соотношения могут описывать компактность объекта и фактор формы Малиновского, а в комбинации с более сложными геометрическими признаками — вытянутость и особенности контура.

    По научному, то что мы хотим получить, называется — выделение связных областей. Существует множество решений данной задачи — основанные на шаблонах [1], некоторые рекурсивные реализации (например листинг 1 из [2]) и др.

    void Labeling(BIT* img[], int* labels[])
    {
     int L = 1;
     for(int y = 0; y < H; ++y)
      for(int x = 0; x < W; ++x)
       Fill(img, labels, x, y, L++);
    }
    void Fill(BIT* img[], int* labels[], int x, int y, int L)
    {
      if( (labels[x][y] == 0) & (img[x][y] == 1) )
      {
        labels[x][y] = L;
        if( x > 0 )
            Fill(img, labels, x – 1, y, L);
        if( x < W — 1 )
            Fill(img, labels, x + 1, y, L);
        if( y > 0 )
            Fill(img, labels, x, y — 1, L);
        if( y < H — 1 )
            Fill(img, labels, x, y + 1, L);
      }
    }

    * This source code was highlighted with Source Code Highlighter.
    Листинг 1 — рекурсивная реализация алгоритма разметки связных областей.

    Это все нам не подходит.
    Область пикселей будем называть связной, если для каждого пикселя из этой области существует сосед из этой же области. Про виды связности я приводил иллюстрацию в этом топике(Рис.1,2). Наш алгоритм будет четырех связным, хотя его без особых умственных усилий можно переделать и для восьми связного варианта.

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


    Идея данного алгоритма основана на использовании уголка — ABC-маски, которая показана на рисунке 3.
    image

    Рисунок 3 — ABC маска и направление последовательного сканирования изображения.

    Проход по изображению данной маской осуществляется слева-направо и сверху-вниз. Считается, что за границей изображения объектов нет, поэтому, если туда попадают B или C — это требует дополнительной проверки при сканировании. На рисунке 4 изображены 5 возможных позиций маски на изображении.
    Рассмотрим их.
    image

    Рисунок 4 — Пять возможных позиций ABC-маски
    1. Позиция под номером 0, когда не размечены все три компонента маски — в этом случае мы просто пропускаем пиксель.
    2. Позиция под номером 1, когда помечен только элемент A — в этом случае мы говорим о создании нового объекта — новый номер.
    3. Позиция под номером 2, когда помечен элемент элемент B — в этом случае мы помечаем текущий пиксель A меткой, расположенной в B.
    4. Позиция под номером 3, когда помечен элемент элемент С — в этом случае мы помечаем текущий пиксель A меткой, расположенной в С.
    5. Позиция под номером 4, тогда мы говорим о том, что метки (номера объектов) B и C связаны — то есть эквивалентны и пиксель A может быть помечен либо как B либо как C. В некоторых реализация составляют граф эквивалентности таких меток, с последующим его разбором, однако на мой взгляд в это нет необходимости. Мы будем делать так — в том случае, если B не равно C то перенумеруем все уже обработанные пиксели помеченные как С в метку B. Но об этом в самом конце.

    Ближе к реализации.
    В начале создадим тестовые данные аналогичные изображенным на рисунке 2, а именно матрицу Image состоящую из нулей и единичек.

    Image = [0 0 0 0 0 0 0 1 0 0 0;
         0 0 0 1 1 0 0 1 0 0 0;
         0 0 0 1 1 0 0 1 1 0 0;
         0 1 1 1 1 0 0 1 1 1 1;
         0 1 1 1 1 0 0 0 0 0 0;
         0 1 1 1 0 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 0 0;
         0 0 0 0 0 0 1 1 1 0 0;
         0 1 1 0 0 1 1 1 1 0 0;
         0 1 1 0 0 1 1 1 1 1 1;
         0 1 1 1 0 1 1 0 1 1 0;
         0 0 0 0 0 0 0 0 0 0 0]
    Matplot(Image*255) // Посмотрим на нашу матрицу как на картинке
    [m,n]=size(Image); // Узнаем горизонтальный и вертикальные размеры матрицы
    km = 0; kn = 0; // Они нам еще пригодятся
    cur = 1; // Переменная для подсчета объектов


    * This source code was highlighted with Source Code Highlighter.
    Листинг 2 — инициализация исходных данных.

    А затем пройдемся по изображению выполняя серию простых и очевидных проверок. ABC-маска, которой мы проходим по картинке, проиллюстрирована на рисунке 3. Возможные комбинации, которые проверяются этой маской изображены на рисунке 4.

    // Цикл по пикселям изображения
    for i = 1:1:m
      for j = 1:1:n
        kn = j — 1;
        if kn <= 0 then
          kn = 1;
          B = 0;
        else
          B = Image(i,kn); // Смотри рисунок 3 в статье
        end
        km = i — 1;
        if km <= 0 then
          km = 1;
          C = 0;
        else
          C = Image(km,j); // Смотри рисунок 3 в статье
        end
        A = Image(i,j); // Смотри рисунок 3 в статье
        if A == 0 then // Если в текущем пикселе пусто — то ничего не делаем
        elseif B == 0 & C == 0 then // Если вокруг нашего пикселя пусто а он не пуст — то это повод подумать о новом объекте
            cur = cur + 1;
            Image(i,j) = cur;
        elseif B ~=0 & C == 0 then
            Image(i,j) = B;
        elseif B == 0 & C ~= 0 then
            Image(i,j) = C;
        elseif B ~= 0 & C ~= 0 then       
            if B == C then
              Image(i,j) = B;
            else
              Image(i,j) = B;
              Image(Image == C) = B;
            end
        end
      end
    end


    * This source code was highlighted with Source Code Highlighter.
    Листинг 3 — последовательное сканирование пикселей изображения и разметка связных областей.

    Результатом выполнения этого скрипта будет размеченная матрица:
    image

    Рисунок 5 — Результат выполнения листингов 2 и 3 — объектам присвоены уникальные номера.

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

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

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

    Спасибо всем прочитавшим, жду ваших отзывов.

    Ссылки и использованные источники


    [1] Алгоритмы / Подсчёт объектов на изображении ссылка (видимо эта статья когда-то была на Хабре, но найти не удалось)
    [2] Лекция: Анализ информации содержащейся в изображении ссылка на ppt файл
    Share post
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 44

      +1
      Это все нам не подходит.
      почему?
        0
        Может быть большая глубина рекурсии?
          +2
          можно заменить на поиск в ширину.
          0
          Ответить на это в комментарии — сложно, слишком спорная тема. У меня все субъективно: мне не нравится рекурсия, особенно, в тех случаях, когда ее можно заменить изящным сканирующим алгоритмом. Если у вас есть какие-то доводы за другие подходы — прошу изложить.
          0
          Когда мне это когда-то давно нужно было решить такую задачу на курсовой, я решил её за 5 минут путём заполнения (у вас это Fill). А вам почему-то не подходит. Почему?
            0
            Я не знаю, как можно реализовать подсчет периметра в рекурсивной версии, когда в сканирующей версии это делается очень просто вместе с остальными геометрическими характеристиками — тема второй части.
              +1
              В первом приближении можно использовать периметр bounding box бинарного объекта :)

              Кстати, а ваш метод с масками работает для объектов, содержащих внутренние незаполненные области?
                0
                bounding box — тогда отличить квадрат от круга не получится, а фактор формы с точно вычисленным периметром и площадью — это позволит сделать. но об этом позже.
                Да конечно, если внутри есть незаполненная область — то все так же будет работать. Для уверенности — перепроверил сейчас. Тут и число Эйлера можно вычислить вроде.
            +1
            Разметка связных областей делается с помощью bwlabel (есть и для Scilab).

            Вычисление характеристик (длина, площадь) с помощью встроенных функций (если есть) или в матричном виде, аккумулируя результаты с помощью vl_binsum, histc или подобных.

            Всяко будет красивее и быстрее, чем писать циклы в Matlab/Scilab.

            Вот например на днях нужно было посчитать длину общей границы между всеми парами областей: (чуть более сложная задача, чем посчитать длины границ объектов):
            function adj = pairwiseBoundaryLengths(map)
                % map содержит разметку изображения по областям
            
                h = size(map, 1);
                w = size(map, 2);
                
                n = max(map(:));
                
                % вертикальные и горизонтальные границы
                % фокус - записываем для каждой границы номера граничащих областей
                % но разворачиваем пару чисел в одномерный индекс, чтобы проще было потом суммировать
                vertical = (map(1 : h - 1, :) - 1) * n + map(2 : h, :);
                horizontal = (map(:, 1 : w - 1) - 1) * n + map(:, 2 : w);
                
                adj = zeros(n);
                adj(:) = vl_binsum(adj(:), 1, vertical(:));
                adj(:) = vl_binsum(adj(:), 1, horizontal(:));
                adj = adj + adj';
                adj = adj .* (ones(n) - eye(n));
            end
            

            Не понадобился ни один цикл.
              0
              А ещё, у вас операция слияния областей «Image(Image == C) = B;» может оказаться дороговатой.

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

              Потому и предоставляется готовая функция bwlabel, написанная на C.
                0
                В пункте 5-ом, описания позиции маски, говорится как раз о перенумерации Image(Image == C) = B. При реальной рализации это эфективно делать с использованием указатлей — см. задачку в самом конце статьи.
                  0
                  > матрица указателей на специальную структуру
                  «специальная структура» — это система непересекающихся множеств, или что-то другое?
                    0
                    да, массив не пересекающихся множеств индексов, если два индекса в одном множестве — области эквивалентны. И при работе алгоритма выставляются указатели, а множества так же могут объединятся.
                      0
                      У меня получилось сделать однопроходный алгоритм (C++) без использования структур/множеств. Храню массив int* вместо int, и в случае, когда текущий пиксель объединяет 2 объекта (последний else у вас в коде), делаю так:

                      int **m_labels;
                      
                      ...
                      
                      } else {
                          if (*m_labels[bIndex] <= *m_labels[cIndex])
                              *m_labels[cIndex] = *m_labels[bIndex];
                          else
                              *m_labels[bIndex] = *m_labels[cIndex];
                          m_labels[index] = m_labels[bIndex];
                      }
                      

                      т.е. изменяю значение указателя с бОльим номером объекта на значение объекта с меньшим номером, с которым происходит объединение. Таким образом остается один (общий) объект с меньшим номером. int* в качестве ячейки позволяет изменить номер объекта для всех пикселей этого объекта без дополнительного прохода по массиву рисунка.

                      Конечно, для такой реализации необходимо иметь матрицу указателей, но в моем случае массив рисунка для этих нужд я использовать не мог.
                        0
                        Пока еще не понял будет ли это работать для всех коллизий, однако идея очень интересная.
                        Вот к примеру пусть разметка такова, что возникают последовательно эквивалентности:
                        2 эквивалентно 3, затем: 3 эквивалентно 5.
                        1. По адресу в указателе на 2 ставится 3.
                        2. По адресу указывающему на на 3 ставится 5.
                        Таким образом коллизия не разрешена и перенумерация выполнена не полностью: осталось и 3 и 5. Когда должно остаться только 5.
                        Или я где-то заблуждаюсь?

                          0
                          1. Соединяются объекты 2 и 3. Объект 3 станет объектом 2.
                          2. Соединяются объекты 3 и 5. Объект 5 станет объектом 2, т.к. объект 3 — уже имеет номер 2 (переименован в предыдущем шаге)

                          Код работает, я сейчас использую этот подход в реальном проекте (кстати, спасибо вам за отличную статью, очень своевременно для меня).

                          Ксати, что характерно, если не делать проверку

                          } else {
                              if (*m_labels[bIndex] <= *m_labels[cIndex])
                                  *m_labels[cIndex] = *m_labels[bIndex];
                              else
                                  *m_labels[bIndex] = *m_labels[cIndex];
                              m_labels[index] = m_labels[bIndex];
                          }
                          

                          а делать просто

                          } else {
                              *m_labels[cIndex] = *m_labels[bIndex];
                              m_labels[index] = m_labels[bIndex];
                          }
                          

                          то выйдет криво (будет куча объектов вместо одного).
                            0
                            Хм, действительно, сейчас подумал — все верно. Только не как доказать это не могу для себя. У вас есть какие-то рассуждения? Поделитесь?
                              +1
                              Да, я могу объяснить. Вечером нарисую картинку, как всё происходит.
                              0
                              А если эквивалентность вида: 2 эквивалентно 3, а 3 эквивалентно 1?
                              1. Соединяются объекты 2 и 3. Объект 3 станет объектом 2. Теперь указателей на тройку нет вообще.
                              2. Соединяются объекты 3 и 1. Объект 3 станет объектом 1. А так как указатели на тройку уже потеряны, то эта эквивалентность не обработается получается?

                                +1
                                Да, вы правы. На новых рисунках отчетливо видны необъединенные объекты. Видимо, придется делать двухпроходным алгоритмом со списком множеств эквивалентных объектов
                                  0
                                  Действительно жаль, что не нашлось более изящного способа.
                                  0
                                  Можете это поподробнее объяснить (или картинку добавить если не очень сложно)? Никак понять не могу.
                                  Я сейчас как раз делаю так же с матрицей указателей. Ошибка возникает, а причина до конца неясна.
                                    0
                                    Пусть изображение представляется матрицей, где в каждой ячейке — указатель на номер объекта, которому принадлежит эта ячейка.
                                    Все номера хранятся в отдельном массиве — выделенном предварительно и с некоторым запасом, так как неизвестно сколько объектов на картинке.

                                    Самое простое и нормально работающее:

                                    Возникла коллизия 2 эквивалентно 3, мы не трогаем указатели в матрице, мы в массиве номеров, на втором месте вместо двойки ставим 3.
                                    Возникла коллизия 3 эквивалентно 1, мы просматриваем массив номеров — и везде, где в нем стоит 3, ставим 1.

                                    То есть обобщая эти два тезиса имеем:
                                    При возникновении эквивалентности C1 эквивалентно C2 мы проходим весь массив номеров и там где C1 ставим C2.

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

                                    В эти выходные я опубликую вторую часть статьи, с которой уже порядком затянул — обстоятельства. Там расскажу про улучшение этого способа, но оно на поверхности.
                                      0
                                      Да. Так и делаю, скорость действительно высокая.
                                      А разве возможен случай, когда две эквивалентности (это ведь позиция 4?) ещё и с разными метками возникает подряд?

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

                                          +1
                                          Спасибо за ответы. Жду выходных. :-)

                                          У меня не было курса связанного с обработкой изображений, а в ВКР этим занимаюсь, вот и маюсь от непривычной специфики.
                                        0
                                        Пришел к аналогичной идее после комментария skb7 с методом присвоения ссылок.
                                        Большое спасибо автору и skb7 за идеи и направления.
                          0
                          " длину общей границы между всеми парами областей" — у меня взорвался мозг, если две области не связанные как между ними есть общая граница? Расскажите подробнее пожалуйста.
                            +1
                            Ну у меня просто разбиение было не на черное/белое, а по цветам (связная красная область, связная синяя область и т.п.), то есть разные области могли соприкасаться и не было т.н. фона (т.е. клеток, не принадлежащих никакой области). Это более общий случай. Тогда у каждой пары областей можно рассматривать длину общей границы. Если области не связаны, то длина общей границы равна нулю.

                            А в данном случае, когда есть несвязанные области и фон — посчитать матрично периметр всех областей ещё проще.
                              +1
                              На всякий случай поясню что будет в вашем случае:
                              * найти площади областей — напрямую использовать vl_binsum/histc
                              * найти периметры — обнулить все внутренние пиксели области, после чего задача сведется к нахождению площадей.

                              Матрично (без циклов) удалить внутренние пиксели области несложно.
                                0
                                Действительно, но что-то не соображу, как без циклов удалить внутренние пиксели, мне сразу представился последовательный проход и маска.
                                  +1
                                  Пусть map содержит 0 — там где фон, и числа больше 0, таи где области. Пусть также на границах массива не будет областей (если это не так, то можно дописать к нему нулевые строки/столбцы).

                                  Тогда:
                                  [h, w] = size(map);
                                  
                                  % логический массив, содержащий 1 для внутренних элементов
                                  inner = map(2:h-1, 3:w)~=0 & map(2:h-1,1:w-2)~=0 & map(3:h, 2:w-1)~=0 & map(1:h-2, 2:w-1)~=0;
                                  
                                  map = map(2:h-1, 2:w-1);
                                  map(inner) = 0;
                                  


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

                                  Для конкретных задач можно использовать функции erode или сразу bwperim.
                                    0
                                    Афигеть просто, здорово, выражаю огромную признательность за фишки которые вы рассказали в комментах!
                            –1
                            данную задачу можно решить существенно быстрее и сильно сэкономив память.
                            сначала обходим объект по внешнему контуру (можно сделать за время порядка количества внешних точек, что обычно не больше 4*X*Y, если интересно, могу написать как делается), а потом заполняем соответствующим индексом все внутренние точки (если это нужно).

                            остается только уметь находить начальные точки для обхода, без предварительного квадратичного (X*Y) прохода по изображению это сделать сложно, но можно совместить этот процесс с считыванием изображения.
                              0
                              Нет, это не так. О какой памяти вообще говорится ?, дополнительной памяти не требуется, и экономить ее не получится, тем более существенно.
                              «уметь находить начальные точки для обхода» — в этом проблема.
                              Более того, после того как контур найден, его нужно убрать с изображения, чтобы снова не найти — дополнительные расходы.
                              И еще, ваш вариант совершенно не подходит если изображение поступает последовательно, придется дождаться полного получения картинки, в моем варианте — достаточно получить две строки.
                                0
                                мой комментарий имеет смысл в предположении, что картинка в памяти доступна только для чтения.
                                собственно если это честный набор бит Ч/Б, то ее эффективнее представлять как набор бит, а если нет — то после работы алгоритма она не должна портиться.
                                контуры не надо убирать с изображения, достаточно просто списка пропусков точек для каждого столбца.

                                и да, если картинка поступает последовательно с построением контуров будут сложности.
                                  0
                                  «достаточно просто списка пропусков точек» — действительно, но все равно как-то криво. Так и не понял преимуществ предлагаемого подхода.
                                  0
                                  и + скорость. мы когда занимались подобной задачей пользовали подобный (или в точности такой? уже не припомню) метод, время работы на картинках 3200x2400 измерялась в секундах. с внешними контурами — десятые доли секунды (реально сильно зависело от картинки, от 5 миллисекунд до секунды).
                                    0
                                    Я пробовал реализовывать описанный вами алгоритм и отказался от него.
                                    Этапы: Поиск стартовой позиции, трассировка контуров, определение внутренних точек.
                                    Самый затратный этап был у меня — определение принадлежности точки к полигону — он сделал строго по алгоритму из книги Кормена.
                                    У меня получилось медленнее чем в предлагаемом в статье варианте.
                                    И я не верю, что что можно сделать существенно быстрее.
                                    Хотя все имеет место в своих приложениях.
                                      +1
                                      завтра попробую написать топик про то как мы это делали, и в частности проблема принадлежности точки к контуру была решена.

                                      хотя, ради справедливости, в финальном релизе мы плюнули на оба варианта и реализовали «самый тупой» (алгоритм заливки) на CUDA. со скоростью обработки в 20Гбайт/с результирующее время даже не мерили.
                                        0
                                        И про CUDA тоже напишите пожалуйста, очень интересно. Спасибо за дискуссию.
                                0
                                Спустя 5 лет все-таки задам вопрос, вы говорите, что алгоритм однопроходный и может работать, когда изображение поступает последовательно пиксель за пикселем (надо две строчки). А как в этом случае реализовавать операцию слияния двух объектов, если B != 0, С != 0 и B != C.
                                Image(Image == C) = B?

                                Мне приходит в голову только опять обойти изображени только теперь снизу вверх, но также справа налево.

                                Only users with full accounts can post comments. Log in, please.