Расстояние от точки до кривой Безье(сплайна)

Общие вопросы программирования, алгоритмы и т.п.

Модератор: Модераторы

Расстояние от точки до кривой Безье(сплайна)

Сообщение Pavia » 05.01.2020 20:44:51

Подскажите алгоритм нахождения минимального расстояние от точки до кривой Безье(сплайна)? Для простоты сплайн 3 порядка на 4 точках.

Собственно отдельно интересно как это решено у zub в его редакторе.
Аватара пользователя
Pavia
постоялец
 
Сообщения: 290
Зарегистрирован: 07.01.2011 12:46:51

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение Дож » 06.01.2020 06:13:20

Обозначим четыре точки сплайна за p0=(x0,y0),p1=(x1,y1),p2=(x2,y2),p3=(x3,y3), а заданную точку за p=(x,y). Если 0<=t<=1 -- параметр спалйна, то вектор от заданной точки до точки t на сплайне задаётся формулой

f1.gif
f1.gif (1.23 КБ) Просмотров: 10567


Задача сводится к минимизации значения длины этого вектора. Вычислим длину вектора и сгруппируем их в многочлен по t, получится многочлен 6ой степени:

f2.gif
f2.gif (1.38 КБ) Просмотров: 10567


Ki вычисляются по некоторым простым (без ветвлений) формулам от изначальных точек. K6>0 всегда (если это действительно невырожденный кубический сплайн). Там, где многочлен достигает минимума, его производная равна нулю. Вычислим производную:

f3.gif
f3.gif (1.48 КБ) Просмотров: 10567


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

Как вариант, подойдёт этот (численный) алгоритм:
https://habr.com/ru/post/303342/

Корни дадут от 1 до 5 кандидатов ближайшей точки на сплайне, плюс к ним нужно добавить точки t=0 и t=1, и остаётся просто перебором найти ту, на которой f(t) минимально.
Аватара пользователя
Дож
энтузиаст
 
Сообщения: 899
Зарегистрирован: 12.10.2008 16:14:47

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение Vadim » 06.01.2020 09:23:33

Пора писать статью по методам оптимизации... :D
Vadim
долгожитель
 
Сообщения: 4112
Зарегистрирован: 05.10.2006 08:52:59
Откуда: Красноярск

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение Pavia » 06.01.2020 12:29:25

Как вариант, подойдёт этот (численный) алгоритм:

Таким смысла не у там дихотомия. Да и корни только действительные.
Я думаю взять алгоритм от сюда
http://mathworld.wolfram.com/PolynomialRoots.html
Изображение
Изображение
Правда не уверен что у них верное описание.

Я вот думаю для Безье известно что она почти точно представляется 16 линиями с шагом t=1/16 может так код быстрее будет?
Аватара пользователя
Pavia
постоялец
 
Сообщения: 290
Зарегистрирован: 07.01.2011 12:46:51

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение Alex2013 » 06.01.2020 12:49:49

В виде части "мозгового штурма"...
Почему нельзя просто прогнать все видимые точки сплайна на предмет расстояния до точки? Минимальное и будет искомым "расстоянием до сплайна" . :idea: Оптимизация методом постепенного приближения тоже понятна .
Последний раз редактировалось Alex2013 07.01.2020 03:50:40, всего редактировалось 2 раз(а).
Alex2013
долгожитель
 
Сообщения: 2922
Зарегистрирован: 03.04.2013 11:59:44

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение zub » 06.01.2020 15:43:14

У меня сплайны поддерживаются "пока" только формально. В кавычках потому что нет ничего более постоянного чем чтото временное))
безье в зкаде нет, есть нурбс. Но он только рисуется (с помощью glu) и редактируется по контролным точкам, привязок на него и операций типа разрезать - нет
В свое время интересовался этим, но так и не разобрался с математикой - забросил((. Мне для зкада нужна длина сплайна, ближайшая точка на сплайне и параметрическая точка на сплайне p(t) - это для базовой поддержки. Для резки\склейки соответственно неплохо иметь резку и склейку)) Буду рад любой помощи
zub
долгожитель
 
Сообщения: 2884
Зарегистрирован: 14.11.2005 23:51:26

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение olegy123 » 09.01.2020 10:01:41

лучше искать по английски, "spline Near point"
https://borjaportugal.com/2018/01/13/fi ... to-spline/
olegy123
долгожитель
 
Сообщения: 1643
Зарегистрирован: 25.02.2016 12:10:20

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение Mirage » 18.01.2020 00:29:37

zub писал(а):В свое время интересовался этим, но так и не разобрался с математикой - забросил((. Мне для зкада нужна длина сплайна, ближайшая точка на сплайне и параметрическая точка на сплайне p(t) - это для базовой поддержки. Для резки\склейки соответственно неплохо иметь резку и склейку)) Буду рад любой помощи


Точка считается по формуле P = At^3 + Bt^2 + Ct + D.
Ближайшая точка брутфорсом проще всего. Думаю достаточно будет.
Длина им же.
Резка/склейка специальным алгоритмом. Довольно простым насколько помню, по крайней мере для безье.
Mirage
энтузиаст
 
