在C++中实现最近点分治算法

Question

我正在尝试实施分而治之的最近点算法。尽其所能，但我的脑袋快要爆炸了，因为我的代码似乎（随机地）给出了错误的答案。我使用 stl 编写了一个随机数生成器用于测试目的，我不断遇到的错误是，每隔几次运行算法 returns 一对明显比最近的一对更远（相隔 1 个距离单位，这我手动输入）。

请原谅全局变量，但这是我第 3 次重写它，我只是觉得它更容易使用。 Pastebin link 对于那些喜欢在屏幕上看到更多内容的人：http://pastebin.com/93dtj81z

[编辑] 不正确的值似乎来自 BruteCP 函数...我认为...这让我很头疼...

#include <vector>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <ctime>
#include <random>

using namespace std;
using point = pair<int, int>;
double MAX = 1000000000.0;
double GLOBAL_minDist = MAX;
pair<point, point> GLOBAL_nearestPoints;

bool Xcmp ( const point &c1, const point &c2 ) {
  return ( c1.first < c2.first );
}

bool Ycmp ( const point &c1, const point &c2 ) {
  return ( c1.second < c2.second );
}

inline ostream& operator<< ( ostream& os, const point& p ) {
  return os << p.first << " " << p.second << "\n";
}

inline ostream& operator<< ( ostream& os, const vector<point> &points ) {
  for( auto itr = points.begin(); itr < points.end(); itr++ ) {
    os << *itr;
  }
  return os;
}

inline ostream& operator<< ( ostream& os, const pair<point, point> nearestPair ) {
  return os << static_cast<int> (nearestPair.second.first) << " " << static_cast<int> (nearestPair.second.second) << "\n"
            << static_cast<int> (nearestPair.first.first) << " " << static_cast<int> (nearestPair.first.second);
}

inline double distance( const point a, const point b ) {
  return sqrt( pow(( a.first - b.first ), 2 ) + pow( a.second - b.second, 2 ));
}

void bruteCP( const vector<point> &Xs ) {
  for( auto it = Xs.begin(); it < Xs.end() - 1; it++ ) {
    for( auto it2 = it + 1; it2 < Xs.end(); it2++ ) {
      double minDist = distance( *it, *it2 );
      if( minDist < GLOBAL_minDist ) {
        cout << minDist << "\n";
        GLOBAL_minDist = minDist;
        GLOBAL_nearestPoints = pair<point, point> ( *it, *it2 );
      }
    }
  }
}

void divConCP( const vector<point>& Xs, const vector<point>& Ys ) {
  int Xsize = Xs.size();
  if( Xsize <= 3 ) { bruteCP( Xs ); return; }

  int mid =  Xsize / 2;
  int median = Xs[mid].first;

  vector<point> leftYs;
  copy_if( Ys.begin(), Ys.end(), back_inserter(leftYs), [median](const point& point) 
          {return point.first <= median;} );
  vector<point>leftXs;
  leftXs.insert( leftXs.end(), Xs.begin(), Xs.begin() + mid );
  divConCP( leftXs, leftYs );

  vector<point> rightYs, rightXs;
  copy_if( Ys.begin(), Ys.end(), back_inserter(leftYs), [median](const point& point) 
          {return point.first > median;} );
  rightXs.insert( rightXs.end(), Xs.begin() + mid, Xs.end() );
  divConCP( rightXs, rightYs );

  vector<point> strip;
  copy_if( Ys.begin(), Ys.end(), back_inserter(strip), [median, GLOBAL_minDist](const point& point) 
          {return abs(point.first - median) < GLOBAL_minDist;} );

  //vector<point> uniqStrip;
  //unique_copy( strip.begin(), strip.end(), uniqStrip.begin() );

  for( auto itr = strip.begin(); itr < strip.end(); itr++ ) {
    for( auto itr2 = itr + 1; itr2 < strip.end() && (*itr2).second < GLOBAL_minDist; itr2++ ) {
      double minDist = distance( *itr, *itr2 );
      if( minDist < GLOBAL_minDist) {
        //cout << minDist << "\n";
        //cout << *itr << " " << *itr2 << "\n";
        GLOBAL_minDist = minDist;
        GLOBAL_nearestPoints = pair<point, point> ( *itr, *itr2 );
      }
    }
  }
}

