Schnittpunkt Linie vs. Kreis

[Namenloser]

Ich war vorhin auf der Suche nach einem Code, mit dem sich die Schnittpunkte eines Kreises und einer Linie ausrechnen lassen. Da ich für Delphi nichst gefunden habe, habe ich folgende Funktion geschrieben:

Delphi-Quellcode:

function LineHitsCircle(LineX1,LineY1,LineX2,LineY2,
  CircleX,CircleY,CircleRadius: double; out OutPos: Double): boolean;
var
  x1, x2, x: double;
  d: double;
  d_helper: double;
  m,n: double;

  tmp: double;
begin
  LineX1 := LineX1 - CircleX;
  LineY1 := LineY1 - CircleY;
  LineX2 := LineX2 - CircleX;
  LineY2 := LineY2 - CircleY;

  if (LineX2-LineX1=0) then
  begin
    if (LineY2-LineY1=0) then
    begin
      OutPos := 0.0;
      result := sqr(LineX1)+sqr(LineY1) <= sqr(CircleRadius);
      exit;
    end
    else
    begin
      tmp := LineY1;
      LineY1 := LineX1;
      LineX1 := tmp;
      
      tmp := LineY2;
      LineY2 := LineX2;
      LineX2 := tmp;
    end;
  end;

  m := (LineY2-LineY1)/(LineX2-LineX1);
  n := LineY1-m*LineX1;
  d_helper := 4*sqr(m)*sqr(n) - 4*(1+sqr(m))*(sqr(n)+2*n-sqr(CircleRadius));

  if d_helper < 0 then
  begin
    OutPos := -1;
    result := False;
  end
  else
  begin
    d := sqrt(d_helper);
    x1 := (-2*m*n - d)/(2+2*sqr(m));
    x2 := (-2*m*n + d)/(2+2*sqr(m));

    if (x1>=LineX1) and ((x1 - LineX1)<=(x2-LineX1))  then
      x := x1
    else
      x := x2;

    OutPos :=  abs((x-LineX1)/(LineX2-LineX1));

    result := (OutPos >= 0.0) and (OutPos <= 1.0);
  end;
end;

Sicher noch optimierbar, und gut möglich, dass es mit Vektoren irgendwie eleganter geht (das war bisher immer so :mrgreen:). Die Funktion gibt keinen Punkt zurück, sondern die Stelle auf der Strecke (Gleitkommazahl von 0 (Startpunkt) bis 1 (Endpunkt)), weil das für meine Zwecke praktischer ist. Man könnte aber auch sehr einfach den Punkt ausrechnen (über die lineare Funktion, m und n sind ja gegeben). Wenn es zwei Schnittpunkte gibt, wählt die Funktion den, der am nächsten beim Startpunkt liegt.
Wenn man die If-Abfragen am Ende weglässt, hat man eine Funktion, die die Schnittpunkte einer geraden und eines Kreises berechnet.

Ich hoffe, dass eine solche Funktion nicht schon gibt in der CodeLibrary existiert. Über die Suche habe ich wie gesagt nichts gefunden.

[edit]Bugfix bei senkrechter Linie[/edit]

[mr_emre_d]

Ein bisschen Theorie:

Bekannte Werte:
Mittelpunkt ( -> Mx, My )
Radius ( -> r )
Linie ( f )

-----------------

k: (x-Mx)²+(y-My)² = r²

f(x) = y = kx + d

-----------------

Gleichsetzen und ausrechnen

EDIT: Ich rechne mal mit Notepad ein Beispiel aus...

MfG

[Namenloser]

Zitat:

Zitat von mr_emre_d

Ein bisschen Theorie:

Bekannte Werte:
Mittelpunkt ( -> Mx, My )
Radius ( -> r )
Linie ( f )

-----------------

k: (x-Mx)²+(y-My)² = r²

f(x) = y = kx + d

-----------------

Gleichsetzen und ausrechnen

EDIT: Ich rechne mal mit Notepad ein Beispiel aus...

MfG

Ich verstehe nicht ganz, was du mir sagen willst :gruebel: Oder sitze ich schon wieder zu lange vor der Kiste?

[mr_emre_d]

xD

Du willst ja die Schnittpunkte ausrechnen oder ? Ich zeig dir nur, wie es mathematisch möglich ist.
Wie weit du das in delphi implementierst ist dir überlassen ;)

MfG

Weiter gehts:

Code:

k: (x-Mx)²+(y-My)² = r²

