Алгоритм поиска изолиний

[Ed_Em]

Вопрос: есть ли в мире открытые реализации алгоритмов поиска изолиний на изображении (на С)?
Сам я набрел на алгоритм шагающих квадратов и попробовал реализовать его, однако в случае резких градиентов интенсивности изображения некоторые изолинии получаются разорванными. Да и алгоритм этот работает последовательно.
Изолинии хранятся у меня в виде связанных списков: для набора изолиний создается массив, каждый член которого - список изофот с одним и тем же уровнем интенсивности; каждая изофота в свою очередь является списком точек - координат ее узлов.
Здесь мне посоветовали распараллелить задачу и разбить ее на два этапа: сначала в каждом подызображении искать изолинии, а затем проходиться по всем изолиниям и сшивать их. Если поступать так, то для N кусков изолиний на M подызображениях придется сделать в среднем еще и $M\cdot(\frac{N}{M})^2\cdot 8$ переборов отрезков. Это, конечно, не так уж и много, но при "сшивании" изолиний придется по сути перегонять довольно большие объемы информации.
Второй способ я, конечно, попробую реализовать, но все-таки хочется узнать: может этот алгоритм уже где-то реализован и мне достаточно просто подключить какую-нибудь библиотеку?

[Circiter]

2Ed_Em
Я так вообще впервые слышу о такой задаче... :) Изолинии это что? Связные кривые, вдоль которых интенсивность пикселей постоянна?

Цитата:

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

Почему они (изолинии) должны быть непрерывными без наложения дополнительных условий на изображение?

Вообще, изолиний у поля их же континуум. :) Во-первых, подразумевается, что сама интенсивность (не поле, а скаляр, e.g. на промежутке $[0;\ 255]$ :) ) с некоторым шагом должна сканироваться, а потом изолинии с одним и тем же уровнем интенсивности как-то должны пространственно разноситься, опять же с шагом (нашли одну, удалили с картинки сумму Минковского этой кривой с небольшим шаром и дальше пошли)...

Вот как, например, должна анализироваться ровно залитая область изображения? Или в ней изолинии не должны детектироваться вовсе?

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

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

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

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

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

Это так, чистая фантазия. Без понимания условия задачи можно много чего нафантазировать, так что извиняйте...

[Ed_Em]

Цитата:

Изолинии это что? Связные кривые, вдоль которых интенсивность пикселей постоянна?

Да, это изофоты.
Изофоты не могут быть разорванными по вполне объективным причинам: представьте себе карту местности, где изолиниями обозначены высоты. Они не всегда замкнуты (если идут за пределы карты, либо если обозначают вершину хребта), но всегда непрерывны. Так же и с изофотами.

Цитата:

Вот как, например, должна анализироваться ровно залитая область изображения? Или в ней изолинии не должны детектироваться вовсе?

Да, на такой области не будет ни одной изолинии.

Цитата:

у меня вообще даже в саму постановку исходной задачи никак въехать не получается

Вкратце объясню: как я уже говорил выше, изофоты (изолинии) - это линии с одинаковой интенсивностью. Т.к. изображение у нас дискретное как по координатам, так и по линиям, задача поиска изофот нетривиальна. Простейшим способом является метод шагающих квадратов: изображение дискретизуется по интенсивности в соответствии с формулой получения изолиний (в простейшем случае равноотстоящих изолиний "новая" интенсивность вычисляется как $I_{new}=\text{\rm floor}(\frac{N\cdot(I_{old}-I_{min})}{I_{max}-I_{min}})$, где N - желаемое количество изолиний для изображения с интенсивностью $I\in[I_{min}, I_{max}]$.
Далее строится специальная маска на основании распределения интенсивностей в соседних пикселях, после чего все пиксели маски сканируются и, в случае обнаружения перепада интенсивности в одной из точек изображения, на основе значения в соответствующих пикселях вспомогательной матрицы происходит поиск соседних пикселей на предмет продолжения изолинии. Для каждого пикселя изображения методом линейной аппроксимации вычисляется положение опорной точки изолинии. В итоге получаем связанный набор точек $M_i(x_i, y_i)$ - изофоту.

[PAV]

Поскольку тема скорее алгоритмическая, чем программерская, то я ее переношу из подфорума "Программирование" в корень CS.

Хорошо, начнем с простейшего подхода: выбирается некоторый порог яркости, по которому изображение бинаризуется. Затем строим границы полученных областей на черно-белом изображении, и считаем их требуемыми изолиниями. Какие недостатки у такого подхода?

[Ed_Em]

PAV в сообщении #489997 писал(а):

Поскольку тема скорее алгоритмическая, чем программерская, то я ее переношу из подфорума "Программирование" в корень CS.

Хорошо, начнем с простейшего подхода: выбирается некоторый порог яркости, по которому изображение бинаризуется. Затем строим границы полученных областей на черно-белом изображении, и считаем их требуемыми изолиниями. Какие недостатки у такого подхода?

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

[PAV]

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

[Ed_Em]

Цитата:

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

И в итоге мы получим метод "шагающих квадратов" :)

