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

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

Классический алгоритм вычисления расстояния Хаусдорфа (brute force).

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

Так как это расстояние направленное, результат может отличаться от (рис 1). Поэтому симметричное расстояние Хаусдорфа обычно считают так:

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

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

directedHausdorff(A, B):
    maxMin = 0
    for a in A:
        minDist = infinity
        for b in B:
            minDist = min(minDist, distance(a, b))
        maxMin = max(maxMin, minDist)
    return maxMin

Алгоритм вычисления расстояния Хаусдорфа с ранним выходом (brute force with early break).

В случае, если нам нужно точное максимальное расстояние, можно применить известную оптимизацию раннего выхода. Если на итерации текущий минимум для уже стал меньше или равен текущему максимуму , значит, эта точка уже не сможет увеличить итоговый максимум , и можно выйти из внутреннего цикла. Эта оптимизация не меняет результат — расстояние остается точным. Мы просто не досматриваем те варианты, которые уже гарантированно не повлияют на максимум. (см. рис 3)

Как видно из визуализации, ранний выход дает существенное ускорение. Но для нашей задачи этого было недостаточно.

directedHausdorffEarlyBreak(A, B):
    maxMin = 0
    for a in A:
        minDist = infinity
        for b in B:
            d = distance(a, b)
            minDist = min(minDist, d)
            if minDist <= maxMin:
                break
        maxMin = max(maxMin, minDist)
    return maxMin

Алгоритм вычисления расстояния Хаусдорфа с перемешиванием и ранним выходом (shuffling with early break)

Идея состоит в том, что если изменить порядок обхода точек контуров и , например случайно их перемешать, то вероятность раннего выхода может возрасти. Мы быстрее можем наткнуться на точку, для которой текущий уже меньше или равен , и выйти из внутреннего цикла. Это особенно заметно на последних итерациях, когда для получения меньшего расстояния нужно проходить практически весь цикл. Оптимизация работает особенно хорошо для сильно отличающихся фигур. Но требует дополнительного выделения памяти (для перемешанных контуров) и вычислительных затрат на само перемешивание (shuffling). А также делает время работы менее предсказуемым, так как результат по числу проверок зависит от случайного порядка обхода. Похожая идея используется, например, в реализации directed_hausdorff из SciPy, где порядок точек случайно перемешивается для повышения эффективности раннего выхода. (см. рисунок 4)

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

Алгоритм вычисления расстояния Хаусдорфа с использованием структуры KD-tree.

Еще один интересный способ ускорить вычисление расстояния Хаусдорфа — использовать структуру KD-tree. Идея такая: вместо того чтобы для каждой точки перебирать все точки контура , можно заранее построить spatial index по точкам . Spatial index — это структура данных для быстрого поиска по координатам. В случае KD-tree точки рекурсивно делятся то по , то по , образуя дерево. После этого ближайшего соседа для точки можно искать не полным перебором, а проходом по дереву с отсечением областей, где ближайшей точки уже не может быть. Построение дерева имеет сложность плюс выделение памяти . Тогда в хороших случаях поиск ближайшего соседа работает существенно быстрее полного перебора и часто близок к . А общая сложность для вычисления направленного расстояния Хаусдорфа (построение + поиск) - . С точки зрения асимптотической оценки это выглядит значительно лучше, чем worst case в предыдущих рассмотреных вариантах. (см. рисунок 5)

Но здесь есть два важных момента:

• Поиск ближайшей точки в KD‑tree не имеет такой же гарантии , как поиск по ключу в обычном дереве. Для построения KD‑tree мы делим пространство: сначала по , потом по , потом снова по и так далее. Для поиска ближайшей точки мы идём в наиболее перспективную ветку, но потом иногда обязаны проверить и соседние ветки: вдруг там есть точка ближе. В хороших данных дерево отсекает много областей, и поиск близок к . Но в плохих случаях отсечь почти ничего не получается: например, точки лежат неудобно, дерево несбалансировано, запрос близок к границам разбиений, или ближайший сосед может оказаться «по другую сторону» многих split‑линий. Тогда алгоритм вынужден проверить много узлов, в худшем случае — почти все точки. Контурные точки могут быть не самым удобным случаем для KD‑tree: они лежат на тонкой кривой в 2D‑пространстве, а не равномерно заполняют область. Из‑за этого ближайшие по расстоянию точки могут оказываться по разные стороны split‑линий, и для точного поиска приходится проверять больше соседних веток.

• Для ускорения поиска в  KD-tree практические реализации часто используют approximate nearest neighbor search.  В этом случае алгоритм не проверяет все потенциальные ветки ради скорости.  В таком случае ближайший сосед может быть найден не точно, а приближенно. Для многих задач это приемлемо, но для поиска дефектов мне был нужен exact Hausdorff distance.

  Практическая реализация, к сожалению, подтвердила подозрения  – построение дерева несет слишком много накладных расходов,а упрощение поиска ближайшего соседа (использование approximate nearest neighbor search)  влияет на точность, что для поиска дефекта критично.

   Как итог для нашей задачи, результат оказался значительно хуже shuffling + early break.

