Geometria obliczeniowa na płaszczyźnie
3.3. Przynależność punktu do figury
Rozwi¸azanie wielu zadań wymaga umiejętności sprawdzania, Literatura [WDA] - 35.1
[ASD] - 8.2 czy pewne punkty należ¸a do różnych obiektów geometrycznych,
ta-kich jak odcinki, proste, wielok¸aty czy koła. W aktualnym rozdziale przedstawione zostan¸a metody umożliwiaj¸ace dokonywanie takich testów. Analizę problemu rozpoczniemy od bardzo prostego
przy-padku sprawdzania przynależności punktu do prostok¸ata o bokach równoległych do osi układu współrzędnych.
W celu sprawdzenia czy punkt A = (x, y) znajduje się we wnętrzu prostok¸ata wyznac-zonego przez dwa jego przeciwległe wierzchołki B = (x1, y1) oraz C = (x2, y2), wystarczy sprawdzić, czy odcięta punktu leży pomiędzy odciętymi wierzchołków prostok¸ata oraz czy podobna zależność jest zachowana przez rzędne. Nierówności, które musz¸a zostać spełnione s¸a następuj¸ace:
min(x1, x2) < x < max(x1, x2), min(y1, y2) < y < max(y1, y2)
W przypadku sprawdzania, czy punkt należy do wnętrza prostok¸ata lub leży na jego obwodzie, należy zamienić ostre nierówności z poprzednich warunków na nieostre:
min(x1, x2) ¬ x ¬ max(x1, x2), min(y1, y2) ¬ y ¬ max(y1, y2)
Do weryfikowania tych dwóch rodzajów przynależności dostępne s¸a makraPointInRectoraz PointInsideRect. Oba makra przyjmuj¸a jako parametry dwa punkty wyznaczaj¸ace pros-tok¸at oraz punkt, dla którego badana jest przynależność. W przypadku, gdy punkt należy do prostok¸ata, makra przyjmuj¸a wartość prawda. Makra te zostały przedstawione na listingach 3.14 oraz 3.15.
Listing 3.14: Implementacja makraPointInsideRect 1 #define PointInsideRect(p1,p2,p3) (min(p1.x, p2.x) < p3.x && \
2 min(p1.y, p2.y) < p3.y && max(p1.x, p2.x) > p3.x && max(p1.y, p2.y) > p3.y)
Listing 3.15: Implementacja makraPointInRect 1 #define PointInRect(p1,p2,p3) (min(p1.x, p2.x) <= p3.x && \
2 min(p1.y, p2.y) <= p3.y && max(p1.x, p2.x) >= p3.x && max(p1.y, p2.y) >= p3.y)
Analogicznymi makrami s¸a PointInsideSegment oraz PointInSegment. Ich zadaniem jest sprawdzanie, czy punkt podany jako trzeci parametr leży wewn¸atrz b¸adź na odcinku wyz-naczonym przez punkty stanowi¸ace pierwszy oraz drugi parametr.
Implementacja tych dwóch makr wykorzystuje zarówno poprzednie makra do sprawdza-nia przynależności punktu do prostok¸ata, jak również iloczyn wektorowy. Najpierw jest
X Y
(-8, 6)
(-3, -8) (-3, 0)
(3, -4) (-2, -6)
(4, 8)
(-8, 2)
(-7, 4)
(7, 7)
Rysunek 3.3: Przykładowy zbiór punktów prezentuj¸acy różne przypadki przynależności punktów do obiektów geometrycznych — odcinka, prostok¸ata i koła
sprawdzane, czy testowany punkt leży na prostej wyznaczonej przez odcinek. Jeśli tak, to musi on jeszcze należeć do prostok¸ata wyznaczonego przez końce odcinka. Implementacja tych makr przedstawiona jest na listingach 3.16 i 3.17.
Listing 3.16: Implementacja makraPointInsideSegment
1 #define PointInsideSegment(p1,p2,l) (Det(p1,p2,l)==0 && PointInsideRect(p1,p2,l))
Listing 3.17: Implementacja makraPointInSegment
1 #define PointInSegment(p1,p2,l) (Det(p1,p2,l)==0 && PointInRect(p1,p2,l))
Ostatnimi z serii prostych i krótkich makr, s¸a PointInsideCircle oraz PointInCircle.
Służ¸a one odpowiednio do sprawdzania, czy punkt określony przez trzeci parametr makra należy do wnętrza b¸adź wnętrza lub obwodu koła o środku w punkcie określonym przez pier-wszy parametr makra oraz promieniu zdefiniowanym przez drugi parametr. Implementacje tych makr zostały umieszczone odpowiednio na listingach 3.18 i 3.19.
Listing 3.18: Implementacja makraPointInsideCircle
1 #define PointInsideCircle(c,r,p) (sqr(c.x-p.x)+sqr(c.y-p.y)+EPS < sqr(r))
Listing 3.19: Implementacja makraPointInCircle 1 #define PointInCircle(c,r,p) (sqr(c.x-p.x)+sqr(c.y-p.y)-EPS < sqr(r))
Listing 3.20: Wynik działania makr służ¸acych do sprawdzania przynależności punktów do różnych obiektów geometrycznych na przykładzie odcinka o końcach w punktach (−2, −6) i (−8, 6), okręgu o środku w punkcie (−3, 2) i promieniu 5 oraz prostok¸ata wyznaczonego przez punkty (−3, −8) i (7, 7) (patrz rysunek 3.3).
Rectangle Segment Circle
Inside In Inside In Inside In
(-8, 2) 0 0 0 0 0 1
(-7, 4) 0 0 1 1 1 1
Listing 3.20: (c.d. listingu z poprzedniej strony)
(-8, 6) 0 0 0 1 0 0
(-3, 0) 0 1 0 0 1 1
(3, -4) 1 1 0 0 0 0
(4, 8) 0 0 0 0 0 0
Listing 3.21: Kod źródłowy programu użytego do wyznaczenia wyniku z listingu 3.20. Pełny kod źródłowy programu znajduje się w pliku pointin inside.cpp
01 intmain() {
02 POINT r1,r2,c,s1,s2,p;
03 int r;
// Wczytaj współrzędne wierzchołków wyznaczających prostokąt
04 cin >> r1.x >> r1.y >> r2.x >> r2.y;
// Wczytaj współrzędne punktów wyznaczających odcinek
05 cin >> s1.x >> s1.y >> s2.x >> s2.y;
// Wczytaj środek oraz promień okręgu
06 cin >> c.x >> c.y >> r;
07 cout << "\t\t\t Rectangle\t\t Segment\t\t Circle" << endl;
08 cout << "\t\t\tInside\t In\t Inside\t In\t Inside\t In" << endl;
// Dla wszystkich punktów wyznacz wynik działania poszczególnych makr 09 while(cin >> p.x >> p.y) {
10 cout << p << "\t\t" << PointInsideRect(r1, r2, p) << "\t\t" <<
11 PointInRect(r1, r2, p) << "\t\t" <<
12 PointInsideSegment(s1, s2, p) << "\t\t" <<
13 PointInSegment(s1, s2, p) << "\t\t" <<
14 PointInsideCircle(c, r, p) << "\t\t" <<
15 PointInCircle(c, r, p) << endl;
16 }
17 return 0;
18}
Bardziej interesuj¸acym zagadnieniem jest sprawdzanie przynależności punktu do wielok¸ata.
Zajmiemy się tym problemem przy założeniu, że wszystkie punkty stanowi¸ace wierzchoł-ki wielok¸ata maj¸a całkowitoliczbowe współrzędne. W celu rozwi¸azania tego problemu, is-totn¸a obserwacj¸a jest to, że dowolna półprosta o pocz¸atku w punkcie należ¸acym do wnętrza wielok¸ata przecina jego obwód nieparzyst¸a liczbę razy. Sytuacja odwrotna występuje w przy-padku punktów nienależ¸acych do wielok¸ata — wtedy liczba przecięć jest parzysta. Zatem w celu sprawdzenia przynależności punktu do wielok¸ata można sprawdzić, z iloma boka-mi wielok¸ata przecina się dowolna półprosta o pocz¸atku w tym punkcie. Wybór jakiejkol-wiek półprostej niesie ze sob¸a jednak wiele problemów. Jak przedstawiono na rysunku 3.4, prosta taka może przecinać obwód wielok¸ata na różne sposoby, co wi¸aże się z dodatkowymi przypadkami do rozpatrzenia, znacznie komplikuj¸acymi implementację. W przypadku przed-stawionym na rysunku 3.4.a, punkt przecięcia półprostej z obwodem wielok¸ata powinien być liczony podwójnie, podczas gdy sytuacja z rysunku 3.4.b — pojedynczo.
Widać, że odpowiedni wybór półprostej jest bardzo ważny — najlepiej byłoby wybrać j¸a w ten sposób, aby nie przechodziła przez żaden wierzchołek wielok¸ata. Okazuje się, że
(a) (b) (c)
Rysunek 3.4: Różne przypadki przecinania obwodu wielok¸ata przez półprost¸a.
wyboru takiego można dokonać w bardzo prosty sposób — dla testowanego punktu (x, y), półprosta przechodz¸aca przez punkt (INF, y+1), nie przecina żadnego wierzchołka wielok¸ata, co wynika bezpośrednio z faktu, iż wszystkie punkty maj¸a całkowitoliczbowe współrzędne.
Implementacja algorytmu wykorzystuj¸acego ten fakt została zrealizowana w funkcjach boolPointInsidePol(vector<POINT>&, POINT)oraz boolPointInPol(vector<POINT>&, POINT)przedstawionych na listingach 3.23 i 3.24 — przyjmuj¸a one jako parametry odpowied-nio wektor wierzchołków wielok¸ata oraz punkt, dla którego przeprowadzany jest test. Pier-wsza z tych funkcji zwraca prawdę, jeśli punkt należy do wnętrza wielok¸ata, natomiast druga zwraca prawdę, jeśli punkt należy do wnętrza lub obwodu wielok¸ata. Czas ich działania jest liniowy — dla każdego boku wielok¸ata funkcje sprawdzaj¸a, czy przecina on półprost¸a o pocz¸atku w testowanym punkcie. Obie przedstawione poniżej funkcje wykorzystuj¸a funkcję bool SegmentCross(POINT&, POINT&, POINT&, POINT&), która przyjmuje jako parame-try dwa odcinki reprezentowane przez pary wierzchołków oraz zwara prawdę, jeśli odcinki te przecinaj¸a się. Jej implementacja została przedstawiona na listingu 3.22.
Listing 3.22: Implementacja funkcji boolSegmentCross(POINT&, POINT&, POINT&, POINT&) 1 inline bool SegmentCross(POINT & p1, POINT & p2, POINT & l1, POINT & l2) { 2 return sgn(Det(p1, p2, l1)) * sgn(Det(p1, p2, l2)) == -1
3 && sgn(Det(l1, l2, p1)) * sgn(Det(l1, l2, p2)) == -1;
4 }
Listing 3.23: Implementacja funkcji boolPointInsidePol(vector<POINT>&, POINT) 1 bool PointInsidePol(vector<POINT> &l, POINT p) {
2 int v = 0, s = SIZE(l);
3 POINT d(INF, p.y + 1);
// Jeśli punkt leży na jednym z boków wielokąta, to nie należy do wnętrza // wielokąta
4 REP(x, s) if(PointInSegment(l[x], l[(x + 1) % s], p)) return false;
// Wyznacz liczbę przecięć obwodu wielokąta z półprostą (p -> d) 5 REP(x, s) v += SegmentCross(p, d, l[x], l[(x + 1) % s]);
// Jeśli półprosta przecina obwód nieparzystą liczbę razy, to punkt należy do // wielokąta
6 return v & 1;
7 }
Listing 3.24: Implementacja funkcji bool PointInPol(vector<POINT>&, POINT) 1 bool PointInPol(vector<POINT> &l, POINT p) {
2 intv = 0, s = SIZE(l);
3 POINT d(INF, p.y + 1);
// Jeśli punkt leży na jednym z boków wielokąta, to zwróć prawdę 4 REP(x, s) if(PointInSegment(l[x], l[(x + 1) % s], p)) return true;
// Wyznacz liczbę przecięć obwodu wielokąta z półprostą (p -> d) 5 REP(x, s) v += SegmentCross(p, d, l[x], l[(x + 1) % s]);
// Jeśli półprosta przecina obwód nieparzystą liczbę razy, to punkt należy do // wielokąta
6 return v & 1;
7 }
W przypadku rozpatrywania przynależności punktu do wielok¸ata wypukłego, wcześniej omów-ione funkcje jak najbardziej działaj¸a poprawnie, jednak nie s¸a optymalne. W takim przy-padku, istnieje możliwość sprawdzenia przynależności punktu w czasie logarytmicznym ze względu na liczbę wierzchołków wielok¸ata.
Zasada działania takiego algorytmu opiera się na metodzie wyszukiwania binarnego — na pocz¸atku algorytm wybiera sobie punkt należ¸acy do wielok¸ata wypukłego (może to być przykładowo jego wierzchołek). Następnie wnętrze wielok¸ata dzielone jest na trójk¸aty (tak jak przedstawiono to na rysunku 3.2.a) przy użyciu rodziny półprostych P1, P2, . . .Jeśli punkt należy do wielok¸ata, to musi znajdować się w jednym z powstałych trójk¸atów. Przy użyciu iloczynu wektorowego, za pomoc¸a którego sprawdzane jest położenie punktu względem kole-jno wybieranych prostych, przeprowadzane jest wyszukiwanie binarne, w celu wyznaczenia pary s¸asiednich półprostych Pk i Pk+1, pomiędzy którymi leży testowany punkt. Wówczas jeśli znajduje się on w trójk¸acie wyznaczonym przez półproste Pk, Pk+1 oraz odpowiedni bok wielok¸ata, to należy on również do wielok¸ata. Realizacja tej metody jest przedstaw-iona na listingach 3.25 oraz 3.26 w postaci funkcji bool PointInsideConvexPol(vector
<POINT>&, POINT)oraz boolPointInConvexPol(vector <POINT>&, POINT).
Przyjmu-j¸a one jako parametry listę wierzchołków wielok¸ata wypukłego oraz punkt, dla którego prze-prowadzany jest test.
Listing 3.25: Implementacja funkcji bool PointInsideConvexPol(vector<POINT>&, POINT) 01 boolPointInsideConvexPol(vector<POINT> &l, POINT p) {
02 int a = 1, b = SIZE(l) - 1, c;
// Jeśli odcinek (l[0] -> l[a]) leży na prawo od odcinka (l[0] -> l[b]) to // następuje zamiana
03 if (Det(l[0], l[a], l[b]) > 0) swap(a, b);
// Jeśli punkt p nie leży po prawej stronie prostej (l[0] -> l[a]) lub po // lewej stronie (l[0] -> l[b]) to nie należy do wielokąta
04 if (Det(l[0], l[a], p) >= 0 || Det(l[0], l[b], p) <= 0) return false;
// Wyszukiwanie binarne wycinka płaszczyzny zawierającego punkt p 05 while (abs(a - b) > 1) {
06 c = (a + b) / 2;
07 if (Det(l[0], l[c], p) > 0) b = c;
08 else a = c;
09 }
// Jeśli punkt p leży w trójkącie (l[0],l[a],l[b]), to należy do wielokąta
X Y
(-7, 5)
(-3, 4)
(-3, -4)
(1, 8)
(6, 1)
(5, -2) (7, 4)
(3, -4)
Rysunek 3.5: Przykładowy wielok¸at wypukły oraz zbiór trzech punktów (3, −4), (6, 1) i (−3, 4), przed-stawiaj¸acych różne przypadki przynależności punktu do wielok¸ata.
Listing 3.25: (c.d. listingu z poprzedniej strony) 10 returnDet(l[a], l[b], p) < 0;
11}
Listing 3.26: Implementacja funkcji boolPointInConvexPol(vector<POINT>&, POINT) 01 bool PointInConvexPol(vector<POINT> &l, POINT p) {
02 inta = 1, b = SIZE(l) - 1, c;
// Jeśli odcinek (l[0] -> l[a]) leży na prawo od odcinka (l[0] -> l[b]) to // następuje zamiana
03 if (Det(l[0], l[a], l[b]) > 0) swap(a, b);
// Jeśli punkt p leży po lewej stronie prostej (l[0] -> l[a]) lub po // prawej stronie (l[0] -> l[b]) to nie należy do wielokąta
04 if (Det(l[0], l[a], p) > 0 || Det(l[0], l[b], p) < 0) return false;
// Wyszukiwanie binarne wycinka płaszczyzny zawierającego punkt p 05 while(abs(a - b) > 1) {
06 c = (a + b) / 2;
07 if (Det(l[0], l[c], p) > 0) b = c;
08 else a = c;
09 }
// Jeśli punkt p leży w trójkącie (l[0],l[a],l[b]), to należy do wielokąta 10 returnDet(l[a], l[b], p) <= 0;
11}
Listing 3.27: Przykład wykorzystania funkcji sprawdzaj¸acych przynależności punktów do wielok¸ata na podstawie zbioru punktów z rysunku 3.5
Polygon Convex Polygon
Inside In Inside In
(-3, 4) 1 1 1 1
(1, 8) 0 1 0 1
(6, 1) 0 1 0 1
(3, -4) 0 0 0 0
Listing 3.28: Kod źródłowy programu użytego do wyznaczenia wyniku z listingu 3.27. Pełny kod źródłowy programu znajduje się w pliku pointincpol.cpp
01 intmain() {
02 vector<POINT> pol;
03 POINT p;
04 int n;
// Wczytaj liczbę wierzchołków wielokąta
05 cin >> n;
// Wczytaj kolejne wierzchołki wielokąta oraz dodaj je do wektora 06 REP(x, n) {
07 cin >> p.x >> p.y;
08 pol.PB(p);
09 }
10 cout << "\t\t\t\t Polygon\t\t\t\t Convex Polygon" << endl;
11 cout << "\t\t Inside\t\t In\t\t Inside\t\t In" << endl;
// Dla kolejnych punktów wyznacz ich przynależność do wielokąta przy użyciu // poszczególnych funkcji
12 while(cin >> p.x >> p.y) cout << p << "\t\t" <<
13 PointInsidePol(pol, p) << "\t\t\t" <<
14 PointInPol(pol, p) << "\t\t\t" <<
15 PointInsideConvexPol(pol, p) << "\t\t\t" <<
16 PointInConvexPol(pol, p) << endl;
17 return 0;
18}