Polygon - Mittelpunkt und Schwerepunkt?
-
Wenn in einer Formel
steht, dann ist nicht
gemeint. Das hast du gründlich missverstanden. Aber ich kann dir nicht im Rahmen eines Forenbeitrags die Indexschreibweise in der Mathematik beibringen. Gemeint ist das "(i+1)'te x", aber das erklärt es nicht wirklich, oder? Zumal das darauf hindeutet, dass du dann bloß den nächsten Teil missverstehen würdest. Du gehst noch zur Schule und hast so etwas noch nie gesehen?
-
So könnte die korrekte Umsetzung in Python aussehen:
import dataclasses import itertools @dataclasses.dataclass class Point2D: x: float y: float def polygon_centroid(polygon: list[Point2D]): # Voraussetzung der Formel: Punkt N ist wieder Punkt 0 polygon = polygon + [polygon[0]] area = 1 / 2 * sum( (point1.x * point2.y - point2.x * point1.y) for point1, point2 in itertools.pairwise(polygon) ) print(area) centroid_x = 1 / (6 * area) * sum( (point1.x + point2.x) * (point1.x * point2.y - point2.x * point1.y) for point1, point2 in itertools.pairwise(polygon) ) centroid_y = 1 / (6 * area) * sum( (point1.y + point2.y) * (point1.x * point2.y - point2.x * point1.y) for point1, point2 in itertools.pairwise(polygon) ) return Point2D(centroid_x, centroid_y) polygon = [Point2D(10, 10), Point2D(30, 20), Point2D(15, 15)] print(polygon_centroid(polygon))Es ist leider ein bisschen schwer zu verstehen, was da genau passiert, weil es halt einfach keine ganz einfache Formel ist. Man geht durch ein Eckpunktpaar nacheinander (in meinem Programm dann jeweils
point1undpoint2), verwurstet dann deren x- und y-Koordinate zu einer Zahl, und summiert das über all diese Eckpunktpaare. Verständlicher kann ich das leider nicht erklären, außer ich könnte es dir an einer Tafel aufzeichnen. Ich hoffe, es wird wenigstens halbwegs klar, was mit der Schreibweise auf Wikipedia (oder wo immer du das nachgeschlagen hast), gemeint ist, indem du das Programm und die Formel vergleichst.PS: Das Ergebnis für dein Polygon ist
Da das ein Dreieck ist, hätte man das natürlich auch einfacher haben können, weil beim Dreieck ist der Schwerpunkt tatsächlich der Mittelwert der x- und y- Koordinaten (jeweilig). Zeigt aber, dass das Programm funktioniert.
-
@SeppJ sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Wenn in einer Formel
steht, dann ist nicht
gemeint. Das hast du gründlich missverstanden. Aber ich kann dir nicht im Rahmen eines Forenbeitrags die Indexschreibweise in der Mathematik beibringen. Zumal das darauf hindeutet, dass du dann bloß den nächsten Teil missverstehen würdest. Du gehst noch zur Schule und hast so etwas noch nie gesehen?
Pardon! Jetzt hab ich's:
#include <vector> #include <iostream> namespace std { vector<double> calculateCenterOfGravity(vector<vector<double>> poly) { vector<double> result; double sum1 = 0; double sum2 = 0; double area = 0; for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double s1 = (x0 + x1) * (x0 * y1 - x1 * y0); double s2 = (y0 + y1) * (x0 * y1 - x1 * y0); double s3 = x0 * y1 - x1 * y0; sum1 += s1 / 6.0; sum2 += s2 / 6.0; area += s3 * 0.5; } result.push_back(sum1 / area); result.push_back(sum2 / area); return result; } } // namespace std int main(int argc, char const *argv[]) { std::vector<std::vector<double>> poly; poly.push_back(std::vector<double>{10, 10}); poly.push_back(std::vector<double>{30, 20}); poly.push_back(std::vector<double>{15, 15}); auto r = calculateCenterOfGravity(poly); std::cout << "x:" << r[0] << " y:" << r[1] << "\n"; return 0; }$ g++ Polygone.cpp -o out $ ./out x:18.3333 y:15Das scheint ja nun richtig zu sein.
@SeppJ sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Du gehst noch zur Schule und hast so etwas noch nie gesehen?
Ich gehe zur Schule.
Verstehe aber nicht, was das damit zu tun hat.
-
@SeppJ sagte in Polygon - Mittelpunkt und Schwerepunkt?:
was mit der Schreibweise auf Wikipedia (oder wo immer du das nachgeschlagen hast), gemeint ist, indem du das Programm und die Formel vergleichst.
-
Ja, das sieht (beim Diagonallesen) so aus, als würde es passen.
@TaucherOne sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Ich gehe zur Schule.
Verstehe aber nicht, was das damit zu tun hat.Weil das in der Mathematik der Oberstufe (aber eventuell auch nur im Leistungskurs) drankommen würde, und dir offenbar noch nicht bekannt war. Und wenn du noch zur Schule gehst, sollte ich bei den Erklärungen nicht zu viel der typischen Ausdrucksweise der Hochschulmathematik benutzen.
Dann verstehe ich auch, wieso dir nach unseren ersten Antworten offenbar nicht so schnell klar wurde, wieso es im allgemeinen keinen Punkt gibt, der von allen Ecken gleich weit entfernt ist, oder wieso der Schwerpunkt nochmal ganz was anderes wäre. Das war von der Erklärung her wahrscheinlich zu knapp für die Schule, wo man mit so etwas gerade erst anfängt und kein intuitives Gefühl dafür hat.
-
@SeppJ Sorry, ich kannte diese Schreibweise einfach noch nicht:
X = SUM[(Xi + Xi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A Y = SUM[(Yi + Yi+1) * (Xi * Yi+1 - Xi+1 * Yi)] / 6 / A
-
Eben. Das kommt in der 11? Oder gar der 12? Keine Ahnung, bin schon lange nicht mehr in einer Schule gewesen
Dann brauchst du das aber auch noch nicht kennen. Und selbst dann ist das hier nicht gerade die üblichste Schreibweise.
-
Der Vollständigkeit wegen hier noch ein Monte-Carlo-Ansatz... (probabilistisches Verfahren)
#include <vector> #include <iostream> namespace std { vector<double> calculateCenterOfGravity(vector<vector<double>> poly) { double sum1 = 0; double sum2 = 0; double area = 0; for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double s1 = (x0 + x1) * (x0 * y1 - x1 * y0); double s2 = (y0 + y1) * (x0 * y1 - x1 * y0); double s3 = x0 * y1 - x1 * y0; sum1 += s1 / 6.0; sum2 += s2 / 6.0; area += s3 * 0.5; } vector<double> result; result.push_back(sum1 / area); result.push_back(sum2 / area); return result; } vector<double> calculateCenterOfGravityMonteCarlo(vector<vector<double>> poly, const unsigned int seed) { srand(seed); double maxx = -(double)(unsigned int)(-1), minx = (unsigned int)(-1), maxy = -(double)(unsigned int)(-1), miny = (unsigned int)(-1); for (auto p : poly) { if (p[0] > maxx) { maxx = p[0]; } if (p[0] < minx) { minx = p[0]; } if (p[1] > maxy) { maxy = p[1]; } if (p[1] < miny) { miny = p[1]; } } double x, y; double currentBest = (unsigned int)(-1); for (size_t i = 0; i < 5000; i++) { double rx = (double)rand() / RAND_MAX * (maxx - minx) + minx; double ry = (double)rand() / RAND_MAX * (maxy - miny) + miny; double rbest = -(double)(unsigned int)(-1); for (auto p : poly) { double d = (p[0] - rx) * (p[0] - rx) + (p[1] - ry) * (p[1] - ry); if (d > rbest) { rbest = d; } } if (rbest < currentBest) { x = rx; y = ry; currentBest = rbest; cout << " Found: " << x << " " << y << "\n"; } } vector<double> result; result.push_back(x); result.push_back(y); return result; } } // namespace std int main(int argc, char const *argv[]) { std::vector<std::vector<double>> poly; poly.push_back(std::vector<double>{10, 10}); poly.push_back(std::vector<double>{30, 20}); poly.push_back(std::vector<double>{5, 5}); auto r = calculateCenterOfGravity(poly); std::cout << "x:" << r[0] << " y:" << r[1] << "\n"; r = calculateCenterOfGravityMonteCarlo(poly, 7); std::cout << "x:" << r[0] << " y:" << r[1] << "\n"; return 0; }Seltsamerweise ist das bei mir aber ziemlich ungenau... Ich habe da
x:17.2598 y:12.9228anstattx:15 y:11.6667raus. Ich kann mir das nicht erklären.
-
@TaucherOne sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Seltsamerweise ist das bei mir aber ziemlich ungenau... Ich habe da
x:17.2598 y:12.9228anstattx:15 y:11.6667raus. Ich kann mir das nicht erklären.Ist ja schonmal die richtige Grössenordnung
... ich kenne den Algo nicht und hab ihn mir nicht genau angesehen, aber vielleicht konvergiert der einfach nicht sonderlich gut (sofern korrekt implementiert). Hast du mal mit der Anzahl der Iterationen gespielt? Ich würde so ins Blaue vermuten, dass das die 5000in Zeile 57 ist. Wirds mit10000oder100000eventuell besser?
-
@Finnegan sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Dadurch ergibt sich das Problem mit dem hohen Rang des zentralen Vertex erst gar nicht, da jedes Dreieck nur einmal in das Ergebnis einfliesst.
Edit: Achja, und wenn alle Dreiecke die gleiche Fläche haben, dann zerfällt das tatsächlich zu (p_0 + p_1 + ... p_n) / n, d.h. für eine gute Triangulierung erhält man selbst mit der falschen Methode eine "ziemlich richtiges" Ergebnis... aber eben nur unter ganz bestimmten VoraussetzungenDanke für die Info, ist ja sehr einfach.
Gute Triangulierungen sehe ich selten. Ich bin aber auch im Bereich der digitalen Geländemodelle unterwegs und so mancher Geländeschnitt oder so manche Bruchkante verändert die Vermaschung öfters negativ. In einem extremen Fall hatte ich eine Vermaschung wo auf 1% der Fläche 99% der Dreiecke waren und das völlig korrekt!
-
@Finnegan Habe den Fehler, glaube ich, entdeckt ...
Wenn man sich noch einmal die Definition des Schwerpunkts ansieht, dann ist die Entfernung zu allen SEITEN maximal. Ich habe aber getestet, ob die Entfernung zu allen Ecken (Punkten) maximal ist ... Das wäre demnach nicht der Schwerpunkt (sondern stattdessen ein anderer Punkt).
Hier noch einmal geändert:
#include <vector> #include <iostream> #include <complex> namespace std { vector<double> calculateCenterOfGravity(vector<vector<double>> poly) { double sum1 = 0; double sum2 = 0; double area = 0; for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double s1 = (x0 + x1) * (x0 * y1 - x1 * y0); double s2 = (y0 + y1) * (x0 * y1 - x1 * y0); double s3 = x0 * y1 - x1 * y0; sum1 += s1 / 6.0; sum2 += s2 / 6.0; area += s3 * 0.5; } vector<double> result; result.push_back(sum1 / area); result.push_back(sum2 / area); return result; } double pDistance(double x, double y, double x1, double y1, double x2, double y2) { double A = x - x1; double B = y - y1; double C = x2 - x1; double D = y2 - y1; double dot = A * C + B * D; double len_sq = C * C + D * D; double param = -1; if (len_sq != 0) { param = dot / len_sq; } double xx, yy; if (param < 0) { xx = x1; yy = y1; } else if (param > 1) { xx = x2; yy = y2; } else { xx = x1 + param * C; yy = y1 + param * D; } double dx = x - xx; double dy = y - yy; return sqrt(dx * dx + dy * dy); } vector<double> calculateCenterOfGravityMonteCarlo(vector<vector<double>> poly, const unsigned int seed, double delta) { srand(seed); double maxx = -(double)(unsigned int)(-1), minx = (unsigned int)(-1), maxy = -(double)(unsigned int)(-1), miny = (unsigned int)(-1); for (auto p : poly) { if (p[0] > maxx) { maxx = p[0]; } if (p[0] < minx) { minx = p[0]; } if (p[1] > maxy) { maxy = p[1]; } if (p[1] < miny) { miny = p[1]; } } double x, y; double currentBest = -(double)(unsigned int)(-1); for (size_t i = 0; i < 5000; i++) { a: double rx = (double)rand() / (double)RAND_MAX * (maxx - minx) + minx; double ry = (double)rand() / (double)RAND_MAX * (maxy - miny) + miny; double rBest = (unsigned int)(-1); for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double d = pDistance(rx, ry, x0, y0, x1, y1); if (d > delta) { i++; goto a; } if (d < rBest) { rBest = d; } } if (rBest > currentBest) { x = rx; y = ry; currentBest = rBest; cout << " Found: " << x << " " << y << " " << rBest << "\n"; } } vector<double> result; result.push_back(x); result.push_back(y); return result; } } // namespace std int main(int argc, char const *argv[]) { std::vector<std::vector<double>> poly; poly.push_back(std::vector<double>{10, 10}); poly.push_back(std::vector<double>{30, 20}); poly.push_back(std::vector<double>{5, 5}); auto r = calculateCenterOfGravity(poly); std::cout << "x:" << r[0] << " y:" << r[1] << "\n"; r = calculateCenterOfGravityMonteCarlo(poly, 7, 0.62); std::cout << "x:" << r[0] << " y:" << r[1] << "\n"; return 0; }Das
gotoist natürlich unschön ... Aber eigentlich fehlt dort noch ein Test, ob die zwei jeweils zufällig gewählten Punkte innerhalb des Polygons liegen. Deshalb habe ichdelta(magischer Wert) undgotoeingefügt.Ich denke, es könnte auch noch dahin optimiert werden, dass es schneller konvergiert (wenn man max verringert oder min erhöht).
BTW:
deltasollte möglichst klein, aber nicht zu klein gewählt werden ... das ist ein bisschen Quacksalberei.
-
Der Schwerpunkt eines Polygons liegt nicht notwendigerweise in seinem Inneren.
-
@SeppJ sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Der Schwerpunkt eines Polygons liegt nicht notwendigerweise in seinem Inneren.
Danke für den Hinweis!
Edit: Mein Beispiel war auch nicht sonderlich toll ...
#include <vector> #include <iostream> #include <complex> namespace std { vector<double> calculateCenterOfGravity(vector<vector<double>> poly) { double sum1 = 0; double sum2 = 0; double area = 0; for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double s1 = (x0 + x1) * (x0 * y1 - x1 * y0); double s2 = (y0 + y1) * (x0 * y1 - x1 * y0); double s3 = x0 * y1 - x1 * y0; sum1 += s1 / 6.0; sum2 += s2 / 6.0; area += s3 * 0.5; } vector<double> result; result.push_back(sum1 / area); result.push_back(sum2 / area); return result; } double pDistance(double x, double y, double x1, double y1, double x2, double y2) { double A = x - x1; double B = y - y1; double C = x2 - x1; double D = y2 - y1; double dot = A * C + B * D; double len_sq = C * C + D * D; double param = -1; if (len_sq != 0) { param = dot / len_sq; } double xx, yy; if (param < 0) { xx = x1; yy = y1; } else if (param > 1) { xx = x2; yy = y2; } else { xx = x1 + param * C; yy = y1 + param * D; } double dx = x - xx; double dy = y - yy; return sqrt(dx * dx + dy * dy); } vector<double> calculateCenterOfGravityMonteCarlo(vector<vector<double>> poly, const unsigned int seed, double delta) { const size_t NUM_ITERATIONS = 50000; srand(seed); double maxx = -(double)(unsigned int)(-1), minx = (unsigned int)(-1), maxy = -(double)(unsigned int)(-1), miny = (unsigned int)(-1); for (auto p : poly) { if (p[0] > maxx) { maxx = p[0]; } if (p[0] < minx) { minx = p[0]; } if (p[1] > maxy) { maxy = p[1]; } if (p[1] < miny) { miny = p[1]; } } double x, y; double currentBest = -(double)(unsigned int)(-1); for (size_t i = 0; i < NUM_ITERATIONS; i++) { double rx = (double)rand() / (double)RAND_MAX * (maxx - minx) + minx; double ry = (double)rand() / (double)RAND_MAX * (maxy - miny) + miny; double rBest = (unsigned int)(-1); for (size_t i = 0; i < poly.size() - 1; i++) { auto p1 = poly[i]; auto p2 = poly[i + 1]; double x0 = p1[0]; double x1 = p2[0]; double y0 = p1[1]; double y1 = p2[1]; double d = pDistance(rx, ry, x0, y0, x1, y1); if (d > delta) { goto a; } if (d < rBest) { rBest = d; } } if (rBest > currentBest) { x = rx; y = ry; currentBest = rBest; cout << " Found: " << x << " " << y << " " << currentBest << "\n"; } a:; } vector<double> result; result.push_back(x); result.push_back(y); return result; } } // namespace std int main(int argc, char const *argv[]) { { // Quadrat: std::vector<std::vector<double>> poly; poly.push_back(std::vector<double>{1, 1}); poly.push_back(std::vector<double>{2, 1}); poly.push_back(std::vector<double>{2, 2}); poly.push_back(std::vector<double>{1, 2}); poly.push_back(std::vector<double>{1, 1}); // schlißende Ecke vergessen auto r1 = calculateCenterOfGravity(poly); auto r2 = calculateCenterOfGravityMonteCarlo(poly, 7, 100.0); // "delta" eigentlich unnötig std::cout << "x:" << r1[0] << " y:" << r1[1] << "\n"; std::cout << "x:" << r2[0] << " y:" << r2[1] << "\n"; } { // Konvexes Polygon: std::vector<std::vector<double>> poly; poly.push_back(std::vector<double>{3, 1}); poly.push_back(std::vector<double>{5, 2}); poly.push_back(std::vector<double>{4, 3.5}); poly.push_back(std::vector<double>{2, 4}); poly.push_back(std::vector<double>{1, 2.5}); poly.push_back(std::vector<double>{3, 1}); // schlißende Ecke vergessen auto r1 = calculateCenterOfGravity(poly); auto r2 = calculateCenterOfGravityMonteCarlo(poly, 7, 100.0); // "delta" eigentlich unnötig std::cout << "x:" << r1[0] << " y:" << r1[1] << "\n"; std::cout << "x:" << r2[0] << " y:" << r2[1] << "\n"; } return 0; }Jetzt sollte es doch zumindest vom Prinzip her richtig sein.
-
Hier wäre eine Optimierung für eine schnellere Konvergenz ...
if (rBest > currentBest) { x = rx; y = ry; currentBest = rBest; cout << " Found: " << x << " " << y << " " << currentBest << "\n"; maxx -= (maxx - minx) * 0.1; minx += (maxx - minx) * 0.1; maxy -= (maxy - miny) * 0.1; miny += (maxy - miny) * 0.1; }Ich kann die Korrektheit aber nicht zeigen.
Sorry, so:
if (rBest > currentBest) { x = rx; y = ry; currentBest = rBest; cout << " Found: " << x << " " << y << " " << currentBest << "\n"; const double FACTOR = 0.2; const double f1 = (maxx - minx) * FACTOR; const double f2 = (maxy - miny) * FACTOR; maxx -= f1; minx += f1; maxy -= f2; miny += f2; }
-
@TaucherOne sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Ich kann die Korrektheit aber nicht zeigen.
Ich bin mir auch nicht sicher, ob da die Korrektheit überhaupt garantiert ist. Da gibt es bestimmt Beispiele, wo dein Algorithmus fröhlich für immer in lokalen Minima verschwinden wird.
-
Da haste wahr.
Obwohl das bei konvexen Polygonen eigentlich nicht der Fall sein sollte ...BTW, zusammen mit der Optimierung geht es fast schon in Richtung Simulated Annealing.
Aber eigentlich ist das alles Overkill, man kann ja die Formel für den Schwerpunkt nutzen. ...
-
@SeppJ sagte in Polygon - Mittelpunkt und Schwerepunkt?:
Der Schwerpunkt eines Polygons liegt nicht notwendigerweise in seinem Inneren.
Im Allgemeinen nicht, aber für konvexe Polygone schon. Wenn ich den Algorithmus richtig verstanden habe, liegen alle Zwischenschritte wie auch das Endergebnis zwischen zwei Punkten des Polygons. Aus der Definition eines konvexen Polygons folgt dann, dass die alle im Inneren* liegen.
* Inklusive Rand. Bin grad etwas verwirrt - gehören die Randpunkte dazu, wenn man vom "Inneren" spricht? Vielleicht bin ich auch einfach gerade etwas zu pedantisch
-
https://stackoverflow.com/questions/47004208/should-point-on-the-edge-of-polygon-be-inside-polygon
The simple answer to the question is that a point on the edge is neither inside, nor outside of polygon, it is on the boundary of the polygon - the third option. [...]
Antwort: Weder noch.
(Solche Fälle hat man aber auch selten)
-
@TaucherOne sagte in Polygon - Mittelpunkt und Schwerepunkt?:
https://stackoverflow.com/questions/47004208/should-point-on-the-edge-of-polygon-be-inside-polygon
The simple answer to the question is that a point on the edge is neither inside, nor outside of polygon, it is on the boundary of the polygon - the third option. [...]
Antwort: Weder noch.
(Solche Fälle hat man aber auch selten)
Der Rand ist eh so dünn, dass er gar nicht vorhanden ist. Den kann man getrost ignorieren, wenn man kein Mathematiker ist
Oder um Stanisław Lem zu zitieren:
Nachdem er den Drahttelegraphen erfunden hatte, zog er schliesslich einen Draht so fein aus, dass gar keiner mehr da war. Und so entstand der drahtlose Telegraph.
-
Ich hab es jetzt noch einmal schöner gemacht:
#include <vector> #include <iostream> #include <complex> #include <functional> using namespace std; void forEachPolyEdge(vector<vector<double>> poly, function<void(vector<double>)> f) { for (auto &p : poly) { f(p); } } void forEachPolySite(vector<vector<double>> poly, function<void(vector<vector<double>>)> f) { for (size_t i = 0; i < poly.size() - 1; i++) { f(vector<vector<double>>{poly[i], poly[i + 1]}); } } double getMaxDbl() { return (unsigned int)(-1); } double getMinDbl() { return -(double)(unsigned int)(-1); } double pointLineDistance(double x, double y, double x1, double y1, double x2, double y2) { double A = x - x1; double B = y - y1; double C = x2 - x1; double D = y2 - y1; double dot = A * C + B * D; double len_sq = C * C + D * D; double param = -1; if (len_sq != 0) { param = dot / len_sq; } double xx, yy; if (param < 0) { xx = x1; yy = y1; } else if (param > 1) { xx = x2; yy = y2; } else { xx = x1 + param * C; yy = y1 + param * D; } double dx = x - xx; double dy = y - yy; return sqrt(dx * dx + dy * dy); } vector<double> calculateCenterOfGravityProbabilistic(vector<vector<double>> poly) { // init random srand(time(NULL)); // bounding box double maxx = getMinDbl(), minx = getMaxDbl(), maxy = getMinDbl(), miny = getMaxDbl(); forEachPolyEdge(poly, [&](vector<double> p) { if (p[0] > maxx) { maxx = p[0]; } if (p[0] < minx) { minx = p[0]; } if (p[1] > maxy) { maxy = p[1]; } if (p[1] < miny) { miny = p[1]; } }); // Monte Carlo double x, y; double currentBest = getMinDbl(); for (;;) { double rx = (double)rand() / (double)RAND_MAX * (maxx - minx) + minx; double ry = (double)rand() / (double)RAND_MAX * (maxy - miny) + miny; double rBest = getMaxDbl(); forEachPolySite(poly, [&](vector<vector<double>> ps) { double x1 = ps[0][0]; double x2 = ps[1][0]; double y1 = ps[0][1]; double y2 = ps[1][1]; double dist = pointLineDistance(rx, ry, x1, y1, x2, y2); if (dist < rBest) { rBest = dist; } }); if (rBest > currentBest) { if (rBest - 0.001 <= currentBest) { break; } x = rx; y = ry; currentBest = rBest; cout << " Found: " << x << " " << y << " " << currentBest << "\n"; const double FACTOR = 0.15; double f1 = (maxx - minx) * FACTOR; double f2 = (maxy - miny) * FACTOR; maxx -= f1; minx += f1; maxy -= f2; miny += f2; } a:; } // prepare results vector<double> result; result.push_back(x); result.push_back(y); return result; } vector<double> calculateCenterOfGravityFallback(vector<vector<double>> poly) { double sum1 = 0; double sum2 = 0; double area = 0; forEachPolySite(poly, [&](vector<vector<double>> ps) { double x1 = ps[0][0]; double x2 = ps[1][0]; double y1 = ps[0][1]; double y2 = ps[1][1]; double s1 = (x1 + x2) * (x1 * y2 - x2 * y1); double s2 = (y1 + y2) * (x1 * y2 - x2 * y1); double s3 = x1 * y2 - x2 * y1; sum1 += s1 / 6.0; sum2 += s2 / 6.0; area += s3 * 0.5; }); vector<double> result; result.push_back(sum1 / area); result.push_back(sum2 / area); return result; } int main(int argc, char const *argv[]) { { // Quadrat: vector<vector<double>> poly; poly.push_back(vector<double>{1, 1}); poly.push_back(vector<double>{2, 1}); poly.push_back(vector<double>{2, 2}); poly.push_back(vector<double>{1, 2}); poly.push_back(vector<double>{1, 1}); auto r1 = calculateCenterOfGravityFallback(poly); auto r2 = calculateCenterOfGravityProbabilistic(poly); cout << "x:" << r1[0] << " y:" << r1[1] << "\n"; cout << "x:" << r2[0] << " y:" << r2[1] << "\n"; } { // Konvexes Polygon: vector<std::vector<double>> poly; poly.push_back(vector<double>{3, 1}); poly.push_back(vector<double>{5, 2}); poly.push_back(vector<double>{4, 3.5}); poly.push_back(vector<double>{2, 4}); poly.push_back(vector<double>{1, 2.5}); poly.push_back(vector<double>{3, 1}); auto r1 = calculateCenterOfGravityFallback(poly); auto r2 = calculateCenterOfGravityProbabilistic(poly); cout << "x:" << r1[0] << " y:" << r1[1] << "\n"; cout << "x:" << r2[0] << " y:" << r2[1] << "\n"; } return 0; }Leider hat die Auto-Formatierung von VSCode etwas gezickt. Na ja, aber so ist es "speicherungswert".
