Поиск кротчайшего расстояния между точкой и отрезком

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

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

Уравнение прямой y=kx+b
условие перпедикулярности k1=-1/k2


Вот а в палне реализации даже не знаю с чего начинать........
Синтаксис ASM-a относительно знаю, но практика практически нулевая - не было времени практиковаться - работа всё занимала ...
HELP!!!!!!!!!!!!!!!! Скоро сдавать!!!

 

Вроде, всё просто:

(x1,y1) - координаты начала отрезка.
(x2,y2) - координаты конца отрезка.
(x3,y3) - координаты точки.

Строим уравнение прямой отрезка (y = kx+b):
k = dy/dx = (y2-y1)/(x2-x1)
b = y-kx = y2-kx2 = y1-kx1

Строим уравнение перпендикуляра (y = Kx+B):
K = -1/k
B = y3-Kx3

Находим точку пересечения (X,Y):
kX+b = KX+B
kX-KX = B-b
X = (B-b)/(k-K)
Y = kX+b = KX+B

Вычисляем расстояния:
d1 : (x3,y3) : (X,Y)
d2 : (x3,y3) : (x1,y1)
d3 : (x3,y3) : (x2,y2)

Чтобы найти dn нужно решить прямоугольный треугольник, т.е. сначала найти угол по формуле тангенса ( |dx|/|dy| ), потом по акрсинусу или арккосинусу вычислить гипотенузу, которая и является расстоянием dn.

Принадлежит ли точка пересечения отрезку?
if(|X-x1| > |x2-x1|){ d1 = +INF }

Возвращаем min(d1,d2,d3)

В плане реализации нужно почитать про FPU, т.к. целочисленной арифметикой тут не обойтись. Если лень разбираться с FPU, то можно поступить проще:
Пишем всё на C
В конфигурации проекта на вижуале указываем, чтобы инлайнились интринсеки (ftol(), fsin() и т.д.) Для этого есть ключик в настройках (см. msdn)
Включаем максимальную оптимизацию
Включаем создание листинга на ассемблере.

Получаем готовое решение на ассемблере.

 

Quantum

А почему нельзя просто по теореме Пифагора sqrt((y1-y2)^2+(x1-x2)^2) ?

Не пойму я что-то логику Во-первых зачем вообще d1 = +INF (Добавлено: это понял ), во-вторых что делает этот if ?

 

Stiver

fsqrt, вроде, тормозит сильно...

Чтобы min(d1,d2,d3) не возвращало d1, если точка не принадлежит отрезку.

Просто проверяет принадлежит ли (X,Y) заданному отрезку. По крайней мере так задумано.

 

Quantum

Хм, если x1=2, x2=0, X=3, то вроде не сработает. Я бы проверял if(|X-x1|+|X-x2|>|x1-x2|) Наверняка можно и как-нибудь красивее
Добавлено: например if((x1-X)*(X-x2)<0)

 

Метод который работает при любом расположении точек:

A=(xa,ya) ;начало отрезка
B=(xb,yb) ;конец отрезка
C=(xc,yc) ;точка

P=(C-A)*(B-A) ;скалярное произведение двух векторов
L=(B-A)*(B-A)
Если 0>=P то R=|C-A| ;длина вектора
Если P>=L то R=|C-B|
Если 0<P<L то R=|C-A-P/L*(B-A)|

 

Это на каком компиляторе так моно?
Сам я С пишу на Turbo C++ 3.0, а С++ на Borland-e 6.

Black_mirror, можешь по-подробней объяснить свой алгоритм, насколько я знаю...
скалярное произведение это ещё и умножение на Cos угла между ними...

 

Quantum

Латентность fsqrt примерно равна латентности деления

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

Код (Text):

