#ifndef __POISSON_DISK_SAMPLING__
#define __POISSON_DISK_SAMPLING__
// ========================================================================== //
// Forward Declare
class BBox;

// Local libraries
#include "Matrix/Matrix.hpp"

// External libraries
#include <GEL/CGLA/CGLA.h>

// Standard libraries
#include <map>
#include <random>
#include <vector>

// ========================================================================== //
// Start of the class

/** Poisson Disk Sampling - Bridson 2007.
 *
 * This Poisson Disk Sampling code was written based on the algorithm
 * developed by Robert Bridson in 2007, as presented in his article "Fast
 * Poisson Disk Sampling in Arbitrary Dimensions"
 *
 * This class can be used to construct a full sampling, or used to maintain the
 * underlying lookup grid.
 */
class PDS {

  private:
    // ====================================================================== //
    // Data Fields

    double                   pds_r;        // PDS radius.
    int                      search_size;  // Search distance for neighbours.
    double                   bbox_diag;    // Diagonal size of the bbox.
    double                   element_size; // Size of each individual element.
    Matrix<double>           bbox;         // Bounding box for the domain.
    std::vector<int>         pds_dim;      // Size of each dimension of the grid
    std::vector<CGLA::Vec3d> points;       // List of points
    std::map<std::vector<int>, int> grid;  // Lookup grid for the sampling

  public:
    // ====================================================================== //
    // Constructors

    /** Default constructor. */
    PDS();
    /** Construct a new PDS object with a given sampling radius. */
    PDS(double radius, BBox &bb);
    /** Construct a new PDS object with a given sampling radius. */
    PDS(double radius, Matrix<double> bbox);
    /** Copy constructor. */
    PDS(const PDS &original);
    /** Move constructor*/
    PDS(PDS &&original);
    /** Destroy the PDS object. */
    ~PDS();

    // ====================================================================== //
    // Operator overloads

    /** Copy assignment operator. */
    PDS &operator=(PDS &original);
    /** Move assignment operator. */
    PDS &operator=(PDS &&original);

    // ====================================================================== //
    // Global operations

    /** Resetting the pds grid.
     *
     * This will reset the sampling grid and the list of points. It will retain
     * the sampling radius, bounding box etc.
     */
    void reset();

    /** Generate full sampling. */
    void generate();
    /** Generate full sampling from origin. */
    void generate(CGLA::Vec3d p0);
    /** Generate full sampling from origin with specified random device. */
    void generate(CGLA::Vec3d p0, std::mt19937 random_device);

    /** Get the radius used for the grid. */
    double get_radius();

    /** Get list of all points. */
    std::vector<CGLA::Vec3d> get_points();

    // ====================================================================== //
    // Local operations

    /** Add point to the pds grid. */
    void add(CGLA::Vec3d p);
    /** Add point to the pds grid with specified index. */
    void add(int idx, CGLA::Vec3d p);
    /** Add point to the pds grid with specified index. */
    void add(int idx, std::vector<double> p);

    /** Remove point to the pds grid. */
    void remove(CGLA::Vec3d p);

    /** Get index stored at position. */
    int get(CGLA::Vec3d p);
    /** Get index stored at position. */
    int get(std::vector<double> p);

    /** Get a list of all neighbouring points. */
    std::vector<int> neighbours(CGLA::Vec3d p);
    /** Get a list of all neighbouring points. */
    std::vector<int> neighbours(std::vector<double> p);
    /** Get a list of all neighbouring points within radius r. */
    std::vector<int> neighbours(CGLA::Vec3d p, double r);

  private:
    // ====================================================================== //
    // Field Access

    /** Compute the grid index for a given point. */
    std::vector<int> index(CGLA::Vec3d p);
    /** Compute the gird index for a given point. */
    std::vector<int> index(std::vector<double> p);

    /** Generate neighbour list from grid index. */
    std::vector<int> neighbours(std::vector<int> p_idx, int r);

    // ====================================================================== //
    // Boolean operations

    /** Check if point is inside the bounding box. */
    bool inside_bbox(CGLA::Vec3d p);
};

#endif /* __STREAM_SURFACES__ */