Jakie są proste metody adaptacyjnego próbkowania funkcji 2D?

Mam dwuwymiarową funkcję której wartości chciałbym próbkować. Ta funkcja jest bardzo droga do obliczenia i ma złożony kształt, więc muszę znaleźć sposób, aby uzyskać jak najwięcej informacji o jej kształcie, używając najmniejszej liczby punktów próbkowania.$f(x,y)$

Jakie są dobre metody, aby to zrobić?

Co mam do tej pory

• Zaczynam od istniejącego zestawu punktów, w którym już obliczyłem wartość funkcji (może to być kwadratowa sieć punktów lub coś innego).

• Następnie obliczam triangulację tych punktów według Delaunaya.

• Jeśli dwa sąsiednie punkty w triangulacji Delaunaya są wystarczająco daleko ( ), a wartość funkcji różni się w nich wystarczająco ( ), wówczas wstawiam między nimi nowy punkt. Robię to dla każdej sąsiedniej pary punktów.$> \Delta X$$> \Delta f$

Co jest nie tak z tą metodą?

Cóż, działa stosunkowo dobrze, ale na funkcjach podobnych do tej nie jest idealny, ponieważ punkty próbkowe mają tendencję do „przeskakiwania” grzbietu i nawet nie zauważenia, że ​​tam jest.

Daje takie wyniki (jeśli rozdzielczość początkowej siatki punktów jest wystarczająco przybliżona):

Powyższy wykres pokazuje punkty, w których obliczana jest wartość funkcji (faktycznie komórki Voronoi wokół nich).

Powyższy wykres pokazuje interpolację liniową wygenerowaną z tych samych punktów i porównuje ją z wbudowaną metodą próbkowania Mathematiki (dla mniej więcej tej samej rozdzielczości początkowej).

Jak to poprawić?

Myślę, że głównym problemem jest to, że moja metoda decyduje, czy dodać punkt uściślenia, czy nie na podstawie gradientu.

Lepiej byłoby wziąć pod uwagę krzywiznę lub przynajmniej drugą pochodną podczas dodawania punktów uściślenia.

Pytanie

Jaki jest bardzo prosty do wdrożenia sposób uwzględnienia drugiej pochodnej lub krzywizny, gdy położenie moich punktów w ogóle nie jest ograniczone? (Niekoniecznie mam kwadratową sieć punktów początkowych, najlepiej powinien być ogólny).

Lub jakie są inne proste sposoby optymalnego obliczenia pozycji punktów uściślenia?

Zamierzam to zaimplementować w Mathematica, ale to pytanie dotyczy głównie metody. W przypadku bitu „łatwego do wdrożenia” liczy się jednak to, że używam Mathematiki (tj. Było to do tej pory łatwe, ponieważ ma pakiet do triangulacji Delaunaya)

Do jakiego problemu praktycznego się odnoszę

Obliczam diagram fazowy. Ma złożony kształt. W jednym regionie jego wartość wynosi 0, w innym regionie między 0 a 1. Nastąpił gwałtowny skok między dwoma regionami (jest nieciągły). W regionie, w którym funkcja jest większa od zera, występują zarówno pewne płynne zmiany, jak i kilka nieciągłości.

Wartość funkcji jest obliczana na podstawie symulacji Monte Carlo, więc czasami można się spodziewać niepoprawnej wartości funkcji lub szumu (jest to bardzo rzadkie, ale dzieje się tak w przypadku dużej liczby punktów, np. Gdy stan ustalony nie jest osiągnięty z powodu jakiś czynnik losowy)

_{Ja zapytałem to już na Mathematica.SE ale nie mogę się połączyć z nim dlatego, że wciąż jest w prywatnej beta. To pytanie dotyczy metody, a nie implementacji.}

Odpowiedz na @suki

Czy sugerujesz taki podział, np. Umieszczenie nowego punktu pośrodku trójkątów?

Obawiam się tutaj, że wydaje się, że wymaga specjalnego traktowania na krawędziach regionu, w przeciwnym razie da bardzo długie i bardzo cienkie trójkąty, jak pokazano powyżej. Poprawiłeś to?

AKTUALIZACJA

Problem, który pojawia się zarówno w metodzie, którą opisuję, jak i w sugestii @ suki, aby umieścić podział na podstawie trójkątów i umieścić punkty podziału wewnątrz trójkąta, polega na tym, że gdy występują nieciągłości (jak w moim problemie), ponowne obliczenie triangulacji Delaunaya po kroku może spowodować zmianę trójkątów i być może pojawienie się dużych trójkątów, które mają różne wartości funkcji w trzech wierzchołkach.

Oto dwa przykłady:

Pierwszy pokazuje wynik końcowy przy próbkowaniu wokół prostej nieciągłości. Drugi pokazuje rozkład punktu próbkowania dla podobnego przypadku.

Jakie są proste sposoby, aby tego uniknąć? Obecnie po prostu dzielę te egdy, które znikają po ponownej aranżacji, ale wydaje mi się, że to hack i należy to zrobić ostrożnie, ponieważ w przypadku siatek symetrycznych (jak kwadratowa siatka) istnieje kilka ważnych triangulacji Delaunaya, stąd krawędzie mogą się zmienić losowo po retriangulacji.

Jakiś czas temu pracowałem nad podobnym problemem.

Myślę, że główna różnica między naszymi implementacjami polega na tym, że decydowałem, gdzie dodać punkty na podstawie trójkątów, a nie krawędzi. Wybieram również nowe punkty wewnątrz trójkątów zamiast na krawędziach.

Mam wrażenie, że dodanie punktów w trójkątach poprawiłoby jego wydajność, powodując niewielki wzrost średniej odległości od starych punktów do nowych.

