Алгоритм Рамера–Дугласа–Пекера

Отличная штука при векторизации растра. Мы нашли контур ромба, но он содержит тысячу точек. Хотя там должно быть всего четыре точки. Четыре, Карл! Вот для того, чтобы точек стало четыре, и может пригодиться алгоритм Рамера-Дугласа-Пекера.

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

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

Здесь будет описан алгоритм с небольшими дополнениями. Дело в том, что классический алгоритм не заточен на замкнутые ломаные, полигоны. А хотелось бы. Поэтому появился ряд нюансов и дополнений, о которых речь пойдёт ниже. Это в основе своей тот же алгоритм, но более пригодный к практическому использованию. Не только в Delphi. Мне он был нужен в распознавании образов. То есть прошёл практическую апробацию.

Содержание скрыть

Описание алгоритма

Алгоритм в картинках

Нюансы реализации

Нюанс первый: Полигон (параметр AClosed)

Нюанс второй: Две ломаные (AStrictOnePolyline)

Нюанс третий: Математика

Формула

Оптимизация формулы

Оптимизация вычислений

Нюанс четвёртый: Лишняя точка (ACheckFirst)

Рекурсивный алгоритм Рамера-Дугласа-Пекера

Подготовительная часть

Рекурсивная часть

Ключевые особенности рекурсивной версии

Преимущества

Недостатки

Итеративный алгоритм Рамера-Дугласа-Пекера

Зачем нужен итеративный подход?

Суть итеративного подхода

Реализация сегмента и стека

Реализация итеративного алгоритма RDP

Ключевые особенности итеративной версии

Преимущества

Недостатки

Итеративный алгоритм RDP: Турбо-режим

Заключение

Скачать

Описание алгоритма

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

Алгоритм в картинках

Итак, у нас есть некая ломаная:

Массив результирующей ломаной пуст. Стартуем алгоритм.

A и B это начальная и конечная точки ломаной. Рассматриваем весь диапазон точек между A и B. Точку A помещаем в результирующий массив.

На диапазоне AB находим максимально удалённую точку С, которая отстоит от отрезка AB на расстояние, больше указанного. Значит, точка С — это новая точка. У нас образовалось два новых диапазона — AC и CB. Точку С пока не добавляем в массив, потому что в процессе прогона по AC вдруг окажутся ещё точки.

Рассматриваем диапазон AC. На нём ничего выдающегося не оказалось.

Добавляем C в массив точек новой ломаной и рассматриваем диапазон CB. На диапазоне CB оказалась выдающаяся точка D, которая делит диапазон на два — CD и DB.

На диапазоне CD обнаружилась выдающаяся точка E, которая снова делит диапазон на CE и ED. По прежнему ничего не добавляем в массив итоговой ломаной. Потому что не знаем, как сложится ситуация дальше, в новых диапазонах.

На диапазоне CE не нашлось ничего выдающегося. Точки m и n находятся на слишком малом расстоянии от отрезка CE, чтобы считаться выдающимися.

Поэтому дописываем в массив точку E и переходим к рассмотрению диапазона ED.

На диапазоне ED не нашлось ничего выдающегося, поэтому дописываем D. И переходим к рассмотрению диапазона DB.

На диапазоне DB не нашлось ничего выдающегося, поэтому дописываем B.

Нюансы реализации

В теории существует две разновидности алгоритма — рекурсивная и итеративная.

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

Для начала определим типы:

type
  // Динамический массив вещественных 2D-точек
  TPointFDynArray = TArray<TPointF>;
  // Список вещественных 2D-точек
  TPointFList = TList<TPointF>;

type

// Динамический массив вещественных 2D-точек

TPointFDynArray = TArray<TPointF>;

// Список вещественных 2D-точек

TPointFList = TList<TPointF>;

Константы:

const
  // Пустая точка
  EmptyPointF: TPointF = (X: MaxSingle; Y: MaxSingle);

const

// Пустая точка

EmptyPointF: TPointF = (X: MaxSingle; Y: MaxSingle);

И функции:

// Равны ли точки
function EqualPts(const P1, P2: TPointF; Eps: Single = 1e-6): Boolean;
begin
  Result := P1.EqualsTo(P2, Eps);
end;

// Пустая ли точка
function IsEmptyPt(const P: TPointF): Boolean;
begin
  Result := EqualPts(P, EmptyPointF);
end;

// Равны ли точки

function EqualPts(const P1, P2: TPointF; Eps: Single = 1e-6): Boolean;

begin

Result := P1.EqualsTo(P2, Eps);

end;

// Пустая ли точка

function IsEmptyPt(const P: TPointF): Boolean;

begin

Result := EqualPts(P, EmptyPointF);

end;

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

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

Нюанс первый: Полигон (параметр AClosed)

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

Это изображение имеет пустой атрибут alt; его имя файла - p11.jpg

A и B — это начало и конец ломаной, которая состоит из 823 точек. Всё находит изумительно, но вместо 5 вершин, я хочу видеть 4. Это же явный ромб.

Поэтому у нас появляется параметр AClosed, который определяет, с каким типом кривой работаем — замкнутой или незамкнутой.

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

    // Если требуется закрытая ломаная, совмещаем начало и конец
    if AClosed and not EqualPts(List[Len - 1], List[0]) then
    begin
      SetLength(List, Len + 1);
      List[Len] := List[0];
      Len := Len + 1;
    end;

// Если требуется закрытая ломаная, совмещаем начало и конец

if AClosed and not EqualPts(List[Len - 1], List[0]) then

begin

SetLength(List, Len + 1);

List[Len] := List[0];

Len := Len + 1;

end;

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

  // Коррекция длины массива
  if AClosed and (Len > 0) and
    EqualPts(ALines[0], ALines[Len-1]) then 
      Dec(Len);
  SetLength(ALines, Len);

// Коррекция длины массива

if AClosed and (Len > 0) and

EqualPts(ALines[0], ALines[Len-1]) then

Dec(Len);

SetLength(ALines, Len);

ALines: TPointFDynArray — это результирующий массив точек упрощённой ломаной.

Нюанс второй: Две ломаные (AStrictOnePolyline)

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

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

Пока мы ищем максимально удалённую точку, можем заодно удалить дубликаты точек.

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

Поэтому у нас появляется параметр AStrictOnePolyline, который отвечает за следующее. При установленном AStrictOnePolyline незамкнутая ломаная будет обработана классическим алгоритмом. Если у нас выставлен AClosed, действие этого параметра отменяется, деление на две ломаные будет происходить в любом случае.

    // Подготовка списка точек
    SetLength(List, AList.Count);
    Len := 0;
    P1 := AList[0];
    List[0] := P1;
    LMx := 0;
    MaxIdx := 0;

    // Поиск максимально удалённой вершины и удаление дублей
    for i := 1 to AList.Count - 1 do
    begin
      P3 := AList[i];
      if EqualPts(P3, List[Len]) then
        Continue;
      Inc(Len);
      List[Len] := P3;
      if AStrictOnePolyline and not AClosed then
        continue;
      dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);
      if dL > LMx then
      begin
        MaxIdx := Len;
        LMx := dL;
      end;
    end;

    if AStrictOnePolyline and not AClosed then
      MaxIdx := Len;

    Len := Len + 1;
    SetLength(List, Len);

// Подготовка списка точек

SetLength(List, AList.Count);

Len := 0;

P1 := AList[0];

List[0] := P1;

LMx := 0;

MaxIdx := 0;

// Поиск максимально удалённой вершины и удаление дублей

for i := 1 to AList.Count - 1 do

begin

P3 := AList[i];

if EqualPts(P3, List[Len]) then

Continue;

Inc(Len);

List[Len] := P3;

if AStrictOnePolyline and not AClosed then

continue;

dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);

if dL > LMx then

begin

MaxIdx := Len;

LMx := dL;

end;

