znajdowanie najwiekszego pustego trojkata w zbiorze

0

Witam, mam taki problem ze mam zbior punktow w otoczce wypuklej i musze znalezc najwiekszy trojkat na tych punktach ktory nie zawiera zadnego innego punktu w srodku
tutaj ilustracja:
user image
http://i.imagebanana.com/img/rlb6cdbw/Zaznaczenie_003.png
aktualnie robie to naiwnie sprawdzajac wszystkie kombinacje 4 punktow, ale zlozonosc O(n^4) jest zbyt duzafor ( i = 0; i < N ; i++ )
for ( j = 0; j < N; j++ )
for ( k = 0; k < N; k++ )
{
for (m = 0; m < N; m++)
{
if(m == k || m == j || m == i)
continue ;
if(jestWewnatrz(x[m],y[m],x[k],y[k],x[i],y[i],x[j],y[j]))
break ;
}
if(m == N)
{
++F;
PoleTrojkata = poleTrojkata(x[i],y[i],x[j],y[j],x[k],y[k]) ;
if(PoleTrojkata > maxPole)
{
a1 = point(x[j],y[j]) ;
b1 = point(x[i],y[i]) ;
c1 = point(x[k],y[k]) ;
maxPole = PoleTrojkata ;
}
}

            } 

Ma ktos pomysl jak mozna to przyspieszyc?
Dziekuje za kazdy pomysl.
0
barillaaaa napisał(a):

Witam, mam taki problem ze mam zbior punktow w otoczce wypuklej i musze znalezc najwiekszy trojkat na tych punktach ktory nie zawiera zadnego innego punktu w srodku

Ciekawe zadanie, skąd ono pochodzi?

barillaaaa napisał(a):

for ( i = 0; i < N ; i++ )
for ( j = 0; j < N; j++ )
for ( k = 0; k < N; k++ )

for ( i = 0; i < N ; i++ )
        for ( j = i+1; j < N; j++ )
                for ( k = j+1; k < N; k++ )

A do wewnętrznej pętli zastosuj strukturę zwaną KD-Tree.

0

Dzieki za pomoc, zadanie dostałem do rozwiązania na analizie algorytmów, wymyślone przez prowadzącego.

0
barillaaaa napisał(a):

Dzieki za pomoc, zadanie dostałem do rozwiązania na analizie algorytmów, wymyślone przez prowadzącego.

Myślę dalej na tym zadaniem, pewnie ktoś myślał przede mną i odkrywam żarówkę... ale nie wiem co wpisać w google.
Intuicja podpowiada mi, że da się to zadanie szybko policzyć, ale udowodnić ciężko.

Wydaje się że nie da się uniknąć pierwszej pary pętli, czyli szukamy
każdej pary punktów:

for( i=0 ; i<N-1 ; i++ ) {
for( j=i+1 ; j<N ; j++ ) {
  p1 = points[i];
  p2 = points[j];
}}

Potem szukamy punktu p3 leżącego blisko odcinka |p1,p2| Dalej tworzymy
dwie półproste p1,p3 i p2,p3. Punkty leżące tak jakby za punktem p3 i
pomiędzy półprostymi, na pewno nie będą tworzyły pustego trójkąta z
punktami p1 i p2. Czyli dla tych punktów nie trzeba już badać.

Generalnie każdy punkt p3 będzie odrzucał dużą część punktów z którymi
na pewno nie da się utworzyć pustego trójkąta. Jeśli dużo odrzuca, to
do sprawdzenia pozostaje mało. Pozostaje tylko pytanie, czy da się
zaprojektować strukturę danych która uwzględnia tę korzyść.

3

Chciałem dzisiaj pomyśleć nad tym zadaniem, ale niestety nie miałem zbyt dużo czasu, a jako że wygląda na to że @mariotti myśli równolegle ze mną to muszę się pośpieszyć i wrzucić na razie nieskończone efekty tego myślenia ;].

Chciałem się pochwalić jak znajdę jakiś dobry algorytm, ale na razie zdążyłem tylko naiwne praktycznie brute-force w \Theta(n^3) :/. No ale konkurencja ściga, więc odpisuję teraz.

Algorytm naiwny - sprawdzanie wszystkich możliwych kombinacji i naiwne sprawdzanie czy zawierają punkty:
T(n) \in \Theta(n^4)
(etap: gotowe do implementacji)
To już sam zrobiłeś, nie ma się co nad tym rozwodzić ;).

