Подтвердить что ты не робот

Эффективное сшивание строк в очень больших изображениях

Фон

Я работаю с очень большими наборами данных со спутников Radar Synthetic Aperture. Их можно рассматривать как изображения с высокой яркостью с высоким динамическим диапазоном порядка 10 тыс. Пикселей на стороне.

В последнее время я разрабатываю приложения одномасштабного варианта метод алгоритма определения масштаба шкалы Lindeberg для определения линейных характеристик в изображении SAR. Это усовершенствование использования направленных фильтров или использование преобразования Hough Transform, которые ранее использовались, потому что они менее дорогостоящи по сравнению с другими. (Я представлю некоторые недавние результаты в JURSE 2011 в апреле, и я могу загрузить препринт, если это будет полезно).

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

struct ridge_t { unsigned char top, left, bottom, right };
int rows, cols;
struct ridge_t *ridges;  /* An array of rows*cols ridge entries */

Запись в ridges содержит сегмент хребта, если ровно два из top, left, right и bottom имеют значения в диапазоне 0 - 128. Предположим, что у меня есть:

ridge_t entry;
entry.top = 25; entry.left = 255; entry.bottom = 255; entry.right = 76;

Затем я могу найти начало сегмента хребта (x1, y1) и end (x2, y2):

float x1, y1, x2, y2;
x1 = (float) col + (float) entry.top / 128.0;
y1 = (float) row;
x2 = (float) col + 1;
y2 = (float) row + (float) entry.right / 128.0;

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

Rendered ridge segments

Каждая из этих длинных кривых выводится из серии крошечных сегментов хребта.

Это тривиально, чтобы определить, связаны ли два соседних местоположения, которые содержат сегменты хребта. Если у меня есть ridge1 at (x, y) и ridge2 at (x + 1, y), то они являются частями одной и той же строки, если 0 <= ridge1.right <= 128 и ridge2.left= ridge1.right.

Проблема

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

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

Что читатели предлагают в качестве возможного подхода? Кажется, что-то вроде того, что кто-то понял бы эффективный способ сделать в прошлом...

4b9b3361

Ответ 1

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

Запрещенный алгоритм непересекающихся множеств

Я предполагаю наличие операции сравнения и замены с двумя указателями по вашему выбору архитектуры процессора. Это доступно для архитектуры x86 и x64 как минимум.

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

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

Союз(), однако, должен измениться:

function Union(x, y)
  redo:
    x = Find(x)
    y = Find(y)
    if x == y
        return
    xSnap = AtomicRead(x) -- read both rank and pointer atomically
    ySnap = AtomicRead(y) -- this operation may be done using a CAS
    if (xSnap.parent != x || ySnap.parent != y)
        goto redo
    -- Ensure x has lower rank (meaning y will be the new root)
    if (xSnap.rank > ySnap.rank)
        swap(xSnap, ySnap)
        swap(x, y)
    -- if same rank, use pointer value as a fallback sort
    else if (xSnap.rank == ySnap.rank && x > y)
        swap(xSnap, ySnap)
        swap(x, y)
    yNew = ySnap
    yNew.rank = max(yNew.rank, xSnap.rank + 1)
    xNew = xSnap
    xNew.parent = y
    if (!CAS(y, ySnap, yNew))
      goto redo
    if (!CAS(x, xSnap, xNew))
      goto redo
    return

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

  • Во-первых, до завершения один из двух корней всегда будет иметь родителя, указывающего на другой. Поэтому, если цикл отсутствует, слияние будет успешным.
  • Во-вторых, ранг всегда увеличивается. После сравнения порядка x и y мы знаем, что x имеет меньший ранг, чем y во время моментального снимка. Для того, чтобы петля сформировалась, другой поток должен был сначала увеличить ранг x, а затем объединить x и y. Однако в CAS, который пишет x родительский указатель, мы проверяем, что ранг не изменился; поэтому y-ранг должен оставаться больше x.

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

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

И теперь к приложению к вашей проблеме...

Предположения

Я делаю предположение, что сегменты хребта могут пересекаться только на своих конечных точках. Если это не так, вам нужно каким-то образом изменить фазу 1.

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

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

Этап 1: идентификация локальной сети

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

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

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

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

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

Как только сектор полностью обработан, все данные в массиве пиксельных позиций отбрасываются; однако данные, соответствующие пикселям на краю сектора, копируются в новый массив и сохраняются для следующей фазы.