if AStrictOnePolyline and not AClosed then

MaxIdx := Len;

Len := Len + 1;

SetLength(List, Len);

Где AList: TPointFList — исходные данные, а List: TPointFDynArray — это массив, куда мы переписали обработанные данные.

Нюанс третий: Математика

Алгоритм заключается в том, что на очередном шаге берётся участок ломаной между некими индексами i₁, i₂. Индексам соответствуют точки P1 и P2. Далее, просматриваем все точки между этими индексами и находим такую точку, которая максимально удалена от отрезка (P1, P2) на расстояние, больше заданного Epsilon. Таким образом, P1 и P2 — это неизменяемые точки отрезка в этом шаге, а P3 — перебираемая (i₁+1 .. i₂-1) точка участка, от которой ищется расстояние dL до отрезка.

Формула

Для нахождения расстояния от точки до отрезка, берём формулу высоты треугольника:

h = 2 * S / L, где:

S — площадь треугольника,
L — длина стороны треугольника, на которую опущена высота.

Площадь треугольника считаем по формуле:

S = |(X₁-X₃)*(Y₂-Y₃) — (X₂-X₃)(Y₁-Y₃)| / 2

Оптимизация формулы

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

(X₁-X₃)*(Y₂-Y₃) — (X₂-X₃)*(Y₁-Y₃) =

Раскроем скобки:

X₁Y₂ — X₃Y₂ — X₁Y₃ + X₃Y₃ — X₂Y₁ + X₃Y₁ + X₂Y₃ — X₃Y₃ =

Сгруппируем:

(X₁Y₂ — X₂Y₁ ) + (X₃Y₃ — X₃Y₃) + X₃ * (Y₁ — Y₂) + Y₃ * (X₂ — X₁) =

Сократим (X₃Y₃ — X₃Y₃). И у блока X₃*(…) поменяем знак, чтобы в скобках получить (Y₂ — Y₁):

(X₁Y₂ — X₂Y₁ ) + Y₃ * (X₂ — X₁) — X₃ * (Y₂ — Y₁) =

Выделим неизменяемое:

Y₃ * dX — X₃ * dY + A, где

A = (X₁Y₂ — X₂Y₁)
dX = (X₂ — X₁)
dY = (Y₂ — Y₁)

Таким образом, площадь треугольника:

S = |(X₁-X₃)*(Y₂-Y₃) — (X₂-X₃)(Y₁-Y₃)| / 2 = |Y₃ * dX — X₃ * dY + A| / 2

Находим высоту:

h = 2 * S / L = 2 * (|Y₃ * dX — X₃ * dY + A| / 2) / L =

|Y₃ * dX — X₃ * dY + A| / L

Оптимизация вычислений

Длина отрезка (P1, P2), или L из формулы выше, равна корню суммы квадратов dX и dY. Операция извлечения квадратного корня одна из самых дорогих операций с вещественным числом. Поэтому мы будем работать с квадратом длины отрезка (основание треугольника).

  dx := P2.x - P1.x;
  dy := P2.y - P1.y;
  // Квадрат длины основания
  L := Sqr(dx)+Sqr(dy);

dx := P2.x - P1.x;

dy := P2.y - P1.y;

// Квадрат длины основания

L := Sqr(dx)+Sqr(dy);

При нахождении расстояния от точки P3 мы не будем брать модуль |Y₃ * dX — X₃ * dY + A|. Мы будем брать квадрат.

      // Вершина, от которой считаем высоту
      P3 := APoints[i];
      // Квадрат высоты треугольника:
      // h = 2 * S/L = 2*(1/2*|Y*DX-X*DY+A|)/L
      dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

// Вершина, от которой считаем высоту

P3 := APoints[i];

// Квадрат высоты треугольника:

// h = 2 * S/L = 2*(1/2*|Y*DX-X*DY+A|)/L

dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

Чтобы сравнение с Epsilon было корректным, будем сравнивать с его квадратом. В самом начале, запоминаем квадрат отклонения.

EpsilonSqr := Sqr(AEpsilon);

1	EpsilonSqr := Sqr(AEpsilon);

В дальнейших вычислениях будем учитывать только его:

  for i := StartIdx + 1 to EndIdx - 1 do
  begin
    P3 := List[i];
    dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;
    if dL > dMax then
    begin
      dMax := dL;
      maxIndex := i;
    end;
  end;

  // Если максимальное отклонение превышает epsilon, обрабатываем дальше
  if (dMax > EpsilonSqr) and (maxIndex > StartIdx) and (maxIndex < EndIdx) then
  begin
    // Обработка сегментов (StartIdx, maxIndex) и (maxIndex, EndIdx)
  end;
end;

for i := StartIdx + 1 to EndIdx - 1 do

begin

P3 := List[i];

dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

if dL > dMax then

begin

dMax := dL;

maxIndex := i;

end;

// Если максимальное отклонение превышает epsilon, обрабатываем дальше

if (dMax > EpsilonSqr) and (maxIndex > StartIdx) and (maxIndex < EndIdx) then

begin

// Обработка сегментов (StartIdx, maxIndex) и (maxIndex, EndIdx)

end;

Нюанс четвёртый: Лишняя точка (ACheckFirst)

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

Это изображение имеет пустой атрибут alt; его имя файла - p12.jpg

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

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

// Проверить начальную точку и, если необходимо, удалить
function CheckFirstPoint(AEpsilonSqr: Single;
  var ALines: TPointFDynArray): Boolean;
var L, dx, dy: Single; P1, P2, P3: TPointF;
begin
  if Length(ALines) <= 2 then
    exit(False);
  P3 := ALines[0];
  P1 := ALines[1];
  P2 := ALines[High(ALines)];
  dx := P2.x - P1.x;
  dy := P2.y - P1.y;
  L := Sqr(dx) + Sqr(dy);
  if IsZero(L) then
    exit(False);
  L := Sqr(P3.Y*dx - P3.X*dy + P1.X*P2.Y - P2.X*P1.Y)/L;
  Result := L < AEpsilonSqr;
  if Result then
    Delete(ALines, 0, 1);
end;

// Проверить начальную точку и, если необходимо, удалить

function CheckFirstPoint(AEpsilonSqr: Single;

var ALines: TPointFDynArray): Boolean;

var L, dx, dy: Single; P1, P2, P3: TPointF;

begin

if Length(ALines) <= 2 then

exit(False);

P3 := ALines[0];

P1 := ALines[1];

P2 := ALines[High(ALines)];

dx := P2.x - P1.x;

dy := P2.y - P1.y;

L := Sqr(dx) + Sqr(dy);

if IsZero(L) then

exit(False);

L := Sqr(P3.Y*dx - P3.X*dy + P1.X*P2.Y - P2.X*P1.Y)/L;

Result := L < AEpsilonSqr;

if Result then

Delete(ALines, 0, 1);

end;

Дополним список параметров аргументом ACheckFirst. И если он равен TRUE, в самом конце нашей запускающей процедуры вызовем CheckFirstPoint.

Рекурсивный алгоритм Рамера-Дугласа-Пекера

Рекурсивная версия алгоритма RDP состоит из двух ключевых частей:

Подготовительная часть

В подготовительной части происходит подготовка данных и запуск рекурсии:

AList: TPointFList — исходные данные, список точек;
AEpsilon: Single — точность упрощения. Величина, превышение которой означает появление новой выдающейся точки;
AClosed, ACheckFirst, AStrictOnePolyline — описаны выше;
out ALines: TPointFDynArray — результат: массив с точками упрощённой линии.

// Подготовительная часть: подготовить данные и вызвать рекурсию
function RamerDouglasPeuckerRecourse(AList: TPointFList;
  AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;
  out ALines: TPointFDynArray): Boolean;
var
  List: TPointFDynArray;
  Len, MaxIdx, i: Integer;
  LMx, dL: Single;
  P, Start: TPointF;
  EpsilonSqr: Single;