Algorytm naiwny + jakaś drzewowa struktura danych - np. odpowiednie KD-drzewo jak @mariotti sugerował
T(n) \in \Theta(n^3 log n)
(etap: gotowe do implementacji)
Trochę kombinowania z napisaniem i testowaniem tego drzewa. W dodatku obawiam się że takie rozwiązanie wiele zysku nie przyniesie (fakt, złożoność lepsza, ale przy mniejszych n będzie odczuwalne mocne zwiększenie się niepozornej i zawsze pomijanej stałej - a większych n nie będzie, no bo przy złożoności n^3 i tak większe zbiory danych są nieosiągalne).

Algorytm dopasowywania punktu do odcinka
T(n) \in \Theta(n^3 log n)
(etap: prawie gotowe do implementacji)
Rozważanie wszystkich możliwych n^2 odcinków i znalezienie do każdego z nich pasującego punktu. Chyba o to chodzi @mariotti -emu w ostatnim poście ;].
Rozważamy po kolei wszystkie n^2 par punktów (A, B). Sortujemy najpierw wszystkie punkty po odległości od prostej AB, i jeśli będziemy rozważać punkty od najbliższych do najdalszych to jesteśmy w stanie w czasie stałym sprawdzać czy kolejny punkt tworzy większy legalny trójkąt, oraz w czasie log(n) dodawać ten punkt do wykorzystywanej pomocniczej struktury danych.
Nad tym trochę myślałem, ale takie podejście nie ma przyszłości - poniżej n<sup>3 log n sie nie da zejśc z powodu n</sup>2 sortowań n punktów.
Nie będę się nad tym więcej rozwodził bo mam lepszy algorytm, tylko się chwalę że coś tam jeszcze wymyśliłem ;].

Algorytm naiwny + naiwna struktuta danych - policzenie jakimś algorytmem dynamicznym ilości punktów wewnątrz dla każdego możliwego trójkąta
T(n) \in \Theta(n^3)
(etap: 'wydaje mi się')
Wydaje mi się że takie coś by było osiągalne w czasie stałym dla każdego trójkąta przez odpowiednie przetworzenie ich wstępnie i rozważanie w odpowiedniej kolejnosci 'wstępująco' - jak pisałem, 'wydaje mi sie', nie myślałem nad tym za dużo.

Algorytm naiwny + zamiatanie.
T(n) \in \Theta(n^3)
(etap: 'zaimplementowany')
Jak na razie najlepsze do czego doszedłem, a doszedłem do niego przed chwilą z herbatą i ołówkiem w ręce ;].

O co chodzi - jakaś skrócona lista kroków:

  • sortujemy punkty po współrzędnej X (kwestia umowna oczywiście).
  • przeglądamy po kolei punkty (będziemy rozważać trójkąty znajdujące się na prawo):
    • bierzemy kolejny punkt, nazwijmy go A
    • sortujemy kątowo wszystkie punkty z prawej strony (te z lewej możemy równie dobrze wywalać przy przechodzeniu na przykład) względem punktu A
    • no i full naiwność (stąd niestety czas n^3) - bierzemy po kolei wszystkie punkty, w kolejności z sortowania kątowego, nazwijmy kolejny C
      • bierzemy po kolei wszystkie punkty po B w kolejności z sortowania kątowego, nazwijmy kolejny taki punkt C
        • możemy, znając najgorszy możliwy odcinający kąt, sprawdzić w czasie stałym czy ABC zawiera jakiś trójkąt.

Daje to narodziny zaskakująco prostemu algorytmowi (mniej niż 50 linijek faktycznego kodu, żadnych skomplikowanych struktur danych)...

Jako że talk is cheap...

(ponad dwie godziny później. Jak zwykle zamotałem się przy znakach cross produktu oraz pół godziny walczyłem z prostym IO w pythonie. No ale wygląda na to że się udało)

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <algorithm>
using namespace std;

struct point {
    point() { }
    point(int x, int y) :x(x), y(y) { }
    int x, y;
};

struct cmp_by_x {
    bool operator()(const point &a, const point &b) const {
        return a.x < b.x;
    }
};

int cross(const point &a, const point &b, const point &c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}

struct cmp_by_angle {
    point base;
    cmp_by_angle(point base) :base(base) { }
    bool operator()(const point &a, const point &b) const {
        return cross(base, a, b) > 0;
    }
};

bool on_plane(point &a, point &b, point &c) { return cross(a, b, c) < 0; }
double area(point &a, point &b, point &c) { return abs(cross(a, b, c) / 2.0); }
point rnd_pt() { return point(rand() % 30, rand() % 30); }

