Autorzy Dmitry Anisimov, David Bommes, Kai Hormann, and Pierre Alliez
Pakiet 2D Generalized Barycentric Coordinates oferuje wydajną i solidną implementację dwuwymiarowych zamkniętych uogólnionych współrzędnych barycentrycznych zdefiniowanych dla prostych dwuwymiarowych wielokątów. Jeśli wymagane są współrzędne w odniesieniu do punktów rozproszonych zamiast wielokąta, proszę odnieść się do współrzędnych naturalnego sąsiada z pakietu 2D i interpolacji funkcji powierzchniowych.
W szczególności, pakiet zawiera implementację współrzędnych Wachspressa, wartości średniej i dyskretnych współrzędnych harmonicznych oraz zapewnia kilka dodatkowych funkcji do obliczania współrzędnych barycentrycznych w odniesieniu do odcinków (współrzędne odcinka) i trójkątów (współrzędne trójkąta). Sekcja Theory of 2D Generalized Barycentric Coordinates zawiera krótkie wprowadzenie do tematu współrzędnych barycentrycznych.
Każda klasa, która oblicza współrzędne barycentryczne, jest parametryzowana przez klasę cech. Ta klasa cech określa typy i prymitywy geometryczne, które są używane w obliczeniach, i musi być modelem koncepcji BarycentricTraits_2.
Głównym punktem wejścia do komponentu jest iterator po wierzchołkach wielokąta. Wierzchołki wielokąta muszą być uporządkowane zgodnie z ruchem wskazówek zegara lub przeciwnie do ruchu wskazówek zegara i mogą być dowolnego typu. Jednak wewnętrznie klasy używają typu CGAL::Point_2, dlatego należy dostarczyć odpowiednią klasę traits, która konwertuje typ podany przez użytkownika na CGAL::Point_2. Ten sam argument odnosi się do punktów zapytania.
Współrzędne średniej wartości są najbardziej ogólnymi współrzędnymi w tym pakiecie, ponieważ pozwalają na podanie dowolnego prostego wielokąta jako danych wejściowych. Współrzędne Wachspressa i dyskretne współrzędne harmoniczne są z definicji ograniczone do ściśle wypukłych wielokątów. Współrzędne segmentu przyjmują jako dane wejściowe dowolny niezdegenerowany segment, a współrzędne trójkąta pozwalają na dowolny niezdegenerowany trójkąt.
Współrzędne segmentu i trójkąta mogą być obliczane przy użyciu funkcji globalnej lub przez utworzenie odpowiedniej klasy. Wszystkie inne uogólnione współrzędne mogą być obliczone przez utworzenie instancji klasy CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2 sparametryzowanej przez odpowiedni typ współrzędnych, który musi być modelem koncepcji BarycentricCoordinates_2.
Dowolny punkt na płaszczyźnie może być potraktowany jako punkt zapytania. Nie zalecamy jednak używania współrzędnych Wachspressa i dyskretnych współrzędnych harmonicznych z punktami zapytania znajdującymi się poza zamknięciem wielokąta, ponieważ w niektórych z tych punktów współrzędne te nie są dobrze zdefiniowane, jak wyjaśniono w rozdziale Degeneracje i przypadki specjalne.
Po zainicjowaniu dla pewnego wielokąta, współrzędne mogą być obliczane wielokrotnie dla różnych punktów zapytania w odniesieniu do wszystkich wierzchołków podanego wielokąta. Szczegółowy interfejs znajduje się w Reference Manual.
Wynikiem obliczeń jest zbiór wartości współrzędnych w bieżącym punkcie zapytania w odniesieniu do wszystkich wierzchołków wielokąta. Dane wyjściowe mogą być przechowywane w dowolnym kontenerze zapewniającym odpowiedni iterator wyjściowy. Dodatkowo wszystkie klasy zwracają wskaźnik do ostatniego przechowywanego elementu oraz status obliczeń (Boolean true lub false).
Współrzędne segmentu
Jest to prosty przykład pokazujący użycie funkcji globalnej CGAL::Barycentric_coordinates::compute_segment_coordinates_2(). Obliczamy współrzędne w trzech zielonych punktach wzdłuż odcinka oraz w dwóch niebieskich punktach poza tym odcinkiem, ale wzdłuż jego linii pomocniczej. Używamy dokładnego jądra i zwracamy współrzędne jako tablicę dwóch wartości. Ponownie, symetria punktów zapytania pomaga nam rozpoznać błędy, które mogły wystąpić podczas obliczeń.
Rysunek 96.1 Wzorzec punktów przykładu.
Plik Barycentric_coordinates_2/Segment_coordinates_example.cpp
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Segment_coordinates_2.h>
// Alias przestrzeni nazw.
namespace BC = CGAL::Barycentric_coordinates;
// Kilka wygodnych typedefów.
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef std::array<Scalar,2> Pair;
using std::cout; using std::endl; using std::string;
int main()
{
// Skonstruuj odcinek.
const Point first_vertex(0, Scalar(2)/Scalar(5));
const Point second_vertex(2, Scalar(2)/Scalar(5));
// Instantiate three interior and two exterior query points.
const Point query_points = { Point(Scalar(2) /Scalar(5), Scalar(2)/Scalar(5)), // wewnętrzne punkty zapytania
Point(1 , Scalar(2)/Scalar(5)),
Point(Scalar(8) /Scalar(5), Scalar(2)/Scalar(5)),
Point(Scalar(-1)/Scalar(5), Scalar(2)/Scalar(5)), // zewnętrzne punkty zapytania
Point(Scalar(11)/Scalar(5), Scalar(2)/Scalar(5))
};
// Oblicz współrzędne odcinka dla wszystkich zdefiniowanych punktów.
// Używamy funkcji globalnej i zwracamy współrzędne odcinka przechowywane w tablicy typu std::array<FT,2>.
cout << endl << „Obliczone współrzędne segmentu: ” << endl << endl;
for(int i = 0; i < 5; ++i) {
const Pair pair = BC::compute_segment_coordinates_2(first_vertex, second_vertex, query_points, Kernel());
// Wyprowadź obie współrzędne dla każdego punktu.
cout << „Para współrzędnych # ” << i + 1 << ” = (” << para << „, ” << para << „);” << endl;
}
cout << endl;
return EXIT_SUCCESS;
}
Współrzędne trójkąta
W tym przykładzie pokazujemy, jak używać klasy CGAL::Barycentric_coordinates::Triangle_coordinates_2 z jądrem Simple_cartesian dla typu double. Obliczamy współrzędne dla trzech zbiorów punktów: wewnętrznego (zielony), brzegowego (czerwony) i zewnętrznego (niebieski). Zauważ, że niektóre z wartości współrzędnych dla punktów zewnętrznych są ujemne. Używamy standardowego kontenera typu std::vector i std::insert_iterator, aby uzyskać dostęp do wynikowych wartości współrzędnych i je przechowywać.
Rysunek 96.2 Przykładowy układ punktów.
Plik Barycentric_coordinates_2/Triangle_coordinates_example.cpp
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Barycentric_coordinates_2/Triangle_coordinates_2.h>
// Kilka wygodnych typedefów.
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef std::vector<Scalar> Scalar_vector;
typedef CGAL::Barycentric_coordinates::Triangle_coordinates_2<Kernel> Triangle_coordinates;
using std::cout; using std::endl; using std::string;
int main()
{
// Skonstruuj trójkąt.
const Point first_vertex(0.0f, 0.0f);
const Point second_vertex(2.0f, 0.5f);
const Point third_vertex(1.0f, 2.0f);
// Tworzymy wektor std::vector do przechowywania współrzędnych.
Scalar_vector coordinates;
// Instantiate the class Triangle_coordinates_2 for the triangle defined above.
Triangle_coordinates triangle_coordinates(first_vertex, second_vertex, third_vertex);
// Wypisz kilka informacji o trójkącie i współrzędnych.
triangle_coordinates.print_information();
// Zinstytucjonalizuj niektóre wewnętrzne, brzegowe i zewnętrzne punkty zapytania, dla których obliczamy współrzędne.
const int number_of_query_points = 18;
const Point query_points = { Point(0.5f , 0.5f ), Point(1.0f, 0.5f ), Point(1.0f , 0.75f), Point(1.0f , 1.0f), // interior query points
Point(1.0f , 1.25f), Point(1.0f , 1.5f ), Point(0.75f , 1.0f ), Point(1.25f , 1.0f), Point(1.5f , 0.75f),
Point(1.0f , 0.25f), Point(0.5f, 1.0f ), Point(1.5f , 1.25f), Point(1.0f , 2.0f), Point(2.0f, 0.5f ), // boundary query points
Point(0.25f, 1.0f ), Point(0.5f, 1.75f), Point(1.5f , 1.75f), Point(1.75f, 1.5f) // zewnętrzne punkty zapytania
};
// Zarezerwuj pamięć do przechowywania współrzędnych trójkątów dla 18 punktów zapytania.
coordinates.reserve(number_of_query_points * 3);
// Oblicz współrzędne trójkąta dla tych punktów.
cout << endl << „Obliczone współrzędne trójkąta: ” << endl << endl;
for(int i = 0; i < number_of_query_points; ++i) {
triangle_coordinates(query_points, std::inserter(coordinates, coordinates.end()));
// Wyprowadź współrzędne dla każdego punktu.
cout << „Punkt ” << i + 1 << „: „;
for(int j = 0; j < 3; ++j)
cout << „współrzędna ” << j + 1 << ” = ” << współrzędne << „; „;
cout << endl << endl;
}
return EXIT_SUCCESS;
}
Współrzędne Wachspressa
W poniższym przykładzie tworzymy 1000 losowych punktów, następnie bierzemy wypukły kadłub tego zbioru punktów jako nasz wielokąt i obliczamy współrzędne Wachspressa we wszystkich zdefiniowanych punktach. Jako klasy cech używamy jądra Simple_cartesian z typem double, a uzyskane wartości współrzędnych przechowujemy w kontenerze typu std::vector. Wyjściowym iteratorem jest std::back_insert_iterator.
Plik Barycentric_coordinates_2/Wachspress_coordinates_example.cpp
#include <CGAL/convex_hull_2.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/point_generators_2.h>
#include <CGAL/Barycentric_coordinates_2/Wachspress_2.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
// Kilka wygodnych typedefów.
typedef CGAL::Simple_cartesian<double> Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point> Point_vector;
typedef CGAL::Creator_uniform_2<double, Point> Creator;
typedef CGAL::Barycentric_coordinates::Wachspress_2<Kernel> Wachspress;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Wachspress, Kernel> Wachspress_coordinates;
using std::cout; using std::endl; using std::string;
int main()
{
// Wybierz, ile losowych punktów chcemy wygenerować.
const int number_of_points = 1000;
// Utwórz wektory do przechowywania wygenerowanych punktów i wierzchołków wielokąta wypukłego.
Point_vector points, vertices;
// Wygeneruj zbiór losowych punktów.
CGAL::Random_points_in_square_2<Point,Creator> point_generator(1.0);
std::copy_n(point_generator, number_of_points, std::back_inserter(points));
// Znajdź wypukły kadłub wygenerowanego zbioru punktów.
// Ten wypukły kadłub daje wierzchołki wielokąta wypukłego, który zawiera wszystkie wygenerowane punkty.
CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(vertices));
const size_t number_of_vertices = vertices.size();
// Instantiate the class with Wachspress coordinates for the convex polygon defined above.
Wachspress_coordinates wachspress_coordinates(vertices.begin(), vertices.end());
// Wypisz kilka informacji o wielokącie i współrzędnych.
wachspress_coordinates.print_information();
// Oblicz współrzędne Wachspressa dla wszystkich losowo zdefiniowanych punktów.
cout << endl << „Obliczone współrzędne Wachspressa: ” << endl << endl;
for(int i = 0; i < number_of_points; ++i) {
// Compute coordinates.
Scalar_vector coordinates;
coordinates.reserve(number_of_vertices);
wachspress_coordinates(points, std::back_inserter(coordinates));
// Wyprowadź obliczone współrzędne.
cout << „Punkt ” << i + 1 << „: ” << endl;
for(int j = 0; j < int(number_of_vertices); ++j) cout << „Współrzędne ” << j + 1 << ” = ” << współrzędne << „; ” << endl;
cout << endl;
}
return EXIT_SUCCESS;
}
Dyskretne współrzędne harmoniczne
W tym przykładzie obliczamy dyskretne współrzędne harmoniczne dla zbioru punktów zielonych (wewnętrznych), czerwonych (brzegowych) i niebieskich (zewnętrznych) względem kwadratu jednostkowego. Pokazujemy również jak można określić położenie punktu zapytania za pomocą dodatkowych parametrów funkcji. Zastosowane jądro jest dokładne, a jako kontener wyjściowy wykorzystujemy kontener typu std::vector. Ponieważ wszystkie punkty są symetryczne, łatwo jest debugować poprawność otrzymanych wartości współrzędnych. Wyjściowy iterator ma postać std::back_insert_iterator.
Rysunek 96.3 Przykładowy układ punktów.
Plik Barycentric_coordinates_2/Discrete_harmonic_coordinates_example.cpp
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Discrete_harmonic_2.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
// Kilka wygodnych typedefów.
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point> Point_vector;
typedef std::back_insert_iterator<Scalar_vector> Vector_insert_iterator;
typedef boost::optional<Vector_insert_iterator> Output_type;
typedef CGAL::Barycentric_coordinates::Discrete_harmonic_2<Kernel> Discrete_harmonic;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Discrete_harmonic, Kernel> Discrete_harmonic_coordinates;
using std::cout; using std::endl; using std::string;
int main()
{
// Konstruujemy kwadrat jednostkowy.
const int number_of_vertices = 4;
Point_vector vertices(number_of_vertices);
vertices = Point(0, 0); vertices = Point(1, 0); vertices = Point(1, 1); vertices = Point(0, 1);
// Utwórz wektor std::vector do przechowywania współrzędnych.
Scalar_vector coordinates;
// Instantiate the class with discrete harmonic coordinates for the unit square defined above.
Discrete_harmonic_coordinates discrete_harmonic_coordinates(vertices.begin(), vertices.end());
// Wypisz kilka informacji o wielokącie i współrzędnych.
discrete_harmonic_coordinates.print_information();
// Instantiate the center point of the unit square.
const Point center(Scalar(1)/Scalar(2), Scalar(1)/Scalar(2));
// Oblicz współrzędne harmoniczne dyskretne dla punktu środkowego.
// Użyj parametru query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE.
Output_type result = discrete_harmonic_coordinates(center, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE);
// Instantiate other 4 interior points.
const int number_of_interior_points = 4;
const Point interior_points = { Point(Scalar(1)/Scalar(5), Scalar(1)/Scalar(5)) ,
Point(Scalar(4)/Scalar(5), Scalar(1)/Scalar(5)) ,
Point(Scalar(4)/Scalar(5), Scalar(4)/Scalar(5)) ,
Point(Scalar(1)/Scalar(5), Scalar(4)/Scalar(5)) };
// Oblicz dyskretne współrzędne harmoniczne dla tych punktów i przechowuj je w tym samym wektorze „współrzędnych”, co poprzednio.
for(int i = 0; i < number_of_interior_points; ++i)
result = discrete_harmonic_coordinates(interior_points, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE);
// Instantiate 2 boundary points on the second and last edges.
const Point second_edge(1, Scalar(4)/Scalar(5));
const Point last_edge(0, Scalar(4)/Scalar(5));
// Oblicz dyskretne współrzędne harmoniczne dla tych 2 punktów.
// Użyj parametru query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDARY.
result = discrete_harmonic_coordinates(second_edge , std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDARY);
result = discrete_harmonic_coordinates(last_edge , std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_BOUNDARY);
// Instantiate 2 inne punkty graniczne na pierwszej i trzeciej krawędzi.
const Point first_edge(Scalar(1)/Scalar(2), 0);
const Point third_edge(Scalar(1)/Scalar(2), 1);
// Oblicz dyskretne współrzędne harmoniczne używając indeksu odpowiedniej krawędzi.
// Nie zapomnij, że liczenie indeksów zaczyna się od zera.
result = discrete_harmonic_coordinates.compute_on_edge(first_edge, 0, std::back_inserter(coordinates));
result = discrete_harmonic_coordinates.compute_on_edge(third_edge, 2, std::back_inserter(coordinates));
// Oblicz współrzędne harmoniczne dyskretne dla punktów w pierwszym i trzecim wierzchołku kwadratu jednostkowego.
result = discrete_harmonic_coordinates.compute_on_vertex(0, std::back_inserter(coordinates));
result = discrete_harmonic_coordinates.compute_on_vertex(2, std::back_inserter(coordinates));
// Instantiate points at the second and fourth vertex of the unit square.
const Point second_vertex(1, 0);
const Point fourth_vertex(0, 1);
// Oblicz dyskretne współrzędne harmoniczne dla tych punktów.
// Użyj parametru query_point_location = CGAL::Barycentric_coordinates::ON_VERTEX.
result = discrete_harmonic_coordinates(second_vertex, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_VERTEX);
result = discrete_harmonic_coordinates(fourth_vertex, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_VERTEX);
// Instantiate 2 points outside the unit square – one from the left and one from the right.
const Point left_most(Scalar(-1)/Scalar(2), Scalar(1)/Scalar(2));
const Point right_most(Scalar(3)/Scalar(2), Scalar(1)/Scalar(2));
// Oblicz dyskretne współrzędne harmoniczne dla tych 2 punktów.
// Użyj parametru query_point_location = CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE.
result = discrete_harmonic_coordinates(left_most , std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE);
result = discrete_harmonic_coordinates(right_most, std::back_inserter(coordinates), CGAL::Barycentric_coordinates::ON_UNBOUNDED_SIDE);
// Wyprowadź obliczone wartości współrzędnych.
cout << endl << „Dokładne dyskretne współrzędne harmoniczne dla wszystkich zdefiniowanych punktów: ” << endl << endl;
const size_t number_of_query_points = coordinates.size();
for(int index = 0; index < int(number_of_query_points); ++index) {
cout << „Coordinate ” << index % number_of_vertices + 1 << ” = ” << coordinates << ” „;
if((index + 1) % number_of_vertices == 0) cout << endl;
if((index + 13) % (4 * number_of_vertices) == 0) cout << endl;
}
// Zwróć status ostatniego obliczenia.
const string status = (result ? „SUCCESS.” : „FAILURE.”);
cout << endl << „Status ostatniego obliczenia: ” << status << endl << endl;
return EXIT_SUCCESS;
}
Współrzędne wartości średniej
To jest przykład, który pokazuje, jak obliczyć współrzędne wartości średniej dla zbioru zielonych punktów w wielokącie w kształcie gwiazdy. Zauważmy, że ten typ współrzędnych jest dobrze zdefiniowany dla takiego wklęsłego wielokąta, podczas gdy współrzędne Wachspressa i dyskretne współrzędne harmoniczne nie są. Jednakże, może on dawać ujemne wartości współrzędnych dla punktów znajdujących się poza jądrem wielokąta (pokazanych na czerwono). Używamy nieścisłego typu danych, kontenera wyjściowego typu std::vector oraz iteratora wyjściowego typu std::back_insert_iterator do obliczania, dostępu i przechowywania wynikowych wartości współrzędnych. Pokazujemy również, jak wybrać różne algorytmy do obliczania uogólnionych współrzędnych barycentrycznych (jeden jest bardziej precyzyjny, a drugi szybszy).
Rysunek 96.4 Wzór punktu przykładu.
Plik Barycentric_coordinates_2/Mean_value_coordinates_example.cpp
#include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
// Kilka wygodnych typedefów.
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point> Point_vector;
typedef std::back_insert_iterator<Scalar_vector> Vector_insert_iterator;
typedef boost::optional<Vector_insert_iterator> Output_type;
typedef CGAL::Barycentric_coordinates::Mean_value_2<Kernel> Mean_value;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Mean_value, Kernel> Mean_value_coordinates;
using std::cout; using std::endl; using std::string;
int main()
{
// Skonstruuj wielokąt w kształcie gwiazdy.
const int number_of_vertices = 10;
Point_vector vertices(number_of_vertices);
vertices = Point(0.0, 0.0); vertices = Point(0.1, -0.8); vertices = Point(0.3, 0.0); vertices = Point(0.6, -0.5); vertices = Point(0.6 , 0.1);
vertices = Point(1.1, 0.6); vertices = Point(0.3, 0.2); vertices = Point(0.1, 0.8); vertices = Point(0.1, 0.2); vertices = Point(-0.7, 0.0);
// Create an std::vector to store coordinates.
Scalar_vector coordinates;
// Instantiate the class with mean value coordinates for the polygon defined above.
Mean_value_coordinates mean_value_coordinates(vertices.begin(), vertices.end());
// Wypisz kilka informacji o wielokącie i współrzędnych.
mean_value_coordinates.print_information();
// Zinstytucjonalizuj niektóre punkty wewnętrzne w wielokącie.
const int number_of_interior_points = 8;
const Point interior_points = { Point(0.12, -0.45), Point(0.55, -0.3), Point(0.9 , 0.45),
Point(0.15, 0.35), Point(-0.4, 0.04), Point(0.11, 0.11),
Point(0.28, 0.12), // jedyny punkt w jądrze wielokąta w kształcie gwiazdy
Point(0.55, 0.11)
};
// Oblicz współrzędne wartości średniej dla wszystkich zdefiniowanych punktów wewnętrznych.
// Przyspieszamy obliczenia stosując algorytm O(n) wywoływany z parametrem CGAL::Barycentric_coordinates::FAST.
// Domyślnym jest CGAL::Barycentric_coordinates::PRECISE.
const CGAL::Barycentric_coordinates::Type_of_algorithm type_of_algorithm = CGAL::Barycentric_coordinates::FAST;
// Przyspieszamy również obliczenia poprzez użycie parametru query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE.
const CGAL::Barycentric_coordinates::Query_point_location query_point_location = CGAL::Barycentric_coordinates::ON_BOUNDED_SIDE;
for(int i = 0; i < number_of_interior_points; ++i) {
const Output_type result = mean_value_coordinates(interior_points, std::back_inserter(coordinates), query_point_location, type_of_algorithm);
// Wyprowadź współrzędne dla każdego punktu.
const string status = (result ? „SUCCESS.” : „FAILURE.”);
cout << endl << „Dla punktu ” << i + 1 << ” status obliczenia: ” << status << endl;
for(int j = 0; j < liczba_wierzchołków; ++j)
cout << „Współrzędne ” << j + 1 << ” = ” << współrzędne << endl;
}
// Jeśli potrzebujemy tylko nienormalizowanych wag dla jakiegoś punktu (weźmy ostatni), możemy je obliczyć w następujący sposób.
// Instantiate an std::vector to store weights.
Scalar_vector weights;
// Oblicz wagę wartości średniej.
const int last_point_index = number_of_interior_points – 1;
const Output_type result = mean_value_coordinates.compute_weights(interior_points, std::back_inserter(weights));
// Oblicz ich sumę.
Scalar mv_denominator = Scalar(0);
for(int j = 0; j < number_of_vertices; ++j) mv_denominator += weights;
// Odwróć tę sumę.
const Scalar mv_inverted_denominator = Scalar(1) / mv_denominator;
// Wyprowadzenie wag wartości średniej.
const string status = (result ? „SUCCESS.” : „FAILURE.”);
cout << endl << „Status obliczeń wag dla punktu ” << ostatni_punkt_index + 1 << „: ” << status << endl;
for(int j = 0; j < liczba_wag; ++j)
cout << „Waga ” << j + 1 << ” = ” << wagi << endl;
// Teraz, jeśli znormalizujemy wagi, odzyskamy wartości współrzędnych wartości średniej dla ostatniego punktu obliczonego wcześniej.
cout << endl << „Po normalizacji, dla punktu ” << last_point_index + 1 << ” współrzędne wartości średniej wynoszą ” << endl;
for(int j = 0; j < liczba_wierzchołków; ++j)
cout << „Współrzędne ” << j + 1 << ” = ” << wagi * mv_odwrócony_denominator << endl;
cout << endl;
return EXIT_SUCCESS;
}
Height Interpolation for Terrain Modeling
To jest zaawansowany przykład, który pokazuje, jak używać uogólnionych współrzędnych barycentrycznych do interpolacji wysokości z zastosowaniem do modelowania terenu. Pokazuje on również, jak używać nie-domyślnej klasy cech z naszym pakietem zamiast klasy Kernel cech. Załóżmy, że znamy granicę trójwymiarowego kawałka terenu, który może być reprezentowany jako wielokąt z kilkoma trójwymiarowymi wierzchołkami, gdzie trzeci wymiar daje odpowiednią wysokość. Zadanie polega na propagacji wysokości ze znanych nam przykładowych punktów na granicy do wnętrza wielokąta. Daje to przybliżone oszacowanie powierzchni terenu w tym regionie.
Rysunek 96.5 Wielokąt 2D o 50 wierzchołkach reprezentujący fragment terenu z częściami wypukłymi i wklęsłymi. Wysokość nie jest pokazana.
W tym przykładzie rzutujemy trójwymiarowy wielokąt ortogonalnie na płaszczyznę dwuwymiarową przy użyciu klasy CGAL::Projection_traits_xy_3, triangulujemy jego wnętrze przy użyciu klasy CGAL::Delaunay_mesher_2 i obliczamy współrzędne wartości średniej dla wszystkich uzyskanych punktów względem wszystkich wierzchołków wielokąta. Na koniec interpolujemy dane wysokościowe od granicy wielokąta do jego wnętrza przy użyciu obliczonych współrzędnych oraz globalnej funkcji interpolacji z pakietu 2D and Surface Function Interpolation.
Plik Barycentric_coordinates_2/Terrain_height_modeling.cpp
#include <CGAL/Delaunay_mesher_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/Projection_traits_xy_3.h>
#include <CGAL/interpolation_functions.h>
#include <CGAL/Delaunay_mesh_face_base_2.h>
#include <CGAL/Delaunay_mesh_size_criteria_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Barycentric_coordinates_2/Mean_value_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Generalized_barycentric_coordinates_2.h>
// Some convenient typedefs.
// General.
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Projection_traits_xy_3<Kernel> Projection;
typedef Projection::FT Scalar;
typedef Projection::Point_2 Point;
typedef std::vector<Scalar> Scalar_vector;
typedef std::vector<Point> Point_vector;
// Współrzędne związane.
typedef CGAL::Barycentric_coordinates::Mean_value_2<Projection> Mean_value;
typedef CGAL::Barycentric_coordinates::Generalized_barycentric_coordinates_2<Mean_value, Projection> Mean_value_coordinates;
/// Związane z triangulacją.
typedef CGAL::Delaunay_mesh_face_base_2<Projection> Face_base;
typedef CGAL::Triangulation_vertex_base_2<Projection> Vertex_base;
typedef CGAL::Triangulation_data_structure_2<Vertex_base, Face_base> TDS;
typedef CGAL::Constrained_Delaunay_triangulation_2<Projection, TDS> CDT;
typedef CGAL::Delaunay_mesh_size_criteria_2<CDT> Criteria;
typedef CGAL::Delaunay_mesher_2<CDT, Criteria> Mesher;
typedef CDT::Finite_vertices_iterator Vertex_iterator;
typedef CDT::Vertex_handle Vertex_handle;
// Związane z interpolacją.
typedef CGAL::Interpolation_traits_2<Projection> Interpolation_traits;
typedef CGAL::Data_access< std::map<Point, Scalar, Projection::Less_xy_2 > > Value_access;
// STD.
using std::cout; using std::endl; using std::string;
int main()
{
// Konstrukcja wielokąta ograniczającego fragment trójwymiarowego terenu.
// Zauważ, że współrzędna z każdego wierzchołka reprezentuje funkcję wysokości.
// Rzutowanie w 2D jest wykonywane automatycznie przez klasę cech Projection.
const int number_of_vertices = 50;
Point_vector vertices(number_of_vertices);
vertices = Point(0.03, 0.05, 0.000); vertices = Point(0.07, 0.04, 10.00); vertices = Point(0.10, 0.04, 20.00);
vertices = Point(0.14, 0.04, 30.00); vertices = Point(0.17, 0.07, 40.00); vertices = Point(0.19, 0.09, 50.00);
vertices = Point(0.22, 0.11, 60.00); vertices = Point(0.25, 0.11, 70.00); vertices = Point(0.27, 0.10, 80.00);
vertices = Point(0.30, 0.07, 90.00); vertices = Point(0.31, 0.04, 100.0); vertices = Point(0.34, 0.03, 110.0);
vertices = Point(0.37, 0.02, 120.0); vertices = Point(0.40, 0.03, 130.0); vertices = Point(0.42, 0.04, 140.0);
vertices = Point(0.44, 0.07, 150.0); vertices = Point(0.45, 0.10, 160.0); vertices = Point(0.46, 0.13, 170.0);
vertices = Point(0.46, 0.19, 180.0); vertices = Point(0.47, 0.26, 190.0); vertices = Point(0.47, 0.31, 200.0);
vertices = Point(0.47, 0.35, 210.0); vertices = Point(0.45, 0.37, 220.0); vertices = Point(0.41, 0.38, 230.0);
vertices = Point(0.38, 0.37, 240.0); vertices = Point(0.35, 0.36, 250.0); vertices = Point(0.32, 0.35, 260.0);
vertices = Point(0.30, 0.37, 270.0); vertices = Point(0.28, 0.39, 280.0); vertices = Point(0.25, 0.40, 290.0);
vertices = Point(0.23, 0.39, 300.0); vertices = Point(0.21, 0.37, 310.0); vertices = Point(0.21, 0.34, 320.0);
vertices = Point(0.23, 0.32, 330.0); vertices = Point(0.24, 0.29, 340.0); vertices = Point(0.27, 0.24, 350.0);
vertices = Point(0.29, 0.21, 360.0); vertices = Point(0.29, 0.18, 370.0); vertices = Point(0.26, 0.16, 380.0);
vertices = Point(0.24, 0.17, 390.0); vertices = Point(0.23, 0.19, 400.0); vertices = Point(0.24, 0.22, 410.0);
vertices = Point(0.24, 0.25, 420.0); vertices = Point(0.21, 0.26, 430.0); vertices = Point(0.17, 0.26, 440.0);
vertices = Point(0.12, 0.24, 450.0); vertices = Point(0.07, 0.20, 460.0); vertices = Point(0.03, 0.15, 470.0);
vertices = Point(0.01, 0.10, 480.0); vertices = Point(0.02, 0.07, 490.0);
// Mesh this polygon.
// Utwórz ograniczoną triangulację Delaunay.
CDT cdt;
std::vector<Vertex_handle> vertex_handle(number_of_vertices);
// Wstaw wierzchołki wielokąta jako nasz początkowy zestaw punktów.
for(int i = 0; i < number_of_vertices; ++i) vertex_handle = cdt.insert(vertices);
// Wstaw ograniczenia – krawędzie wielokąta – aby zazębić tylko wnętrze wielokąta.
for(int i = 0; i < number_of_vertices; ++i) cdt.insert_constraint(vertex_handle, vertex_handle);
Mesher mesher(cdt);
// Ustaw kryteria jak meshować.
mesher.set_criteria(Criteria(0.01, 0.01));
// Mesh the polygon.
mesher.refine_mesh();
// Oblicz współrzędne wartości średniej i użyj ich do interpolacji danych od granicy wielokąta do jego wnętrza.
// Skojarz każdy punkt z odpowiadającą mu wartością funkcji i współrzędnymi.
std::map<Point, Scalar, Projection::Less_xy_2> point_function_value;
std::vector< std::pair<Point, Scalar> > point_coordinates(number_of_vertices);
for(int i = 0; i < number_of_vertices; ++i)
point_function_value.insert(std::make_pair(vertices, vertices.z()));
// Utwórz instancję klasy ze współrzędnymi wartości średniej.
Mean_value_coordinates mean_value_coordinates(vertices.begin(), vertices.end());
// Przechowuj tutaj wszystkie nowe punkty wewnętrzne z interpolowanymi danymi.
std::vector<Point> points(cdt.number_of_vertices()));
cout << endl << „Wynik interpolacji wysokości: ” << endl << endl;
// Oblicz współrzędne i interpoluj dane brzegowe do wnętrza wielokąta.
int index = 0;
for(Vertex_iterator vertex_iterator = cdt.finite_vertices_begin(); vertex_iterator != cdt.finite_vertices_end(); ++vertex_iterator) {
Scalar_vector coordinates;
const Point &point = vertex_iterator->point();
mean_value_coordinates(point, std::back_inserter(coordinates));
for(int j = 0; j < number_of_vertices; ++j)
point_coordinates = std::make_pair(vertices, coordinates);
Scalar f = CGAL::linear_interpolation(point_coordinates.begin(), point_coordinates.end(), Scalar(1), Value_access(point_function_value));
points = Point(point.x(), point.y(), f);
cout << „Interpolowana wysokość o indeksie ” << index << ” wynosi ” << f << „;” << endl;
++index;
}
cout << endl;
return EXIT_SUCCESS;
}
W wyniku otrzymujemy gładką funkcję wewnątrz wielokąta, która aproksymuje bazową powierzchnię terenu.
Rysunek 96.6 Interpolowane dane. Kolorowy pasek przedstawia odpowiednią wysokość.
Zwyrodnienia i przypadki specjalne
Współrzędne odcinka
Współrzędne odcinka mogą być obliczane dokładnie, jeżeli wybrany zostanie dokładny typ danych. Sam odcinek, względem którego obliczamy współrzędne, musi być niezdegenerowany. Jeśli oba te warunki są spełnione, to obliczenie nigdy nie kończy się niepowodzeniem. Jednakże, aby obliczyć współrzędne, użytkownik musi być pewien, że punkt zapytania znajduje się dokładnie na prostej wspierającej odcinek. Ponieważ w wielu zastosowaniach tak nie jest i punkt zapytania może leżeć bardzo blisko, ale nie dokładnie na tej prostej, klasa jest w stanie poradzić sobie również z taką sytuacją.
Rysunek 96.7 Rzut ortogonalny \(p’\) wektora \(p\) (zielony) na wektor \(q\) (czerwony).
Załóżmy, że pewien punkt zapytania nie leży dokładnie na prostej \(L\), ale jest od niej oddalony o pewną odległość \(d\), jak pokazano na powyższym rysunku. Jeśli chcemy obliczyć współrzędną barycentryczną odcinka \(b_1(v)\) w odniesieniu do wierzchołka \(v_1\), najpierw znajdujemy rzut ortogonalny \(p’\) wektora \(p\) na wektor \(q\), a następnie normalizujemy go przez długość \(q\). W ten sposób otrzymujemy współrzędną barycentryczną odcinka \(b_1(v’) = b_1(v)\), jeśli punkt \(v’) leży dokładnie na prostej.
Ostrzeżenie: nie należy nadużywać opisanej powyżej funkcji, ponieważ nie daje ona poprawnych współrzędnych barycentrycznych odcinka dla punktu \(v’), lecz dla \(v’\). Co więcej, współrzędne barycentryczne odcinka dla punktu \(v), który nie leży dokładnie na prostej \(L), nie istnieją. Jeśli jednak niezerowa odległość jest spowodowana niestabilnością numeryczną przy obliczaniu położenia punktu lub innym problemem, który powoduje, że punkt nie leży dokładnie na prostej, to ostateczne współrzędne odcinka będą, przynajmniej w przybliżeniu, poprawne.
W przypadku nieścisłych typów danych, wynikowe wartości współrzędnych są poprawne do precyzji wybranego typu.
Współrzędne trójkąta
Współrzędne te mogą być obliczone dokładnie, jeśli wybrany jest dokładny typ danych, dla dowolnego punktu zapytania na płaszczyźnie i w odniesieniu do dowolnego niezdegenerowanego trójkąta. Nie są obsługiwane żadne specjalne przypadki. Obliczenia zawsze dają poprawny wynik. Pojęcie poprawności zależy od precyzji użytego typu danych. Zauważ, że dla punktów zewnętrznych niektóre wartości współrzędnych będą ujemne.
Współrzędne Wachspressa
Współrzędne Wachspressa są dobrze zdefiniowane w zamknięciu dowolnego ściśle wypukłego wielokąta. Dlatego dla każdego punktu zapytania z zamknięcia wielokąta z dokładnym typem danych, współrzędne te są obliczane dokładnie i nie oczekuje się żadnego fałszywego wyniku. Dla niedokładnych typów danych, precyzja obliczeń wynika z zastosowanego algorytmu i wybranego typu danych. W następnym paragrafie omówimy dwa dostępne algorytmy do obliczania współrzędnych Wachspress. Jeden z nich to CGAL::Barycentric_coordinates::PRECISE, drugi to CGAL::Barycentric_coordinates::FAST.
Rysunek 96.8 Notacja dla współrzędnych Wachspress.
Aby obliczyć wagi Wachspressa, podążamy i korzystamy ze wzoru
\(w_i = \frac{C_i}{A_{i-1}A_i}\)
z \(i = 1\ kropki n\), gdzie \(n\) jest liczbą wierzchołków wielokąta. Aby obliczyć współrzędne, normalizujemy te wagi,
\(b_i = \frac{w_i}{W^{wp}}}\qquad\) z \(\qquad W^{wp} = \sum_{j=1}^n w_j.\)
Ta formuła staje się niestabilna, gdy zbliżamy się do granicy wielokąta ( \ ok. 1.0e-10 \) i bliżej). Aby rozwiązać ten problem, modyfikujemy wagi \(w_i\) jako
\(\bar{w}_i = C_i\prod_{j\not=i-1,i} A_j\).
Po normalizacji jak wyżej, daje nam to precyzyjny algorytm do obliczania współrzędnych Wachspress, ale tylko z \(O(n^2)\) wydajnością. Szybki algorytm \(O(n)\) używa standardowych wag \(w_i\). Zauważ, że matematycznie ta modyfikacja nie zmienia współrzędnych.
Wiadomo, że dla ściśle wypukłych wielokątów zbiór zerowy mianownika współrzędnych Wachspressa ( ^(W^{wp} = 0~)) jest krzywą, która (w wielu przypadkach) leży dość daleko od wielokąta. Mówiąc dokładniej, interpoluje ona punkty przecięcia kontynuacji krawędzi wielokąta. Dlatego obliczenie współrzędnych Wachspressa poza wielokątem jest możliwe tylko w punktach, które nie należą do tej krzywej.
Rysunek 96.9 Zbiór zerowy (czerwony) mianownika współrzędnych Wachspressa ^{wp}} dla nieregularnego sześciokąta foremnego.
Ostrzeżenie: nie zalecamy stosowania współrzędnych Wachspressa dla punktów zewnętrznych!
Dyskretne współrzędne harmoniczne
Dyskretne współrzędne harmoniczne mają takie same wymagania jak współrzędne Wachspressa. Są one dobrze zdefiniowane w zamknięciu dowolnego ściśle wypukłego wielokąta i, jeśli wybrany jest dokładny typ danych, są one obliczane dokładnie. Jednak, w przeciwieństwie do funkcji bazowych Wachspressa, współrzędne te nie muszą być dodatnie. W szczególności, waga w_i jest dodatnia wtedy i tylko wtedy, gdy \(\alpha+\beta < \pi\) (patrz poniższy rysunek dla notacji). W przypadku nieścisłych typów danych precyzja obliczeń wynika z zastosowanego algorytmu i wybranego typu danych. Ponownie, opisujemy dwa algorytmy obliczania współrzędnych: jeden jest precyzyjny, a drugi szybki.
Rysunek 96.10 Notacja dla dyskretnych współrzędnych harmonicznych.
Aby obliczyć dyskretne wagi harmoniczne, podążamy i korzystamy ze wzoru
\(w_i = \frac{r_{i+1}^2A_{i-1}-r_i^2B_i+r_{i-1}^2A_i}{A_{i-1}A_i}})
z \(i = 1\ n\) gdzie \(n\) jest liczbą wierzchołków wielokąta. Aby obliczyć współrzędne, normalizujemy te wagi,
\(b_i = \frac{w_i}{W^{dh}}}\qquad\) z \(\qquad W^{dh} = \sum_{j=1}^n w_j.\)
Ta formuła staje się niestabilna, gdy zbliżamy się do granicy wielokąta ( \) i bliżej). Aby rozwiązać ten problem, podobnie jak w poprzednim podrozdziale, modyfikujemy wagi \(w_i\) jako
\(\(\bar{w}_i = (r_{i+1}^2A_{i-1}-r_i^2B_i+r_{i-1}^2A_i)\prod_{j\not=i-1,i} A_j}).
Po normalizacji jak wyżej, daje to dokładny algorytm do obliczania dyskretnych współrzędnych harmonicznych, ale tylko z wydajnością \(O(n^2)\). Szybki algorytm \(O(n)\) używa standardowych wag \(w_i\). Ponownie, matematycznie ta modyfikacja nie zmienia współrzędnych.
Ostrzeżenie: jeśli chodzi o współrzędne Wachspress, nie zalecamy używania dyskretnych współrzędnych harmonicznych dla punktów zewnętrznych, ponieważ krzywa \(W^{dh} = 0\) może mieć kilka składowych, a jedna z nich interpoluje wierzchołki wielokąta. Jeśli jednak jesteś pewien, że punkt zapytania nie należy do tej krzywej, możesz obliczyć współrzędne, jak pokazano w tym przykładzie.
Współrzędne wartości średniej
W przeciwieństwie do poprzednich współrzędnych, współrzędnych wartości średniej nie można obliczyć dokładnie z powodu nieuniknionej operacji pierwiastka kwadratowego. Chociaż, jeśli używany jest dokładny typ danych, domyślna precyzja obliczeń zależy tylko od dwóch funkcji CGAL: CGAL::to_double() i CGAL::sqrt(). Z drugiej strony, współrzędne wartości średniej są dobrze zdefiniowane wszędzie na płaszczyźnie dla dowolnego prostego wielokąta. Ponadto, jeśli Twoja klasa cech udostępnia bardziej precyzyjną wersję funkcji pierwiastka kwadratowego, ostateczna precyzja obliczeń z dokładnymi typami danych będzie zależała tylko od precyzji tej funkcji.
Rysunek 96.11 Notacja dla współrzędnych wartości średniej.
Dla tych współrzędnych mamy również dwa algorytmy: jeden jest precyzyjny, a drugi szybki. Pierwszy z nich działa wszędzie na płaszczyźnie, a precyzja obliczeń zależy tylko od wybranego typu danych, z uwzględnieniem powyższych uwag. Algorytm ten opiera się na następującej formule wagowej z
\(w_i = \sigma_i\bar{w}_i\qquad\) gdzie \(\quad\{w}_i = (r_{i-1}r_{i+1}-d_{i-1}d_{i+1})^{1/2}prod_{j\not= i-1,i}(r_jr_{j+1} + d_jd_{j+1})^{1/2}) gdzie \(\quad r_i = \|d_i\. Zasadniczo, waga ta jest zawsze dodatnia na lewo od czerwonej krzywej liniowej i jest ujemna na prawo od tej krzywej, poruszając się w kierunku przeciwnym do ruchu wskazówek zegara.
Rysunek 96.12 Znaki wagi wartości średniej \ w_i\ w zależności od regionu względem wielokąta wypukłego \ i wielokąta wklęsłego \.
Po normalizacji tych wag tak jak poprzednio
\(b_i = \frac{w_i}{W^{mv}}}quad\) z \(\quad W^{mv} = \sum_{j=1}^n w_j\)
otrzymujemy dokładny algorytm \(O(n^2)\). Szybki algorytm O(n) oblicza wagi \(w_i\) używając pseudokodu stąd. Wagi te
\(w_i = \frac{t_{i-1} + t_i}{r_i}} \qquad\) z \(\quad t_i = \frac{tekst{det}(d_i, d_{i+1})}{r_ir_{i+1} + d_id_{i+1}})
są również znormalizowane. Zauważ, że są one niestabilne, jeśli punkt zapytania jest bliżej niż ∗ (∗ ok. 1.0e-10} granicy wielokąta, podobnie jak w przypadku współrzędnych Wachspressa i dyskretnych współrzędnych harmonicznych.
Wydajność
Oprócz najważniejszego wymogu stawianego współrzędnym barycentrycznym, aby były jak najdokładniejsze, bardzo ważne jest, aby były jak najszybsze w obliczaniu. Współrzędne te są używane w wielu aplikacjach, gdzie muszą być obliczane dla milionów punktów, a więc wykorzystanie współrzędnych w czasie rzeczywistym jest kluczowe. Pisząc kod, staraliśmy się spełnić to ważne wymaganie, a w tym rozdziale przedstawiamy kilka wyników dotyczących czasów obliczeń zaimplementowanych współrzędnych.
Struktura testu szybkości, który przeprowadziliśmy dla wszystkich funkcji, polega na obliczeniu wartości współrzędnych (lub wag) w >= 1 milionie ściśle wewnętrznych punktów względem pewnego wielokąta (lub trójkąta, lub odcinka). W każdej iteracji pętli tworzymy punkt zapytania, przekazujemy go do funkcji i obliczamy wszystkie związane z nim współrzędne. Pętlę tę wykonujemy 10 razy z rzędu, a czas przedstawiony na wykresie w skali logarytmicznej na końcu rozdziału jest średnią arytmetyczną wszystkich prób.
Typowy przykład tego testu wydajności dla współrzędnych trójkąta z ograniczoną liczbą punktów zapytania można znaleźć poniżej. Przykład ten ilustruje również jak skonstruować iterator i przekazać go do klasy. W tym przykładzie tworzymy iterator, który zapisuje wartości współrzędnych dla każdego nowego punktu zapytania nad wartościami współrzędnych poprzedniego punktu w standardowej tablicy C++ o stałym rozmiarze, dzięki czemu pamięć jest alokowana tylko raz.
Plik Barycentric_coordinates_2/Triangle_coordinates_speed_test.cpp
#include <CGAL/Real_timer.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Barycentric_coordinates_2/Triangle_coordinates_2.h>
/// Skonstruuj iterator, który pobiera jako dane wejściowe bieżący typ danych i wskaźnik do pierwszego elementu standardowej tablicy C++.
template<typename Scalar>
class overwrite_iterator
{
private:
Scalar* pointer;
public:
explicit overwrite_iterator(Scalar* new_pointer) : pointer(new_pointer) { }
// Istnieją tylko dwie operacje, które musimy przeciążyć, aby móc korzystać z klasy Triangle_coordinates_2.
// Ta operacja ma na celu zwrócenie bieżącej wartości współrzędnych.
inline Scalar& operator* () { return *pointer; }
// Ta operacja ma na celu zwiększenie indeksu współrzędnej.
inline void operator++ () { ++pointer; }
};
// Kilka wygodnych typedefów.
typedef CGAL::Real_timer Timer;
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef Kernel::FT Scalar;
typedef Kernel::Point_2 Point;
typedef overwrite_iterator<Scalar> Overwrite_iterator;
typedef CGAL::Barycentric_coordinates::Triangle_coordinates_2<Kernel> Triangle_coordinates;
using std::cout; using std::endl; using std::string;
int main()
{
// Liczba współrzędnych x i y razem daje liczbę punktów.
const int number_of_x_coordinates = 100000;
const int number_of_y_coordinates = 1000;
// Liczba przebiegów do obliczenia średniej arytmetycznej czasu.
const int number_of_runs = 10;
// Oblicz jednolitą wielkość kroku wzdłuż kierunków x i y, aby zmienić współrzędne.
const Scalar zero = Scalar(0);
const Scalar one = Scalar(1);
const Scalar two = Scalar(2);
const Scalar x_step = one / Scalar(number_of_x_coordinates);
const Scalar y_step = one / Scalar(number_of_y_coordinates);
// Utwórz trójkąt prosty z niewielkim przesunięciem od zera.
const Point first_vertex(zero – x_step, zero – x_step);
const Point second_vertex(two + y_step, zero – x_step);
const Point third_vertex(zero – x_step, two + y_step);
// Instantiate the class Triangle_coordinates_2 for the right triangle defined above.
Triangle_coordinates triangle_coordinates(first_vertex, second_vertex, third_vertex);
// Utwórz instancję standardowej tablicy C++ do przechowywania współrzędnych.
// Ma ona stały rozmiar = 3 = liczba wierzchołków.
Scalar coordinates = {0, 0, 0};
// Przekaż wskaźnik do pierwszego elementu tablicy ze współrzędnymi, aby je nadpisać.
Overwrite_iterator it( &(współrzędne) );
// Utwórz timer.
Timer time_to_compute;
double time = 0.0;
for(int i = 0; i < number_of_runs; ++i) { // Liczba przebiegów
time_to_compute.start(); // Zegar startowy
for(Scalar x = zero; x <= jeden; x += x_step) {
for(Scalar y = zero; y <= jeden; y += y_step)
triangle_coordinates(Point(x, y), it); // Oblicz 3 wartości współrzędnych dla każdego wygenerowanego punktu
}
time_to_compute.stop(); // Zatrzymaj zegar
time += time_to_compute.time(); // Dodaj czas bieżącego przebiegu do całego czasu
time_to_compute.reset(); // Zresetuj czas
}
// Oblicz średnią arytmetyczną wszystkich przebiegów.
const double mean_time = time / number_of_runs;
// Wyprowadź wynikowy czas.
cout.precision(10);
cout << endl << „CPU time to compute triangle coordinates for ”
<< number_of_x_coordinates * number_of_y_coordinates << ” points = ” << mean_time << ” seconds.”;
cout << endl << endl;
return EXIT_SUCCESS;
}
Czas obliczania współrzędnych zależy od wielu czynników, takich jak alokacja pamięci, jądro wejściowe, kontener wyjściowy, liczba punktów itp. W naszych testach wykorzystaliśmy najbardziej standardowe funkcje C++ i CGAL z minimalną alokacją pamięci. Dlatego też, przedstawiony czas końcowy jest średnim czasem, jakiego można się spodziewać bez głębokiej optymalizacji, ale wciąż z wydajną alokacją pamięci. Oznacza to również, że może się on różnić w zależności od wykorzystania pakietu.
Do wszystkich testów wykorzystaliśmy MacBooka Pro 2011 z procesorem Intel Core i7 2 GHz (2 rdzenie) i 8 GB pamięci DDR3 1333 MHz. Zainstalowanym systemem operacyjnym był OS X 10.9 Maverick. Do skompilowania zestawu testów prędkości wykorzystaliśmy kompilator Clang 5.0 64bit. Uzyskane timingi można znaleźć na poniższym rysunku.
Rysunek 96.13 Czas w sekundach na obliczenie wartości współrzędnych dla wielokąta o rzednących wierzchołkach w 1 milionie punktów za pomocą szybkich algorytmów (przerywany) i powolnych (stały) dla współrzędnych Wachspressa (niebieski), dyskretnej harmonicznej (czerwony) i wartości średniej (zielony).
Z powyższego rysunku łatwo zauważyć, że algorytm \(O(n^2)\) jest tak samo szybki jak algorytm \(O(n)\), jeśli mamy wielokąt z małą liczbą wierzchołków. Jednak wraz ze wzrostem liczby wierzchołków, zgodnie z oczekiwaniami, algorytm liniowy przewyższa algorytm kwadratowy. Jedną z przyczyn takiego zachowania jest to, że dla małej liczby wierzchołków mnożenie elementów wewnątrz algorytmu ∗(O(n^2)∗) za pomocą szybkich algorytmów ∗(O(n)∗) (przerywane) i powolnych algorytmów ∗(O(n^2)∗) (stałe) zajmuje prawie tyle samo czasu, co odpowiednie dzielenie w algorytmie ∗(O(n)∗). Dla wielokąta z wieloma wierzchołkami to mnożenie jest znacznie wolniejsze.
Szczegóły implementacji
Ogólny projekt pakietu został opracowany w 2013 roku przez Dmitry’ego Anisimova i Davida Bommesa z wieloma przydatnymi komentarzami Kaia Hormanna i Pierre’a Allieza. Pakiet składa się z 6 klas, 2 enumeracji i jednej przestrzeni nazw. Odpowiednie iteratory są używane do zapewnienia efektywnego dostępu do danych i przekazania ich do jednego z generycznych algorytmów obliczania współrzędnych. Raz zainicjowane dla wielokąta (trójkąta, odcinka) współrzędne mogą być obliczane wielokrotnie dla różnych punktów zapytania w odniesieniu do wszystkich wierzchołków danego wielokąta (trójkąta, odcinka). Wszystkie klasy są w pełni szablonowe i mają prosty i podobny wygląd. W szczególności, stosujemy tę samą konwencję nazewnictwa dla wszystkich funkcji. Jednak liczba funkcji może się różnić w zależności od klasy.
Zaimplementowane algorytmy obliczania współrzędnych nie zależą od konkretnego jądra, a wszystkie współrzędne mogą być obliczone dokładnie, jeśli używane jest dokładne jądro, z wyjątkiem współrzędnych wartości średniej. Te ostatnie współrzędne wymagają operacji pierwiastka kwadratowego, co powoduje nieco gorszą precyzję w przypadku dokładnych typów danych ze względu na czasową konwersję na typ zmiennoprzecinkowy. Obliczone współrzędne mogą być przechowywane w dowolnym kontenerze, jeśli podany jest odpowiedni iterator wyjściowy.
Warto zauważyć, że klasa CGAL::Barycentric_coordinates::Segment_coordinates_2 służy do obliczania uogólnionych współrzędnych barycentrycznych wzdłuż granicy wielokąta. Można więc zastosować sztuczkę ze współrzędnymi odcinka z rozdziału Degeneracje i przypadki specjalne, jeśli jest się przekonanym, że punkt musi leżeć dokładnie na granicy wielokąta, ale z powodu pewnych niestabilności numerycznych nie leży.
Pakiet jest zaimplementowany w taki sposób, że później, jeśli zajdzie taka potrzeba, można do niego łatwo dodać inne dwuwymiarowe uogólnione współrzędne barycentryczne.
Teoria dwuwymiarowych uogólnionych współrzędnych barycentrycznych
W 1827 roku niemiecki matematyk i astronom teoretyczny August Ferdinand Möbius (1790-1868) zaproponował metodę znajdowania współrzędnych punktu na płaszczyźnie w odniesieniu do wierzchołków trójkąta. Współrzędne te nazywane są współrzędnymi barycentrycznymi trójkąta (czasami współrzędnymi powierzchniowymi) i są one szeroko stosowane w różnych aplikacjach. Niektóre z tych zastosowań to interpolacja liniowa na trójkącie i test inkluzji trójkąta. Pierwszy z nich jest używany do tak zwanego cieniowania, a drugi pojawia się w kroku rasteryzacji, gdy obraz w formacie grafiki wektorowej musi być przekształcony w obraz rastrowy.
Trójkątne współrzędne barycentryczne mają wiele ważnych właściwości, w tym stałą i liniową precyzję, własność Lagrange’a i dodatniość wewnątrz trójkąta. Te właściwości sprawiają, że współrzędne te są unikalnym narzędziem w wielu dziedzinach nauki. Jeśli ograniczymy współrzędne trójkąta do jednej z krawędzi trójkąta i jego linii pomocniczej, otrzymamy współrzędne barycentryczne względem odcinka i nazwiemy je współrzędnymi segmentu.
Pokażmy kilka wykresów dla współrzędnych opisanych powyżej. Aby wykreślić współrzędne odcinka, bierzemy prostą y = 0,4 i definiujemy odcinek y na tej prostej. Następnie pobieramy próbkę tego odcinka i obliczamy współrzędne odcinka dla wszystkich punktów próbki. Jeśli wykreślimy funkcję współrzędnych odcinka we wszystkich zdefiniowanych punktach względem wierzchołka, otrzymamy niebieską linię przedstawioną na poniższym rysunku. Rośnie ona od zera przy wierzchołku \(v_0\) do jednego przy wierzchołku \(v_1\).
Rysunek 96.14 Współrzędne odcinka (niebieski) dla wszystkich punktów odcinka (zielony) względem wierzchołka \(v_1 = (2,0,\ 0,4)\).
Jeżeli chcemy wykreślić współrzędne trójkątów, postępujemy podobnie. Bierzemy trójkąt \(\) na płaszczyźnie i próbkujemy jego wnętrze i granicę pewną liczbą punktów. Po pobraniu próbki, wykreślamy jedną z funkcji współrzędnych trójkąta (tutaj w odniesieniu do trzeciego wierzchołka trójkąta) we wszystkich zdefiniowanych punktach próbki. Podobnie, możemy wykreślić funkcję współrzędnych w odniesieniu do pierwszego lub drugiego wierzchołka. Otrzymana funkcja jest liniowa (pokazana na rysunku poniżej), która rośnie od zera wzdłuż pierwszej krawędzi do jedynki przy wybranym wierzchołku \(v_2\).
Rysunek 96.15 Współrzędne trójkąta względem \(v_2 = (1,0,\ 2,0)\). Kolorowy pasek wskazuje zakres wartości dla wybranej współrzędnej.
Ponieważ wiele zastosowań wymaga pracy z bardziej złożonymi planarnymi figurami geometrycznymi niż odcinki i trójkąty, naturalnym wydaje się zbadanie uogólnionej wersji współrzędnych trójkątów względem dowolnych wielokątów. Pierwsza próba została podjęta w 1975 roku przez E. L. Wachspressa, a wynikające z niej uogólnione współrzędne barycentryczne są obecnie nazywane współrzędnymi Wachspressa. Współrzędne te są dobrze zdefiniowane dla dowolnych ściśle wypukłych wielokątów i mają wszystkie własności współrzędnych trójkątnych. Niestety, nie są one dobrze zdefiniowane dla słabo wypukłych i wklęsłych wielokątów.
Analogicznie do poprzednich przypadków, chcemy wykreślić współrzędne Wachspressa i zobaczyć, jak one wyglądają. Wybierzmy nieregularny sześciokąt, lekko go obróćmy i przesuńmy jeden z jego wierzchołków w kierunku prostej przechodzącej przez jego dwa sąsiednie wierzchołki. Próbkujemy wnętrze i granicę tego wielokąta tak jak poprzednio i wykreślamy funkcję współrzędnych względem wierzchołka, który przesunęliśmy we wszystkich punktach próbki. Widzimy, że otrzymujemy gładką funkcję, która jest liniowa wzdłuż wszystkich krawędzi i rośnie od zera do jedynki, na co wskazuje kolorowy pasek.
Rysunek 96.16 Funkcja współrzędnych Wachspressa względem wskazanego wierzchołka o wartościach od zera do jednego, jak wskazuje kolorowy pasek.
Inny rodzaj uogólnionych współrzędnych barycentrycznych sięga Pinkall i Polthier w 1993 roku oraz Eck et al. w 1995 roku w kontekście parametryzacji siatki trójkątów. Nazywane są one dyskretnymi współrzędnymi harmonicznymi. Współrzędne te są dobrze zdefiniowane, podobnie jak współrzędne Wachspressa, dla dowolnych ściśle wypukłych wielokątów i dziedziczą wszystkie własności współrzędnych trójkątnych poza dodatniością wewnątrz wielokąta, gdyż dla niektórych wielokątów mogą przyjmować wartości ujemne. Inną ciekawą własnością tych funkcji współrzędnych jest to, że pokrywają się one ze współrzędnymi Wachspressa dla dowolnego wielokąta, którego wierzchołki leżą na wspólnym okręgu.
Aby wykreślić dyskretne współrzędne harmoniczne bierzemy ten sam wielokąt co dla współrzędnych Wachspressa i wykreślamy funkcję współrzędnych względem tego samego wierzchołka. Ponownie, otrzymujemy gładką funkcję, która jest liniowa wzdłuż wszystkich krawędzi i rośnie od zera do jednego. Izolinie na wykresie pokazują różnicę między współrzędnymi harmonicznymi dyskretnymi i Wachspressa dla wybranego wielokąta i wierzchołka.
Rysunek 96.17 Funkcja współrzędnych harmonicznych dyskretnych względem wskazanego wierzchołka o wartościach od zera do jeden, jak wskazuje kolorowy pasek.
Ostatnim rodzajem uogólnionych współrzędnych barycentrycznych, które omawiamy, są współrzędne wartości średniej zaproponowane przez M. Floatera w 2003 roku. W oparciu o twierdzenie o średniej wartości, współrzędne te, w odróżnieniu od współrzędnych Wachspressa i dyskretnych współrzędnych harmonicznych, są dobrze zdefiniowane dla dowolnych wielokątów prostych, dziedziczą wszystkie własności współrzędnych trójkątnych dla dowolnego wielokąta wypukłego i brakuje im jedynie własności pozytywności dla ogólnych wielokątów wklęsłych. Hormann i Floater udowadniają, że współrzędne te są dodatnie wewnątrz jądra wielokąta gwiaździstego. Są one również dodatnie w zamknięciu dowolnego czworokąta. Podobnie jak dyskretne wagi harmoniczne, wagi wartości średniej są często używane w kontekście parametryzacji siatki trójkątów.
Aby pokazać szczególne zachowanie współrzędnych wartości średniej w zastosowaniu do wielokątów wklęsłych, bierzemy wielokąt gwiaździsty o dziesięciu wierzchołkach, próbkujemy jego wnętrze i granicę, i wykreślamy funkcję współrzędnych względem czwartego wierzchołka \(v_3\). Jak wskazuje kolorowy pasek, otrzymana funkcja rośnie od wartości lekko ujemnej do jedynki w wybranym wierzchołku. Jest ona również gładka wewnątrz wielokąta i liniowa wzdłuż wszystkich krawędzi.
Rysunek 96.18 Współrzędne wartości średniej względem wierzchołka \(v_3\). Kolorowy pasek wskazuje zakres wartości dla wybranej funkcji współrzędnych.
Ciekawy fakt: wszystkie współrzędne omawiane w tym rozdziale i zaimplementowane w pakiecie pochodzą z jednej i tej samej rodziny uogólnionych współrzędnych barycentrycznych o nazwie 3-punktowa rodzina współrzędnych .
Podziękowania