1. type
2.   pPointDouble = ^TPointDouble;
3.   TPointDouble = record X,Y:double; end;
4.  
5. //PSDistance - минимальное расстояние от точки до отрезка
6. {Result = 0 - точка проецируется на отрезок: Dist = расст.от P до прямой SP1-SP2
7.           1 - точка P лежит "левее"  SP1: Dist:=PPDistance(P,SP1)
8.           2 - точка P лежит "правее" SP2: Dist:=PPDistance(P,SP2)
9. function PSDistance(const P,SP1,SP2:TPointDouble;var Distance:double):integer;
10. begin
11.   Result:=PSDistanceEx(P,SP1,SP2,Distance,Nil)
12. end;
13.  
14. //PSDistanceEx - тоже с расчетом точки проекции PP
15. если PP <> Nil, то дополнительно расчитывает PP
16. Result    PP
17.   0       проекция точки P на отрезок SP1-SP2
18.   1       SP1
19.   2       SP2 }
20. function PSDistanceEx(const P,SP1,SP2:TPointDouble;var Distance:double;
21.                       PP:pPointDouble):integer;
22. var
23.   dX21,dY21,dX1,dY1,dX2,dY2,aSide,d,L2:double;
24.   d,L2:extended;
25. begin
26.   dX21:=SP2.X-SP1.X;
27.   dY21:=SP2.Y-SP1.Y;
28.   dX1:=P.X-SP1.X;
29.   dY1:=P.Y-SP1.Y;
30.   //если угол между векторами (SP1->SP2) и (SP1->P)  > +-90 градусов
31.   //=> P лежит "левее" точки SP1 => мин.Dist = расст.до точки SP1
32.   aSide:=dX21*dX1+dY21*dY1;
33.   if aSide <= 0 then
34.   begin
35.     if aSide = 0 then
36.       Result:=0
37.     else
38.       Result:=1;
39.     Distance:=sqrt(dX1*dX1+dY1*dY1);
40.     if PP <> Nil then
41.       PP^:=SP1;
42.   end
43.   else
44.   begin
45.     dX2:=P.X-SP2.X;
46.     dY2:=P.Y-SP2.Y;
47.     //если угол между векторами (SP2->SP1) и (SP2->P)  > +-90 градусов
48.     //=> P лежит "правее" точки SP2 => мин.Dist = расст.до точки SP2
49.     aSide:=dX21*dX2+dY21*dY2;
50.     if aSide >= 0 then  //= ((-dX21)*dX2+(-dY21)*dY2 < 0)
51.     begin
52.       if aSide = 0 then
53.         Result:=0
54.       else
55.         Result:=2;
56.       Distance:=sqrt(dX2*dX2+dY2*dY2);
57.       if PP <> Nil then
58.         PP^:=SP2;
59.     end
60.     else
61.     //оба угла <= 90гр. => точка проецируется на отрезок (SP1-SP2)
62.     //=> мин.расст = расстоянию до прямой SP1-SP2
63.     begin
64.       d:=dY21*dX1-dX21*dY1;       //ненормированное расстояние до отрезка со знаком
65.       L2:=dX21*dX21+dY21*dY21;    //квадрат длины отрезка
66.       Distance:=abs(d)/sqrt(L2);  //расстояние
67.       if PP <> Nil then
68.       begin
69.         d:=d/L2;
70.         PP^.X:=P.X-dY21*d;
71.         PP^.Y:=P.Y+dX21*d;
72.       end;
73.       Result:=0;
74.     end;
75.   end;
76. end;

PS: Опять я долго копался - Black_mirror уже указал путь истинный ))

 

Формула из курса аналитической геометрии:
расстояние от точки x0, y0 до прямой с каноническим уравнением Ax+By+C=0 r = abs(A*x0 + B*y0 + C) / sqrt(A^2+B^2+C^2)

 

Ух ты!!! немного поспал - и всё стало ясно...
Black_mirror, leo упростили в несколько раз весь алгоритм....
Теперь ещё больший соблазн возникает оформить всё это в Asm-e (конечно использование FPU немного настораживает, но всё-таки)
Буду втыкать...

 

triam

Ну тогда немного поясню только третий пункт:
R=|C-A-P/L*(B-A)|=|C-(A+P/L*(B-A))|
(A+P/L*(B-A)) - точка на отрезке, в которую проецируется точка C.

 