int main() {
  int n, x, y;
  vector<point> Xs, Ys;
/*
  cin >> n;
  for( int i = 0; i < n; i++ ) {
    cin >> x;
    cin >> y;
   // x = -i;
   // y = -i;

    point xy( x, y );
    Xs.push_back( xy );
    Ys.push_back( xy );
  }
*/
    // DEBUG //

  n = 100000;
  srand(time(0));
  std::default_random_engine gen(time(0));
  std::uniform_int_distribution<int> dis(-20000000, 20000000);
  for( int i = 0; i < n - 2; i++ ) {
    x = dis( gen );
    y = dis( gen );
    //x = i;
    //y = i;
    point xy( x, y );
    Xs.push_back( xy );
    Ys.push_back( xy );
  }

    Xs.push_back( point( 20001, 20001 ));
    Ys.push_back( point( 20001, 20001 ));
    Xs.push_back( point( 20000, 20001 ));
    Ys.push_back( point( 20000, 20001 ));

  // DEBUG //

  sort( Xs.begin(), Xs.end(), Xcmp );
  sort( Ys.begin(), Ys.end(), Ycmp );

  divConCP( Xs, Ys );
  //bruteCP( Xs );
  cout << GLOBAL_minDist << "\n";
  cout << GLOBAL_nearestPoints << "\n";

}

Answer 1

我已经通过重写程序解决了这个问题。这次我不是创建按 y 坐标排序的向量，将其切碎并将相关部分向下传递到递归调用中。相反，我从一个空向量开始，它一直向下传递到递归的基本情况，在那里它被填充、排序并向上传递，在那里它与其他也以这种方式填充的向量合并。然后，从这个向量中选择相关点来制作我们需要测试的条带，以寻找位于给定中位数的最近点。

最后一个问题，虽然 - 我一直在 std::merge 中遇到我无法理解的分段错误，所以我编写了自己的合并函数，它运行良好......关于这些错误在哪里的任何想法来自（哪里？要重现它们，您只需替换：

Ys = mergee(leftYs, rightYs);

//merge(leftYs.begin(), leftYs.end(), rightYs.begin(), rightYs.end(), Ys.begin());

与：

//Ys = mergee(leftYs, rightYs);

merge(leftYs.begin(), leftYs.end(), rightYs.begin(), rightYs.end(), Ys.begin());

代码，带有用于测试的随机数据生成器： http://pastebin.com/yQCNh2J9

如果投反对票的人至少留下评论解释原因，我将不胜感激。

#include <vector>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <ctime>
#include <random>

using namespace std;
using point = pair<int, int>;

bool Xcmp ( const point &c1, const point &c2 ) {
  return ( c1.first < c2.first );
}

bool Ycmp ( const point &c1, const point &c2 ) {
  return ( c1.second < c2.second );
}

inline ostream& operator<< ( ostream& os, const point& p ) {
  return os << p.first << " " << p.second << "\n";
}

inline ostream& operator<< ( ostream& os, const vector<point> &points ) {
  for( auto itr = points.begin(); itr < points.end(); itr++ ) {
    os << *itr;
  }
  return os;
}

inline ostream& operator<< ( ostream& os, const pair<point, point> nearestPair ) {
  return os << nearestPair.second.first << " " << nearestPair.second.second << "\n"
            << nearestPair.first.first << " " <<  nearestPair.first.second;
}

inline double distance( const point a, const point b ) {
  return sqrt( pow(( a.first - b.first ), 2 ) + pow( a.second - b.second, 2 ));
}

pair<pair<point, point>, double> bruteCP( const vector<point> &Xs, const int l, const int r ) {
  pair<point, point> nearestPair;
  double tempDist = 1000000000.0;
  for( int i = l; i < r; i++ ) {
    for( int j = i + 1; j <= r; j++ ) {
      double minDist = distance( Xs[i], Xs[j] );
      if( minDist < tempDist ) {
        tempDist = minDist;
        nearestPair = pair<point, point> ( Xs[i], Xs[j] );
      }
    }
  }
  return pair<pair<point, point>, double> (nearestPair, tempDist);
}

