Geometria obliczeniowa na płaszczyźnie
3.4. Punkty przecięcia
W niektórych sytuacjach stwierdzenie samego faktu przecina- Literatura [WDA] - 35.1 nia się obiektów geometrycznych jest niewystarczaj¸ace —
potrzeb-ne okazuje się również wyliczenie współrzędnych punktów przecię-cia. W niniejszym rozdziale zajmiemy się tego typu problemami.
Przedstawimy sposób wyznaczania punktów przecięcia prostych, odcinków, prostej z okręgiem oraz dwóch okręgów. Rozpoczniemy od przedstawienia problemu wyznaczania punktu prze-cięcia dwóch prostych.
Jeżeli przeanalizujemy zbiór punktów przecięcia prostych P1 → P2 oraz L1 → L2 to za-uważymy, że może zachodzić jeden z trzech przypadków: proste mog¸a być do siebie równoległe i się w ogóle nie przecinać, mog¸a się pokrywać czyli mieć nieskończenie wiele punktów prze-cięcia, albo mieć tylko jeden punkt wspólny. Pierwsze dwa przypadki w prosty sposób można zbadać, sprawdzaj¸ac wartość iloczynu wektorowego (P1 7→ P2) × (L1 7→ L2) — jeśli jest on równy 0, to znaczy, że proste s¸a równoległe. W takim przypadku trzeba jeszcze sprawdzić, czy się pokrywaj¸a, co można w łatwy sposób zweryfikować choćby poprzez zbadanie przy-należności punktu L1 do prostej P1→ P2.
Jeśli proste nie s¸a równoległe, to istnieje dokładnie jeden punkt przecięcia, który moż-na wyzmoż-naczyć, rozwi¸azuj¸ac odpowiedni układ rówmoż-nań. Zauważmy, że dowolny punkt (x, y) należ¸acy do prostej P1 → P2 można wyrazić w postaci (b ∗ P1.x+ (1 − b) ∗ P2.x, b∗ P1.y+ (1 − b) ∗ P2.y) dla pewnej rzeczywistej liczby b. Ponieważ poszukiwany punkt należy do obu prostych P1 → P2 oraz L1 → L2, to uzyskujemy w ten sposób układ równań z czterema
niewiadomymi x, y, b i c:
x= b ∗ x1+ (1 − b) ∗ x2
y= b ∗ y1+ (1 − b) ∗ y2
x= c ∗ x3+ (1 − c) ∗ x4
y = c ∗ y3+ (1 − c) ∗ y4
Wyznaczaj¸ac z niego niewiadome x oraz y, otrzymujemy:
x= x2∗ (x3∗ y4− x4∗ y3+ (x4− x3) ∗ y1) + x1∗ (x4∗ y3− x3∗ y4+ (x3− x4) ∗ y2) (x1− x2) ∗ (y3− y4) − (x3− x4) ∗ (y1− y2)
y= y2∗ (x3∗ y4− x4∗ y3) + y1∗ (x4∗ y3− x3∗ y4) + x2∗ y1∗ (y4− y3) + x1∗ y2∗ (y3− y4) (x1− x2) ∗ (y3− y4) − (x3− x4) ∗ (y1− y2)
Funkcj¸a realizuj¸ac¸a powyższ¸a metodę jest intLineCrossPoint(POINTD, POINTD, POINTD, POINTD, POINTD&), która jako parametry przyjmuje dwie proste wyznaczone przez pary punktów oraz dodatkowy punkt przekazywany przez referencję. Zwracana przez funkcję wartość może być następuj¸aca:
0 — proste nie przecinaj¸a się
1 — proste przecinaj¸a się w jednym punkcie
2 — proste pokrywaj¸a się (istnieje nieskończenie wiele punktów przecięcia)
W przypadku, gdy funkcja zwraca wartość 1, przekazany przez referencję punkt jest ustawiany na punkt przecięcia prostych. Implementacja tej funkcji przedstawiona jest na listingu 3.29.
Listing 3.29: Implementacja funkcji int LineCrossPoint(POINTD, POINTD, POINTD, POINTD, POINTD&)
1 int LineCrossPoint(POINTD p1, POINTD p2, POINTD l1, POINTD l2, POINTD & prz) { // Iloczyn wektorowy (p1 -> p2) i (l1 -> l2)
2 double s, t = (p1.x - p2.x) * (l1.y - l2.y) - (p1.y - p2.y) * (l1.x - l2.x);
// Iloczyn wektorowy (l2 -> p2) i (l1 -> l2)
3 s = (l2.x - p2.x) * (l1.y - l2.y) - (l2.y - p2.y) * (l1.x - l2.x);
// Jeśli proste są równoległe (t == 0), to istnieje nieskończenie wiele // punktów wspólnych wtw gdy proste się pokrywają (iloczyn wektorowy s == 0) 4 if (IsZero(t)) return IsZero(s) ? 2 : 0;
5 s = s / t;
// Istnieje jeden punkt wspólny - wyznacz jego współrzędne 6 prz.x = s * p1.x + (1 - s) * p2.x;
7 prz.y = s * p1.y + (1 - s) * p2.y;
8 return 1;
9 }
Problem wyznaczania punktów przecięcia dwóch odcinków z teoretycznego punktu widzenia jest stosunkowo podobny, jednak istniej¸a tu pewne różnice zwi¸azane ze sposobem reprezen-tacji wspólnego obszaru. Gdy odcinki częściowo się pokrywaj¸a, ich część wspólna nie jest
cał¸a prost¸a je zawieraj¸ac¸a. W takim przypadku, obszar stanowi¸acy ich przecięcie można reprezentować również przy użyciu odcinka. Jeśli odcinki się przecinaj¸a i nie s¸a równoległe, to sposób wyznaczania punktu przecięcia jest podobny do tego, wykorzystanego w przy-padku dwóch prostych. Funkcja bool SegmentCrossPoint(POINT, POINT, POINT, POINT, vector<POINTD>&)przyjmuje jako parametry dwa odcinki reprezentowane przez pary punk-tów oraz referencję na wektor punkpunk-tów v. Zwracanym wynikiem jest fałsz, jeżeli odcinki się nie przecinaj¸a, a prawda wpp. Jeśli część wspólna odcinków składa się tylko z jednego punk-tu, to jest on umieszczany w wektorze v. W przypadku, gdy odcinki pokrywaj¸a się, wektor v wypełniany jest dwoma punktami — końcami odcinka reprezentuj¸acego wspólny obszar.
Implementacja omawianej funkcji znajduje się na listingu 3.30.
Listing 3.30: Implementacja funkcji bool SegmentCrossPoint(POINT, POINT, POINT, POINT, vector<POINTD>&)
01 inline boolSegmentCrossPoint(POINT p1, POINT p2, POINT l1, POINT l2, 02 vector<POINTD> &r) {
03 r.clear();
04 int w1 = sgn(Det(p1, p2, l1)), w2 = sgn(Det(p1, p2, l2)), 05 v1 = sgn(Det(l1, l2, p1)), v2 = sgn(Det(l1, l2, p2));
// Jeśli punkty l1 i l2 znajdują się po tej samej stronie prostej p1 -> p2 // lub p1 i p2 znajdują się po tej samej stronie prostej l1 -> l2,
// to odcinki nie przecinają się
06 if (w1*w2 > 0 || v1*v2 > 0) return false;
// Jeśli punkty l1 i l2 leżą na prostej p1 -> p2, to odcinki // l1 -> l2 i p1 -> p2 są współliniowe
07 if (!w1 && !w2) {
// Zamiana par punktów reprezentujących odcinki, tak aby pierwsze punkty // w parach były mniejsze w porządku po współrzędnych (x,y)
08 if (OrdXY(&p2, &p1)) swap(p1, p2);
09 if (OrdXY(&l2, &l1)) swap(l1, l2);
// Jeśli odcinki są rozłączne, to nie ma punktów przecięcia 10 if (OrdXY(&p2, &l1) || OrdXY(&l2,&p1)) return false;
// Wyznacz krańcowe punkty wspólne
11 if (p2 == l1) r.PB(POINTD(p2.x, p2.y));
12 else if(p1 == l2) r.PB(POINTD(l2.x, l2.y));
13 else {
14 r.PB(OrdXY(&p1, &l1) ? l1 : p1);
15 r.PB(OrdXY(&p2, &l2) ? p2 : l2);
16 }
17 }
// Jeśli jeden z odcinków jest zdegenerowany, to jest on punktem przecięcia 18 else if (l1 == l2) r.PB(l1);
19 else if (p1 == p2) r.PB(p2);
20 else {
// Wyznacz punkt przecięcia
21 doublet = double(LL(l2.x - p2.x) * LL(l1.y - l2.y) - LL(l2.y - p2.y) * 22 LL(l1.x - l2.x)) / double(LL(p1.x - p2.x) * LL(l1.y - l2.y)
-23 LL(p1.y - p2.y) * LL(l1.x - l2.x));
24 r.PB(POINTD(t * p1.x + (1.0 - t) * p2.x, t * p1.y + (1.0 - t) * p2.y));
25 }
X Y
(-6, -3)
(5, 4) (-5, 4)
(-6, 1)
(3, 4) (4, 8)
Rysunek 3.6: Możliwe przypadki przecinania się zbioru odcinków oraz prostych je zawieraj¸acych
Listing 3.30: (c.d. listingu z poprzedniej strony) 26 return true;
27}
Listing 3.31: Wynik działania funkcji boolSegmentCrossPoint(POINT, POINT, POINT, POINT, vector<POINTD>&) oraz int LineCrossPoint(POINTD, POINTD, POINTD, POINTD, POINTD&) dla zbioru odcinków i prostych z rysunku 3.6
(4, 8) - (-6, 1) oraz (-5, 4) - (3, 4) Punkt przeciecia prostych: 1 (-1.71429, 4) Punkt przeciecia odcinkow: 1 (-1.71429, 4) (-6, -3) - (5, 4) oraz (-5, 4) - (3, 4) Punkt przeciecia prostych: 1 (5, 4) Punkt przeciecia odcinkow: 0
(-6, -3) - (5, 4) oraz (4, 8) - (-6, 1) Punkt przeciecia prostych: 1 (-68.8571, -43) Punkt przeciecia odcinkow: 0
Listing 3.32: Kod źródłowy programu użytego do wyznaczenia wyniku z listingu 3.31. Pełny kod źródłowy programu znajduje się w pliku linesegcross.cpp
01 int main() {
02 vector<POINT> b, e;
03 vector<POINTD> l;
04 intres;
05 POINT p1, p2;
06 POINTD p;
// Wczytaj wszystkie pary punktów wyznaczające proste i odcinki 07 while(cin >> p1.x >> p1.y >> p2.x >> p2.y) {
08 b.PB(p1); e.PB(p2);
09 }
// Dla każdej pary punktów wykonaj funkcje LineCrossPoint oraz SegmentCrossPoint 10 REP(x, SIZE(b)) REP(y, x) {
11
Listing 3.32: (c.d. listingu z poprzedniej strony)
12 cout << b[x] << " - " << e[x] <<
13 " oraz " << b[y] << " - " << e[y] << endl;
14 cout << "Punkt przeciecia prostych: " <<
15 (res = LineCrossPoint(b[x], e[x], b[y], e[y], p));
16 if (res == 1) cout << " " << p;
17 cout << endl;
18 cout << "Punkt przeciecia odcinkow: " <<
19 SegmentCrossPoint(b[x], e[x], b[y], e[y], l);
20 FOREACH(it, l) cout << " " << (*it);
21 cout << endl;
22 }
23 return 0;
24}
Kolejn¸a funkcj¸a, służ¸ac¸a do wyznaczania punktów przecięcia obiektów geometrycznych, jest vector<POINTD> LineCircleCross (POINTD, double, POINTD, POINTD). Wyznacza ona punkty przecięcia prostej z okręgiem. Przyjmowane przez tę funkcję parametry, to odpowied-nio punkt stanowi¸acy środek okręgu, jego promień oraz dwa punkty wyznaczaj¸ace prost¸a.
Wynikiem działania funkcji jest wektor punktów przecięcia. W zależności od wzajemnego położenia okręgu i prostej, wektor ten zawiera zero, jeden lub dwa elementy. W celu wyz-naczenia poszukiwanych punktów przecięcia między okręgiem o środku w punkcie P i promie-niu r oraz prostej wyznaczonej przez dwa punkty L1 oraz L2, należy rozwi¸azać następuj¸acy układ równań:
x= b ∗ L1.x+ (1 − b) ∗ L2.x y= b ∗ L1.y+ (1 − b) ∗ L2.y (x − P.x)2+ (y − P.y)2 = r2
Działanie omawianej funkcji sprowadza się w zasadzie do wyznaczenia rozwi¸azań tego układu równań.
Listing 3.33: Implementacja funkcji vector<POINTD> LineCircleCross(POINTD, double, POINTD, POINTD)
// Funkcja wyznacza punkty przecięcia okręgu i prostej
01vector<POINTD> LineCircleCross(POINTD p, doublecr, POINTD p1, POINTD p2){
02 double a = sqr(p1.x) + sqr(p1.y) + sqr(p2.x) + sqr(p2.y) -03 2.0 * (p1.x * p2.x + p1.y * p2.y);
04 double b = 2.0 * (p.x * (p2.x - p1.x) + p.y * (p2.y - p1.y) + 05 p1.x * p2.x + p1.y * p2.y - sqr(p2.x)-sqr(p2.y));
06 double c = -sqr(cr) + sqr(p2.x) + sqr(p2.y) + sqr(p.x) + 07 sqr(p.y) - 2.0 * (p.x * p2.x + p.y * p2.y);
08 double d = b * b - 4.0 * a * c;
09 vector<POINTD> r;
// Jeśli nie istnieje rozwiązanie równania kwadratowego, // to brak punktów przecięcia
10 if (d < -EPS) return r;
11 double t = -b / (2.0 * a), e = sqrt(abs(d)) / (2.0 * a);
12 if (IsZero(d))
Listing 3.33: (c.d. listingu z poprzedniej strony) // Istnieje tylko jeden punkt przecięcia...
13 r.PB(POINTD(t * p1.x + (1.0 - t) * p2.x, t * p1.y + (1.0 - t) * p2.y));
14 else {
// Istnieją dwa punkty przecięcia
15 r.PB(POINTD((t + e) * p1.x + (1.0 - t - e) * p2.x, 16 (t + e) * p1.y + (1.0 - t - e) * p2.y));
17 r.PB(POINTD((t - e) * p1.x + (1.0 - t + e) * p2.x, 18 (t - e) * p1.y + (1.0 - t + e) * p2.y));
19 }
20 returnr;
21}
Ostatni¸a prezentowan¸a funkcj¸a, służ¸ac¸a do wyznaczania punktów przecięcia, będzievector
<POINTD> CirclesCross(POINTD, double, POINTD, double). Jako parametry funkcja ta
przyjmuje środki oraz promienie dwóch okręgów (P1, r1, P2, r2), natomiast jako wynik zwraca wektor punktów przecięcia. Wyznaczenie punktów przecięcia dwóch okręgów sprowadza się jak i poprzednio do rozwi¸azania odpowiedniego układu równań.
( (x − P1.x)2+ (y − P1.y)2= r12 (x − P2.x)2+ (y − P2.y)2= r22
Listing 3.34: Implementacja funkcjivector<POINTD> CirclesCross(POINTD, double,POINTD, double)
01vector<POINTD> CirclesCross(POINTD c1, double c1r, POINTD c2, double c2r) { 02 vector<POINTD> r;
03 c2.x -= c1.x;
04 c2.y -= c1.y;
// Jeśli okręgi są współśrodkowe, to nie wyznaczamy punktów przecięcia 05 if (IsZero(c2.x) && IsZero(c2.y)) return r;
06 doubleA = (-sqr(c2r) + sqr(c1r) + sqr(c2.x) + sqr(c2.y)) / 2.0;
// Jeśli środki okręgów mają tą samą współrzędną y...
07 if (IsZero(c2.y)) { 08 double x = A / c2.x;
09 double y2 = sqr(c1r) - sqr(x);
10 if (y2 < -EPS) returnr;
// Jeśli okręgi są styczne...
11 if (IsZero(y2)) r.PB(POINTD(c1.x + x, c1.y));
12 else {
// Okręgi przecinają się
13 r.PB(POINTD(c1.x + x, c1.y + sqrt(y2)));
14 r.PB(POINTD(c1.x + x, c1.y - sqrt(y2)));
15 }
16 return r;
17 }
18 doublea = sqr(c2.x) + sqr(c2.y);
19 doubleb = -2.0 * A * c2.x;
20 doublec = A * A - sqr(c1r) * sqr(c2.y);
X Y
(8, -3) (-8, -3)
(-5, 0) (2, 0)
(4, 5)
Rysunek 3.7: Różne możliwości przecinania się okręgów oraz prostych
Listing 3.34: (c.d. listingu z poprzedniej strony) 21 double d = b * b - 4.0 * a * c;
22 if (d < -EPS) return r;
23 double x = -b / (2.0 * a);
// Jeśli okręgi są styczne...
24 if (IsZero(d)) r.PB(POINTD(c1.x + x, c1.y + (A - c2.x * x) / c2.y));
25 else {
// Okręgi przecinają się
26 doublee = sqrt(d) / (2.0 * a);
27 r.PB(POINTD(c1.x + x + e, c1.y + (A - c2.x * (x + e)) / c2.y));
28 r.PB(POINTD(c1.x + x - e, c1.y + (A - c2.x * (x - e)) / c2.y));
29 }
30 return r;
31}
Listing 3.35: Wynik działania funkcji vector<POINTD> CirclesCross(POINTD, dou-ble, POINTD, double) oraz vector<POINTD> LineCircleCross(POINTD, double, POINTD, POINTD)na przykładzie rysunku 3.7
Przeciecie prostej (-8, -3) - (8, -3) i okregu:
(-5, 0, 3) - (-5, -3)
(2, 0, 4) - (-0.645751, -3) (4.64575, -3) (4, 5, 3) - brak punktow przeciecia Przeciecie okregow:
(2, 0, 4) oraz (-5, 0, 3) - (-2, 0)
(4, 5, 3) oraz (-5, 0, 3) - brak punktow przeciecia
(4, 5, 3) oraz (2, 0, 4) - (5.28141, 2.28744) (1.20135, 3.91946)
Listing 3.36: Kod źródłowy programu użytego do wyznaczenia wyniku z listingu 3.35. Pełny kod źródłowy programu znajduje się w pliku linecirclecross.cpp
01 intmain() {
02 vector<POINT> cCen;
03 VI cR;
04 vector<POINTD> res;
05 POINT l1, l2, p; intr;
06
Listing 3.36: (c.d. listingu z poprzedniej strony) // Wczytaj współrzędne punktów wyznaczających prostą
07 cin >> l1.x >> l1.y >> l2.x >> l2.y;
// Wczytaj opisy kolejnych okręgów 08 while(cin >> p.x >> p.y >> r) { 09 cR.PB(r); cCen.PB(p);
10 }
// Wyznacz punkty przecięcia prostej i okręgu
11 cout << "Przeciecie prostej " << l1 << " - " << l2 <<
12 " i okregu:" << endl;
13 REP(x, SIZE(cCen)) {
14 res = LineCircleCross(cCen[x], cR[x], l1, l2);
15 cout << "\t(" << cCen[x].x << ", " << cCen[x].y <<
16 ", " << cR[x] << ") - ";
17 if (!SIZE(res)) cout << "brak punktow przeciecia";
18 FOREACH(it, res) cout << " " << (*it);
19 cout << endl;
20 }
// Wyznacz punkty przecięcia dwóch okręgów
21 cout << "Przeciecie okregow:" << endl;
22 REP(x, SIZE(cCen)) REP(y, x) {
23 res = CirclesCross(cCen[x], cR[x], cCen[y], cR[y]);
24 cout << "\t(" << cCen[x].x << ", " << cCen[x].y << ", " << cR[x] <<
25 ") oraz (" <<
26 cCen[y].x << ", " << cCen[y].y << ", " << cR[y] << ") - ";
27 if (!SIZE(res)) cout << "brak punktow przeciecia";
28 FOREACH(it, res) cout << " " << *(it);
29 cout << endl;
30 }
31 return0;
32}