Schnittpunkt Linie vs. Kreis

[NamenLozer]

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 ). 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

[NamenLozer]

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 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

[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.

[NamenLozer]

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]

m := 1E99 sollte reichen