Цитата:

А в чем именно видится грубость данного подхода?

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

[Pavia]

Толи я чего-то не понимаю. Толи я не понимаю ничего.

Цитата:

причем через некоторые точки может проходить несколько изолиний)

Точки у которых высота варьируется в природе не встречаются. Или я неправ?

А зачем вы градиент считает?

[Ed_Em]

Цитата:

Точки у которых высота варьируется в природе не встречаются. Или я неправ?

Да нет же. Смотрите: например, мы ищем изолинии с шагом 100, начиная с нуля. Пусть в одной точке изображения интенсивность равна 50, а в соседней - 130. Следовательно, между этими точками проходит одна изолиния с номером 1. Пусть теперь где-нибудь в одной точке интенсивность равна 420, а в соседней - 680. Следовательно, между этими точками проходят две изолинии: с номером 4 и номером 5.
Конечно, после линейной интерполяции узлы будут в разных точках, но через данную точку маски (квадрат 2х2 точки) будет проходить две изолинии.

Цитата:

А зачем вы градиент считает?

Это я вам как пример привел. Естественно, я градиент не считаю (это очень ресурсоемкая операция).

[Circiter]

Я кажется немного въехал в тему. И попробовал набросать в поле радактирования сообщения C-образный псевдокод, реализующий ваш метод шагающих квадратов.

Запуск выглядит как-то так:

Используется синтаксис C

void Test(matrix *Image, int Threshold)
{
    contour *Isolines;

    /* Make a regular covering mesh, which consist of triangular cells. */
    grid *Grid=TriangularGrid(Image);

    ExtractFeatures(Image, Grid, Threshold);
    Interpolate(Grid, Image);
    Isolines=TraceIsolines(Grid);

    ...
}
 

Функция TriangularGrid строит треугольную сетку на изображении, она лучше чем прямоугольная (и вариантов расположения отрезков в ячейках меньше и неопределенностей никаких не возникает, да и интерполяция по ним проще наверное будет). Основная логика заключена в функциях ExtractFeatures и TraceIsolines, о них ниже. Функция Interpolate необязательна и её я пока не написал; все что она делает -- сдвигает уже построенные куски изолиний в каждой ячейке сетки так, чтобы они проходили ближе к настощим изолиниям.

Функция ExtractFeatures может выглядеть примерно так:

Используется синтаксис C

void ExtractFeatures(matrix *Image, grid *Grid, int Threshold)
{
    do /* Go through all the cells. */
    {
        /* Get coordinates of the vertices A, B and C.
           N.B., each cell holds only the indices of its vertices. */
        vector A=Grid->Vertices[Grid->Cell->A],
            B=Grid->Vertices[Grid->Cell->B],
            C=Grid->Vertices[Grid->Cell->C];

        /* Extract actual data from the image and threshold it,
           assigning a 3-bits code to the current cell. */
        Grid->Cell->Code=Image->Data[A.x][A.y]>Threshold;
        (Grid->Cell->Code<<=1)|=(Image->Data[B.x][B.y]>Threshold);
        (Grid->Cell->Code<<=1)|=(Image->Data[C.x][C.y]>Threshold);
    }
    while(Grid->Cell = Grid->Cell->Next);
}
 

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

Теперь давайте разберемся с TraceIsolines:

Используется синтаксис C

/* Lookup table for six predefines edges;
   actually there are seven cells for technical reasons. */
edge PredefinedEdges[7]={...};

