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/