Сообщения: 881
Зарегистрирован: 06.05.2005 20:29:07
Откуда: Russia

Re: Расстояние от точки до кривой Безье(сплайна)

Сообщение MylnikovDm » 02.03.2020 15:50:42

Держите готовые проверенные функции для вычисления точек кривой Безье с заданным количеством шагов для 2 и 3 порядка.
Код: Выделить всё
interface

type

  TDblPoint = packed record
    x, y: double;
  end;
  PDblPoint = ^TDblPoint;

//расчёт точек для кривой Безье
//если количество точек aStep указано равным 0, либо оно меньше количества точек в aPoints
//то количество точек определяется из параметров массива aPoints
//первая и последняя точка включаются в состав точек и указанное количество
//то есть, если aStep = 5, то будет добавлено 3 точки внутри кривой, плюс первая и
//последняя опорные точки
procedure Bezie2Points(aP1, aP2, aP3: TPoint; var aPoints: array of TPoint; aStep: integer = 0); overload;
procedure Bezie2Points(aP1, aP2, aP3: TDblPoint; var aPoints: array of TDblPoint; aStep: integer = 0); overload;

procedure Bezie3Points(aP1, aP2, aP3, aP4: TPoint; var aPoints: array of TPoint; aStep: integer = 0); overload;
procedure Bezie3Points(aP1, aP2, aP3, aP4: TDblPoint; var aPoints: array of TDblPoint; aStep: integer = 0); overload;

implementation

procedure Bezie3Points(aP1, aP2, aP3, aP4: TPoint; var aPoints: array of TPoint; aStep: integer = 0);
const
  MsInt = 4;
var
  i, SegmCount: integer;
  pl1, pl2, pl3, tp1, tp2, tp3, bp1, bp2: TPoint;
begin
  SegmCount := Length(aPoints);
  if (aStep = 0) or (aStep > SegmCount) then begin
    aStep := SegmCount;
    if aStep < 2 then exit;
  end; // if

  aStep -= 2;
  SegmCount := aStep + 1;

  aPoints[0].X := aP1.x;
  aPoints[0].Y := aP1.Y;

  if aStep < 1 then begin
    aPoints[1].X := aP4.x;
    aPoints[1].Y := aP4.Y;
    Exit;
  end; // if
  aPoints[SegmCount].X := aP4.x;
  aPoints[SegmCount].Y := aP4.Y;

  aP1.x := aP1.x shl MsInt;
  aP1.y := aP1.y shl MsInt;
  aP2.x := aP2.x shl MsInt;
  aP2.y := aP2.y shl MsInt;
  aP3.x := aP3.x shl MsInt;
  aP3.y := aP3.y shl MsInt;
  aP4.x := aP4.x shl MsInt;
  aP4.y := aP4.y shl MsInt;

  pl1.x := aP2.x - aP1.x;
  pl1.y := aP2.y - aP1.y;

  pl2.x := aP3.x - aP2.x;
  pl2.y := aP3.y - aP2.y;

  pl3.x := aP4.x - aP3.x;
  pl3.y := aP4.y - aP3.y;

  for i := 1 to aStep do begin
    tp1.x := aP1.x + ((pl1.x * i) div SegmCount);
    tp1.y := aP1.y + ((pl1.y * i) div SegmCount);
    tp2.x := aP2.x + ((pl2.x * i) div SegmCount);
    tp2.y := aP2.y + ((pl2.y * i) div SegmCount);
    tp3.x := aP3.x + ((pl3.x * i) div SegmCount);
    tp3.y := aP3.y + ((pl3.y * i) div SegmCount);

    bp1.x := tp1.x + (((tp2.x - tp1.x) * i) div SegmCount);
    bp1.y := tp1.Y + (((tp2.y - tp1.y) * i) div SegmCount);
    bp2.x := tp2.x + (((tp3.x - tp2.x) * i) div SegmCount);
    bp2.y := tp2.Y + (((tp3.y - tp2.y) * i) div SegmCount);


    aPoints[i].X := (bp1.x + (((bp2.x - bp1.x) * i) div SegmCount)) shr MsInt;
    aPoints[i].Y := (bp1.Y + (((bp2.y - bp1.y) * i) div SegmCount)) shr MsInt;
  end; // for
end;

procedure Bezie3Points(aP1, aP2, aP3, aP4: TDblPoint; var aPoints: array of TDblPoint; aStep: integer = 0);
var
  i, SegmCount: integer;
  pl1, pl2, pl3, tp1, tp2, tp3, bp1, bp2: TDblPoint;
