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

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

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

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

Лично мне алгоритм понадобился в своё время для векторизации замкнутых контуров, которые получаются после обработки образов. Поэтому, разговор будет с акцентом в замкнутые ломаные, полигоны.

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

Запуск

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

Скачать

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

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

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

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

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

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.

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

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

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

Вначале, небольшие вспомогательные типы и функции.

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

type
  // Динамический массив вещественных 2D-точек
  TPointFDynArray = array of TPointF;

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

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

const

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

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

type

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

TPointFDynArray = array of TPointF;

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

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

begin

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

end;

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

function IsEmptyPts(P: TPointF): Boolean;

begin

Result := EqualPts(P, EmptyPointF);

end;

Рекурсивная часть будет называться, как в Вики — DouglasPeucker.

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

$h_{a}={\frac {2S}{a}},$ где $S$ — площадь треугольника, $a$ — длина стороны треугольника, на которую опущена высота.

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

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

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

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

// Рекурсивный алгоритм
procedure DouglasPeucker(const APoints: TPointFDynArray;
  AStartIdx, AEndIdx: Integer; AEpsilon: 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];
  // Площадь треугольника:
  // S = 1/2 * |(X1-X)(Y2-Y)-(X2-X)(Y1-Y)| =
  //   = 1/2 * |Y*(X2-X1)-X*(Y2-Y1)+X1*Y2-X2*Y1|
  dx := P2.x - P1.x; // DX
  dy := P2.y - P1.y; // DY
  // Длина основания
  L :=  Hypot(dx, dy); // sqrt(dx*dx + dy*dy);
  // Неизменная часть формулы
  A := P1.X*P2.Y - P2.X*p1.Y; // X1*Y2-X2*Y1
  // Если концы отрезка не совпадают
  if not IsZero(L) then
    for i := AStartIdx+1 to AEndIdx-1 do
    begin
      // Вершина, от которой считаем высоту
      P3 := APoints[i];
      // Высота треугольника:
      // h = 2 * S/L = 2*(1/2*|Y*DX-X*DY+A|)/L
      dL := Abs(P3.Y*dx - P3.X*dy + A)/L;
      if dL > dMax then
      begin
        Idx := i;
        dMax := dL;
      end;
    end;
  // Нашлась выдающаяся точка
  if (dMax > AEpsilon) and (Idx>-1) then
  begin
    // Погнали рекурсию в диапазоне
    // от начала до этой вершины
    DouglasPeucker(APoints, AStartIdx, Idx,
      AEpsilon, AResLength, ARes);
    // Запоминаем выдающуюся точку
    ARes[AResLength] := APoints[Idx];
    // Инкремент длины массива результата
    Inc(AResLength);
    // Погнали рекурсию в диапазоне
    // от этой вершины до конца
    DouglasPeucker(APoints, Idx, AEndIdx,
      AEpsilon, AResLength, ARes);
  end;
end;

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

procedure DouglasPeucker(const APoints: TPointFDynArray;

AStartIdx, AEndIdx: Integer; AEpsilon: 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];

// Площадь треугольника:

// S = 1/2 * |(X1-X)(Y2-Y)-(X2-X)(Y1-Y)| =

// = 1/2 * |Y*(X2-X1)-X*(Y2-Y1)+X1*Y2-X2*Y1|

dx := P2.x - P1.x; // DX

dy := P2.y - P1.y; // DY

// Длина основания

L := Hypot(dx, dy); // sqrt(dx*dx + dy*dy);

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

A := P1.X*P2.Y - P2.X*p1.Y; // X1*Y2-X2*Y1

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

if not IsZero(L) then

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

begin

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

P3 := APoints[i];

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

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

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

if dL > dMax then

begin

Idx := i;

dMax := dL;

end;

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

if (dMax > AEpsilon) and (Idx>-1) then

begin

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

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

DouglasPeucker(APoints, AStartIdx, Idx,

AEpsilon, AResLength, ARes);

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

ARes[AResLength] := APoints[Idx];

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

Inc(AResLength);

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

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

DouglasPeucker(APoints, Idx, AEndIdx,

AEpsilon, AResLength, ARes);

end;

Запуск алгоритма

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

Так исторически сложилось, что у меня очень часто для хранения точек используется список:

TList<TPointF>

1	TList<TPointF>

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

Подготовка

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

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

В конце будет листинг полностью. Сейчас фрагментарно. Если нужна информация полностью, легче перейти непосредственно к листингу.