ОФФ Блин я в своё время в институте решал такую задачку на Паскале. Только у меня в условии была не прямая а треугольник и произвольно расположенная точка и надо было найти расстояние до ближайшей стороны треугольника ://

 

Black_mirror, да да - это я понял...

 

Ну, вот и я проснулся

triam

Вижуал: 6й, 7й и т.д. - все поддерживают эти опции.

leo

Хмм, точно. Я почему-то считал, что эта инструкция тормозит больше всех.

Stiver

Ага, переоптимизировал я вчера эту формулу - вот и не работает.

 

Quantum, Ладно надо будет попробовать с вижуалом...

Но вот поделюсь до чего я дошёл:

Стараниями Black_mirror мы имеем выражение P=(C-A)*(B-A)
т.е. еслт более подробно расписать получается:
P=(cx-ax)(bx-ax)+(cy-ay)(by-ay) // в результате должны получить по идее число, потому как скалярное проиведение...

Но учитывая тот факт, что мы иммем дело с FPU, мне кажется, это выражение прийдётся немного переписать а именно

подогнать под стандарт "обратной польской нотаци" оно же RPN.... Про этого зверя я знаю пока мало, но знаю, что

операторы указываются после своих аргументов...
следовательно по-моему выражение выше можно переписать так:

cx ax - bx ax -* cy ay - by ay -* +

И вот вытекающий из этого и моей больной фантазии фрагмент кода:

Код (Text):

1.    .model small
2.    .stack 100h
3.    .code
4.  
5. start  proc near
6.  
7.       xc  dw ?
8.       yc  dw ?
9.       xa  dw ?
10.       xb  dw ?
11.       ya  dw ?
12.       yb  dw ?
13.  
14.    
15.  fld yc      ; priemnik pomechaem v st(0)
16.  fsub ya     ;  priemnik - istochnik
17.  fld1 yb     ; priemnik pomechaem v st(1)
18.  fsub ya     ;  priemnik - istochnik
19.  fmul        ; umnozaem
20.  fld2 xc      ; priemnik pomechaem v st(2)
21.  fsub xa     ;  priemnik - istochnik
22.  fld3 xb     ; priemnik pomechaem v st(3)
23.  fsub xa     ;  priemnik - istochnik
24.  fmul        ; umnozaem
25.  fadd         ; skladivaem
26.  ret

Чувствую наламерил жёстко, но надо же с чего-то начинать )
Жду ваших наставлений!!!

 

triam
В обратной польской нотации выражение записано правильно, а вот код не совсем.
Инструкций fld2,fld3 не существует. Инструкция fld1 загружает на вершину стека FPU-регистров единицу. Инструкция fld перемещает st(i) в st(i+1), а в st(0) загружает указанный операнд. Загружать можно также один из регистров сопроцессора: fld st(i).
Так как чисел с плавающей точкой размером в одно слово не существует (float занимает двойное слово), то нужно либо dw заменить на dd, либо использовать инструкции fild и fisub.
При объявлении чисел с плавающей точкой, точку нужно писать обязательно.

А вот и сам код, в комментариях состояние стека сопроцессора, слева st0:

Код (Text):

