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

Фон

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

Недавно я разрабатывал приложения одномасштабного варианта алгоритма Линдеберга для обнаружения гребней в масштабе и пространстве. для обнаружения линейных объектов на изображении SAR. Это улучшение по сравнению с использованием направленных фильтров или методов преобразования Хафа, которые использовались ранее, потому что они менее затратны в вычислительном отношении, чем любой из них. (Я представлю некоторые последние результаты на 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) и конец (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;

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

Визуализированные сегменты гребня

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

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

Проблема

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

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

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


person Peter T.B. Brett    schedule 29.01.2011    source источник
comment
+1 за такой хорошо написанный вопрос.   -  person    schedule 29.01.2011


Ответы (4)


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

Алгоритм непересекающихся множеств без блокировки

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

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

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

Однако Union() должен измениться:

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 больше не является корнем (в этом случае ранг не имеет значения), либо ранг 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.

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

person bdonlan    schedule 30.01.2011
comment
Это выглядит как очень интересный подход; в настоящее время работаю над этим, чтобы убедиться, что я понимаю, как все должно работать. - person Peter T.B. Brett; 31.01.2011
comment
Разумно ли использовать атомарные встроенные модули GCC (gcc.gnu .org/onlinedocs/gcc-4.1.2/gcc/Atomic-Builtins.html) для операции сравнения и замены? - person Peter T.B. Brett; 31.01.2011
comment
... не думаю: GCC разрешит любой целочисленный скаляр или тип указателя длиной 1, 2, 4 или 8 байтов. - person Peter T.B. Brett; 31.01.2011
comment
В системе x86 вы можете сделать это, упаковав два поля указателя в файл unsigned long long. Однако на x64 вам нужно будет использовать встроенную сборку для вызова инструкции cmpxchg16b. - person bdonlan; 01.02.2011
comment
Я должен указать - другой альтернативой является блокировка корней с использованием запасного бита в структуре в качестве флага блокировки и фьютексов для блокировки. Вам просто нужно убедиться, что заблокирован согласованный порядок (порядок указателя подойдет для этого). Однако это подразумевает больше атомарных операций, чем вращающийся алгоритм слияния без блокировки. - person bdonlan; 01.02.2011
comment
Оказывается, GCC превращает __sync_bool_compare_and_swap() со 128-битными аргументами (кстати, скрытый тип __int128_t полезен) в вызов __sync_bool_compare_and_swap_16(), для которого GCC не предоставляет реализации! На данный момент я запускаю этот код только на x86-64, поэтому я просто предоставил свою собственную реализацию на основе CMPXCHG16B. - person Peter T.B. Brett; 07.02.2011
comment
Я решил выполнить шаги 1 и 2 одновременно, а затем на месте свернуть леса в двусвязные списки. Честно говоря, результаты чертовски удивительны; на моей рабочей станции сшивание строк с четырьмя процессами на 35-мегапиксельном изображении занимает меньше секунды. Спасибо за отличный ответ! - person Peter T.B. Brett; 07.02.2011
comment
@Peter, рад слышать, что это сработало :) Одна вещь, на которую следует обратить внимание, если вы объединяете шаги 1 и 2, заключается в том, что операции объединения на шаге 1 должны будут использовать алгоритм, который я перечислил, - вы не можете использовать неатомарный -операция объединения в этом случае, так как вы больше не можете гарантировать, что ваш набор является локальным для текущего потока - person bdonlan; 07.02.2011
comment
Привет, bdonlan, не могли бы вы написать мне по электронной почте [email protected]? Спасибо. - person Peter T.B. Brett; 20.06.2011

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

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

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

Проходы

person David Titarenco    schedule 29.01.2011
comment
Привет, Дэвид! В литературе использовались два метода для извлечения линейных признаков из изображений SAR: наборы направленных фильтров и преобразование Хафа. Исследование, которое я провожу, предназначено для того, чтобы попытаться улучшить точность метода направленного фильтра, избегая при этом дорогостоящего шага оптимизации в преобразовании Хафа. Что меня раздражает, так это тот факт, что если вы начинаете с заданной структуры гребня, то легко пройтись по изображению, найдя все соединенные сегменты гребня, но вы не можете определить, идет ли другой поток по той же линии с другого направления. - person Peter T.B. Brett; 29.01.2011
comment
@Peter: ваша проблема нетривиальна, я предоставил предварительное решение в своем посте выше, но я думаю, что это может быть слишком наивный подход. - person David Titarenco; 29.01.2011
comment
Привет Дэвид, спасибо за это. Это определенно подход, который я рассматривал. Больше всего меня беспокоит необходимость отслеживать линии внутри рабочих единиц отдельно от линий между рабочими единицами, но об этом определенно стоит подумать! Я был бы признателен за ваше мнение о потенциальном подходе, который я только что добавил в качестве дополнительного ответа на этот вопрос... - person Peter T.B. Brett; 29.01.2011

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

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

person Martin Beckett    schedule 29.01.2011
comment
Поскольку я знаю расположение гребней с субпиксельной точностью и могу легко проверить, соединены ли два соседних сегмента гребня, алгоритмы с потерями, такие как расширение/соседи/размывание, мне не нужны. Действительно, основная причина использования методов обнаружения гребней состоит в том, чтобы обойти недостатки традиционных методов поиска линий/OCR. И да, я понимаю, что то, о чем я спрашиваю, гораздо сложнее. ;-) - person Peter T.B. Brett; 29.01.2011
comment
Хорошо - мой единственный опыт - создание контуров из нарезки TIN, где вы получаете много несвязанных сегментов прямой линии, но порядок перепутан. - person Martin Beckett; 30.01.2011

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

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_STARTs, и один раз, чтобы выполнить окончательную переупаковку, что приводит к неэффективности вычислений. Тем не менее, это по-прежнему O (N) с точки зрения хранения и времени вычислений, что является хорошим знаком...

person Peter T.B. Brett    schedule 29.01.2011
comment
Что произойдет, если линии пересекутся? - person bdonlan; 29.01.2011
comment
Ах, шарики. Хороший улов. Кроме того, я думаю, требуется обнаружение петли. - person Peter T.B. Brett; 29.01.2011