Итак, инициализируем массив List и ищем максимально удалённую точку в нём с индексом Idx. Исходные данные ломаной хранятся в списке AList.

  // Инициализируем индексы, длину будущего массива
  Idx := 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;
    // Находим текущий максимум
    H := sqr(p.X-Start.X) + sqr(p.Y-Start.Y);
    if H > LMx then
    begin
      Idx := Len;
      LMx := H;
    end;
  end;
  // Высчитываем конечную длину массива
  Len := Len+1;
  // Выставляем длину массива данных
  SetLength(List, Len);

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

Idx := 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;

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

H := sqr(p.X-Start.X) + sqr(p.Y-Start.Y);

if H > LMx then

begin

Idx := Len;

LMx := H;

end;

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

Len := Len+1;

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

SetLength(List, Len);

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

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

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

После инициализации массива допишем следующий фрагмент, совмещающий начало и конец ломаной.

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

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

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

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

begin

Len := Len+1;

SetLength(List, Len);

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

end;

Запуск

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

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

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

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

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

SetLength(ALines, Len);

ALines[0] := List[0];

Len := 1;

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

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

DouglasPeucker(List, 0, Idx, AEpsilon, Len, ALines);

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

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

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

if Len < Length(ALines) then

begin

ALines[Len] := List[Idx];

Inc(Len);

end;

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

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

if Idx < High(List) then

begin