int main() {
    vector<point> pts(20); generate(pts.begin(), pts.end(), rnd_pt);
    sort(pts.begin(), pts.end(), cmp_by_x());
    double best = 0;
    for (size_t i = 0; i < pts.size(); i++) {
        vector<point> ang(pts.begin() + i + 1, pts.end());
        sort(ang.begin(), ang.end(), cmp_by_angle(pts[i]));
        for (size_t b = 0, cut = b; b < ang.size(); b++) {
            for (size_t c = b + 1; c < ang.size(); c++) {
                if (cut == b || !on_plane(ang[b], ang[cut], ang[c])) {
                    cut = c; best = max(best, area(pts[i], ang[b], ang[c]));
                }
            }
        }
    } printf("%.1f\n", best);
}

I... to tyle. Czas działania to okrągłe n^3 (bo czasy sortowania sa dominowane przez przeglądanie wszystkich wierzchołków).

Napisałem też prosty (trywialny wręcz) skrypt w pythonie do wisualizowania punktów, żeby mieć pewność że pole (a ściśle, punkty na brzegach) są wyznaczane dobrze. Jeśli ktoś jest zainteresowany, wrzucę tutaj.

Zakładam że punkt wewnątrz trójkąta to nie punkt na jego krawędzi - jeśli tak nie jest, wystarczy zmienić < na <=.

No to tyle do tego algorytmu.

Algorytm mniej naiwny + zamiatanie.
T(n) \in \Theta(n^2 log n)
(etap: 'wydaje mi się')
Wydaje mi się, że sprawdzanie trójkątów wychodzących z punktu można zrobić lepiej niż w n<sup>2. A ściśle, w n * log n (za pomocą jakiegoś podejścia dynamicznego), co dawałoby czas n</sup>2 log n - to już powoli zaczyna przypominać praktyczną złożoność (chociaż i tak bardzo sporą). Jak na razie mi się wydaje, spróbuję jutro nad tym pomyśleć w wolnym czasie.

Marzenia?
T(n) \in \Theta(n^2)
(etap: 'wydaje mi się')
Wydaje mi się, że ten problem da się rozwiązać jeszcze wydajniej, co najmniej n^2 (a może jeszcze lepiej?). To tylko intuicja, nawet jeśli się da być może nie wpadnę na to rozwiązanie ;].

0

Dziekuje bardzo za pomoc, to rozwiazanie z O(n^3) jest calkiem spoko, ja myslalem nad budowaniem grafu i minimalnym drzewem rozpinajacym w ogole w inna strone kombinowalem. Jeszcze raz dzięki i jakbyście mieli jeszcze jakieś ciekawe pomysły to na pewno się przydadzą.

0

Hm, jednak jest jeden problem, sam algorytm wydaje się słuszny ale dla niektórych danych metoda naiwna daje inne lepsze wyniki np dla zestawu danych 150 punktow.
http://i.imagebanana.com/img/q1r8ckx6/Zaznaczenie_004.png

http://4programmers.net/Pastebin/2176
(Wrzuciłem wejście na pastebin - msm)

0

A dla większości danych działa?
Szczerze mówiąc nie porównałem wyników z bruteforcem (jedynie sprawdzałem czy trójkąt jest poprawnie wyznaczony i wygląda na największy), ale jeśli w większości przypadków działa to

  • Twój program odrzuca czy akceptuje trójkąty które mają punkty leżące dokładnie na krawędzi trójkąta? - to by był dobry powód na dawanie innych wyników.
  • Twój program działa tylko na intach czy używa floatów do testowania punktów (błędy zaokrągleń, niestety)

Porównam z jakimś bruteforce za chwilę i powiem co mi wyszło.

http://i.imagebanana.com/img/q1r8ckx6/Zaznaczenie_004.png
Ten po prawej wygląda lepiej ;]

0

Wszystko sie dzieje na calkowitych, pole licze w ten sam sposob i punkty lezace na prostej nie naleza do trojkata w obu programach. Dziwne bo jak patrze na ten rysunek to ten algorytm ktory wymysliles powinien znalezc ten wiekszy trojkat. Ale jeszcze popatrze na to wszystko od poczatku bo mam straszny balagan w kodzie, pomieszane funkcje opengl od rysowania z samym algorytmem takze musze to uporzadkowac;p

0
barillaaaa napisał(a):