begin
  // С пустыми данными не работаем
  Result := (AList<>nil) and (AList.Count>0);
  if not Result then exit;
  EpsilonSqr := Sqr(AEpsilon);

  // Инициализируем индексы, длину будущего массива
  MaxIdx := 0; // Индекс равен стартовому
  LMx := 0; // Максимальная длина
  Len := 0; // Длина массива
  // Пока длина массива равна исходной
  SetLength(List, AList.Count);
  // Поиск максимально удалённой будет происходить
  // от первой точки ломаной
  Start := AList[0];
  // Фиксируем в исходном массиве стартовую точку
  List[0] := Start;
  // Поиск максимально удалённой вершины
  // Заодно удаление дублей
  for i := 1 to AList.Count-1 do
  begin
    P := AList[i];
    // Совпадает с предыдущей, ничего не делаем
    if EqualPts(P, List[Len]) then
      continue;
    Inc(Len);
    List[Len] := P;
    // Если это не полигон, на две части не делим
    // Следовательно, искать самую дальнюю не надо
    if AStrictOnePolyline and not AClosed then
      continue;
    // Находим текущий максимум
    dL := Sqr(p.X-Start.X) + Sqr(p.Y-Start.Y);
    if dL > LMx then
    begin
      MaxIdx := Len;
      LMx := dL;
    end;
  end;

  // Если это не полигон, максимальный индекс - конец ломаной
  if AStrictOnePolyline and not AClosed then
    MaxIdx := Len;

  // Высчитываем конечную длину массива
  Len := Len+1;
  // Выставляем длину массива данных
  SetLength(List, Len);

  // Если требуется закрытая ломаная, совмещаем
  // начало и конец массива
  if AClosed and not EqualPts(List[Len-1], List[0]) then
  begin
    Len := Len+1;
    SetLength(List, Len);
    List[Len-1] := List[0];
  end;

  // Инициализируем результат
  SetLength(ALines, Len);
  ALines[0] := List[0];
  Len := 1;

  // Пошла рекурсия в диапазоне от старта
  // до максимально удалённой (либо конечной)
  DouglasPeucker(List, 0, MaxIdx, EpsilonSqr, Len, ALines);
  // Дописываем последнюю точку диапазона
  // Которая одновременно будет являться начальной
  // для возможной следующей ломаной
  if Len < Length(ALines) then
  begin
    ALines[Len] := List[MaxIdx];
    Inc(Len);
  end;
  // Пошла рекурсия в диапазоне от
  // максимально удалённой до конца
  if MaxIdx < High(List) then
  begin
    DouglasPeucker(List, MaxIdx, High(List),
      EpsilonSqr, Len, ALines);
    // Дописываем последнюю точку диапазона
    if Len < Length(ALines) then
    begin
      ALines[Len] := List[High(List)];
      Inc(Len);
    end;
  end;

  // Коррекция длины массива
  if AClosed and (Len > 0) and
    EqualPts(ALines[0], ALines[Len-1]) then Dec(Len);
  SetLength(ALines, Len);

  // Проверить начальную точку
  if AClosed and ACheckFirst then
    while CheckFirstPoint(EpsilonSqr, ALines) do;
end;

100

101

102

103

104

105

106

107

108

// Подготовительная часть: подготовить данные и вызвать рекурсию

function RamerDouglasPeuckerRecourse(AList: TPointFList;

AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;

out ALines: TPointFDynArray): Boolean;

var

List: TPointFDynArray;

Len, MaxIdx, i: Integer;

LMx, dL: Single;

P, Start: TPointF;

EpsilonSqr: Single;

begin

// С пустыми данными не работаем

Result := (AList<>nil) and (AList.Count>0);

if not Result then exit;

EpsilonSqr := Sqr(AEpsilon);

// Инициализируем индексы, длину будущего массива

MaxIdx := 0; // Индекс равен стартовому

LMx := 0; // Максимальная длина

Len := 0; // Длина массива

// Пока длина массива равна исходной

SetLength(List, AList.Count);

// Поиск максимально удалённой будет происходить

// от первой точки ломаной

Start := AList[0];

// Фиксируем в исходном массиве стартовую точку

List[0] := Start;

// Поиск максимально удалённой вершины

// Заодно удаление дублей

for i := 1 to AList.Count-1 do

begin

P := AList[i];

// Совпадает с предыдущей, ничего не делаем

if EqualPts(P, List[Len]) then

continue;

Inc(Len);

List[Len] := P;

// Если это не полигон, на две части не делим

// Следовательно, искать самую дальнюю не надо

if AStrictOnePolyline and not AClosed then

continue;

// Находим текущий максимум

dL := Sqr(p.X-Start.X) + Sqr(p.Y-Start.Y);

if dL > LMx then

begin

MaxIdx := Len;

LMx := dL;

end;

// Если это не полигон, максимальный индекс - конец ломаной

if AStrictOnePolyline and not AClosed then

MaxIdx := Len;

// Высчитываем конечную длину массива

Len := Len+1;

// Выставляем длину массива данных

SetLength(List, Len);

// Если требуется закрытая ломаная, совмещаем

// начало и конец массива

if AClosed and not EqualPts(List[Len-1], List[0]) then

begin

Len := Len+1;

SetLength(List, Len);

List[Len-1] := List[0];

end;

// Инициализируем результат

SetLength(ALines, Len);

ALines[0] := List[0];

Len := 1;

// Пошла рекурсия в диапазоне от старта

// до максимально удалённой (либо конечной)

DouglasPeucker(List, 0, MaxIdx, EpsilonSqr, Len, ALines);

// Дописываем последнюю точку диапазона

// Которая одновременно будет являться начальной

// для возможной следующей ломаной

if Len < Length(ALines) then

begin

ALines[Len] := List[MaxIdx];

Inc(Len);

end;

// Пошла рекурсия в диапазоне от

// максимально удалённой до конца

if MaxIdx < High(List) then

begin

DouglasPeucker(List, MaxIdx, High(List),

EpsilonSqr, Len, ALines);

// Дописываем последнюю точку диапазона

if Len < Length(ALines) then

begin

ALines[Len] := List[High(List)];

Inc(Len);

end;

// Коррекция длины массива

if AClosed and (Len > 0) and

EqualPts(ALines[0], ALines[Len-1]) then Dec(Len);

SetLength(ALines, Len);

// Проверить начальную точку

if AClosed and ACheckFirst then

while CheckFirstPoint(EpsilonSqr, ALines) do;

end;

Рекурсивная часть

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

APoints: TPointFDynArray — Массив вещественных точек исходной ломаной;
AStartIdx: Integer — Индекс точки в APoints, означающий начало текущего рассматриваемого диапазона;
AEndIdx: Integer — Индекс точки в APoints, означающий конец текущего рассматриваемого диапазона;
AEpsilon: Single — Величина, превышение которой означает появление новой выдающейся точки. Все точки, отстоящие от отрезка, заданного точками AStartIdx и AEndIdx, и меньшие этой величины, считаются принадлежащими этому отрезку;
var AResLength: Integer — Текущая длина массива точек результирующей ломаной;
var ARes: TPointFDynArray — Массив точек результирующей ломаной.

// Рекурсивная часть алгоритм
procedure DouglasPeucker(const APoints: TPointFDynArray;
  AStartIdx, AEndIdx: Integer; AEpsilonSqr: Single;
  var AResLength: Integer; var ARes: TPointFDynArray);
var
  P1, P2, P3: TPointF;
  dx, dy, L, A, dL, dMax: Double;
  i, Idx: Integer;
