Пересечение отрезков

[Kuguar]

Добрый день!
Уважаемые специалисты, помогите найти решение следующей задачи:
Необходимо найти координаты пересечения двух отрезков на поверхности земли.
Есть географические координаты начала и конца отрезков (широта, долгота).
Нужна математика (формула расчета).
Графический вариант (построения на карте) не подходит в ввиду большого количества отрезков, которые надо обсчитать.

[bingeomap]

А топологические функции в каком нибудь ГИС программе не подходит?

[Kuguar]

А топологические функции в каком нибудь ГИС программе не подходит?

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

[Игорь Белов]

Предлагаю такой алгоритм.

Для первичного отрезка находим длину D₀ и параметры C и S:

Код: Выделить всё

∆X₀ = X₀₂ − X₀₁
∆Y₀ = Y₀₂ − Y₀₁
D₀² = ∆X₀² + ∆Y₀²
C = ∆X₀ / D₀
S = ∆Y₀ / D₀

Последовательно для всех вторичных отрезков находим параметр Ξ:

Код: Выделить всё

ξ₁ = C (X₁ − X₀₁) + S (Y₁ − Y₀₁)
η₁ = −S (X₁ − X₀₁) + C (Y₁ − Y₀₁)
ξ₂ = C (X₂ − X₀₁) + S (Y₂ − Y₀₁)
η₂ = −S (X₂ − X₀₁) + C (Y₂ − Y₀₁)
Ξ = ξ₁ − η₁ (ξ₂ − ξ₁) / (η₂ − η₁)

Это есть расстояние (со знаком) от начала первичного отрезка (X₀₁, Y₀₁) до точки пересечения со вторичным отрезком. Если оно меньше ноля или больше длины D₀, то пересечение прямых лежит за пределами первичного отрезка.

Если Вам действительно зачем-то нужны координаты точки пересечения двух прямых, то вот они:

Код: Выделить всё

X = X₀₁ + C Ξ
Y = Y₀₁ + S Ξ

Осторожно! Это решение для плоскости.

[Kuguar]

Спасибо большое за интересный вариант расчета!
я думал попробовать воспользоваться теоремой прямых (тоже для плоскости),
НО что брать за X и Y , если координаты заданы в формате ГГ ММ СС.
Надо сначала перевести в прямоугольные или как ?
Не подскажите?

[Игорь Белов]

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

[Kuguar]

Спасибо огромное!
попробую - отпишусь.

[Kuguar]

to: ErnieBoyd
Спасибо еще раз за методику! Считает в пределах нужной точности.
Похоже, это уже мои уши торчат (много лишних срабатываний- находит точки пересечения там, где их не должно быть)

[Franklin1967]

Kuguar, в принципе на такие вопросы я не отвечаю на этом форуме, но для Вас сделаю исключение.
Вы пишите что один из модулей моей программы значит Вам нужен код а не методика.
Попробую дать и то и другое.
Методика
1. грубая на малых расстояниях
- Во заданным координатам двух точек находите линию на плоскости переведя из геодезической системы координат в какую то прямоугольную. Для простоты советую систему прямоугольную Меркатора, известную как "прямоугольная секта им.Ген.Штаба" или "UTM" если преобразования производится с эллипсоида WGS84
- полученные пару прямоугольных координат "пропускаете" через другую функцию , функцию пересечения отрезков
2. точная, на больших расстояниях.
- разбиваете начало и окончание отрезков на более мелкие части в массивы отрезков.
- Каждый кусочек преобразовываете в прямоугольную систему координат (согласно п1 предыдущей методике).
- проверяете пересечение и ищите точку.
Второй метод ищет пересечение с учётом кривизны Вашего отрезка в прямоугольной системе координат.
Теперь код
Я пишу в Delphi поэтому коды представляю на Паскале. Думаю разберётесь по ходу.
1. Функция преобразования координат в геодезической системе в координаты в проекции Меркатора

Код: Выделить всё

// =============================================================================
// ГАУСА-КРЮГЕРА
// на вход lp - координата в радианах, FKrassovsky -параметры эллипсоида (ниже по тексту)
// результат - координата в метрах
function GaussForward(lp : T3dPoint; FKrassovsky : TEllipsoidProperty): T3dPoint;
var N, X       : extended;
    l_ro,l_ro2 : extended;
    T2, n2     : extended;