f(x) = y = kx + d 

Konkretes Beispiel:

Mittelpunkt = ( 5 | 6 )
Radius = 5
f(x) = 2x + 3
--:

(x-5)²+(y-6)² = 25
x²-10x+25 + y²-12y+36 = 25
x²-10x + y²-12y = -36
--------------------------

x²-10x + y²-12y = -36
   (2x+3)² -> 4x²+12x+9
      -12*(2x+3) -> -24x-36

x²-10x + 4x²+12x+9 -24x-36 = -36

x²-10x -12x + 4x² + 9 =
 0


5x² - 22x + 9 = 0   |   / 5

x² - 4,4x + 1,8 = 0

( Quadratisch Gleichung auflösen: -p/2 +- SQRT( (-p/2)² - q )   |   p = -4.4   q = 1.8 )
   
   x1,2 = 2,2 +- SQRT( 4,84 - 1.8 )
   x1 = 2.2 + 1,74 -> 3.94
   x2 = 2.2 - 1,74 -> 0.46

Punkt1( 3.94 | y )
Punkt1( 0.46 | y )

y = 2x + 3 ->
Punkt1 y = 2*3.94 + 3 = 10.88
Punkt2 y = 2*0.46 + 3 = 3.92

Schnittpunkte
   Punkt1( 3.94 | 10.88 )
   Punkt2( 0.46 | 3.92 )

Graph (Kreis + F) im Anhang :)

EDIT:
Wohoo .. ich dacht mir du brauchst Hilfe - bin auf den Thread hier durch die Startseite gestoßen und habe nur deinen ersten Satz gelesen -> wusste damit nicht, dass es ein Vorschlag für die Code-Lib ist :S )

EDIT2:

Delphi-Quellcode:

type
  TPointF = record
    X, Y: Single;
  end;

  T2Points = Array[0..1] of TPointF;

...

function Intersection_Circle_Line( cRadius, cMX, cMY: Single; lKX: Single;
  const lD: Single = 0.0 ): T2Points;
{
  CIRCLE:
    cRadius   = Radius of the circle
    cMX / cMY = Center Coordinates
  LINE:
    lKX       = INCREASE
    lD        = Intersection with y-axis
}
var
  xX, X,
  C, T: Single;
begin
  xX := 1 + lKX * lKX;
  X  := -cMX*2 + lKX*lD*2 + (-cMY*2)*lKX;
  C  := cMX*cMX + lD*lD + (-cMY*2)*lD + cMY*cMY - cRadius*cRadius;
  if xX > 1 then
  begin
    X := X / xX;
    C := C / xX;
  end;          
  if (X*X/4) < C then
    Exit;
  t := SQRT( X*X/4 - C );;
  Result[0].X := -X/2 + t;
  Result[0].Y := lKX * Result[0].X + lD;
  Result[1].X := -X/2 - t;
  Result[1].Y := lKX * Result[1].X + lD;
end;

Mir war langweilig :P

[Medium]

Und da es schon angesprochen wurde, hier das ganze noch vektoriell in Pseudocode (um den ganzen Krams zu sparen den man betreiben muss um Vektoroperationen durch Funktionen zu ersetzen). Im Grunde läuft es auf eine pq- bzw. Mitternachtsformel hinaus.

Delphi-Quellcode:

// M : Mittelpunkt des Kreises
// r : Kreisradius
// A : Stützpunkt der zu schneidenden Strecke/Gerade
// V : Richtungsvektor der Strecke/Gerade (Wenn diese aus 2 Punkten A und B gebildet ist, ist dies B-A)
// U : Hilfsvektor
// p, q, d, x1, x2 : Hilfsvariablen (Skalare, keine Vektoren)

  U := A-M;
  p := U*V;
  q := U^2 - r^2;
  d := p^2 - q;

  if d < 0 then 
    Kein_Schnittpunkt
  else
  begin
    d  := sqrt(d / |V|);
    x1 := -p + d;
    x2 := -p - d;
  end;

  if (x2>x1) and (x1>0) then
    Schnittpunkt := A + x1*V
  else
  if (x1>x2) and (x2>0) then
    Schnittpunkt := A + x2*V
  else
    Kein_Schnittpunkt_Oder_Schnitt_Hinter_A;

