Schnittpunkt zweier Strecken ermitteln

[Namenloser]

Hallo Community,

nachdem ich schon länger nach einer Funktion, die den Schnittpunkt von zwei Strecken errechnet, gesucht habe, bin ich heute bei Google fündig geworden und habe einen ActionScript-Code dafür gefunden. Diesen habe ich dann noch schnell in Delphi überetzt, et voilà: Es funktioniert tatsächlich

Die Funktion IntersectLines ist einmal überladen. Die erste Funktion erwartet Fließkommazahlen als Parameter, die zwiete arbeitet mit Points, wobei p1 und p2 die Endpuntke der ersten, und p2 und p3 die Endpunkte der zweiten Strecke markieren. Die Funktion gibt False zurück, wenn sich die Strecken nicht schneiden, andernfalls steht in output der Schnittpunkt.

Es gibt allerdings eine Division durch 0, wenn die Strecken exakt aufeinander liegen, man sollte diesen Fall also vorher im Programm abfangen.

Delphi-Quellcode:

type
  tvec2d = class
    x,
    y: extended;
    constructor create (ax,ay: extended);
    function plus(v: tvec2d): tvec2d;
    function minus(v: tvec2d): tvec2d;
    procedure multiply(f: extended);
  end;
  constructor tvec2d.create(ax, ay: extended);
  begin
    x := ax;
    y := ay;
  end;
  function tvec2d.minus(v: tvec2d): tvec2d;
  begin
    result := tvec2d.create(x-v.x,y-v.y);
  end;
  procedure tvec2d.multiply(f: extended);
  begin
    x := x*f;
    y := y*f;
  end;
  function tvec2d.plus(v: tvec2d): tvec2d;
  begin
    result := tvec2d.create(x+v.x,y+v.y);
  end;
function intersectlines(x1,y1,x2,y2,x3,y3,x4,y4: extended;
  var output: tpoint): boolean; overload;
var
  v0,v,v1: tvec2d; // strecke 1  [p1,p2,richtung]
  u0,u,u1: tvec2d; // strecke 2  [p1,p2,richtung]
  fdiv,
  t,s: extended;
  p: tvec2d;
begin
  result := false;
  // Strecke 1
  v0 := tvec2d.create(x1,y1);
  v := tvec2d.create(x2,y2);
  v1 := v.minus(v0);
  // Strecke 2
  u0 := tvec2d.create(x3,y3);
  u := tvec2d.create(x4,y4);
  u1 := u.minus(u0);

  fdiv := v1.x * u1.y - v1.y * u1.x;
  t := -(v0.x * u1.y - v0.y * u1.x - u0.x * u1.y + u0.y * u1.x) / fdiv;
  s := -(v0.x * v1.y - v0.y * v1.x + v1.x * u0.y - v1.y * u0.x) / fdiv;

  v1.multiply(t);
  // Punkt wo sie sich schneiden:
  p := v0.plus(v1);
  if ((t>0) AND (s>0)) AND ((t<1) AND (s<1)) then
  begin
    output := point(round(p.x),round(p.y));
    result := true;
  end;

  freeandnil(v0);
  freeandnil(v);
  freeandnil(v1);
  freeandnil(u0);
  freeandnil(u);
  freeandnil(u1);
  freeandnil(p);
end;

function intersectlines(p1,p2,p3,p4: TPoint;
  var output: tpoint): boolean; overload;
begin
  result := intersectlines(p1.x,p1.y,p2.X,p2.y,p3.X,p3.Y,p4.x,p4.y,output);
end;

Ich hoffe, das finden auch noch andere nützlich

[shmia]

Zitat von NamenLozer:

Es gibt allerdings eine Division durch 0, wenn die Strecken exakt aufeinander liegen, man sollte diesen Fall also vorher im Programm abfangen.

Das sollte von intersectlines gemanaged werden:

Delphi-Quellcode:

  fdiv := v1.x * u1.y - v1.y * u1.x; 
  if Abs(fdix) < 1e-20 then
     result := True;  
  else
  begin
    t := -(v0.x * u1.y - v0.y * u1.x - u0.x * u1.y + u0.y * u1.x) / fdiv; 
    s := -(v0.x * v1.y - v0.y * v1.x + v1.x * u0.y - v1.y * u0.x) / fdiv; 

    v1.multiply(t); 
    // Punkt wo sie sich schneiden: 
    p := v0.plus(v1); 
    if ((t>0) AND (s>0)) AND ((t<1) AND (s<1)) then 
    begin 
      output := point(round(p.x),round(p.y)); 
      result := true; 
    end; 
  end;

[Namenloser]

Hallo,