begin
  dMax := 0; // Максимальная высота
  Idx := -1; // Индекс пока не найден
  // Концы отрезка
  P1 := APoints[AStartIdx];
  P2 := APoints[AEndIdx];
  dx := P2.x - P1.x;
  dy := P2.y - P1.y;
  // Квадрат длины основания
  L := Sqr(dx)+Sqr(dy);
  // Неизменная часть формулы
  A := P1.X*P2.Y - P2.X*P1.Y;
  // Если концы отрезка не совпадают
  if not IsZero(L) then
    for i := AStartIdx+1 to AEndIdx-1 do
    begin
      // Вершина, от которой считаем высоту
      P3 := APoints[i];
      // Квадрат высоты треугольника:
      dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;
      if dL > dMax then
      begin
        Idx := i;
        dMax := dL;
      end;
    end;
  // Нашлась выдающаяся точка
  if dMax > AEpsilonSqr then
  begin
    // Погнали рекурсию в диапазоне
    // от начала до этой вершины
    DouglasPeucker(APoints, AStartIdx, Idx,
      AEpsilonSqr, AResLength, ARes);
    // Запоминаем выдающуюся точку
    ARes[AResLength] := APoints[Idx];
    // Инкремент длины массива результата
    Inc(AResLength);
    // Погнали рекурсию в диапазоне
    // от этой вершины до конца
    DouglasPeucker(APoints, Idx, AEndIdx,
      AEpsilonSqr, AResLength, ARes);
  end;
end;

// Рекурсивная часть алгоритм

procedure DouglasPeucker(const APoints: TPointFDynArray;

AStartIdx, AEndIdx: Integer; AEpsilonSqr: Single;

var AResLength: Integer; var ARes: TPointFDynArray);

var

P1, P2, P3: TPointF;

dx, dy, L, A, dL, dMax: Double;

i, Idx: Integer;

begin

dMax := 0; // Максимальная высота

Idx := -1; // Индекс пока не найден

// Концы отрезка

P1 := APoints[AStartIdx];

P2 := APoints[AEndIdx];

dx := P2.x - P1.x;

dy := P2.y - P1.y;

// Квадрат длины основания

L := Sqr(dx)+Sqr(dy);

// Неизменная часть формулы

A := P1.X*P2.Y - P2.X*P1.Y;

// Если концы отрезка не совпадают

if not IsZero(L) then

for i := AStartIdx+1 to AEndIdx-1 do

begin

// Вершина, от которой считаем высоту

P3 := APoints[i];

// Квадрат высоты треугольника:

dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

if dL > dMax then

begin

Idx := i;

dMax := dL;

end;

// Нашлась выдающаяся точка

if dMax > AEpsilonSqr then

begin

// Погнали рекурсию в диапазоне

// от начала до этой вершины

DouglasPeucker(APoints, AStartIdx, Idx,

AEpsilonSqr, AResLength, ARes);

// Запоминаем выдающуюся точку

ARes[AResLength] := APoints[Idx];

// Инкремент длины массива результата

Inc(AResLength);

// Погнали рекурсию в диапазоне

// от этой вершины до конца

DouglasPeucker(APoints, Idx, AEndIdx,

AEpsilonSqr, AResLength, ARes);

end;

Ключевые особенности рекурсивной версии

Рекурсивная реализация — это самый понятный способ записать алгоритм RDP, который почти дословно повторяет его описание: «берём отрезок, если он простой — оставляем концы, если нет — находим важную точку и повторяем для обеих половин».

Преимущества

Простота понимания — каждый шаг делает одно и то же, только для меньшего отрезка;
Наглядность — чётко видна логика «разделяй и властвуй»;
Минимальный расход памяти — не нужно хранит промежуточные массивы, всё формируется на лету.

Недостатки

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

Итеративный алгоритм Рамера-Дугласа-Пекера

Зачем нужен итеративный подход?

Хотя рекурсивная реализация алгоритма элегантна и наглядна, она имеет существенное ограничение — риск переполнения стека вызовов при обработке больших наборов данных. Каждый рекурсивный вызов добавляет новый фрейм в стек, и при глубокой вложенности это может привести к исключению Stack Overflow или Out Of Memory или к другой не менее приятной вещице.

К тому же, когда стал делать рекурсивный алгоритм в JS, столкнулся с рядом неприятных вещей, которые привели к сильным тормозам. Не всегда можно повторить рекурсию. Алгоритм, который реализован исключительно в цикле, можно повторить где угодно.

Итеративный алгоритм прост, надёжен и в теории быстрее рекурсивного (но это не так, мы его сделаем таким в самом конце).

Суть итеративного подхода

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

Основные этапы итеративного алгоритма:

Инициализация: Создаем стек и помещаем в него начальный сегмент (0, n-1);
Массив отметок: Создаем булевый массив для отметки точек, которые войдут в результат
Цикл обработки: Пока стек не пуст:
- Извлекаем сегмент (start, end);
- Находим точку с максимальным отклонением (используя ту же математику, что и в рекурсивной версии);
- Если отклонение превышает Epsilon, отмечаем точку и добавляем два новых сегмента в стек;
Формирование результата: Собираем все отмеченные точки

Реализация сегмента и стека

Для сегмента вводим специальный тип:

type
  // Структура для хранения сегментов в стеке
  TSegment = packed record
    StartIdx, EndIdx: Integer;
    // Конструктор работает медленнее
    class function Create(AStartIdx, AEndIdx: Integer): TSegment; static;
  end;

{ TSegment }

class function TSegment.Create(AStartIdx, AEndIdx: Integer): TSegment;
begin
  Result.StartIdx := AStartIdx;
  Result.EndIdx := AEndIdx;
end;

type

// Структура для хранения сегментов в стеке

TSegment = packed record

StartIdx, EndIdx: Integer;

// Конструктор работает медленнее

class function Create(AStartIdx, AEndIdx: Integer): TSegment; static;

end;

{ TSegment }

class function TSegment.Create(AStartIdx, AEndIdx: Integer): TSegment;

begin

Result.StartIdx := AStartIdx;

Result.EndIdx := AEndIdx;

end;

Задача TSegment‘а просто хранить индексы участка ломаной, который надо обработать. Классовая функция Create выполняет ту же функцию, что и конструктор, только быстрее.

Структуры TSegment помещаются в стек, чтобы потом быть извлечёнными и обработанными. В Delphi есть специальный тип для реализации стека — TStack<T>. Но в борьбе за быстродействие пришлось реализовать свой очень минималистический стек.

type
  TSegmentStack = record
  private
    FItems: array of TSegment;
    FCount: Integer;
    FCapacity: Integer;
  public
    procedure Initialize(InitialCapacity: Integer = 128);
    procedure Free;
    function Count: Integer; inline;
    procedure Push(const Segment: TSegment);
    function Pop(out Segment: TSegment): Boolean;
  end;

{ TSegmentStack }

procedure TSegmentStack.Initialize(InitialCapacity: Integer);
begin
  FCapacity := InitialCapacity;
  SetLength(FItems, FCapacity);
  FCount := 0;
end;

procedure TSegmentStack.Free;
begin
  SetLength(FItems, 0);
  FCount := 0;
  FCapacity := 0;
end;

function TSegmentStack.Count: Integer;
begin
  Result := FCount;
end;

procedure TSegmentStack.Push(const Segment: TSegment);
begin
  if FCount = FCapacity then
  begin
    FCapacity := FCapacity * 2;
    SetLength(FItems, FCapacity);
  end;
  FItems[FCount] := Segment;
  Inc(FCount);
end;

function TSegmentStack.Pop(out Segment: TSegment): Boolean;
begin
  Result := FCount > 0;
  if Result then
  begin
    Dec(FCount);
    Segment := FItems[FCount];
  end;
end;

type

TSegmentStack = record

private

FItems: array of TSegment;

FCount: Integer;

FCapacity: Integer;

public

procedure Initialize(InitialCapacity: Integer = 128);

procedure Free;

function Count: Integer; inline;

procedure Push(const Segment: TSegment);

function Pop(out Segment: TSegment): Boolean;

end;

{ TSegmentStack }