1. minR:
2.         fld [xc]
3.         fsub [xa]
4.         fld [yc]
5.         fsub [ya]       ;yc-ya xc-xa
6.  
7.         fld [xb]
8.         fsub [xa]
9.         fld [yb]
10.         fsub [ya]       ;yb-ya xb-xa  yc-ya xc-xa
11.  
12.         fld st1         ;xb-xa  yb-ya xb-xa  yc-ya xc-xa
13.         fmul st0,st4    ;(xb-xa)*(xc-xa)  yb-ya xb-xa  yc-ya xc-xa
14.         fld st1         ;yb-ya (xc-xa)*(xb-xa)  yb-ya xb-xa  yc-ya xc-xa
15.         fmul st0,st4    ;(yb-ya)*(yc-ya) (xc-xa)*(xb-xa)  yb-ya xb-xa  yc-ya xc-xa
16.         faddp           ;P  yb-ya xb-xa  yc-ya xc-xa
17.        
18.         fldz            ;0.0 P  yb-ya xb-xa  yc-ya xc-xa
19.         fcomp           ;P  yb-ya xb-xa  yc-ya xc-xa
20.         fstsw ax
21.         sahf
22.         jb .l1
23.  
24. ;0>=P
25.         fstp st0
26.         fstp st0
27.         fstp st0        ;yc-ya xc-xa
28.         jmp .l3
29.  
30. ;0<P
31.     .l1:                ;P  yb-ya xb-xa  yc-ya xc-xa
32.         fld st2         ;xb-xa  P  yb-ya xb-xa  yc-ya xc-xa
33.         fmul st0,st3    ;(xb-xa)^2  P  yb-ya xb-xa  yc-ya xc-xa
34.         fld st2         ;yb-ya xb-xa  P  yb-ya xb-xa  yc-ya xc-xa
35.         fmul st0,st3    ;(yb-ya)^2 (xb-xa)^2  P  yb-ya xb-xa  yc-ya xc-xa
36.         faddp           ;L P  yb-ya xb-xa  yc-ya xc-xa
37.  
38.         fcom
39.         fstsw ax
40.         sahf
41.         ja .l2
42.  
43. ;L<=P
44.         fstp st0
45.         fstp st0        ;yb-ya xb-xa  yc-ya xc-xa
46.         fsubp st2,st0   ;xb-xa  yc-yb xc-xa
47.         fsubp st2,st0   ;yc-yb xc-xb
48.         jmp .l3
49.  
50. ;0<P<L
51.     .l2:                ;L P  yb-ya xb-xa  yc-ya xc-xa
52.         fdivp           ;P/L  yb-ya xb-xa  yc-ya xc-xa
53.         fmul st2,st0    ;P/L  yb-ya P/L*(xb-xa)  yc-ya xc-xa
54.         fmulp           ;P/L*(yb-ya) P/L*(xb-xa)  yc-ya xc-xa
55.         fsubp st2,st0   ;P/L*(xb-xa)  yc-ya-P/L*(yb-ya) xc-xa
56.         fsubp st2,st0   ;yc-ya-P/L*(yb-ya) xc-xa-P/L*(xb-xa)
57.  
58.     .l3:                ;y x
59.         fld st1         ;x y x
60.         fmulp st2,st0   ;y x*x
61.         fmul st0,st0    ;y*y x*x
62.         faddp           ;y*y+x*x
63.         fsqrt
64.         ret

 

вот относительно этого фрагмента

Код (Text):

1. faddp           ;P  yb-ya xb-xa  yc-ya xc-xa

Здесь не указано никаких операндов, разве обработчик поймёт что мы складываем st1 и st0
и ещё
у меня в букваре (Зубков) написано, после сложения FADDP выталкивает st(0) из стека а сам st0 помечает как пустой...

Далее...

Код (Text):

1. fldz            ;0.0 P  yb-ya xb-xa  yc-ya xc-xa
2. fcomp           ;P  yb-ya xb-xa  yc-ya xc-xa
3. fstsw ax
4.  sahf
5.  jb .l1

После того как сравнили резкльтирующие флаги С3 С2 и С0 переводятся во флаги ZF PF и CF
не совсем понятна строчка

Код (Text):

1. jb .l1

 

triam

faddp == faddp st(1),st(0)

jb - это условный переход, который выполняется, если CF == 1.

 

Код (Text):

1.  
2.                
3.         fstsw ax   ; вот здесь

 

Код (Text):

1.  
2.         fldz            
3.         fcomp        
4.         fstsw ax   ; вот здесь
5.         sahf
6.         jb .l1

не пойму, почему компилятор ругается
вещает: Illegal immediate
в букваре написано в качестве приёмника должен выступать AX