begin
  SegmCount := Length(aPoints);
  if (aStep = 0) or (aStep > SegmCount) then begin
    aStep := SegmCount;
    if aStep < 2 then exit;
  end;

  aStep -= 2;
  SegmCount := aStep + 1;

  aPoints[0].X := aP1.x;
  aPoints[0].Y := aP1.Y;

  if aStep < 1 then begin
    aPoints[1].X := aP4.x;
    aPoints[1].Y := aP4.Y;
    Exit;
  end; // if
  aPoints[SegmCount].X := aP4.x;
  aPoints[SegmCount].Y := aP4.Y;

  pl1.x := aP2.x - aP1.x;
  pl1.y := aP2.y - aP1.y;

  pl2.x := aP3.x - aP2.x;
  pl2.y := aP3.y - aP2.y;

  pl3.x := aP4.x - aP3.x;
  pl3.y := aP4.y - aP3.y;

  for i := 1 to aStep do begin
    tp1.x := aP1.x + ((pl1.x * i) / SegmCount);
    tp1.y := aP1.y + ((pl1.y * i) / SegmCount);
    tp2.x := aP2.x + ((pl2.x * i) / SegmCount);
    tp2.y := aP2.y + ((pl2.y * i) / SegmCount);
    tp3.x := aP3.x + ((pl3.x * i) / SegmCount);
    tp3.y := aP3.y + ((pl3.y * i) / SegmCount);

    bp1.x := tp1.x + (((tp2.x - tp1.x) * i) / SegmCount);
    bp1.y := tp1.Y + (((tp2.y - tp1.y) * i) / SegmCount);
    bp2.x := tp2.x + (((tp3.x - tp2.x) * i) / SegmCount);
    bp2.y := tp2.Y + (((tp3.y - tp2.y) * i) / SegmCount);


    aPoints[i].X := bp1.x + (((bp2.x - bp1.x) * i) / SegmCount);
    aPoints[i].Y := bp1.Y + (((bp2.y - bp1.y) * i) / SegmCount);
  end; // for
end;

procedure Bezie2Points(aP1, aP2, aP3: TPoint; var aPoints: array of TPoint; aStep: integer = 0);
const
  MsInt = 4;
var
  i, SegmCount: integer;
  pl1, pl2, tp1, tp2: TPoint;
begin
  SegmCount := Length(aPoints);
  if (aStep = 0) or (aStep > SegmCount) then begin
    aStep := SegmCount;
    if aStep < 2 then exit;
  end; // if

  aStep -= 2;
  SegmCount := aStep + 1;

  aPoints[0].X := aP1.x;
  aPoints[0].Y := aP1.Y;

  if aStep < 1 then begin
    aPoints[1].X := aP3.x;
    aPoints[1].Y := aP3.Y;
    Exit;
  end; // if
  aPoints[SegmCount].X := aP3.x;
  aPoints[SegmCount].Y := aP3.Y;

  pl1.x := (aP2.x - aP1.x) shl MsInt;
  pl1.y := (aP2.y - aP1.y) shl MsInt;

  pl2.x := (aP3.x - aP2.x) shl MsInt;
  pl2.y := (aP3.y - aP2.y) shl MsInt;

  for i := 1 to aStep do begin
    tp1.x := aP1.x + ((pl1.x * i) div SegmCount);
    tp1.y := aP1.y + ((pl1.y * i) div SegmCount);
    tp2.x := aP2.x + ((pl2.x * i) div SegmCount);
    tp2.y := aP2.y + ((pl2.y * i) div SegmCount);

    aPoints[i].X := (tp1.x + (((tp2.x - tp1.x) * i) div SegmCount)) shr MsInt;
    aPoints[i].Y := (tp1.Y + (((tp2.y - tp1.y) * i) div SegmCount)) shr MsInt;
  end; // for
end;

procedure Bezie2Points(aP1, aP2, aP3: TDblPoint; var aPoints: array of TDblPoint; aStep: integer = 0);
var
  i, SegmCount: integer;
  pl1, pl2, tp1, tp2: TDblPoint;
begin
  SegmCount := Length(aPoints);
  if (aStep = 0) or (aStep > SegmCount) then begin
    aStep := SegmCount;
    if aStep < 2 then exit;
  end;

  aStep -= 2;
  SegmCount := aStep + 1;

  aPoints[0].X := aP1.x;
  aPoints[0].Y := aP1.Y;

  if aStep < 1 then begin
    aPoints[1].X := aP3.x;
    aPoints[1].Y := aP3.Y;
    Exit;
  end; // if
  aPoints[SegmCount].X := aP3.x;
  aPoints[SegmCount].Y := aP3.Y;

  pl1.x := aP2.x - aP1.x;
  pl1.y := aP2.y - aP1.y;

  pl2.x := aP3.x - aP2.x;
  pl2.y := aP3.y - aP2.y;

  for i := 1 to aStep do begin
    tp1.x := aP1.x + ((pl1.x * i) / SegmCount);
    tp1.y := aP1.y + ((pl1.y * i) / SegmCount);
    tp2.x := aP2.x + ((pl2.x * i) / SegmCount);
    tp2.y := aP2.y + ((pl2.y * i) / SegmCount);

    aPoints[i].X := tp1.x + (((tp2.x - tp1.x) * i) / SegmCount);
    aPoints[i].Y := tp1.Y + (((tp2.y - tp1.y) * i) / SegmCount);
  end; // for
end;

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

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

Zub'у для ЗКада тоже может пригодиться. Дарю! :)
MylnikovDm
постоялец
 
Сообщения: 103
Зарегистрирован: 15.02.2007 21:26:10
Откуда: Челябинск


Вернуться в Общее

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и гости: 10

Рейтинг@Mail.ru
cron