Fast Poisson Disk Sampling
Coarse C++ implementation for "Fast Poisson disk sampling in arbitrary dimensions, R. Bridson, ACM SIGGRAPH 2007 Sketches Program".
在公司实在闲的蛋疼，于是做了个快速泊松碟采样的实现。
2D FPDS Code/**
* Copyright (c) Bo Zhou<bo.schwarzstein@gmail.com>
* An 2D implementation of "Fast Poisson Disk Sampling in Arbitrary Dimensions, R. Bridson, ACM SIGGRAPH 2007 Sketches Program".
* Anybody could use this code freely, if you feel it good or helpful, please tell me, thank you very much.
*/
#ifndef FPDS_HPP
#define FPDS_HPP
#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#include <math.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <math.h>
#include <boost/multi_array.hpp>
// Define our basis data structure.
typedef std::pair<int,int> Index2D;
typedef std::pair<float,float> Point2D;
// Use C standard rand() to get uniform random number in (0,1)
inline float URand()
{
return (float)rand() / (float)RAND_MAX;
}
// Calculate distance between 2 points.
inline float LengthSqrt(const Point2D& X, const Point2D& Y)
{
float a = X.first - Y.first;
float b = X.second - Y.second;
return sqrtf( a*a + b*b );
}
float FastPoissonDiskSampling2D(const float r, std::vector<Point2D>& PointList)
{
PointList.clear();
const int N = 2;
const float SqrtN = sqrtf((float)N);
const float R = r*SqrtN; // This R should be used to draw circle.
const float RDivSqrtN = R / SqrtN;
const int K = 30; // Paper suggested.
const float CellSize = (float)N*RDivSqrtN;
const int GridSize = (int)ceilf(1.0f/CellSize);
const int INVALID_IDX = -1;
// Allocate a grid to accelerate.
// Each cell could have one and only one sampler.
boost::multi_array<int,N> Grid( boost::extents[GridSize][GridSize] );
memset( Grid.data(), INVALID_IDX, sizeof(int)*GridSize*GridSize );
std::vector<int> ActiveList;
Point2D P0;
Index2D I0;
P0.first = URand();
P0.second = URand();
I0.first = (int)floorf( P0.first*GridSize );
I0.second = (int)floorf( P0.second*GridSize );
ActiveList.push_back(0);
PointList.push_back(P0);
Grid[I0.first][I0.second] = 0;
int Done = 1;
while( ActiveList.size() != 0 )
{
// Initialize a sampler.
size_t MagicIdx = (size_t)floorf( URand() * ActiveList.size() );
size_t StartIdx = ActiveList[MagicIdx];
Point2D StartPoint = PointList[ StartIdx ];
bool Found = false;
for( size_t i=0; i<K; ++i )
{
// Generate point in radius (R,2R)
float t = URand()*M_PI*2.0;
float r = URand()*R+R;
float X = r*cosf( t );
float Y = r*sinf( t );
// Move to current center
Point2D CurrentPoint;
CurrentPoint.first = X+StartPoint.first;
CurrentPoint.second = Y+StartPoint.second;
// Discard if out of domain
if( CurrentPoint.first < 0.0f || CurrentPoint.first >= 1.0f )
continue;
if( CurrentPoint.second < 0.0f || CurrentPoint.second >= 1.0f )
continue;
// Which cell the test point located in
Index2D TargetCell;
TargetCell.first = (int)floorf( CurrentPoint.first / CellSize );
TargetCell.second = (int)floorf( CurrentPoint.second / CellSize );
if( TargetCell.first < 0 || TargetCell.first >= GridSize )
continue;
if( TargetCell.second < 0 || TargetCell.second >= GridSize )
continue;
if( Grid[TargetCell.first][TargetCell.second] != INVALID_IDX )
continue;
Index2D TCLeft( TargetCell );
Index2D TCRight( TargetCell );
Index2D TCDown( TargetCell );
Index2D TCUp( TargetCell );
Index2D TCLU( TargetCell );
Index2D TCRU( TargetCell );
Index2D TCLD( TargetCell );
Index2D TCRD( TargetCell );
bool A = false, B = false, C = false, D = false;
bool E = false, F = false, G = false, H = false;
TCLeft.first--;
if( TCLeft.first > -1 )
{
int Idx2 = Grid[TCLeft.first][TCLeft.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
A = true;
}
}else
{
A = true;
}
}else
{
A = true;
}
TCRight.first++;
if( TCRight.first < GridSize )
{
int Idx2 = Grid[TCRight.first][TCRight.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
B = true;
}
}else
{
B = true;
}
}else
{
B = true;
}
TCDown.second--;
if( TCDown.second > -1 )
{
int Idx2 = Grid[TCDown.first][TCDown.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
C = true;
}
}else
{
C = true;
}
}else
{
C = true;
}
TCUp.second++;
if( TCUp.second < GridSize )
{
int Idx2 = Grid[TCUp.first][TCUp.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
D = true;
}
}else
{
D = true;
}
}else
{
D = true;
}
// 4 Corner
// Left Up
TCLU.first--;TCLU.second++;
if( TCLU.first > -1 && TCLU.second < GridSize )
{
int Idx2 = Grid[TCLU.first][TCLU.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
E = true;
}
}else
{
E = true;
}
}else
{
E = true;
}
// Right Up
TCRU.first++;TCRU.second++;
if( TCRU.first < GridSize && TCRU.second < GridSize )
{
int Idx2 = Grid[TCRU.first][TCRU.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
F = true;
}
}else
{
F = true;
}
}else
{
F = true;
}
// Left Bottom
TCLD.first--;TCLD.second--;
if( TCLD.first > -1 && TCLD.second > -1 )
{
int Idx2 = Grid[TCLD.first][TCLD.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
G = true;
}
}else
{
G = true;
}
}else
{
G = true;
}
// Right Bottom
TCRD.first++;TCRD.second--;
if( TCRD.first < GridSize && TCRD.second > -1 )
{
int Idx2 = Grid[TCRD.first][TCRD.second];
if( Idx2 >= 0 )
{
if( LengthSqrt( PointList[Idx2],CurrentPoint ) > R )
{
H = true;
}
}else
{
H = true;
}
}else
{
H = true;
}
if( A&B&C&D&E&F&G&H )
{
Grid[TargetCell.first][TargetCell.second] = Done;
PointList.push_back(CurrentPoint);
ActiveList.push_back(Done);
++Done;
Found = true;
break;
}
}
// We have to remove this test sampler from active list.
if( Found == false && ActiveList.size() > 0 )
ActiveList.erase( ActiveList.begin()+MagicIdx );
}
return R;
}
#endif