procedure TSegmentStack.Initialize(InitialCapacity: Integer);

begin

FCapacity := InitialCapacity;

SetLength(FItems, FCapacity);

FCount := 0;

end;

procedure TSegmentStack.Free;

begin

SetLength(FItems, 0);

FCount := 0;

FCapacity := 0;

end;

function TSegmentStack.Count: Integer;

begin

Result := FCount;

end;

procedure TSegmentStack.Push(const Segment: TSegment);

begin

if FCount = FCapacity then

begin

FCapacity := FCapacity * 2;

SetLength(FItems, FCapacity);

end;

FItems[FCount] := Segment;

Inc(FCount);

end;

function TSegmentStack.Pop(out Segment: TSegment): Boolean;

begin

Result := FCount > 0;

if Result then

begin

Dec(FCount);

Segment := FItems[FCount];

end;

В итеративном алгоритме ниже, мы инициализируем стек такой записью:

  // Инициализация
  Stack.Initialize(AList.Count div 4 + 8);

1 2	// Инициализация Stack.Initialize(AList.Count div 4 + 8);

Откуда взялось именно AList.Count div 4 + 8. Это эмпирическая формула, основанная на особенностях алгоритма Рамера-Дугласа-Пекера. Суть инициализации стека — изначально создать массив, которого скорее всего хватит на всю обработку без дополнительного перераспределения памяти.

Почему делим на 4:

Эмпирические наблюдения: Для 1000 точек стек редко превышает 250 элементов;
Бинарное дерево сегментов: Максимально возможных сегментов N-1, но большинство отсекается из-за Epsilon; Типичное количество активных сегментов: N/4.

Зачем добавляем 8:

Минимальный размер для мелких полигонов: Даже для 10 точек: 10/4 + 8 =10 элементов. Без +8 было бы всего 2 элемента;
Избегаем частого реаллока для маленьких массивов;
Запас для начальных сегментов:

// В начале алгоритма кладем 2 сегмента:
Stack.Push(0, MaxIdx);      // 1
Stack.Push(MaxIdx, Len-1);  // 2
// + возможные дополнительные

// В начале алгоритма кладем 2 сегмента:

Stack.Push(0, MaxIdx); // 1

Stack.Push(MaxIdx, Len-1); // 2

// + возможные дополнительные

Реализация итеративного алгоритма RDP

AList: TPointFList — исходные данные, список точек;
AEpsilon: Single — точность упрощения. Величина, превышение которой означает появление новой выдающейся точки;
AClosed, ACheckFirst, AStrictOnePolyline — описаны выше;
out ALines: TPointFDynArray — результат: массив с точками упрощённой линии.

// Итеративный алгоритм Рамера-Дугласа-Пекера
function RamerDouglasPeuckerIterative(AList: TPointFList;
  AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;
  out ALines: TPointFDynArray): Boolean;
var
  List: TPointFDynArray;
  Keep: array of Boolean;
  Stack: TSegmentStack;
  Segment: TSegment;
  i, Len, MaxIdx: Integer;
  P1, P2, P3: TPointF;
  LMx, dx, dy, L, dL, A: Single;
  ResLength: Integer;
  EpsilonSqr: Single;
begin
  // С пустыми данными не работаем
  Result := (AList <> nil) and (AList.Count > 0);
  if not Result then
    exit;

  EpsilonSqr := Sqr(AEpsilon);
  // Инициализация
  Stack.Initialize(AList.Count div 4 + 8);
  try
    // Подготовка списка точек
    SetLength(List, AList.Count);
    Len := 0;
    P1 := AList[0];
    List[0] := P1;
    LMx := 0;
    MaxIdx := 0;

    // Поиск максимально удалённой вершины и удаление дублей
    for i := 1 to AList.Count - 1 do
    begin
      P3 := AList[i];
      if EqualPts(P3, List[Len]) then
        Continue;
      Inc(Len);
      List[Len] := P3;
      if AStrictOnePolyline and not AClosed then
        continue;
      dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);
      if dL > LMx then
      begin
        MaxIdx := Len;
        LMx := dL;
      end;
    end;

    if AStrictOnePolyline and not AClosed then
      MaxIdx := Len;

    Len := Len + 1;
    SetLength(List, Len);

    // Если требуется закрытая ломаная, совмещаем начало и конец
    if AClosed and not EqualPts(List[Len - 1], List[0]) then
    begin
      SetLength(List, Len + 1);
      List[Len] := List[0];
      Len := Len + 1;
    end;

    // Массив для отметки точек, которые нужно сохранить
    SetLength(Keep, Len);
    for i := 1 to Len - 2 do
      Keep[i] := False;

    // Всегда сохраняем первую, последнюю и
    // максимально удалённую точки
    Keep[0] := True;
    Keep[Len-1] := True;
    // Если AStrictOnePolyline and not AClosed,
    // то Keep[MaxIdx] = Keep[Len-1]
    Keep[MaxIdx] := True;

    // Добавляем начальные сегменты в стек
    if MaxIdx > 0 then
      Stack.Push(TSegment.Create(0, MaxIdx));

    if MaxIdx < Len - 1 then
      Stack.Push(TSegment.Create(MaxIdx, Len - 1));

    // Обрабатываем сегменты пока стек не пуст
    while Stack.Count > 0 do
    begin
      Stack.Pop(Segment);
      if Segment.EndIdx - Segment.StartIdx < 2 then
        continue;

      // Находим точку с максимальным отклонением
      LMx := 0;
      MaxIdx := Segment.StartIdx;
      P1 := List[Segment.StartIdx];
      P2 := List[Segment.EndIdx];

      dx := P2.X - P1.X;
      dy := P2.Y - P1.Y;
      L := Sqr(dx) + Sqr(dy);
      if IsZero(L) then
        continue;

      A := P1.X*P2.Y - P2.X*P1.Y;

      for i := Segment.StartIdx + 1 to Segment.EndIdx - 1 do
      begin
        P3 := List[i];
        dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;
        if dL > LMx then
        begin
          LMx := dL;
          MaxIdx := i;
        end;
      end;

      // Если максимальное отклонение превышает epsilon, обрабатываем дальше
      if (LMx > EpsilonSqr) and
        (MaxIdx > Segment.StartIdx) and
        (MaxIdx < Segment.EndIdx)
      then
      begin
        Keep[MaxIdx] := True;
        Stack.Push(TSegment.Create(Segment.StartIdx, MaxIdx));
        Stack.Push(TSegment.Create(MaxIdx, Segment.EndIdx));
      end;
    end;

    // Формируем результирующий массив
    ResLength := 0;
    SetLength(ALines, Len);// Максимально возможный размер

    for i := 0 to Len - 1 do
    begin
      if Keep[i] then
      begin
        ALines[ResLength] := List[i];
        Inc(ResLength);
      end;
    end;

    // Для замкнутой ломаной проверяем,
    // что первая и последняя точки не дублируются
    if AClosed and (ResLength > 1) and
      EqualPts(ALines[0], ALines[ResLength - 1]) then
        // Убираем дублирующую последнюю точку
        Dec(ResLength);

    // Корректируем длину массива
    SetLength(ALines, ResLength);

    // Проверить начальную точку для замкнутой ломаной
    if ACheckFirst and AClosed then
      while CheckFirstPoint(EpsilonSqr, ALines) do;

  finally
    Stack.Free;
  end;
end;

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

// Итеративный алгоритм Рамера-Дугласа-Пекера

function RamerDouglasPeuckerIterative(AList: TPointFList;

AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;

out ALines: TPointFDynArray): Boolean;

var

List: TPointFDynArray;

Keep: array of Boolean;

Stack: TSegmentStack;

Segment: TSegment;

i, Len, MaxIdx: Integer;

P1, P2, P3: TPointF;

LMx, dx, dy, L, dL, A: Single;

ResLength: Integer;

EpsilonSqr: Single;

begin