begin
   X:=FKrassovsky.a*(1-FKrassovsky.es)*(1.0050517739*lp.phi-0.5*0.0050623776*sin(2*lp.phi)+
      0.25*0.0000106245*sin(4*lp.phi)+0.0000000208*sin(6*lp.phi)/6);
   N:=FKrassovsky.a/sqrt(1-FKrassovsky.es*sqr(sin(lp.phi)));
   T2:=sqr(tan(lp.phi));
   n2:=FKrassovsky.es*sqr(cos(lp.phi))/(1-FKrassovsky.es);
   L_ro:=0.017453276125372700167260562868155*(180*lp.lam/pi)*cos(lp.phi);
   l_ro2:=l_ro*l_ro;
   result.X:=l_ro2*(1+l_ro2*((5-T2+9*n2+4*n2*n2)+l_ro2*(61-58*t2+t2*t2+270*n2-330*n2*t2)/30)/12);
   result.X:=X+0.5*result.X*N*tan(lp.phi);
   result.Y:=1+l_ro2*(1-T2-n2+l_ro2*((5-18*t2+t2*t2+14*n2-58*n2*t2+13*n2*n2-64*n2*n2*t2)+
     l_ro2*(61-479*t2+179*t2*t2-t2*t2*t2)/42)/20)/6;
   result.Y:=l_ro*N*result.Y;
end;

где:
 T3DPoint =   packed record
  case byte of
    0: (X,Y,Z   : double);
    1: (phi,lam,height :double);
    2: (B,L,H   : double);
    3: (REC     : array[0..2] of double);
 end;
T3dArray = array of T3DPoint ;


 TEllipsoidProperty = packed record
  eName   : string;
  a       : double;  // major axis or radius if es=0
  b       : double;  // minor axis or radius if es=0
  e       : double;  // eccentricity
  es      : double;  // e ^ 2
  one_es  : double;  // 1 - e^2
  rone_es : double;  // 1/(1 - e^2)
  ra      : double;  // 1/A
  Rf      : double;  // references
  Code    : string;
 end;

в котором параметры для

Код: Выделить всё

FKrassovsky =  (ename:'(CK-42) Krassovsky 1942';a:6378245;b:6356863.01877305;e: 0.081813334; es: 0.00669342162296504;one_es:   0.993306578377035;rone_es: 1.00673852541468;ra: 0.156782939507655E-6;rf: 0.00335232986925913;code:'krass')

2. Пересечение двух отрезков в метрах

Код: Выделить всё

{
 NO_CROSS     = $00; // нет пересечения
 CR_LINE_SS   = $01; // сегмент-сегмент
 CR_LINE_SV   = $02; // сегмент-вершина
 CR_LINE_VV   = $03; // одна вершина совпадает    1----1(2)----2
 CR_LINE_SV2  = $04; // обе вершине на сегментах  2----1---2----1
 CR_LINE_VV2  = $05; // две линии имеют доинаковые координаты  1(2)----2(1)
}
 // ПЕРЕСЕЧЕНИЕ
 // отрезок с отрезком
// Res  - массив пересечений. Может быть пуст, с одим или двумя результатами
function CrossLineToLine(p1,p2, l1,l2 : T3dPoint; var Res : T3dArray; const Decimal : byte = 4): byte;
var V1,V2 : T3dpoint;
    i     : integer;
    minP, maxP  : T3dPoint;
    minL, maxL  : T3dPoint;
    delta : double;