an soetwas hatte ich auch gedacht, allerdings hat man dann das Problem das in output irgendwelche nicht initialisierten Werte stehen. (es sei denn natürlich, man initialisiert sie davor ). Wenn ein Programm den ausgegebenen Wert einfach übernimtm wird das zu Fehlern führen, also müsste man so oder so diesen Fall im Programm gesondert behandeln.

[sirius]

Sowas kann schon nützlich sein (gibts das etwa noch nicht ) wenn man wieder mal 1 Stunde über diesen blöden Schnittpunkt zweier Geraden sitzt. Tausendmal in Mathe gemacht (mit Funktionen und mit Vektoren) und dann will man genau das "nur mal schnell" hinprogrammieren ....

Aber:
1. Müssen für die Vektoren unbedingt gleich Objekte erzeugt werden?
2. Warum unbedingt auf Strecken begrenzt (ich denke Geradenschnittpunkte sind häufiger von Interesse)?

[Namenloser]

1. Weil ich mich möglichst exakt an den Code halten wollte um keine Fehler zu produzieren
2. Weil es dafür bereits tausend Lösungen gibt (glaub sogar in der CodeLib), mich aber speziell Streckenschnittpunkte interessiert haben.

[sirius]

zu 2
Ia der einzige Unterschied (zwischen Gerade und Strecke) ist eben deine If-Abfrage am Ende. Das hätte man dann auch unter jeden anderen code (wenn es ihn denn in der CodeLib gibt) setzen können.

[LoCrux]

Da ich gerade mal wieder hier am stöbern bin...

..diese (folgende) Lösung hatte ich mal vor Ewigkeiten zusammengebastelt. Fragt mich bitte nichtmehr zwecks Dokumentation.
Die Lösung funktioniert trozdem....

Delphi-Quellcode:

interface {Part}

//==[ CONSTANTS ]========//

const
  YES         = TRUE;
  NO          = FALSE;
  cPi         = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
  c2Pi        = 2*cPi;
  cPi90       = cPi/2;
  cPi270      = cPi+cPi90;
  cDegToRad   = (cPi/180);
  cRadToDeg   = (180/cPi);

//==[ TYPE-DEFS ]========//

type

  TSimplePoint = packed record
                   x,y : Double;
                 end;

  TLineLinear   = packed record
                    p1,p2 : TSimplePoint;  // The Building Points;
                    A,B,C : Double;        // 0 = A*x + B*y + C
                  end;

  TLinesBehave  = (lbSame, lbIntersect, lbOrthogonal, lbParallel,
                   lbline1ParaX, lbline1ParaY, lbline2ParaX, lbline2ParaY);

  TLinesBehaves = set of TLinesBehave;

implementation {Part}

//==[ LINE OPERATIONS ]==//

function IsSame(val1,val2:Extended):Boolean; 
begin
  result := IsZero(val1-val2);
end;

function GetMin(values:array of Double;IncZero:Boolean):Double;
var
  s,cnt   : INTEGER;
  AllZero : BOOLEAN;
begin
  s := LOW(values);
  AllZero := NO;
  if IncZero
  then begin
    result := values[s];
    INC(s);
  end
  else begin
    repeat
      AllZero := (s>High(Values));
      if not(AllZero)
      then begin
        result  := values[s];
        INC(s);
      end;
    until Not(IsZero(result)) or AllZero;
  end;
  if not(AllZero)
  then begin
    for cnt := LOW(values)+s to HIGH(values)
    do if (result>values[cnt]) and (IncZero and IsZero(values[cnt]))
       then result:=values[cnt];
  end
  else result := 0;
end;

function Det2x2(p1,p2:TSimplePoint):Double;
begin
  result := (p1.x*p2.y)-(p1.y*p2.x);
end;

function BuildLineLinFromPoints(p1,p2:TSimplePoint):TLineLinear;
var
  min   : Double;
begin
  result.p1 := p1;
  result.p2 := p2;
  result.A  := p1.y-p2.y;
  result.B  := p2.x-p1.x;
  result.C  := DET2x2(p1,p2);
  min := GetMin([ABS(result.A),ABS(result.B),ABS(result.C)],no);
  if not(IsZero(min))
  then begin
    result.A := result.A/min;
    result.B := result.B/min;
    result.C := result.C/min;
  end;
end;