// С пустыми данными не работаем

Result := (AList <> nil) and (AList.Count > 0);

if not Result then

exit;

EpsilonSqr := Sqr(AEpsilon);

// Инициализация

Stack.Initialize(AList.Count div 4 + 8);

try

// Подготовка списка точек

SetLength(List, AList.Count);

Len := 0;

P1 := AList[0];

List[0] := P1;

LMx := 0;

MaxIdx := 0;

// Поиск максимально удалённой вершины и удаление дублей

for i := 1 to AList.Count - 1 do

begin

P3 := AList[i];

if EqualPts(P3, List[Len]) then

Continue;

Inc(Len);

List[Len] := P3;

if AStrictOnePolyline and not AClosed then

continue;

dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);

if dL > LMx then

begin

MaxIdx := Len;

LMx := dL;

end;

if AStrictOnePolyline and not AClosed then

MaxIdx := Len;

Len := Len + 1;

SetLength(List, Len);

// Если требуется закрытая ломаная, совмещаем начало и конец

if AClosed and not EqualPts(List[Len - 1], List[0]) then

begin

SetLength(List, Len + 1);

List[Len] := List[0];

Len := Len + 1;

end;

// Массив для отметки точек, которые нужно сохранить

SetLength(Keep, Len);

for i := 1 to Len - 2 do

Keep[i] := False;

// Всегда сохраняем первую, последнюю и

// максимально удалённую точки

Keep[0] := True;

Keep[Len-1] := True;

// Если AStrictOnePolyline and not AClosed,

// то Keep[MaxIdx] = Keep[Len-1]

Keep[MaxIdx] := True;

// Добавляем начальные сегменты в стек

if MaxIdx > 0 then

Stack.Push(TSegment.Create(0, MaxIdx));

if MaxIdx < Len - 1 then

Stack.Push(TSegment.Create(MaxIdx, Len - 1));

// Обрабатываем сегменты пока стек не пуст

while Stack.Count > 0 do

begin

Stack.Pop(Segment);

if Segment.EndIdx - Segment.StartIdx < 2 then

continue;

// Находим точку с максимальным отклонением

LMx := 0;

MaxIdx := Segment.StartIdx;

P1 := List[Segment.StartIdx];

P2 := List[Segment.EndIdx];

dx := P2.X - P1.X;

dy := P2.Y - P1.Y;

L := Sqr(dx) + Sqr(dy);

if IsZero(L) then

continue;

A := P1.X*P2.Y - P2.X*P1.Y;

for i := Segment.StartIdx + 1 to Segment.EndIdx - 1 do

begin

P3 := List[i];

dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

if dL > LMx then

begin

LMx := dL;

MaxIdx := i;

end;

// Если максимальное отклонение превышает epsilon, обрабатываем дальше

if (LMx > EpsilonSqr) and

(MaxIdx > Segment.StartIdx) and

(MaxIdx < Segment.EndIdx)

then

begin

Keep[MaxIdx] := True;

Stack.Push(TSegment.Create(Segment.StartIdx, MaxIdx));

Stack.Push(TSegment.Create(MaxIdx, Segment.EndIdx));

end;

// Формируем результирующий массив

ResLength := 0;

SetLength(ALines, Len);// Максимально возможный размер

for i := 0 to Len - 1 do

begin

if Keep[i] then

begin

ALines[ResLength] := List[i];

Inc(ResLength);

end;

// Для замкнутой ломаной проверяем,

// что первая и последняя точки не дублируются

if AClosed and (ResLength > 1) and

EqualPts(ALines[0], ALines[ResLength - 1]) then

// Убираем дублирующую последнюю точку

Dec(ResLength);

// Корректируем длину массива

SetLength(ALines, ResLength);

// Проверить начальную точку для замкнутой ломаной

if ACheckFirst and AClosed then

while CheckFirstPoint(EpsilonSqr, ALines) do;

finally

Stack.Free;

end;

Ключевые особенности итеративной версии

Преимущества

Безопасность и надёжность

Нет риска переполнения стека. Контроль над памятью и простота отладки.

// Предсказуемое потребление
Stack.Initialize(AList.Count div 4 + 8);  // Заранее выделяем

1 2	// Предсказуемое потребление Stack.Initialize(AList.Count div 4 + 8); // Заранее выделяем

Возможность оптимизации

Например, большие сегменты обрабатываются первыми.

// Приоритетная обработка (большие сегменты первыми) 
// или другая специфическая экзотика
if (EndIdx - StartIdx) > 100 then
  ProcessLargeSegmentFirst;

// Приоритетная обработка (большие сегменты первыми)

// или другая специфическая экзотика

if (EndIdx - StartIdx) > 100 then

ProcessLargeSegmentFirst;

Параллелизация

Можно разделить на несколько потоков.

// Поток 1: обрабатывает 0..N/2
// Поток 2: обрабатывает N/2..N
// Объединяем результаты

// Поток 1: обрабатывает 0..N/2

// Поток 2: обрабатывает N/2..N

// Объединяем результаты

Измеримость прогресса

// Прогресс обработки
Progress := (TotalSegments - Stack.Count) / TotalSegments;
OnProgress(Progress * 100);  // Callback для UI

// Прогресс обработки

Progress := (TotalSegments - Stack.Count) / TotalSegments;

OnProgress(Progress * 100); // Callback для UI

Оптимизации доступа к памяти

// Линейный доступ к массивам
for i := StartIdx to EndIdx do
  Points[i]...  // Предсказуемо для кэша процессора

// vs рекурсия: прыжки по памяти

// Линейный доступ к массивам

for i := StartIdx to EndIdx do

Points[i]... // Предсказуемо для кэша процессора

// vs рекурсия: прыжки по памяти

Недостатки

Больше служебного кода

// Нужно управлять стеком вручную
Stack.Push(TSegment.Create(Start, Idx));
Stack.Push(TSegment.Create(Idx, End));

// vs рекурсия:
Recurse(Start, Idx);
Recurse(Idx, End);

// Нужно управлять стеком вручную

Stack.Push(TSegment.Create(Start, Idx));

Stack.Push(TSegment.Create(Idx, End));

// vs рекурсия:

Recurse(Start, Idx);

Recurse(Idx, End);

Дополнительные структуры данных

// Требуется массив флагов
Keep: array of Boolean;           // + память
for i := 0 to High(Keep) do       // + время инициализации
  Keep[i] := False;

// И финальный сбор
for i := 0 to High(Points) do     // + дополнительный проход
  if Keep[i] then AddToResult;

// Требуется массив флагов

Keep: array of Boolean; // + память

for i := 0 to High(Keep) do // + время инициализации

Keep[i] := False;

// И финальный сбор

for i := 0 to High(Points) do // + дополнительный проход

if Keep[i] then AddToResult;

Затраты на управление стеком

// Каждый Push/Pop - операции
// Медленнее, чем вызов функции (в некоторых случаях)
Stack.Push(...);  // Проверка capacity, увеличение если нужно
Segment := Stack.Pop;  // Проверка на пустоту

// Каждый Push/Pop - операции

// Медленнее, чем вызов функции (в некоторых случаях)

Stack.Push(...); // Проверка capacity, увеличение если нужно

Segment := Stack.Pop; // Проверка на пустоту

Потенциальные накладные расходы

// - В худшем случае (очень маленький epsilon) стек может расти до O(N)
// - Каждый сегмент разбивается на два
// - Нужно выделять много памяти

// - В худшем случае (очень маленький epsilon) стек может расти до O(N)

// - Каждый сегмент разбивается на два

// - Нужно выделять много памяти

Итеративный алгоритм RDP: Турбо-режим

В теории, итеративный алгоритм Рамера-Дугласа-Пекера должен работать быстрее рекурсивного. Но это в теории. У меня же изначально, на родном TStack<T>, он значительно отставал от рекурсивного. Затем я переделал работу со стеком. Этот вариант представлен выше. Стало быстрее, но рекурсивный всё равно побеждал чаще, чем хотелось бы.