DouglasPeucker(List, Idx, High(List),

AEpsilon, 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 (Len > 0) and

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

SetLength(ALines, Len);

Тут возникает нюанс второй.

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

Пишем дополнительную функцию, которая возвращает TRUE, когда нулевая вершина успешно убрана.

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

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

// и если необходимо - удалить

function CheckFirstPoint(AEpsilon: Single;

var ALines: TPointFDynArray): Boolean;

var

L,dx,dy: Single;

P1,P2,P3: TPointF;

begin

Result := Length(ALines) > 2;

if not Result then exit;

P3 := ALines[0];

P1 := ALines[1];

P2 := ALines[High(ALines)];

dx := P2.x - P1.x;

dy := P2.y - P1.y;

L := Hypot(dx, dy);

if IsZero(L) then exit(False);

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

Result := L < AEpsilon;

if Result then Delete(ALines,0,1);

end;

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

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

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

type
  // Динамический массив вещественных 2D-точек
  TPointFDynArray = array of TPointF;

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

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

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

// Рекурсивный алгоритм
procedure DouglasPeucker(const APoints: TPointFDynArray;
  AStartIdx, AEndIdx: Integer; AEpsilon: 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];
  // Площадь треугольника:
  // S = 1/2 * |(X1-X)(Y2-Y)-(X2-X)(Y1-Y)| =
  //   = 1/2 * |Y*(X2-X1)-X*(Y2-Y1)+X1*Y2-X2*Y1|
  dx := P2.x - P1.x; // DX
  dy := P2.y - P1.y; // DY
  // Длина основания
  L :=  Hypot(dx, dy); // sqrt(dx*dx + dy*dy);
  // Неизменная часть формулы
  A := P1.X*P2.Y - P2.X*p1.Y; // X1*Y2-X2*Y1
  // Если концы отрезка не совпадают
  if not IsZero(L) then
    for i := AStartIdx+1 to AEndIdx-1 do
    begin
      // Вершина, от которой считаем высоту
      P3 := APoints[i];
      // Высота треугольника:
      // h = 2 * S/L = 2*(1/2*|Y*DX-X*DY+A|)/L
      dL := Abs(P3.Y*dx - P3.X*dy + A)/L;
      if dL > dMax then
      begin
        Idx := i;
        dMax := dL;
      end;
    end;
  // Нашлась выдающаяся точка
  if (dMax > AEpsilon) and (Idx>-1) then
  begin
    // Погнали рекурсию в диапазоне
    // от начала до этой вершины
    DouglasPeucker(APoints, AStartIdx, Idx,
      AEpsilon, AResLength, ARes);
    // Запоминаем выдающуюся точку
    ARes[AResLength] := APoints[Idx];
    // Инкремент длины массива результата
    Inc(AResLength);
    // Погнали рекурсию в диапазоне
    // от этой вершины до конца
    DouglasPeucker(APoints, Idx, AEndIdx,
      AEpsilon, AResLength, ARes);
  end;
end;

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

// Подготовить данные и вызвать алгоритм
function RamerDouglasPeucker(AList: TList<TPointF>;
  AEpsilon: Single; AClosed: Boolean; ACheckFirst: Boolean;
  out ALines: TPointFDynArray): Boolean;
var
  List: TPointFDynArray;
  Len,Idx,i: Integer;
  LMx,H: Single;
  P,Start: TPointF;
begin
  // С пустыми данными не работаем
  Result := (AList<>nil) and (AList.Count>0);
  if not Result then exit;
  // Инициализируем индексы, длину будущего массива
  Idx := 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;
    // Находим текущий максимум
    H := sqr(p.X-Start.X) + sqr(p.Y-Start.Y);
    if H > LMx then
    begin
      Idx := Len;
      LMx := H;
    end;
  end;
  // Высчитываем конечную длину массива
  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, Idx, AEpsilon, Len, ALines);
  // Дописываем последнюю точку диапазона
  // Которая одновременно будет являться начальной
  // для возможной следующей ломаной
  if Len < Length(ALines) then
  begin
    ALines[Len] := List[Idx];
    Inc(Len);
  end;
  // Пошла рекурсия в диапазоне от
  // максимально удалённой до конца
  if Idx < High(List) then
  begin
    DouglasPeucker(List, Idx, High(List),
      AEpsilon, 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 ACheckFirst and AClosed then
    while CheckFirstPoint(AEpsilon,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

type

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

TPointFDynArray = array of TPointF;

const

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

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

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

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

begin

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

end;

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

function IsEmptyPts(P: TPointF): Boolean;

begin

Result := EqualPts(P, EmptyPointF);

end;

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

procedure DouglasPeucker(const APoints: TPointFDynArray;

AStartIdx, AEndIdx: Integer; AEpsilon: 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];

// Площадь треугольника:

// S = 1/2 * |(X1-X)(Y2-Y)-(X2-X)(Y1-Y)| =

// = 1/2 * |Y*(X2-X1)-X*(Y2-Y1)+X1*Y2-X2*Y1|

dx := P2.x - P1.x; // DX

dy := P2.y - P1.y; // DY

// Длина основания

L := Hypot(dx, dy); // sqrt(dx*dx + dy*dy);

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

A := P1.X*P2.Y - P2.X*p1.Y; // X1*Y2-X2*Y1

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

if not IsZero(L) then

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

begin

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

P3 := APoints[i];

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

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

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

if dL > dMax then

begin

Idx := i;

dMax := dL;

end;

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

if (dMax > AEpsilon) and (Idx>-1) then

begin

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

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

DouglasPeucker(APoints, AStartIdx, Idx,

AEpsilon, AResLength, ARes);

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

ARes[AResLength] := APoints[Idx];

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

Inc(AResLength);

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

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

DouglasPeucker(APoints, Idx, AEndIdx,

AEpsilon, AResLength, ARes);

end;

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

// и если необходимо - удалить

function CheckFirstPoint(AEpsilon: Single;

var ALines: TPointFDynArray): Boolean;

var

L,dx,dy: Single;

P1,P2,P3: TPointF;

begin

Result := Length(ALines) > 2;

if not Result then exit;

P3 := ALines[0];

P1 := ALines[1];

P2 := ALines[High(ALines)];

dx := P2.x - P1.x;

dy := P2.y - P1.y;

L := Hypot(dx, dy);

if IsZero(L) then exit(False);

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

Result := L < AEpsilon;

if Result then Delete(ALines,0,1);

end;

// Подготовить данные и вызвать алгоритм

function RamerDouglasPeucker(AList: TList<TPointF>;

AEpsilon: Single; AClosed: Boolean; ACheckFirst: Boolean;

out ALines: TPointFDynArray): Boolean;

var

List: TPointFDynArray;

Len,Idx,i: Integer;

LMx,H: Single;

P,Start: TPointF;

begin

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

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

if not Result then exit;

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

Idx := 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;

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

H := sqr(p.X-Start.X) + sqr(p.Y-Start.Y);

if H > LMx then

begin

Idx := Len;

LMx := H;

end;

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

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, Idx, AEpsilon, Len, ALines);

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

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

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

if Len < Length(ALines) then

begin

ALines[Len] := List[Idx];

Inc(Len);

end;

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

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

if Idx < High(List) then

begin

DouglasPeucker(List, Idx, High(List),

AEpsilon, 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 ACheckFirst and AClosed then

while CheckFirstPoint(AEpsilon,ALines) do;

end;

Скачать

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

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

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

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

5 7 голоса

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

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

Старые

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

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

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