Millionen von 3D-Punkten: Wie finden Sie die 10 von ihnen am nächsten an einem bestimmten Punkt?

Ein Punkt in 3-d ist definiert durch (x, y, z). Der Abstand d zwischen zwei beliebigen Punkten (X, Y, Z) und (x, y, z) ist d = Sqrt [(Xx) ^ 2 + (Yy) ^ 2 + (Zz) ^ 2]. Jetzt gibt es eine Million Einträge in einer Datei, jeder Eintrag ist ein Punkt im Raum, in keiner bestimmten Reihenfolge. Zu jedem Punkt (a, b, c) finden Sie die nächsten 10 Punkte. Wie würden Sie die Millionenpunkte speichern und wie würden Sie diese 10 Punkte aus dieser Datenstruktur abrufen?

Millionen Punkte sind eine kleine Zahl. Der einfachste Ansatz funktioniert hier (Code basierend auf KDTree ist langsamer (um nur einen Punkt abzufragen)).

Brute-Force-Ansatz (Zeit ~ 1 Sekunde)

#!/usr/bin/env python import numpy NDIM = 3 # number of dimensions # read points into array a = numpy.fromfile('million_3D_points.txt', sep=' ') a.shape = a.size / NDIM, NDIM point = numpy.random.uniform(0, 100, NDIM) # choose random point print 'point:', point d = ((a-point)**2).sum(axis=1) # compute distances ndx = d.argsort() # indirect sort # print 10 nearest points to the chosen one import pprint pprint.pprint(zip(a[ndx[:10]], d[ndx[:10]])) 

Starte es:

 $ time python nearest.py point: [ 69.06310224 2.23409409 50.41979143] [(array([ 69., 2., 50.]), 0.23500677815852947), (array([ 69., 2., 51.]), 0.39542392750839772), (array([ 69., 3., 50.]), 0.76681859086988302), (array([ 69., 3., 50.]), 0.76681859086988302), (array([ 69., 3., 51.]), 0.9272357402197513), (array([ 70., 2., 50.]), 1.1088022980015722), (array([ 70., 2., 51.]), 1.2692194473514404), (array([ 70., 2., 51.]), 1.2692194473514404), (array([ 70., 3., 51.]), 1.801031260062794), (array([ 69., 1., 51.]), 1.8636121147970444)] real 0m1.122s user 0m1.010s sys 0m0.120s 

Hier ist das Skript, das Millionen von 3D-Punkten generiert:

 #!/usr/bin/env python import random for _ in xrange(10**6): print ' '.join(str(random.randrange(100)) for _ in range(3)) 

Ausgabe:

 $ head million_3D_points.txt 18 56 26 19 35 74 47 43 71 82 63 28 43 82 0 34 40 16 75 85 69 88 58 3 0 63 90 81 78 98 

Sie könnten diesen Code verwenden, um komplexere Datenstrukturen und Algorithmen zu testen (zum Beispiel, ob sie tatsächlich weniger Speicher verbrauchen oder schneller als der oben genannte Ansatz). Es ist erwähnenswert, dass es im Moment die einzige Antwort ist, die funktionierenden Code enthält.

Lösung basierend auf KDTree (Zeit ~ 1,4 Sekunden)

 #!/usr/bin/env python import numpy NDIM = 3 # number of dimensions # read points into array a = numpy.fromfile('million_3D_points.txt', sep=' ') a.shape = a.size / NDIM, NDIM point = [ 69.06310224, 2.23409409, 50.41979143] # use the same point as above print 'point:', point from scipy.spatial import KDTree # find 10 nearest points tree = KDTree(a, leafsize=a.shape[0]+1) distances, ndx = tree.query([point], k=10) # print 10 nearest points to the chosen one print a[ndx] 

Starte es:

 $ time python nearest_kdtree.py point: [69.063102240000006, 2.2340940900000001, 50.419791429999997] [[[ 69. 2. 50.] [ 69. 2. 51.] [ 69. 3. 50.] [ 69. 3. 50.] [ 69. 3. 51.] [ 70. 2. 50.] [ 70. 2. 51.] [ 70. 2. 51.] [ 70. 3. 51.] [ 69. 1. 51.]]] real 0m1.359s user 0m1.280s sys 0m0.080s 