Поэтому, меняем работу с данными внутри алгоритма.

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

type
  PSegment = ^TSegment;
  TSegment = record
    StartIdx, EndIdx: Integer;
  end;

type

PSegment = ^TSegment;

TSegment = record

StartIdx, EndIdx: Integer;

end;

Во-вторых, вся работа со стеком происходит напрямую. Вместо Push использую:

 // Stack.Push(TSegment.Create(0, MaxIdx));
 Stack[SegmentsCount].StartIdx := 0;
 Stack[SegmentsCount].EndIdx := MaxIdx;
 Inc(SegmentsCount);

// Stack.Push(TSegment.Create(0, MaxIdx));

Stack[SegmentsCount].StartIdx := 0;

Stack[SegmentsCount].EndIdx := MaxIdx;

Inc(SegmentsCount);

И вместо Pop:

 // Stack.Pop(Segment);
 Dec(SegmentsCount);
 Segment := @Stack[SegmentsCount];

// Stack.Pop(Segment);

Dec(SegmentsCount);

Segment := @Stack[SegmentsCount];

В-третьих, вся работа с массивами происходит только по указателям. У нас в алгоритме везде циклы от индекса до индекса. Поэтому, получив указатель на нужный элемент массива, мы можем его просто инкрементировать в цикле. Например, как это выглядит для инициализации массива и поиска дублей:

  // Поиск максимально удалённой вершины и удаление дублей
  P3 := @AList.List[0];
  P2 := @List[0];
  for i := 1 to AList.Count - 1 do
  begin
    Inc(P3);
    // Такая запись работает дольше
    // if EqualPts(P3^, P2^ {List[Len]}) then
    if IsZero(P3^.X-P2^.X) and IsZero(P3^.Y-P2^.Y) then
      continue;
    Inc(Len);
    Inc(P2); P2^ := P3^; // List[Len] := P3^;

    if AStrictOnePolyline and not AClosed then
      continue;
    dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);
    if dL > LMx then
    begin
      MaxIdx := Len;
      LMx := dL;
    end;
  end;

// Поиск максимально удалённой вершины и удаление дублей

P3 := @AList.List[0];

P2 := @List[0];

for i := 1 to AList.Count - 1 do

begin

Inc(P3);

// Такая запись работает дольше

// if EqualPts(P3^, P2^ {List[Len]}) then

if IsZero(P3^.X-P2^.X) and IsZero(P3^.Y-P2^.Y) then

continue;

Inc(Len);

Inc(P2); P2^ := P3^; // List[Len] := P3^;

if AStrictOnePolyline and not AClosed then

continue;

dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);

if dL > LMx then

begin

MaxIdx := Len;

LMx := dL;

end;

Листинг RDP-Turbo

function RamerDouglasPeuckerIterativeFast(AList: TPointFList;
  AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;
  out ALines: TPointFDynArray): Boolean;

type
  PSegment = ^TSegment;
  TSegment = record
    StartIdx, EndIdx: Integer;
  end;

var
  List: TPointFDynArray;
  Keep: array of Boolean;
  Stack: array of TSegment;
  SegmentsCount: Integer;
  Segment: PSegment;
  i, Len, MaxIdx: Integer;
  P1, P2, P3: PPointF;
  LMx, dx, dy, L, dL, A: Single;
  ResLength: Integer;
  EpsilonSqr: Single;
  B0: PBoolean;
begin
  // С пустыми данными не работаем
  Result := (AList <> nil) and (AList.Count > 2);
  if not Result then
    exit;

  EpsilonSqr := Sqr(AEpsilon);
  // Инициализация
  SetLength(Stack, AList.Count*2);
  SegmentsCount := 0;

  // Подготовка списка точек (удаление дубликатов)
  SetLength(List, AList.Count);
  Len := 0;
  P1 := @AList.List[0];
  List[0] := P1^;
  LMx := 0;
  MaxIdx := 0;

  // Поиск максимально удалённой вершины и удаление дублей
  P3 := @AList.List[0];
  P2 := @List[0];
  for i := 1 to AList.Count - 1 do
  begin
    Inc(P3);
    // Такая запись работает дольше
    // if EqualPts(P3^, P2^ {List[Len]}) then
    if IsZero(P3^.X-P2^.X) and IsZero(P3^.Y-P2^.Y) then
      continue;
    Inc(Len);
    Inc(P2); P2^ := P3^; // List[Len] := P3^;
    if AStrictOnePolyline and not AClosed then
      continue;
    dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);
    if dL > LMx then
    begin
      MaxIdx := Len;
      LMx := dL;
    end;
  end;
  if AStrictOnePolyline and not AClosed then
    MaxIdx := Len;

  Len := Len + 1;
  // Этого делать не будем, сейчас массив и так больше Len
  // А все обработки массива всё равно анализируют Len
  // SetLength(List, Len);

  // Если требуется закрытая ломаная, совмещаем начало и конец
  if AClosed and not EqualPts(List[Len - 1], List[0]) then
  begin
    SetLength(List, Len + 1);
    List[Len] := List[0];
    Len := Len + 1;
  end;

  // Массив для отметки точек, которые нужно сохранить
  SetLength(Keep, Len);
  // Этого делать тоже не будем, при установке
  // длины массива это происходит автоматически
  //  for i := 1 to Len - 2 do
  //    Keep[i] := False;

  // Всегда сохраняем первую, последнюю и
  // максимально удалённую точки
  Keep[0] := True;
  Keep[Len-1] := True;
  Keep[MaxIdx] := True;

  // Добавляем начальные сегменты в стек
  if MaxIdx > 0 then
  begin
    // Stack.Push(TSegment.Create(0, MaxIdx));
    Stack[SegmentsCount].StartIdx := 0;
    Stack[SegmentsCount].EndIdx := MaxIdx;
    Inc(SegmentsCount);
  end;

  if MaxIdx < Len - 1 then
  begin
    Stack[SegmentsCount].StartIdx := MaxIdx;
    Stack[SegmentsCount].EndIdx := Len - 1;
    Inc(SegmentsCount);
  end;

  // Обрабатываем сегменты пока стек не пуст
  while SegmentsCount > 0 do
  begin
    // Stack.Pop(Segment);
    Dec(SegmentsCount);
    Segment := @Stack[SegmentsCount];

    if Segment.EndIdx - Segment.StartIdx < 2 then
      continue;

    // Находим точку с максимальным отклонением
    LMx := 0;
    MaxIdx := Segment.StartIdx;
    P1 := @List[Segment.StartIdx];
    P2 := @List[Segment.EndIdx];

    dx := P2.X - P1.X;
    dy := P2.Y - P1.Y;
    L := dx*dx + dy*dy;
    if IsZero(L) then
      continue;

    A := P1.X*P2.Y - P2.X*P1.Y;

    P3 := @List[Segment.StartIdx + 1];
    for i := Segment.StartIdx + 1 to Segment.EndIdx - 1 do
    begin
      dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;
      if dL > LMx then
      begin
        LMx := dL;
        MaxIdx := i;
      end;
      Inc(P3);
    end;

    // Если максимальное отклонение превышает epsilon, обрабатываем дальше
    if (LMx > EpsilonSqr) and
      (MaxIdx > Segment.StartIdx) and (MaxIdx < Segment.EndIdx)
    then
    // Ниже ускорение через указатель на сегмент
    // begin
    //   Keep[MaxIdx] := True;
    //   Stack[SegmentsCount+1].StartIdx := MaxIdx;
    //   Stack[SegmentsCount+1].EndIdx := Segment.EndIdx;
    //   Segment.EndIdx := MaxIdx;
    //   Inc(SegmentsCount, 2);
    // end;    
    begin
      Keep[MaxIdx] := True;

      EndIdx := Segment.EndIdx;
      Segment.EndIdx := MaxIdx;

      Inc(Segment);
      Segment.StartIdx := MaxIdx;
      Segment.EndIdx := EndIdx;

      Inc(SegmentsCount, 2);
    end;

  end;

  // Формируем результирующий массив
  ResLength := 0;
  SetLength(ALines, Len); // Максимально возможный размер

  P3 := @List[0];
  P2 := @ALines[0];
  B0 := @Keep[0];
  for i := 0 to Len - 1 do
  begin
    // if Keep[i] then
    if B0^ then
    begin
      P2^ := P3^;
      Inc(P2);
      Inc(ResLength);
    end;
    Inc(P3);
    Inc(B0);
  end;

  // Для замкнутой ломаной проверяем, что первая и последняя точки не дублируются
  if AClosed and (ResLength > 1) and
    EqualPts(ALines[0], ALines[ResLength - 1]) then
      // Убираем дублирующую последнюю точку
      Dec(ResLength);

  // Корректируем длину массива
  SetLength(ALines, ResLength);

  // Проверить начальную точку для замкнутой ломаной
  if AClosed and ACheckFirst then
    while CheckFirstPoint(EpsilonSqr, ALines) do;