vector<point> makeStrip( const vector<point> Ys, int median, double minDist ) {
  vector<point> strip;
  for( auto it = Ys.begin(); it < Ys.end(); it++ )
    if( abs((*it).first - median) < minDist ) {
      strip.push_back( *it );
    }
    return strip;
}

vector<point> mergee(const vector<point> left, const vector<point> right) {
  vector<point> result;

  auto leftIterator = left.begin();
  auto rightIterator = right.begin();

  while(leftIterator != left.end() && rightIterator != right.end()) {
    if((*leftIterator).second < (*rightIterator).second) {
        result.push_back(*leftIterator);
        ++leftIterator;
    }
    else {
        result.push_back(*rightIterator);
        ++rightIterator;
    }
  }
  while(leftIterator != left.end()) {
    result.push_back(*leftIterator);
    ++leftIterator;
  }
  while(rightIterator != right.end()) {
    result.push_back(*rightIterator);
    ++rightIterator;
  }
  return result;
}

pair<pair<point, point>, double>  divConCP( const vector<point> &Xs, vector<point> &Ys, const int l, const int r ) {
  int Xsize = r - l + 1;
  if( Xsize <= 3 ) {
    for( int i = l; i <= r; i++ ) {
      Ys.push_back( Xs[i] );
    }
    sort( Ys.begin(), Ys.end(), Ycmp );
    return bruteCP( Xs, l, r );
  }

  int mid =  l - 1 + (r - l + 1) / 2;
  int median = Xs[mid].first;

  vector<point> leftYs;
  vector<point> rightYs;

  auto leftClosest = divConCP( Xs, leftYs, l, mid );
  auto rightClosest = divConCP( Xs, rightYs, mid + 1, r );

  auto lrClosest = leftClosest.second < rightClosest.second ? leftClosest : rightClosest;

  Ys = mergee(leftYs, rightYs);

  //merge(leftYs.begin(), leftYs.end(), rightYs.begin(), rightYs.end(), Ys.begin());

  auto strip = makeStrip( Ys, median, lrClosest.second );

  for( auto itr = strip.begin(); itr < strip.end(); itr++ ) {
    for( auto itr2 = itr + 1; itr2 < strip.end() && itr2 < itr + 8; itr2++ ) {
      double tempDist = distance( *itr, *itr2 );
      if( tempDist < lrClosest.second ) {
        //cout << minDist << "\n";
        //cout << *itr << " " << *itr2 << "\n";
        lrClosest.second = tempDist;
        lrClosest.first = pair<point, point> ( *itr, *itr2 );
      }
    }
  }
  return lrClosest;
}

int main() {
  int n, x, y, MAX = 1000000000;
  vector<point> Xs, Ys;
  pair<int, int> minPair( -MAX, -MAX );
  Xs.push_back(minPair);
  //Ys.push_back(minPair);

  cin >> n;
  for( int i = 1; i <= n; i++ ) {
    cin >> x;
    cin >> y;
   // x = -i;
   // y = -i;

    point xy( x, y );
    Xs.push_back( xy );
    //Ys.push_back( xy );
  }

    // DEBUG //
/*
   n = 1000000;
   srand(time(0));
   std::default_random_engine gen(time(0));
   std::uniform_int_distribution<int> dis(-20000000, 20000000);
   for( int i = 0; i < n - 2; i++ ) {
     x = dis( gen );
     y = dis( gen );
     //x = i;
     //y = i;
     point xy( x, y );
     Xs.push_back( xy );
     //Ys.push_back( xy );
   }

     Xs.push_back( point( 20001, 20001 ));
     //Ys.push_back( point( 20001, 20001 ));
     Xs.push_back( point( 20000, 20001 ));
     //Ys.push_back( point( 20000, 20001 ));
*/
  // DEBUG //

  sort( Xs.begin(), Xs.end(), Xcmp );
  //sort( Ys.begin(), Ys.end(), Ycmp );

  auto p = divConCP( Xs, Ys, 1, n );
  cout << p.first;
  //bruteCP( Xs );
  return 0;
}

在C++中实现最近点分治算法

Implementing closest points divide and conquer algorithm in C++

c++

c++11

closest-points