function  CheckLineLinBehaviour(ll1,ll2:TLineLinear;var Alpha:Double):TLinesBehaves;
var
  NotReady  : Boolean;
  m1,m2     : Double;

  function DoAlphaLin:Double;
  var
    Divident,
    Divisor   : Double;
  begin
    Divisor  := (ll1.A*ll2.A)-(ll1.B*ll2.B);
    Divident := (ll1.A*ll2.B)-(ll1.B*ll2.A);
    if not(IsZero(Divisor))
    then result := arctan(Divident/Divisor)
    else result := 0;
  end;

  function DoAlphaNor:Double;
  var
    Divident,
    Divisor   : Double;
  begin
    m1       := ll1.A/(-ll1.B);
    m2       := ll2.A/(-ll2.B);
    Divisor  := m2-m1;
    Divident := 1+(m1*m2);
    if not(IsZero(Divisor))
    then result := arctan(Divident/Divisor)
    else result := 0;
  end;

begin
  result := [];
  if IsZero(ll1.A) then result := result+[lbLine1ParaX];
  if IsZero(ll2.A) then result := result+[lbLine2ParaX];
  if IsZero(ll1.B) then result := result+[lbLine1ParaY];
  if IsZero(ll2.B) then result := result+[lbLine2ParaY];

  if (lbLine1ParaX in result) or
     (lbLine1ParaY in result) or
     (lbLine2ParaX in result) or
     (lbLine2ParaY in result)
  then begin
    NotReady := YES;

    if (lbLine1ParaX in result) and NotReady
    then begin
      if (lbLine2ParaX in result)
      then begin
        Alpha    := 0;
        result   := result + [lbParallel];
        NotReady := NO;
      end
      else if (lbLine2ParaY in result)
      then begin
        Alpha    := cPi90;
        result   := result + [lbOrthogonal];
        NotReady := NO;
      end 
      else begin
        Alpha    := DoAlphaLin;
        result   := result + [lbIntersect];
        NotReady := NO;
      end;
    end;

    if (lbLine1ParaY in result) and NotReady
    then begin
      if (lbLine2ParaX in result)
      then begin
        Alpha    := cPI90;
        result   := result + [lbOrthogonal];
        NotReady := NO;
      end
      else if (lbLine2ParaY in result)
      then begin
        Alpha    := 0;
        result   := result + [lbParallel];
        NotReady := NO;
      end 
      else begin
        Alpha    := DoAlphaLin;
        result   := result + [lbIntersect];
        NotReady := NO;
      end;
    end;

    if (lbLine2ParaX in result) and NotReady
    then begin
      if (lbLine1ParaX in result)
      then begin
        Alpha    := 0;
        result   := result + [lbParallel];
        NotReady := NO;
      end
      else if (lbLine1ParaY in result)
      then begin
        Alpha    := cPi90;
        result   := result + [lbOrthogonal];
        NotReady := NO;
      end
      else begin
        Alpha    := DoAlphaLin;
        result   := result + [lbIntersect];
        NotReady := NO;
      end;
    end;

    if (lbLine2ParaY in result) and NotReady
    then begin
      if (lbLine1ParaX in result)
      then begin
        Alpha    := cPi90;
        result   := result + [lbOrthogonal];
        NotReady := NO;
      end
      else if (lbLine1ParaY in result)
      then begin
        Alpha    := 0;
        result   := result + [lbParallel];
        NotReady := NO;
      end 
      else begin
        Alpha    := DoAlphaLin;
        result   := result + [lbIntersect];
        NotReady := NO;
      end;
    end;

  end
  else begin
    Alpha    := DoAlphaNor;
    if IsSame(m1,m2)
    then begin
      alpha  := 0;
      result := result + [lbParallel];
    end
    else if IsSame(m1,-(1/m2))
    then begin
      alpha  := cPI90;
      result := result + [lbOrthogonal];
    end
    else begin
      result   := result + [lbIntersect];
    end;
  end;
end;

function  FindLineIntersection(ll1,ll2:TLineLinear;var intersect:TSimplePoint;var TanAlpha:Double):TLinesBehaves;
var
  D,DX,DY : Double;
  hp1,hp2 : TSimplePoint;
begin
  result := [];
  hp1.x := ll1.A; hp1.y := ll2.A;
  hp2.x := ll1.B; hp2.y := ll2.B;
  D := Det2x2(hp1,hp2);
  if (D<>0)
  then begin
    hp1.x := ll1.C; hp1.y := ll2.C;
    hp2.x := ll1.B; hp2.y := ll2.B;
    DX := Det2x2(hp1,hp2);

    hp1.x := ll1.A; hp1.y := ll2.A;
    hp2.x := ll1.C; hp2.y := ll2.C;
    DY := Det2x2(hp1,hp2);

    intersect.x := -(1/D)*DX;
    intersect.y := -(1/D)*DY;
  end;
  result := CheckLineLinBehaviour(ll1,ll2,TanAlpha);
end;