Идея оптимизации вычисления расстояния Хаусдорфа

Получается интересная ситуация.

Brute force слишком дорогой. Early break помогает, но все еще проверяет точки в фиксированном порядке. Shuffling помогает раннему выходу, но требует дополнительной памяти и накладных расходов. KD-tree ускоряет поиск ближайшего соседа в общем случае, но требует дополнительной структуры данных и тоже не использует порядок точек.

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

Новый алгоритм вычисления расстояния Хаусдорфа с использованием порядка точек и ранним выходом (contour-coherent)

Главная идея алгоритма – для каждой точки начинаем поиск ближайшей точки на контуре не с нулевого индекса, а с .  Где – точка на контуре , найденная на предыдущей итерации для . Расстояние от до это либо текущий , либо меньше, если на предыдущем цикле мы вышли по earlyBreak. Обходим точки на контуре  влево и вправо от , то есть в следующем порядке:

…

Если в одном направлении дошли до границы, идем только в противоположном. Суть – раньше выйти из цикла по earlyBreak. То есть мы не строим дерево, не выделяем дополнительные массивы, не меняем точный результат — просто обходим B в более удачном порядке. (см. рисунок 6)

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

directedHausdorffContourCoherent(A, B):
    maxMin = 0
    startIndex = 0
    for i in 0 .. A.size - 1:
        minDist = infinity
        bestIndex = startIndex
        for j in indicesAround(startIndex, B.size):
            d = distance(A[i], B[j])
            if d < minDist:
                minDist = d
                bestIndex = j
            if minDist <= maxMin:
                break
        maxMin = max(maxMin, minDist)
        startIndex = bestIndex + 1
    return maxMin

Очевидно и просто? Очень! Написала код, протестировала локально на небольших пакетах, результат отличный. На следующий день тестеры удивленно разводят руками – проблема с производительностью испарилась, ускорение более 10 раз!

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

Как проводились тесты

Пакеты изображений

Для тестов использовались бинарные изображения: объект на фоне. Использовались сгенерированные мной изображения, а также стандартные наборы в Computer Vision, например, MPEG7, Kimia99, Kimia216. В таблице ниже приведена часть результатов, чтобы показать общую картину сравнения алгоритмов на разных типах контуров.

Предобработка. Извлечение контура.

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

Вычисление Хаусдорфа

Дальше внутри каждого набора изображений контуры сравнивались “каждый с каждым”. Для всех точных алгоритмов использовались одни и те же извлеченные контуры и одна и та же функция расстояния (в таблице ниже приведены замеры для евклидовой метрики).

Сравниваемые алгоритмы:

• Contour-Coherent (новая оптимизация)

• Shuffling + early break

• brute-force  + early break

• Boost.geometry (скорее для проверки точности)

• kd-tree (построение дерева с использованием OpenCV flann).

В алгоритмы Contour-Coherent, Shuffling + early break  и brute-force + early break был добавлен параметр checks – то есть количество вычислений расстояния между двумя точками. Это не полное время работы, но хороший показатель объема основной алгоритмической работы. Время дополнительно зависит от аллокаций, кэша, shuffle overhead, построения вспомогательных структур и конкретной реализации.

Бенчмарки

Замерялось: медианное время обработки всех пар в пакете (только вычисление расстояния Хаусдорфа), среднее время обработки пакета, медианное время вычисления расстояния Хаусдорфа, общее количество checks на пакете, среднее количество checks для одной пары, среднее расстояние Хаусдорфа (для проверки точности)

Результаты обозначили область применения

Лучший сценарий

Очень близкие по геометрии контуры (поиск дефекта, минимальных отличий). Это самый близкий к исходной задаче случай: небольшие локальные отличия, общий порядок точек сохраняется. Здесь ускорение было от 2,5 до 14 раз и выше по сравнению с shuffling + early break и более чем в 17–77 раз быстрее классического brute force + early break. Для статьи в качестве примера взяты наборы WalkingSilhouette, SmallDeformations, SyntheticDefects, SmallContoursDefects (все пакеты доступны на GitHub).

Хороший сценарий

Похожие контуры, но различия больше, чем в первом случае. Например, стандартные наборы MPEG7/children, MPEG7/chopper, MPEG7/car - тут ускорение от 1.3 до 5 раз по сравнению с shuffle + early break, и в 10-12 раз для brute force + early break.

Cложный, но все ещё подходящий сценарий