Hm, jednak jest jeden problem, sam algorytm wydaje się słuszny ale dla niektórych danych metoda naiwna daje inne lepsze wyniki np dla zestawu danych 150 punktow.

Właśnie mnie się nie wydawał od początku słuszny, ale pomyślałem jest poprawny, a ja go zwyczajnie nie rozumiem.

0

Och, fałszywy alarm, ja po prostu jestem głupi.

for (size_t b = 0, cut = b; b < ang.size(); b++) {

vs.

for (size_t b = 0; b < ang.size(); b++) {
    int cut = b;

Jeśli jeszcze raz będę chciał oszczędzić linijkę kodu kosztem użycia jakiejś konstrukcji w sposób zupełnie do tego nie przeznaczony, niech mnie ktoś doprowadzi do normalności ;)

(oczywiście w zamierzeniu to miało być w pętli poniżej, ale może bym to zauważył bez 20 minut debuggowania czegoś zupełnie innego gdybym to rozpisał...)

for (size_t c = b + 1, cut = b; c < ang.size(); c++) {

Dla tego przykładu wychodzi mi wynik poprawny teraz, powinno być już wszystko OK.

0

No teraz jest super, dzięki naprawdę bardzo mi pomogłeś! [-;

0
barillaaaa napisał(a):

No teraz jest super, dzięki naprawdę bardzo mi pomogłeś! [-;

Możesz podać ile czasu działa dla dużej (1-10tys) ilości punktów i na jakim komputerze?

0

Spróbuję coś wymyślić z O(n*n log n), ale na razie słabo mi idzie, trzeba by było zrobić warunek na to że jakieś punkty nigdy nie mogą zwiększyć pola trójkąta...

Ale jeszcze małe wytłumaczenie algorytmu dla niewiernych (Właśnie mnie się nie wydawał od początku słuszny) ;)
A poważniej to faktycznie rzuciłem jakimś (kiepsko dodatkowo napisanym) kodem (poza tym jak mogłem nazwać zmienne ang, b, c oraz best? Skąd wziąłem cut i on_plane? Nie wiem sam) nie tłumacząc go ani trochę, wypadałoby przynajmniej opisać o co chodzi.
Jakiś (narysowany na kolanie w 2 minuty, więc nie ma zbyt wielkiej jakości - ale takie mniej więcej ilustracje sobie szkicuję) obrazek:

img.png

Punkty na lewo, to punkty dla których rozważyliśmy już wszystkie trójkąty zawierające go. Nimi się już w ogóle nie przejmujemy.

Punkt przez który przechodzi pionowa prosta, to punkt z którego wychodzące trójkąty właśnie analizujemy (pts[i] w kodzie).

Przetwarzamy po kolei punkty B (na obrazku już jeden został przetworzony, i obecnie punktem B jest punkt nr. 2)

Teraz po kolei idziemy punktami na dole - jeden z punktów jest odcięciem które wyznacza półpłaszczyznę za którą nie ma co szukać lepszych punktów (szukany w jednym kroku tylko trójkątów ABC przypominam, czyli muszą zawierać punkt B - a jeśli jakiś punkt C leży za tą półpłaszczyzną, to znaczy że punkt cut który znalazłby się wewnątrz trójkąta ABC. NIezbyt to prosto wytłumaczone, ale na rysunku widać, mam nadzieję).

Czyli punkty 4 i 5 odrzucamy (nie rozważamy trójkątów z nimi). Teraz patrzymy na punkt 6, leży po dobrej stronie prostej B-Cut. Czyli sprawdzamy pole trójkąta ABC (C = 6), oraz zmieniamy cut na 6 (bo teraz punkt 6 jest najbardziej wysuniętym punktem).
W wyniku tego punkt 7 za prostą i jest odrzucany, pozostaje 8 do rozważenia i koniec. Teraz to samo powtórzyć dla B = 3/4/5/6/7 i możemy przejść dalej z prostą zamiatającą (tzn. punktem A).

0

Na procesorze ktory ma 2 tryby pracy (przez wiekszosc czasu pracowal na tym wydajniejszym 2.23Ghz)

Dla 1000 punktow
real 0m4.686s
user 0m4.664s
sys 0m0.004s

Dla 3000 punktow
real 2m1.590s
user 2m1.252s
sys 0m0.060s

Dla 5000 punktow
real 9m18.482s
user 9m16.887s
sys 0m0.340s

(poprawione na prośbę autora - msm)

1 użytkowników online, w tym zalogowanych: 0, gości: 1