begin
 delta := power(10, -(Decimal+1));
 result:=0;  Finalize(res);
 // одна из линий нулевая
 if(Distance(p1,p2)=0) or (Distance(l1,l2)=0) then exit;
 // проверяем вершины с вершинами
 if (Distance(p1,l1)<delta) or (Distance(p1,l2)<delta) then IncArray(Res,[p1]);
 if (Distance(p2,l1)<delta) or (Distance(p2,l2)<delta) then IncArray(Res,[p2]);
 // если вых.массив заполнен хоть одним элементом то ВЫХОД
 case LenGth(Res) of
  1: result:=CR_LINE_VV;
  2: result:=CR_LINE_VV2;
 end;
 if result>0 then exit;

 minP:=Set3dPoint(Min(p1.x,p2.x),Min(p1.Y,p2.Y));
 maxP:=Set3dPoint(Max(p1.x,p2.x),Max(p1.Y,p2.Y));
 minL:=Set3dPoint(Min(l1.x,l2.x),Min(l1.Y,l2.Y));
 maxL:=Set3dPoint(Max(l1.x,l2.x),Max(l1.Y,l2.Y));

 V1:=Set3dPoint(RoundTo(p2.x-p1.X,-Decimal),RoundTo(p2.Y-p1.y,-Decimal),0);
 V2:=Set3dPoint(RoundTo(l2.x-l1.X,-Decimal),RoundTo(l2.Y-l1.y,-Decimal),0);
 // первая вертикальная dp=0
 if (V1.X=0) and (V2.X<>0) then
 begin
  V2.Z:=V2.Y/V2.X; V2.X:=l1.y-l1.x*V2.Z;
  V1:=Set3dPoint(RoundTo(p1.x,-Decimal),RoundTo(V2.Z*p1.x+V2.X,-Decimal),0);
  if ((V1.Y>=minP.y) and (V1.Y<=MaxP.y)) and
     ((V1.X>=minL.x) and (V1.X<=maxL.x)) then IncArray(Res,[v1]);
  result:=CR_LINE_SS*byte(LenGth(Res)>0)+
  byte(((V1.Y=minP.y) and (V1.Y=MaxP.y)) and ((V1.X=minL.x) and (V1.X=maxL.x)));

 end else
 // вторая вертикальная   dl=0
 if (V1.X<>0) and (V2.X=0) then
 begin
  V1.Z:=V1.Y/V1.X; V1.X:=p1.y-p1.x*V1.Z;
  V2:=Set3dPoint(RoundTo(l1.x,-Decimal),RoundTo(V1.Z*l1.x+V1.X,-Decimal),0);
  if ((V2.Y>=minL.y) and (V2.Y<=maxL.y)) and
     ((V2.X>=minP.x) and (V2.X<=MaxP.x)) then IncArray(Res,[v2]);
  result:=CR_LINE_SS*byte(LenGth(Res)>0)+
  byte(((V2.Y=minL.y) and (V2.Y=maxL.y)) and((V2.X=minP.x) and (V2.X<=MaxP.x)));
 end else
 // обе вертикальные и совпадающие
 if ((V1.X=0) and (V2.X=0)) and (p1.X=l1.X) then
 begin

  if (maxL.y>=p2.Y) and (p2.y>=MinL.y) then IncArray(Res,[p2]) else
  if (maxL.y>=p1.Y) and (p1.y>=MinL.y) then IncArray(Res,[p1]);

  if (MaxP.y>=l2.Y) and (l2.y>=MinP.y) then IncArray(Res,[l2]) else
  if (MaxP.y>=p1.Y) and (l1.y>=MinP.y) then IncArray(Res,[L1]);
  result:=CR_LINE_SV*byte(LenGth(Res)=1)+CR_LINE_SV2*byte(LenGth(Res)=2);

 end else
 // обе нормальные
 begin
  V1.Z := MaxDouble;  V2.Z := MaxDouble;
  if v1.x<>0 then
  begin
   V1.Z:=V1.Y/V1.X;
   V1.X:=p1.y-p1.x*V1.Z;
  end else
  V1.X := 0;
  if v2.x<>0 then
  begin
   V2.Z:=V2.Y/V2.X;
   V2.X:=l1.y-l1.x*V2.Z;
  end else
  V2.X := 0;
  // паралельные и совпадающие
  if (RoundTo(V1.Z-V2.Z,-Decimal)+RoundTo(V1.X-V2.X,-Decimal)=0) then
  begin
   if PointInRect(l1,l2,p1, Decimal) then IncArray(Res,[p1]);
   if PointInRect(l1,l2,p2, Decimal) then IncArray(Res,[p2]);
   if PointInRect(p1,p2,l1, Decimal) then IncArray(Res,[l1]);
   if PointInRect(p1,p2,l2, Decimal) then IncArray(Res,[l2]);
   result:=CR_LINE_SV*byte(LenGth(Res)=1)+CR_LINE_SV2*byte(LenGth(Res)=2);
  end else
  // совпадение по ребрам
  begin
   IncArray(Res,[v2]);
   if V1.Z-V2.Z<>0 then
   begin
    Res[0].X:=(V2.X-V1.X)/(V1.Z-V2.Z);
    Res[0].Y:=V1.Z*Res[0].X+V1.X;
    if not PointInRect(l1,l2,res[0], Decimal) or not PointInRect(p1,p2,res[0], Decimal) then
    Finalize(Res);
   end else Finalize(Res);
   result:=CR_LINE_SS*byte(LenGth(Res)>0);
  end;
 end;
end;

Дерзайте =)

+ функция PointInRect пару раз встречается

Код: Выделить всё

Function PointInRect(D0,D1, Pnt : T3dPoint;const Decimal : byte = 4) : boolean;
var FPnt: T3dPoint;
begin
 FPnt:=Set3dPoint(RoundTo(Pnt.X,-Decimal),RoundTo(Pnt.Y,-Decimal),0);
 result:=
  (RoundTo(Max(D0.X,D1.X),-Decimal)>=fPnt.X) and (fPnt.X>=RoundTo(Min(D0.X,D1.X),-Decimal)) and
  (RoundTo(Max(D0.y,D1.y),-Decimal)>=fPnt.y) and (fPnt.y>=RoundTo(Min(D0.Y,D1.Y),-Decimal));
end;

+ дополняю готовым модулем для проверки

[Kuguar]

Уважаемый Franklin1967, ОГРОМНОЕ СПАСИБО за помощь и предоставленный код!