end;

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

function RamerDouglasPeuckerIterativeFast(AList: TPointFList;

AEpsilon: Single; AClosed, ACheckFirst, AStrictOnePolyline: Boolean;

out ALines: TPointFDynArray): Boolean;

type

PSegment = ^TSegment;

TSegment = record

StartIdx, EndIdx: Integer;

end;

var

List: TPointFDynArray;

Keep: array of Boolean;

Stack: array of TSegment;

SegmentsCount: Integer;

Segment: PSegment;

i, Len, MaxIdx: Integer;

P1, P2, P3: PPointF;

LMx, dx, dy, L, dL, A: Single;

ResLength: Integer;

EpsilonSqr: Single;

B0: PBoolean;

begin

// С пустыми данными не работаем

Result := (AList <> nil) and (AList.Count > 2);

if not Result then

exit;

EpsilonSqr := Sqr(AEpsilon);

// Инициализация

SetLength(Stack, AList.Count*2);

SegmentsCount := 0;

// Подготовка списка точек (удаление дубликатов)

SetLength(List, AList.Count);

Len := 0;

P1 := @AList.List[0];

List[0] := P1^;

LMx := 0;

MaxIdx := 0;

// Поиск максимально удалённой вершины и удаление дублей

P3 := @AList.List[0];

P2 := @List[0];

for i := 1 to AList.Count - 1 do

begin

Inc(P3);

// Такая запись работает дольше

// if EqualPts(P3^, P2^ {List[Len]}) then

if IsZero(P3^.X-P2^.X) and IsZero(P3^.Y-P2^.Y) then

continue;

Inc(Len);

Inc(P2); P2^ := P3^; // List[Len] := P3^;

if AStrictOnePolyline and not AClosed then

continue;

dL := Sqr(P3.X-P1.X) + Sqr(P3.Y-P1.Y);

if dL > LMx then

begin

MaxIdx := Len;

LMx := dL;

end;

if AStrictOnePolyline and not AClosed then

MaxIdx := Len;

Len := Len + 1;

// Этого делать не будем, сейчас массив и так больше Len

// А все обработки массива всё равно анализируют Len

// SetLength(List, Len);

// Если требуется закрытая ломаная, совмещаем начало и конец

if AClosed and not EqualPts(List[Len - 1], List[0]) then

begin

SetLength(List, Len + 1);

List[Len] := List[0];

Len := Len + 1;

end;

// Массив для отметки точек, которые нужно сохранить

SetLength(Keep, Len);

// Этого делать тоже не будем, при установке

// длины массива это происходит автоматически

// for i := 1 to Len - 2 do

// Keep[i] := False;

// Всегда сохраняем первую, последнюю и

// максимально удалённую точки

Keep[0] := True;

Keep[Len-1] := True;

Keep[MaxIdx] := True;

// Добавляем начальные сегменты в стек

if MaxIdx > 0 then

begin

// Stack.Push(TSegment.Create(0, MaxIdx));

Stack[SegmentsCount].StartIdx := 0;

Stack[SegmentsCount].EndIdx := MaxIdx;

Inc(SegmentsCount);

end;

if MaxIdx < Len - 1 then

begin

Stack[SegmentsCount].StartIdx := MaxIdx;

Stack[SegmentsCount].EndIdx := Len - 1;

Inc(SegmentsCount);

end;

// Обрабатываем сегменты пока стек не пуст

while SegmentsCount > 0 do

begin

// Stack.Pop(Segment);

Dec(SegmentsCount);

Segment := @Stack[SegmentsCount];

if Segment.EndIdx - Segment.StartIdx < 2 then

continue;

// Находим точку с максимальным отклонением

LMx := 0;

MaxIdx := Segment.StartIdx;

P1 := @List[Segment.StartIdx];

P2 := @List[Segment.EndIdx];

dx := P2.X - P1.X;

dy := P2.Y - P1.Y;

L := dx*dx + dy*dy;

if IsZero(L) then

continue;

A := P1.X*P2.Y - P2.X*P1.Y;

P3 := @List[Segment.StartIdx + 1];

for i := Segment.StartIdx + 1 to Segment.EndIdx - 1 do

begin

dL := Sqr(P3.Y*dx - P3.X*dy + A)/L;

if dL > LMx then

begin

LMx := dL;

MaxIdx := i;

end;

Inc(P3);

end;

// Если максимальное отклонение превышает epsilon, обрабатываем дальше

if (LMx > EpsilonSqr) and

(MaxIdx > Segment.StartIdx) and (MaxIdx < Segment.EndIdx)

then

// Ниже ускорение через указатель на сегмент

// begin

// Keep[MaxIdx] := True;

// Stack[SegmentsCount+1].StartIdx := MaxIdx;

// Stack[SegmentsCount+1].EndIdx := Segment.EndIdx;

// Segment.EndIdx := MaxIdx;

// Inc(SegmentsCount, 2);

// end;

begin

Keep[MaxIdx] := True;

EndIdx := Segment.EndIdx;

Segment.EndIdx := MaxIdx;

Inc(Segment);

Segment.StartIdx := MaxIdx;

Segment.EndIdx := EndIdx;

Inc(SegmentsCount, 2);

end;

// Формируем результирующий массив

ResLength := 0;

SetLength(ALines, Len); // Максимально возможный размер

P3 := @List[0];

P2 := @ALines[0];

B0 := @Keep[0];

for i := 0 to Len - 1 do

begin

// if Keep[i] then

if B0^ then

begin

P2^ := P3^;

Inc(P2);

Inc(ResLength);

end;

Inc(P3);

Inc(B0);

end;

// Для замкнутой ломаной проверяем, что первая и последняя точки не дублируются

if AClosed and (ResLength > 1) and

EqualPts(ALines[0], ALines[ResLength - 1]) then

// Убираем дублирующую последнюю точку

Dec(ResLength);

// Корректируем длину массива

SetLength(ALines, ResLength);

// Проверить начальную точку для замкнутой ломаной

if AClosed and ACheckFirst then

while CheckFirstPoint(EpsilonSqr, ALines) do;

end;

[свернуть]

Вот в таком виде итеративный алгоритм Рамера-Дугласа-Пекера стал быстрее своего рекурсивного родственника. В комментариях оставил примеры, что конкретно ускоряю.

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

Заключение

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

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

Скачать

Друзья, спасибо за внимание!

Надеюсь, материл будет полезен.

Исходник (zip) 107 Кб. Delphi XE 7

Исполняемый файл (zip) 961 Кб.

5 11 голоса

Рейтинг статьи

0 комментариев

Старые

Новые Популярные

Межтекстовые Отзывы

Посмотреть все комментарии