Teilsortierung in C ++ (Zeit ~ 1,1 Sekunden)

 // $ g++ nearest.cc && (time ./a.out < million_3D_points.txt ) #include  #include  #include  #include  // _1 #include  // bind() #include  namespace { typedef double coord_t; typedef boost::tuple point_t; coord_t distance_sq(const point_t& a, const point_t& b) { // or boost::geometry::distance coord_t x = a.get<0>() - b.get<0>(); coord_t y = a.get<1>() - b.get<1>(); coord_t z = a.get<2>() - b.get<2>(); return x*x + y*y + z*z; } } int main() { using namespace std; using namespace boost::lambda; // _1, _2, bind() // read array from stdin vector points; cin.exceptions(ios::badbit); // throw exception on bad input while(cin) { coord_t x,y,z; cin >> x >> y >> z; points.push_back(boost::make_tuple(x,y,z)); } // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout < < "point: " << point << endl; // 1.14s // find 10 nearest points using partial_sort() // Complexity: O(N)*log(m) comparisons (O(N)*log(N) worst case for the implementation) const size_t m = 10; partial_sort(points.begin(), points.begin() + m, points.end(), bind(less(), // compare by distance to the point bind(distance_sq, _1, point), bind(distance_sq, _2, point))); for_each(points.begin(), points.begin() + m, cout < < _1 << "\n"); // 1.16s } 

Starte es:

 g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) point: (69.0631 2.23409 50.4198) (69 2 50) (69 2 51) (69 3 50) (69 3 50) (69 3 51) (70 2 50) (70 2 51) (70 2 51) (70 3 51) (69 1 51) real 0m1.152s user 0m1.140s sys 0m0.010s 

Prioritätswarteschlange in C ++ (Zeit ~ 1,2 Sekunden)

 #include  // make_heap #include  // binary_function<> #include  #include  // boost::begin(), boost::end() #include  // get<>, tuple<>, cout < < namespace { typedef double coord_t; typedef std::tr1::tuple point_t; // calculate distance (squared) between points `a` & `b` coord_t distance_sq(const point_t& a, const point_t& b) { // boost::geometry::distance() squared using std::tr1::get; coord_t x = get<0>(a) - get<0>(b); coord_t y = get<1>(a) - get<1>(b); coord_t z = get<2>(a) - get<2>(b); return x*x + y*y + z*z; } // read from input stream `in` to the point `point_out` std::istream& getpoint(std::istream& in, point_t& point_out) { using std::tr1::get; return (in >> get<0>(point_out) >> get<1>(point_out) >> get<2>(point_out)); } // Adaptable binary predicate that defines whether the first // argument is nearer than the second one to given reference point template class less_distance : public std::binary_function { const T& point; public: less_distance(const T& reference_point) : point(reference_point) {} bool operator () (const T& a, const T& b) const { return distance_sq(a, point) < distance_sq(b, point); } }; } int main() { using namespace std; // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout << "point: " << point << endl; const size_t nneighbours = 10; // number of nearest neighbours to find point_t points[nneighbours+1]; // populate `points` for (size_t i = 0; getpoint(cin, points[i]) && i < nneighbours; ++i) ; less_distance less_distance_point(point); make_heap (boost::begin(points), boost::end(points), less_distance_point); // Complexity: O(N*log(m)) while(getpoint(cin, points[nneighbours])) { // add points[-1] to the heap; O(log(m)) push_heap(boost::begin(points), boost::end(points), less_distance_point); // remove (move to last position) the most distant from the // `point` point; O(log(m)) pop_heap (boost::begin(points), boost::end(points), less_distance_point); } // print results push_heap (boost::begin(points), boost::end(points), less_distance_point); // O(m*log(m)) sort_heap (boost::begin(points), boost::end(points), less_distance_point); for (size_t i = 0; i < nneighbours; ++i) { cout << points[i] << ' ' << distance_sq(points[i], point) << '\n'; } } 

Starte es:

 $ g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) point: (69.0631 2.23409 50.4198) (69 2 50) 0.235007 (69 2 51) 0.395424 (69 3 50) 0.766819 (69 3 50) 0.766819 (69 3 51) 0.927236 (70 2 50) 1.1088 (70 2 51) 1.26922 (70 2 51) 1.26922 (70 3 51) 1.80103 (69 1 51) 1.86361 real 0m1.174s user 0m1.180s sys 0m0.000s 