Контуры похожей, но сложной геометрии. Например, снежинки (набор snowflakes) или другие формы с большим количеством деталей. Если порядок точек сохраняет локальность, contour-coherent подход всё еще хорошо работает. В результате — ускорение в 2,5 раза по сравнению с shuffle + early break и в 170 раз для  brute force + early break (да, в таких случаях, если не менять порядок обхода, получается катастрофа).

Где алгоритм был менее успешным:

• На смешанных или менее согласованных наборах, например MPEG7/phone, Kimia99/Kimia216, shuffling + early break иногда оказывался быстрее contour-coherent: от почти паритета до примерно 1.5x по времени. При этом contour-coherent всё равно оставался быстрее brute force + early break примерно в 3–7 раз.

• Для очень разных контуров, например MPEG7/Mix, shuffling + early break оказался быстрее contour-coherent примерно в 4.3 раза. Это ожидаемое ограничение метода: если ближайшая точка для часто не находится рядом с предыдущей итерации, преимущество contour-coherent обхода уменьшается.

Таблица 1. Медианное время обработки всех пар контуров в пакете, ms

Таблица 2. Среднее количество distance checks на одно симметричное расстояние Хаусдорфа

Почему в таблицах нет Boost.Geometry и KD-tree

• Boost.Geometry. Я использовала Boost в первую очередь как reference implementation для проверки корректности результата. Он считает тот же exact discrete Hausdorff distance, но не возвращает количество проверок distance checks, поэтому его неудобно сравнивать в таблице checks/HD. Кроме того, по времени он оказался заметно медленнее остальных exact-baseline вариантов, поэтому в основной таблице я оставила более информативное сравнение: contour-coherent, shuffling + early break и brute force + early break.

• KD-tree. KD-tree использовался как practical nearest-neighbor baseline. Для реализации использовалась библиотека OpenCV flann. В итоге значение Avg HD  отличался от других алгоритмов, то есть не был строго равным exact baseline. Для задачи поиска дефектов это критично. Кроме того, при all-pairs сравнении дерево приходится строить много раз. Это добавляет накладные расходы и требует дополнительной памяти. В итоге и по времени этот baseline также оказался медленнее алгоритмов, оставленных в основной таблице.

Оба baseline реализованы в проекте на GitHub и могут быть запущены наряду с другими рассмотренными алгоритмами.

Итого

Алгоритм Contour-Coherent не только прошел «проверку боем» на реальных данных в высоконагруженном приложении, но и хорошо себя показал на наборах изображений разного типа.

• На близких и похожих контурах он дает стабильный выигрыш по сравнению с shuffle: от 1.3x до 14x по времени. По сравнению с BF + early break выигрыш еще заметнее: от 6.8x до 171x по времени и от 8.8x до 182x по числу checks.

• На непрофильных данных результат ожидаемо хуже: если контуры сильно отличаются или локальное соответствие между соседними точками нарушается, shuffling + early break может быть быстрее. Но даже в этих случаях contour-coherent оставался быстрее BF + early break.

Выводы

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

2. Contour‑coherent алгоритм использует эту информацию максимально дешево: не строит KD‑tree, не выделяет массивы под перемешивание, не меняет результат и не превращает exact Hausdorff distance в approximation. Он просто начинает поиск ближайшей точки рядом с тем местом, где ближайшая точка была найдена на предыдущей итерации.

3. На протестированных наборах алгоритм во всех случаях оказался быстрее обычного brute force + early break. В наиболее подходящих сценариях — близкие контуры, small deformations, synthetic defects, walking silhouettes — выигрыш составил от 14x до 77x по времени относительно brute force + early break. На Snowflake выигрыш дошел до 171x, потому что фиксированный порядок обхода оказался особенно неудачным.

4. По сравнению с shuffling + early break результат зависит от данных. На близких и согласованных контурах contour‑coherent выигрывал примерно от 2.2x до 14.5x по времени. По числу distance checks выигрыш на таких наборах составил примерно от 2x до 16x, а на SmallContoursDefects — почти 11x. Это как раз те случаи, где порядок точек действительно помогает предсказать, где искать ближайшего кандидата.

5. На смешанных или плохо согласованных наборах shuffling + early break иногда оказывается быстрее. Это ожидаемо: если ближайшая точка для часто не находится рядом с предыдущей итерации, contour‑coherent предположение теряет главное преимущество. Но даже в этих случаях contour‑coherent оставался быстрее обычного brute force + early break примерно в 3–7 раз.

6. Contour‑coherent — это не универсальный алгоритм для любых облаков точек. Но для своего класса задач — похожие упорядоченные контуры, локальные дефекты, небольшие деформации, бинарные маски — contour‑coherent обход дает редкое сочетание: простая реализация, точный результат, отсутствие дополнительной структуры данных и большой практический выигрыш.

Код, тестовые изображения, бенчмарки и визуализации доступны в GitHub-репозитории проекта.