Miljoenen 3D-punten: hoe vind je de 10 die het dichtst bij een bepaald punt liggen?

Een punt in 3-d wordt gedefinieerd door (x,y,z). Afstand d tussen twee willekeurige punten (X,Y,Z) en (x,y,z) is d= Sqrt[(X-x)^2 + (Y-y)^2 + (Z-z)^2].
Nu zijn er een miljoen items in een bestand, elk item is een punt in de ruimte, in geen specifieke volgorde. Gegeven een willekeurig punt (a,b,c) vind je de dichtstbijzijnde 10 punten. Hoe zou je de miljoen punten opslaan en hoe zou je die 10 punten uit die datastructuur halen.


Antwoord 1, autoriteit 100%

Miljoen punten is een klein aantal. De meest rechttoe rechtaan aanpak werkt hier (code gebaseerd op KDTree is langzamer (voor het opvragen van slechts één punt)).

Brute-force nadering (tijd ~1 seconde)

#!/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]]))

Voer het uit:

$ 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

Dit is het script dat miljoen 3D-punten genereert:

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

Uitvoer:

$ 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

Je zou die code kunnen gebruiken om complexere datastructuren en algoritmen te testen (bijvoorbeeld of ze daadwerkelijk minder geheugen of sneller verbruiken dan de eenvoudigste methode hierboven). Het is vermeldenswaard dat dit op dit moment het enige antwoord is dat werkende code bevat.

Oplossing gebaseerd op KDTree(tijd ~1,4 seconden)

#!/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]

Voer het uit:

$ 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

Gedeeltelijke sortering in C++ (tijd ~1.1 seconden)

// $ g++ nearest.cc && (time ./a.out < million_3D_points.txt )
#include <algorithm>
#include <iostream>
#include <vector>
#include <boost/lambda/lambda.hpp>  // _1
#include <boost/lambda/bind.hpp>    // bind()
#include <boost/tuple/tuple_io.hpp>
namespace {
  typedef double coord_t;
  typedef boost::tuple<coord_t,coord_t,coord_t> 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<point_t> 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<coord_t>(), // 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
}

Voer het uit:

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

Prioriteitswachtrij in C++ (tijd ~ 1,2 seconden)

#include <algorithm>           // make_heap
#include <functional>          // binary_function<>
#include <iostream>
#include <boost/range.hpp>     // boost::begin(), boost::end()
#include <boost/tr1/tuple.hpp> // get<>, tuple<>, cout <<
namespace {
  typedef double coord_t;
  typedef std::tr1::tuple<coord_t,coord_t,coord_t> 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 T>
  class less_distance : public std::binary_function<T, T, bool> {
    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<point_t> 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';  
  }
}

Voer het uit:

$ 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

Op lineaire zoeken gebaseerde benadering (tijd ~ 1,15 seconden)

// $ g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt )
#include <algorithm>           // sort
#include <functional>          // binary_function<>
#include <iostream>
#include <boost/foreach.hpp>
#include <boost/range.hpp>     // begin(), end()
#include <boost/tr1/tuple.hpp> // get<>, tuple<>, cout <<
#define foreach BOOST_FOREACH
namespace {
  typedef double coord_t;
  typedef std::tr1::tuple<coord_t,coord_t,coord_t> 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<point_t, point_t, bool> {
    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));
  }
}

Metingen laten zien dat de meeste tijd wordt besteed aan het lezen van de array uit het bestand, de werkelijke berekeningen nemen in orde van grootte minder tijd in beslag.


Antwoord 2, autoriteit 21%

Als de miljoen vermeldingen al in een bestand zitten, is het niet nodig om ze allemaal in een gegevensstructuur in het geheugen te laden. Houd gewoon een array bij met de top tien tot nu toe gevonden punten en scan meer dan de miljoen punten, terwijl je je top tien-lijst bijwerkt terwijl je bezig bent.

Dit is O(n) in het aantal punten.


Antwoord 3, autoriteit 15%

Je zou de punten kunnen opslaan in een k-dimensionale boom(kd-tree) . Kd-trees zijn geoptimaliseerd voor zoekopdrachten naar de dichtstbijzijnde buren (de npunten vinden die zich het dichtst bij een bepaald punt bevinden).


Antwoord 4, autoriteit 10%

Ik denk dat dit een lastige vraag is die test of je niet probeert te overdrijven.

Beschouw het eenvoudigste algoritme dat mensen hierboven al hebben gegeven: houd een tabel bij van de tien beste kandidaten tot nu toe en doorloop alle punten één voor één. Als je een punt dichterbij vindt dan een van de tien beste tot nu toe, vervang het dan. Wat is de complexiteit? Welnu, we moeten elk punt uit het bestand eenmaalbekijken, de afstand berekenen (of het kwadraat van de afstand eigenlijk) en vergelijken met het 10e dichtstbijzijnde punt. Als het beter is, plaatst u het op de juiste plaats in de tabel met 10 beste tot nu toe.

Dus wat is de complexiteit? We kijken één keer naar elk punt, dus het zijn n berekeningen van de afstand en n vergelijkingen. Als het punt beter is, moeten we het op de juiste positie invoegen, dit vereist wat meer vergelijkingen, maar het is een constante factor aangezien de tabel met beste kandidaten een constante grootte van 10 heeft.

