Poisson Disk Sampling

poisson-disk-sampling.cc

#include <cmath>
#include <iostream>
#include <memory>
#include <random>

using namespace std;

#define FLOAT long double

struct Point {
  FLOAT x, y;
  Point() { }
  Point(FLOAT x, FLOAT y) : x(x), y(y) { }
  Point operator -(const Point &other) const {
    return Point(x - other.x, y - other.y);
  }
  FLOAT dist2() {
    return x * x + y * y;
  }
};
ostream &operator <<(ostream &out, const Point &v) {
  return out << '(' << v.x << ", " << v.y << ')';
}

inline size_t imageToGrid(const Point &point, FLOAT cellSize, size_t gridWidth) {
  return (size_t) (point.x / cellSize) + (size_t) (point.y / cellSize) * gridWidth;
}

// References: http://devmag.org.za/2009/05/03/poisson-disk-sampling/
int main() {
  FLOAT width, height, minDist;
  size_t newPointsCount;
  cin >> width >> height >> minDist >> newPointsCount;
  auto cellSize = minDist / M_SQRT2l, minDist2 = minDist * minDist;
  auto gridWidth = (size_t) ceil(width / cellSize), gridHeight = (size_t) ceil(height / cellSize);
  vector<shared_ptr<Point>> grid(gridWidth * gridHeight, nullptr), processList;
  mt19937 rng;
  auto firstPoint = make_shared<Point>(uniform_real_distribution<FLOAT>(0, width)(rng),
                                       uniform_real_distribution<FLOAT>(0, height)(rng));
  processList.push_back(firstPoint);
  cout << *firstPoint << endl;
  grid[imageToGrid(*firstPoint, cellSize, gridWidth)] = firstPoint;
  uniform_real_distribution<FLOAT> radiusDist(minDist, minDist * 2), angleDist(0, M_PIl * 2);
  while (processList.size()) {
    swap(processList[uniform_int_distribution<size_t>(0, processList.size() - 1)(rng)], processList.back());
    auto point = processList.back();
    processList.pop_back();
    for (size_t i = 0; i < newPointsCount; ++i) {
      FLOAT radius = radiusDist(rng), angle = angleDist(rng);
      auto newPoint = make_shared<Point>(point->x + radius * cos(angle), point->y + radius * sin(angle));
      if (newPoint->x < 0 || newPoint->y < 0 || newPoint->x >= width || newPoint->y >= height) continue;
      auto x = (size_t) (newPoint->x / cellSize), y = (size_t) (newPoint->y / cellSize);
      for (long j = -2; j <= 2; ++j) {
        long x0 = x + j;
        if (x0 < 0) continue;
        if (x0 >= gridWidth) break;
        for (long k = -2; k <= 2; ++k) {
          long y0 = y + k;
          if (y0 < 0) continue;
          if (y0 >= gridHeight) break;
          auto cell = grid[x0 + y0 * gridWidth];
          if (cell && (*cell - *newPoint).dist2() <= minDist2) goto next;
        }
      }
      processList.push_back(newPoint);
      cout << *newPoint << endl;
      grid[imageToGrid(*newPoint, cellSize, gridWidth)] = newPoint;
      next:;
    }
  }
}