Der je kleinere Wert von x1 und x2 gibt den ersten Schnittpunkt von A aus in Richtung V gesehen, die Bedingung >0 stellt sicher dass kein Punkt "hinter" der Blickrichtung "V gestützt auf A" gefunden wird. Liegt dieser in ]0|1[, ist der Schnittpunkt zwischen den Punkten A und B bzw. A und A+V.

Das ganze klappt nun auch für beliebige Dimensionen, man muss lediglich entsprechende Vektortypen und die zugehörigen Operationen nehmen. Ich hab das z.B. in der Form für 3D in einem kleinen Raytracer im Kugel-Primitive eingesetzt.

[Dipl Phys Ernst Winter]

Es gibt zwei oder keinen Schnittpunkt.
Tangiert die Linie den Kreis, so zählt man den Schnittpunkt doppelt.

[Namenloser]

Es geht hier um eine Strecke, keine Gerade. Der Code gibt immer den ersten Schnittpunkt zurück. Das ganze lässt sich aber sehr einfach für Geraden und mehrere Schnittpunkte anpassen.

[Dipl Phys Ernst Winter]

Da fand sich bei mir in AnwMath ein diesbezügliches Programm.

Delphi-Quellcode:

type
  TPkt = record x, y: extended end;

// Kreis: M0 Mittelpunkt, r: Radius
// Gerade: P0 Punkt auf der Geraden, m: Steigung
// S1, S2: Schnittpunkte
// Der Funktioswert gibt mit 0..2 die Anzahl der Schnittpunkt zurück
function TForm1.Schnittpunkte(M0, P0: TPkt; r, m: extended;
                                                   var S1, S2: TPkt): integer;
var
  n, a, d, p, q: extended;
begin
  n:= P0.y - m*P0.x - M0.y;
  a:= 1 + m*m;
  p:= 2*(n*m - M0.x);
  q:= M0.x*M0.x + n*n -r*r;
  p:= p/a; q:= q/a; d:= p*p/4 - q;
  if d>eps then begin                  // Diskriminante > 0, zwei Schnittpunkte
    Result:= 2; d:= Sqrt(d); S1.x:= -p/2 - d; S2.x:= -p/2 + d;
    S1.y:= P0.y + m*(S1.x-P0.x); S2.y:= P0.y + m*(S2.x-P0.x) end
  else if d<-eps then Result:= 0       // Diskriminante < 0, kein Schnittpunk
  else begin                           // Diskriminante = 0, Berührungspunkt
    Result:= 1; S1.x:= -p/2; S1.y:= P0.y + m*(S1.x-P0.x);
    S2.x:= S1.x; S2.y:= S1.y end;
end;

[Uwe Raabe]

Zitat:

// Gerade: P0 Punkt auf der Geraden, m: Steigung

Und wie gibt man damit eine Senkrechte an?

[alzaimar]

Delphi-Quellcode:

m := 1E99

sollte reichen

[Namenloser]

Zitat:

Zitat von Uwe Raabe

Zitat:

// Gerade: P0 Punkt auf der Geraden, m: Steigung

Und wie gibt man damit eine Senkrechte an?

Indem man trickst und das Koordinatensystem um 90° dreht (so habe ich es gemacht, s.o.).

[Draos]

Oder man hat mal was in Mathe von Analytik gehört und löst das Problem über Vektoren.

k: [vx-vm]²=r²

g: vx=vn+s*va

mit
vm.. Ortsvektor des Kreismittelpunktes
vn.. Stützvektor der Geraden
va.. Richtungsvektor der Geraden


1. Gleichsetzten der beiden Gleichungen
2. Umformen nach s (ensteht quadratische Gleichung, maximal 2 Werte oder auch keine Lösung falls Wurzelwert<0)
3. Einsetzen von s in g bringt die Ortsvektoren der Punkte

[Nikolas]

Zitat:

Oder man hat mal was in Mathe von Analytik gehört und löst das Problem über Vektoren.

Oder man hat wirklich etwas von Analysis mitbekommen und weiss, das Vektoren in den Bereich der Algebra fallen.

Zitat:

[vx-vm]²=r²

meinst du vielleicht |vx-vm|=r^2 ?

[3_of_8]

Jep, Analytik ist Chemie. :mrgreen:

(Der Ansatz über Vektoren führt aber im Endeffekt auch zum gleichen Ergebnis)

(Und wenn mans ganz genau nimmt, sind reelle Zahlen auch Vektoren, von daher ist der ursprüngliche Ansatz auch über Vektoren :lol: )

[Draos]

1. Analytik war bei uns immer kurzform von Analytische Geometrie.

2. Entweder die Form
k: [vx-vm]²=r² oder die
k: |vx-vm|=r. Deine wäre falsch

[Dipl Phys Ernst Winter]

Zitat:

Und wie gibt man damit eine Senkrechte an?

Vielen Dank für diesen Einwand.

m:= 1e99 würde zwar das Problem lösen, ist aber nicht eingebbar.

Gerade mit Punkt und Steigung war zwar das Einfachste, ist aber ein Fehlgriff. Ich habe das Programm abgeändert auf Gerade durch zwei Punkte.

[3_of_8]

Naja, bei einer senkrechten Gerade/Strecke ist der Fall doch trivial.
Es sei eine Strecke/Gerade [P1P2]/P1P2 gegeben durch die Punkte P1(x1|y1), P2(x2,y2) und ein Kreis k(M,r) mit dem Mittelpunkt M(xm,ym). Es gelte x1=x2.

Delphi-Quellcode:

dx := x1-xm;
disc := r*r-dx*dx;
if disc > 0 then
begin
  setlength(result, 2);
  result[0].X := x1;
  result[0].Y := ym + sqrt(disc);
  result[1].X := x1;
  result[2].Y := ym - sqrt(disc);
end
else if disc = 0 then
begin
  setlength(result, 1);
  result[0].X := x1;
  result[0].Y := ym;
end;

Die Schnittpunkte werden in result zurückgegeben.

(Bei einer Strecke muss man noch überprüfen, ob die gefundenen Schnittpunkte auch wirklich zwischen P1 und P2 liegen, da reicht die Überprüfung ((resultY <= y1) and (resultY >=y2)) or ((resultY >= <1) and (resultY <= y2)) )


Ich habe bisher alle Probleme dieser Art (z.B. Billardsimulation, Dreiecksaufgabe beim BWINF) so gelöst (Betrachtung als Funktion für x1!=x2, ansonsten Sonderfallbehandlung als senkrechte Gerade) und es funktioniert prima.

[Xeridar]

Ich grab mal diesen uralten Beitrag für eine Korrektur aus, da das einer der führenden Treffer in der Suche ist bei passenden Stichworten.

Zitat:

Zitat von Medium (Beitrag 896178)

Delphi-Quellcode:

// M : Mittelpunkt des Kreises
// r : Kreisradius
// A : Stützpunkt der zu schneidenden Strecke/Gerade
// V : Richtungsvektor der Strecke/Gerade (Wenn diese aus 2 Punkten A und B gebildet ist, ist dies B-A)
// U : Hilfsvektor
// p, q, d, x1, x2 : Hilfsvariablen (Skalare, keine Vektoren)

  U := A-M;
  p := U*V;
  q := U^2 - r^2;
  d := p^2 - q;

  if d < 0 then 
    Kein_Schnittpunkt
  else
  begin
    d  := sqrt(d / |V|);
    x1 := -p + d;
    x2 := -p - d;
  end;

  if (x2>x1) and (x1>0) then
    Schnittpunkt := A + x1*V
  else
  if (x1>x2) and (x2>0) then
    Schnittpunkt := A + x2*V
  else
    Kein_Schnittpunkt_Oder_Schnitt_Hinter_A;

Ich war zu faul zum selber rechnen, also hab ich Onkel Google gefragt, bin hier gelandet und dachte mir "hey, schön, fertiger Code zum übernehmen... ist zwar delphi statt java, aber ist ja selbsterklärend".

Ja, das funktioniert allerdings nicht ganz.

Hab jetzt doch selbst gerechnet und möchte das hier eben anmerken:
Das Problem liegt in diesen Zeilen:

Delphi-Quellcode:

p := U*V;
q := U^2 - r^2;
d := p^2 - q;
d := sqrt(d / |V|);

|V| soll wohl die quadrierte euklidische Norm sein, also V.x*V.x + V.y*V.y
Das ist elegant geschrieben, denn die Terme tauchen in der Rechnung wirklich als Nenner bei sowohl p als auch q auf.
So wie es hier geschrieben ist wird |V| aber nicht mit quadriert bei der Berechnung von "p*p - q". Außerdem müsste auch die Berechnung von x1 bzw. x2 dadurch falsch sein, wenn |V|!=1.

Delphi-Quellcode:

p := U*V / |V|;
q := (U^2 - r^2) / |V|;
d := p^2 - q;
d := sqrt(d);

Auf diese Lösung komme ich und das wirft auch ganz andere (und erwartete) Ergebnisse in meinem Programm.