We eindigen met een algoritme dat in lineaire tijd loopt, O(n) in het aantal punten.

Maar bedenk nu wat de ondergrensis voor zo’n algoritme? Als er geen volgorde is in de invoergegevens, moeten wenaar elk punt kijken om te zien of het niet een van de dichtstbijzijnde is. Dus voor zover ik kan zien, is de ondergrens Omega(n) en dus is het bovenstaande algoritme optimaal.


Antwoord 5, autoriteit 6%

Het is niet nodig om de afstand te berekenen. Alleen het kwadraat van de afstand zou aan uw behoeften moeten voldoen. Zou sneller moeten denk ik. Met andere woorden, u kunt het sqrtbit overslaan.


Antwoord 6, autoriteit 4%

Dit is toch geen huiswerkvraag? 😉

Mijn gedachte: herhaal alle punten en plaats ze in een min heap of begrensde prioriteitswachtrij, gecodeerd op afstand tot het doel.


Antwoord 7, autoriteit 4%

Deze vraag test in wezen uw kennis en/of intuïtie van algoritmen voor het partitioneren van ruimten. Ik zou zeggen dat het opslaan van de gegevens in een octreede beste keuze is. Het wordt vaak gebruikt in 3D-engines die precies dit soort problemen oplossen (miljoenen hoekpunten opslaan, raytracing, botsingen vinden, enz.). De opzoektijd zal in het ergste geval (denk ik) in de orde van log(n)zijn.


Antwoord 8, autoriteit 2%

Eenvoudig algoritme:

Sla de punten op als een lijst met tuples en scan over de punten, bereken de afstand en houd een ‘dichtstbijzijnde’ lijst bij.

Creatiever:

Groepeer punten in regio’s (zoals de kubus beschreven door “0,0,0” tot “50,50,50” of “0,0,0” tot “-20,-20,-20”) , zodat u er vanaf het richtpunt in kunt “indexeren”. Controleer in welke kubus het doel ligt en zoek alleen door de punten in die kubus. Als er minder dan 10 punten in die kubus zijn, controleer dan de “naburige” kubussen, enzovoort.

Bij nader inzien is dit geen erg goed algoritme: als je richtpunt dichter bij de muur van een kubus ligt dan 10 punten, dan moet je ook in de aangrenzende kubus zoeken.

Ik zou voor de kd-tree-benadering gaan en de dichtstbijzijnde vinden, dan die dichtstbijzijnde knoop verwijderen (of markeren) en opnieuw zoeken naar de nieuwe dichtstbijzijnde knoop. Spoel en herhaal.


Antwoord 9, autoriteit 2%

Voor elke twee punten P1 (x1, y1, z1) en P2 (x2, y2, z2), als de afstand tussen de punten d is, moet al het volgende waar zijn:

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

Houd de dichtstbijzijnde 10 vast terwijl je je hele set herhaalt, maar houd ook de afstand tot de dichtstbijzijnde 10 vast. Bespaar uzelf een hoop complexiteit door deze drie voorwaarden te gebruiken voordat u de afstand berekent voor elk punt waar u naar kijkt.


Antwoord 10

eigenlijk een combinatie van de eerste twee antwoorden boven mij. aangezien de punten in een bestand staan, is het niet nodig om ze in het geheugen te bewaren. In plaats van een array of een min heap zou ik een max heap gebruiken, omdat je alleen wilt controleren op afstanden die kleiner zijn dan het 10e dichtstbijzijnde punt. Voor een array zou je elke nieuw berekende afstand moeten vergelijken met alle 10 afstanden die je hebt aangehouden. Voor een min heap moet je 3 vergelijkingen maken met elke nieuw berekende afstand. Met een max heap voer je maar 1 vergelijking uit als de nieuw berekende afstand groter is dan de root node.


Antwoord 11

Bereken de afstand voor elk van hen en voer Select(1..10, n) in O(n) tijd uit. Dat zou het naïeve algoritme zijn, denk ik.


Antwoord 12

Deze vraag moet verder worden gedefinieerd.

1)
De beslissing met betrekking tot de algoritmen die gegevens vooraf indexeren, varieert sterk, afhankelijk van of u de volledige gegevens in het geheugen kunt bewaren of niet.

Met kdtrees en octrees hoeft u de gegevens niet in het geheugen te bewaren en de prestaties profiteren daarvan, niet alleen omdat de geheugenvoetafdruk kleiner is, maar simpelweg omdat u niet het hele bestand hoeft te lezen.

Met bruteforce moet je het hele bestand lezen en de afstanden opnieuw berekenen voor elk nieuw punt waarnaar je zoekt.

Toch is dit misschien niet belangrijk voor je.

2)
Een andere factor is hoe vaak je naar een punt moet zoeken.

Zoals aangegeven door JF Sebastian is bruteforce soms zelfs sneller op grote datasets, maar houd er rekening mee dat zijn benchmarks het lezen van de hele dataset van schijf meten (wat niet nodig is als kd-tree of octree eenmaal is gebouwd en ergens geschreven) en dat ze slechts één zoekopdracht meten.

Other episodes