Poisson disk sample

Test in houdini

float r =1;
int k = 30;
int dim = 2;
int width = 200;
int height = 100;
// fill grid with -1
int grid[];        // -1 is no sample
vector grid_pts[]; // -1 is default set(0,0,0)

int active[];    // an array of sample indices
vector active_pts[]; 

vector process_pts[];


// cell size of per voxel
float w = r / sqrt(dim);

int cols = floor(width  / w);
int rows = floor(height / w);


for(int i= 0;i<cols * rows ; i++){
    push(grid, -1);
    push(grid_pts, set(0,0,0));
}
printf("cellsize: %d, cols: %d, rows: %d
",w, cols, rows);

// step1 
// Select the initial sample,x0, randomly chosen uniformlyfrom the domain
float init_x = width / 2.0f;
float init_y = height / 2.0f;
vector init_pos = set(init_x, init_y, 0 ); // init point from center!
int init_i = floor(init_x / w);
int init_j = floor(init_y / w);

grid_pts[init_i + init_j *cols] = init_pos;
grid[ init_i + init_j *cols] = 1;

push(active_pts, init_pos); // Insert it into the background grid
push(active, 0); // push the first point=0 is active!


// ---------------------  step2
int sel_active_seed = 1237;
int push_active_index = 0;

while(len(active)!=0 ){
    // choose a random index from it (say i) from active list
    //printf("length active:%d
", len(active) );
    int rand_index =  floor(rand(sel_active_seed) * float(len(active) ) );
  
    
    //printf("sel rand active index: %d
", rand_index);
    vector pos = active_pts[rand_index];
    
    
    
    int found = false;
    //Generate up to k points chosen uniformly from thespherical annulus between radius r and 2r) aroundxi.
    for(int ns = 0; ns < k; ns++){
        float genkseed = rand( (sel_active_seed + ns )*564  ) * 2 * PI;
        vector genkpos;
        genkpos.x = sin(genkseed);
        genkpos.y = cos(genkseed);
        genkpos.z = 0;
          
        float vlen = rand(ns*135 + sel_active_seed *54);
        vlen = fit(vlen, 0, 1, r, 2*r);
        
        genkpos *= vlen;
        genkpos += pos; // generate our k points
        
        //printf("%f %f 
", genkpos.x, genkpos.y);
        //
        int col = floor(genkpos.x / w);
        int row = floor(genkpos.y / w);
        //printf("col:%d, row:%d
",col,row);
        
        int cond_0 = col >-1 && col < cols;
        int cond_1 = row >-1 && row < rows;
        int cond_2 = grid[col + row * cols] == 1;
  
   
        if(cond_0 && cond_1 && !cond_2){
            // query neibours
            int should_insert = 1;
            
            for(int i= -1; i<=1; i++ ){
                for(int j =-1; j<=1; j++){
                    int neibour_index = (col+i) + (row+j) * cols;
                    int neibour_grid_data = grid[neibour_index];
                    
                    vector neibour_pos = grid_pts[neibour_index];
                    if(neibour_grid_data == 1){
                        float dist = distance(genkpos, neibour_pos);
                        if(dist < r) should_insert = 0;
                    }
                }
            }// end of query near points
            
            if(should_insert) {
                //printf("insert a point
");
                found = 1;
                
                // update the grid
                grid[col + row * cols] = 1;
                grid_pts[col + row * cols] = genkpos;
                
                // update active
                push(active_pts, genkpos);
                push_active_index +=1;
                push(active, push_active_index);
                push(process_pts, genkpos);
                break;
            }
        
        }
        
        
    }
    
    if(!found){
    /*
        printf("{
");
        for(int i: active){
            printf("now active list value: %d
", i);
        }
        printf("}
");
        printf("---now remove id:%d
", rand_index);
        if(removevalue(active,rand_index) ){
            printf("---after remove active, now has length:%d 
", len(active));
        }
        else{
            printf("****remove id:%d faild
", rand_index);
        }*/
        
        pop(active, rand_index);
        pop(active_pts, rand_index);
    }
    
    sel_active_seed += 135;
}

int i=0;
for(vector p: process_pts){
    int pt = addpoint(geoself(), p);
    vector c0 = set(1,0,1);
    vector c1 = set(0,1,1);
    float bias = (sin( float(i*0.002) ) + 1.0) / 2.0f;
    setpointattrib(geoself(), "Cd" , pt , lerp(c0,c1, bias)   ); 
    ++i;
    
}

https://www.jasondavies.com/poisson-disc/

https://www.youtube.com/watch?v=flQgnCUxHlw&t=438s

http://devmag.org.za/2009/05/03/poisson-disk-sampling/