Linear Search -basierter Ansatz (Zeit ~ 1.15 Sekunden)

 // $ g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) #include  // sort #include  // binary_function<> #include  #include  #include  // begin(), end() #include  // get<>, tuple<>, cout < < #define foreach BOOST_FOREACH namespace { typedef double coord_t; typedef std::tr1::tuple point_t; // calculate distance (squared) between points `a` & `b` coord_t distance_sq(const point_t& a, const point_t& b); // read from input stream `in` to the point `point_out` std::istream& getpoint(std::istream& in, point_t& point_out); // Adaptable binary predicate that defines whether the first // argument is nearer than the second one to given reference point class less_distance : public std::binary_function { const point_t& point; public: explicit less_distance(const point_t& reference_point) : point(reference_point) {} bool operator () (const point_t& a, const point_t& b) const { return distance_sq(a, point) < distance_sq(b, point); } }; } int main() { using namespace std; // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout << "point: " << point << endl; less_distance nearer(point); const size_t nneighbours = 10; // number of nearest neighbours to find point_t points[nneighbours]; // populate `points` foreach (point_t& p, points) if (! getpoint(cin, p)) break; // Complexity: O(N*m) point_t current_point; while(cin) { getpoint(cin, current_point); //NOTE: `cin` fails after the last //point, so one can't lift it up to //the while condition // move to the last position the most distant from the // `point` point; O(m) foreach (point_t& p, points) if (nearer(current_point, p)) // found point that is nearer to the `point` //NOTE: could use insert (on sorted sequence) & break instead //of swap but in that case it might be better to use //heap-based algorithm altogether std::swap(current_point, p); } // print results; O(m*log(m)) sort(boost::begin(points), boost::end(points), nearer); foreach (point_t p, points) cout << p << ' ' << distance_sq(p, point) << '\n'; } namespace { coord_t distance_sq(const point_t& a, const point_t& b) { // boost::geometry::distance() squared using std::tr1::get; coord_t x = get<0>(a) - get<0>(b); coord_t y = get<1>(a) - get<1>(b); coord_t z = get<2>(a) - get<2>(b); return x*x + y*y + z*z; } std::istream& getpoint(std::istream& in, point_t& point_out) { using std::tr1::get; return (in >> get<0>(point_out) >> get<1>(point_out) >> get<2>(point_out)); } } 

Messungen zeigen, dass die meiste Zeit damit verbracht wird, ein Array aus der Datei zu lesen. Tatsächliche Berechnungen nehmen eine Größenordnung weniger Zeit in Anspruch.

Wenn sich die Millioneneinträge bereits in einer Datei befinden, müssen sie nicht alle in eine Datenstruktur im Speicher geladen werden. Behalte einfach ein Array mit den Top-Ten-Punkten und scanne über die Millionen Punkte und aktualisiere deine Top-Ten-Liste, während du gehst.

Dies ist O (n) in der Anzahl der Punkte.

Sie könnten die Punkte in einem k-dimensionalen Baum (kd-Baum) speichern. Kd-Bäume sind für Nearest-Neighbor-Suchen optimiert (Finden der n Punkte, die einem gegebenen Punkt am nächsten sind).

Ich denke, das ist eine schwierige Frage, die testet, wenn Sie nicht versuchen, Dinge zu übertreiben.

Betrachten Sie den einfachsten Algorithmus, den die Leute bereits oben angegeben haben: Behalten Sie eine Tabelle mit zehn Best-so-far-Kandidaten und gehen Sie alle Punkte eins nach dem anderen durch. Wenn Sie einen Punkt finden, der näher ist als einer der zehn bisher besten, ersetzen Sie ihn. Was ist die Komplexität? Nun, wir müssen jeden Punkt aus der Datei einmal betrachten , seine Entfernung (oder das Quadrat der tatsächlichen Entfernung) berechnen und mit dem 10. nächsten Punkt vergleichen. Wenn es besser ist, fügen Sie es an der entsprechenden Stelle in der 10 besten-so-weit-Tabelle ein.

Also, was ist die Komplexität? Wir betrachten jeden Punkt einmal, also sind es n Berechnungen der Entfernung und n Vergleiche. Wenn der Punkt besser ist, müssen wir ihn an der richtigen Stelle einfügen, dies erfordert einige weitere Vergleiche, aber es ist ein konstanter Faktor, da die Tabelle der besten Kandidaten eine konstante Größe 10 hat.

Wir enden mit einem Algorithmus, der in linearer Zeit läuft, O (n) in der Anzahl der Punkte.

Aber nun überlegen, was ist die untere Grenze für einen solchen Algorithmus? Wenn in den Eingabedaten keine Reihenfolge vorhanden ist, müssen wir jeden Punkt betrachten, um zu sehen, ob er nicht einer der nächsten ist. Soweit ich sehen kann, ist die untere Grenze Omega (n) und somit ist der obige Algorithmus optimal.

Keine Notwendigkeit, die Entfernung zu berechnen. Nur das Quadrat der Entfernung sollte Ihren Bedürfnissen entsprechen. Sollte schneller sein, denke ich. Mit anderen Worten, Sie können das sqrt Bit überspringen.

Das ist keine Hausaufgabenfrage, oder? 😉