contour *TraceIsolines(grid *Grid)
{
    edge Edge;
    cell *Cell;
    contour *Contour, *Isolines=0;

    /* Enumerate the isolines. */
    while(Cell = Grid->Cell)
    {
        /* Handle empty cell. */
        if(!Cell->Code||Cell->Code==7)
        {
            RemoveCell(Grid, Cell);
            continue;
        }

        Contour=BeginContour();

        /* Trace one isoline. */
        do
        {
            Edge=PredefinedEdges[Cell->Code];

            /* Shift Edge by Position and rotate/reflect
               according to Orientation. */
            TransformEdge(&Edge, Cell->Position, Cell->Orientation);

            AddEdge(&Edge, Contour);

            RemoveCell(Grid, Cell);

            /* Determine an appropriate neighbour of current
               cell using its known position/orientation. */
            Neighbour = ... ;
            Cell = Cell->Neighbours[Neighbour]; /* Go forward. */
        } while(InGrid(Cell))

        /* Construct the linked list of isolines. */
        Contour->Next=Isolines;
        Isolines=Contour;
    }

    return Isolines;
}
 

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

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

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

В общем, интересная задачка...

[Ed_Em]

Circiter, у вас не совсем правильно: ваша ExtractFeatures работает лишь для одного уровня интенсивности, т.е. для построения N изолиний придется N раз сканировать изображение, что нецелесообразно.

По ссылке в самом начале я уже привел то, что у меня имеется - можете посмотреть на код. Сетку я делаю квадратную (так проще), а чтобы было точнее, соседние квадраты у меня накладываются друг на друга (т.е. для изображения NxM я получаю маску размером не (N/2)x(M/2), а (N-1)x(M-1)). Проблема у меня в том, что некоторые линии получаются разорванными.

Ну, а распараллелить на вашем уровне может и получится, но лишь для небольшого количества изолиний. Если, скажем, тип маски будет char (там все равно больше 4х бит не нужно - можно даже в один байт вставлять сведения о двух соседних квадратах), то в лучшем случае размер ее будет (N-1)x(M-1)/2 байт. Для небольшого изображения в 9Мпикс получаем 4.5МБ. Если это вычислять на видеокарте средней мощности, то уже где-то 50 изолиний вызовут переполнение памяти. На CPU, конечно, все ОК (т.к. оперативки обычно довольно много).

[Circiter]

2Ed_Em

Цитата:

Для небольшого изображения в 9Мпикс получаем 4.5МБ

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

Цитата:

Сетку я делаю квадратную (так проще)

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

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

Цитата:

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

Я в самом начале прошел по вашей ссылке, испугался и закрыл вкладку. :) Я поэтому и выбрал треугольную сетку, чтобы код получился попроще...

Цитата:

соседние квадраты у меня накладываются друг на друга

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

Цитата:

для построения N изолиний придется N раз сканировать изображение

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

Да, есть над чем поразмышлять... Вот, например, выяснилось ещё одно требование -- экономичность по памяти. Плюс параллельность. Ужос.

[Ed_Em]

Цитата:

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

Если сетка неперекрываемая, то только треугольная, как у вас. С квадратной такой "финт" не пройдет.

Цитата:

Вот, например, выяснилось ещё одно требование -- экономичность по памяти

Ну, естественно: когда задача рассматривается в абстракции, о ресурсах никто не задумывается, а когда доходит до дела, оказывается, что многие алгоритмы не работают из-за нехватки памяти либо чего-то еще. У меня так с CUDA было, когда я фильтры переносил с CPU на GPU (реализация свертки через БПФ), оказалось, что даже для небольшой картинки 3000x3000 512МБ памяти не хватает.

Ну, а на параллельность данного алгоритма наплевать, если он будет работать быстро.

[Pavia]

Цитата:

небольшой картинки 3000x3000

Небольшая это 320х240 разрешение веб камеры. Большая это Full HD 1920×1080 А у вас просто огромная 9 Мпикселей. Правда есть и гинабитные.

А во вторых алгоритм работает. Просто для таких больших объёмов данных его надо слегка доработать.
Методика "регион интереса" (РИ) позволит легко доработать нужные алгоритмы. Обрабатываешь кусками, а потом результаты сшиваешь.
Если РИ не помогает, то может помочь пирамида.

А вообще я бы присмотрелся к OpenCv + cvBlobsLib

Цитата:

для построения N изолиний придется N раз сканировать изображение, что нецелесообразно.

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

[Ed_Em]

Цитата:

А у вас просто огромная 9 Мпикселей.

Не такая уж и огромная: некоторые FITS'ы и еще больше бывают (широкоугольные камеры).

Цитата:

Так зачем гнаться за оптимизацией надо с начало с ошибками разобраться

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