W każdym razie inną fajną rzeczą w używaniu trójkątów zamiast krawędzi jest to, że daje oszacowanie wektora gradientu, zamiast nachylenia wzdłuż tej konkretnej krawędzi.

W moim kodzie matlab użyłem klasy podstawowej do obsługi większości maszyn, za pomocą kilku abstrakcyjnych metod:

• weight(self) zdecydować o priorytecie, dla którego trójkąty będą dalej dzielone.
• choosePoints(self,npoints = "auto") decydować o nowych punktach do oceny na podstawie masy każdego trójkąta.

Uważam, że ta konfiguracja jest bardzo elastyczna:

• ustawienie funkcji podklasy weight()na obszar trójkąta daje stałą gęstość siatki.
• ustawienie weight()do obliczania średniej wartości funkcji razy powierzchnia trójkąta daje rodzaj quasi-losowego próbkowania prawdopodobieństwa.
• użycie funkcji var(triangle.zs)może zrobić dla funkcji, które mają wyjście binarne, to, co czuję, jest uogólnieniem wyszukiwania dzielenia na więcej niż 1 wymiar.
• używanie area + var(triangle.zs)było dość skuteczne we wprowadzaniu stałej gęstości wszędzie i zwiększonej gęstości na dowolnym zboczu (prawie to, co masz teraz).

Użyłem wariancji wartości z do przybliżenia znaczenia efektów pierwszego rzędu (nachylenie), ponieważ wariancja nigdy nie osiągnie nieskończoności, jak to może zrobić nachylenie.

W ostatnim przykładzie gęstość tła była dobra, ponieważ szukałem nieciągłych plam o wysokiej wartości w przestrzeni o niskiej wartości. Powoli wypełniałby całą siatkę, a kiedy znalazłby kroplę, koncentrowałby się na podążaniu za krawędzią kropelki dookoła ze względu na dużą wagę, którą przykładam do gradientu (i że wypełniał tylko górne ntrójkąty na każdej iteracji). Na koniec mogłem wiedzieć, że nie ma (odpowiednio ukształtowanych) plam (lub dziur w moich plamach) o rozmiarze większym niż wynikowa gęstość siatki tła.

Tak jak ty, otrzymałem kilka złych punktów w moich wynikach, nie były one dla mnie problemem, ponieważ błąd był taki, że jeśli ponownie wybierzesz punkty w pobliżu, prawdopodobnie dadzą poprawną odpowiedź. Skończyło się na tym, że wokół moich złych punktów pojawiły się plamki o zwiększonej gęstości siatki.

Cokolwiek robisz, zawsze zalecam, aby ciężary były powiązane z rozmiarem trójkąta, aby wszystkie pozostałe elementy były równe, najpierw rozbijane były duże trójkąty.

Być może rozwiązaniem jest pójście o krok dalej i zamiast oceny trójkątów na podstawie zawartości tej trójkątnej komórki, oceń na podstawie tego jednego i wszystkich trzech sąsiednich trójkątów.

Będzie zawierał wystarczającą ilość informacji, aby uzyskać oszacowanie pełnej macierzy Hesji. możesz to zrobić, dopasowując co najmniej kwadraty do z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2wszystkich wierzchołków trójkątów będących przedmiotem zainteresowania (najpierw wyśrodkuj układ współrzędnych na trójkącie).

Nie użyłbym bezpośrednio gradientu ani Hesji (tych stałych) bezpośrednio, ponieważ przejdą w nieskończoność z nieciągłością.

Być może błąd sumy kwadratów wartości z w stosunku do płaskiego aproksymacji tych punktów byłby użyteczną miarą tego, jak interesujące są efekty drugiego rzędu.

Zaktualizowano:

Wydaje mi się to rozsądne.

Tak naprawdę nigdy nie udało mi się uzyskać specjalnej osłony krawędzi. Trochę mnie to niepokoiło, ale do tego, co robiłem, wystarczyło zacząć od wielu punktów wokół krawędzi.

bardziej eleganckie byłoby połączenie naszych dwóch podejść, ważących krawędzi i trójkątów. Następnie, jeśli krawędź jest zbyt długa, przeciąć ją na pół ... Podoba mi się sposób, w jaki koncepcja uogólnia się na wyższe wymiary (ale liczby szybko się zwiększają) ...

Ale ponieważ nie oczekujesz, że główna bryła siatki będzie miała trójkąty o wysokim współczynniku kształtu, możesz więc użyć funkcji takiej jak funkcja swobodnej granicy Matlaba, aby znaleźć granicę, a następnie uruchomić ten sam algorytm w jednym mniejszym wymiarze na granicy. Jeśli zrobisz to dobrze, na przykład na kostce, możesz uzyskać taką samą gęstość siatki na krawędziach, ścianach i wewnątrz kostki. Ciekawy.

Jedną z rzeczy, dla których nigdy nie znalazłem dobrego rozwiązania, był fakt, że moja wersja nigdy nie eksplorowałaby poza wypukłym kadłubem początkowego zestawu punktów.

Myślę, że głównym problemem w twojej heurystyce jest to, że rozważasz gradient tylko w jednym wymiarze, a zatem w regionach, w których dfdx jest mały, ale dfdy jest duży (jak to dzieje się w środku twojego przykładu), ominiesz punkty, patrząc w „złym” wymiarze.

Jednym szybkim rozwiązaniem byłoby rozważenie zbiorów czterech punktów, biorąc ich środek ciężkości i przybliżając | dfdx | + | dfdy | używając tych czterech punktów. Inną alternatywą jest wzięcie trzech punktów (tj. Trójkąta) i wzięcie nad nimi maksymalnego gradientu powierzchni.