Mein Gedanke: iteriere über alle Punkte und lege sie in einen Min-Heap oder eine beschränkte Prioritätswarteschlange, die durch die Entfernung vom Ziel bestimmt sind.

Diese Frage testet im Wesentlichen Ihr Wissen und / oder Ihre Intuition von Algorithmen zur Raumpartitionierung . Ich würde sagen, dass das Speichern der Daten in einem Octree die beste Wahl ist. Es wird häufig in 3D-Engines verwendet, die genau diese Art von Problem behandeln (Speichern von Millionen von Vertices, Raytracing, Auffinden von Kollisionen usw.). Die Nachschlagezeit wird im schlimmsten Fall in der Größenordnung von log(n) (glaube ich).

Einfacher Algorithmus:

Speichern Sie die Punkte als eine Liste von Tupeln und scannen Sie über die Punkte, berechnen Sie die Entfernung und halten Sie eine “nächste” Liste.

Kreativer:

Gruppieren Sie Punkte in Regionen (z. B. den Würfel beschrieben von “0,0,0” bis “50,50,50” oder “0,0,0” bis “-20, -20, -20”), also Sie Sie können vom Zielpunkt aus auf sie “indexieren”. Überprüfen Sie, in welchem ​​Würfel das Ziel liegt und durchsuchen Sie nur die Punkte in diesem Würfel. Wenn es weniger als 10 Punkte in diesem Würfel gibt, überprüfen Sie die “benachbarten” Würfel und so weiter.

Bei näherer Betrachtung ist dies kein sehr guter Algorithmus: Wenn Ihr Zielpunkt näher an der Wand eines Würfels ist als 10 Punkte, dann müssen Sie auch in den benachbarten Würfel suchen.

Ich würde mit dem kd-tree Ansatz gehen und den nächsten finden, dann den nächsten Knoten entfernen (oder markieren) und erneut nach dem nächsten Knoten suchen. Spülen und wiederholen.

Für zwei beliebige Punkte P1 (x1, y1, z1) und P2 (x2, y2, z2) muss, wenn der Abstand zwischen den Punkten d ist, alles Folgende wahr sein:

 |x1 - x2| < = d |y1 - y2| <= d |z1 - z2| <= d 

Halten Sie die 10 am nächsten, wenn Sie über Ihren gesamten Satz iterieren, aber halten Sie auch den Abstand zum 10. nächsten. Sparen Sie sich viel Komplexität, indem Sie diese drei Bedingungen verwenden, bevor Sie die Entfernung für jeden Punkt berechnen, den Sie betrachten.

im Grunde eine Kombination der ersten beiden Antwort über mich. Da sich die Punkte in einer Datei befinden, müssen sie nicht gespeichert werden. Statt eines Arrays oder eines Min-Heaps würde ich einen Max-Heap verwenden, da Sie nur nach Distanzen suchen möchten, die kleiner sind als der zehnte Punkt. Für ein Array müssten Sie jede neu berechnete Entfernung mit allen 10 der von Ihnen beibehaltenen Entfernungen vergleichen. Für einen Min-Heap müssen Sie 3 Vergleiche mit jeder neu berechneten Distanz durchführen. Mit einem Max-Heap führen Sie nur einen Vergleich durch, wenn der neu berechnete Abstand größer als der Root-Knoten ist.

Diese Frage bedarf einer weiteren Definition.

1) Die Entscheidung bezüglich der Algorithmen, die Daten vorindizieren, hängt sehr davon ab, ob Sie die gesamten Daten im Speicher halten können oder nicht.

Mit kdtrees und octrees müssen Sie die Daten nicht im Speicher halten, und die Performance profitiert von dieser Tatsache, nicht nur weil der Speicherbedarf geringer ist, sondern weil Sie nicht die gesamte Datei lesen müssen.

Mit Bruteforce müssen Sie die gesamte Datei lesen und Entfernungen für jeden neuen Punkt, den Sie suchen, neu berechnen.

Trotzdem ist das für dich vielleicht nicht wichtig.

2) Ein weiterer Faktor ist, wie oft Sie nach einem Punkt suchen müssen.

Laut JF Sebastian ist Bruteforce manchmal sogar bei großen Datenmengen schneller, aber beachte, dass seine Benchmarks das Lesen des gesamten Datensatzes von der Platte aus messen (was nicht notwendig ist, sobald kd-tree oder octree irgendwo gebaut und geschrieben wurde) und dass sie nur eine Suche messen.

Berechnen Sie die Entfernung für jede von ihnen, und wählen Sie (1..10, n) in O (n) Zeit. Das wäre der naive Algorithmus, denke ich.