Так как итерация по всему изображению O (x * y), а дизъюнктное слияние эффективно O (1), эта операция O (x * y) и требует рабочей памяти O (m + 2 * x * y/k + k ^ 2) = O (x * y/k + k ^ 2), где t - число секторов, k - ширина сектора, а m - число сегментов частичной части сектора ( в зависимости от того, как часто линии пересекают границы, m может значительно различаться, но никогда не будет превышать количество сегментов линии). Память, перенесенная на следующую операцию, равна O (m + 2 * x * y/k) = O (x * y/k)

Этап 2: кросс-секторные слияния

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

Эта операция имеет время работы O (x + y) и потребляет память O (1) (мы должны сохранить память из фазы 1, однако). По завершении реберные массивы могут быть отброшены.

Фаза 3: Сбор строк

Теперь мы выполняем операцию многопоточной карты над всеми выделенными объектами структуры с разбивкой. Сначала пропустите любой объект, который не является корнем (т.е. Где obj.parent!= Obj). Затем, начиная с репрезентативного сегмента линии, мы выходим оттуда и собираем и записываем любую информацию, необходимую о соответствующей строке. Мы уверены, что только один поток смотрит на любую заданную строку за раз, так как пересекающиеся линии оказались бы в одной и той же структуре с непересекающимися сетями.

Это время работы O (m) и использование памяти зависят от того, какую информацию вам нужно собирать по этим сегментам линии.

Резюме

В целом, этот алгоритм должен иметь время работы O (x * y) и использование памяти O (x * y/k + k ^ 2). Корректировка k дает компромисс между использованием переходной памяти в процессах фазы 1 и более длительным использованием памяти для массивов смежности и структур с несвязанной сетью, перенесенных в фазу 2.

Обратите внимание, что я действительно не тестировал производительность этого алгоритма в реальном мире; также возможно, что я упустил проблемы concurrency в алгоритме объединения-слияния без блокировки, описанном выше. Комментарии приветствуются:)

Ответ 2

Вы можете использовать не обобщенную форму Hough Transform. Похоже, что он достигает впечатляющей сложности времени O (N) в массивах сетки N x N (если у вас есть доступ к массивам SIMD ~ 10000x10000, а ваша сетка N x N - обратите внимание: в вашем случае N будет структура хребта или кластер из A x B хребтов, а не пиксель). Нажмите для источника. Более консервативные (не ядро) решения перечисляют сложность как O (kN ^ 2), где k = [-π/2, π]. Source.

Тем не менее, у Hough Transform есть некоторые требования к крутой памяти, и сложность пространства будет O (kN), но если вы прекомпилируете sin() и cos() и предоставляете соответствующие таблицы поиска, она переходит к O ( k + N), который все еще может быть слишком большим, в зависимости от того, насколько велик ваш N, но я не вижу, чтобы вы делали его ниже.

Изменить. Проблема элементов перекрестных потоков/ядра/SIMD/линии процесса нетривиальна. Мой первый импульс подсказывает мне разделить меш на рекурсивные квадраты (зависящие от определенного допуска), проверить непосредственные границы и игнорировать все структуры реберного хребта (вы можете фактически обозначить их как "потенциальные длинные строки" и поделиться ими по всей распределенной системе ); просто делайте работу над всем ВНУТРИ этого конкретного квадроцикла и постепенно двигайтесь наружу. Здесь графическое представление (зеленый - первый проход, красный - второй и т.д.). Тем не менее, моя интуиция подсказывает мне, что это вычислительно дорого.

Passes

Ответ 3

Если хребты разрешены достаточно, чтобы разрывы составляли всего несколько пикселей, то стандартные расширения - найдите соседей - шаги эрозии, которые вы сделали бы для поиска линий /OCR, должны работать.

Соединение более длинных контуров из многих сегментов и знание того, когда создавать шею или когда сделать отдельный остров, намного сложнее

Ответ 4

Хорошо, поэтому, подумав об этом немного дольше, у меня есть предложение, которое кажется слишком простым, чтобы быть эффективным... Я был бы признателен за некоторые отзывы о том, кажется ли это разумным!

1) Так как я могу легко определить, связан ли каждый сегмент хребта ridge_t с нолем, один или два соседних сегмента, я мог бы соответствующим образом покрасить каждый (LINE_NONE, LINE_END или LINE_MID). Это можно легко сделать параллельно, так как нет никаких шансов на состояние гонки.

2) После завершения окраски:

for each `LINE_END` ridge segment X found:
    traverse line until another `LINE_END` ridge segment Y found
    if X is earlier in memory than Y:
        change X to `LINE_START`
    else:
        change Y to `LINE_START`

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

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

Возможно, мне следует подумать, нужно ли собирать статистику, такую ​​как длина строки на шаге 2), чтобы помочь с окончательной переупаковкой...

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

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