poisson.c (6454B)
1 2 #include <stdlib.h> 3 #include <math.h> 4 #include "util.h" 5 #include "common.h" 6 #include "poisson.h" 7 8 9 struct grid_info { 10 double size; 11 int *grid; 12 int *active; 13 int cellw; 14 int cells; 15 float cell_size; 16 struct point *samples; 17 int n_samples; 18 int n_active; 19 }; 20 21 void 22 draw_samples(struct point *samples, double point_dist, const int nsamples, const int size) { 23 const int dimx = size, dimy = size; 24 int x, y; 25 int close = 0; 26 FILE *fp = fopen("poisson.ppm", "wb"); /* b - binary mode */ 27 fprintf(fp, "P6\n%d %d\n255\n", dimx, dimy); 28 for (y = 0; y < dimy; ++y) 29 for (x = 0; x < dimx; ++x) { 30 static unsigned char color[3]; 31 double dx, dy, dist; 32 for (int i = 0; i < nsamples; ++i) { 33 dx = samples[i].x - (double)x; 34 dy = samples[i].y - (double)y; 35 dist = dx * dx + dy * dy; 36 if (dist <= 2) { 37 close = 1; 38 break; 39 } 40 } 41 /* if (fmod(x, point_dist/sqrt(2)) <= 1 || */ 42 /* fmod(y, point_dist/sqrt(2)) <= 1) { */ 43 /* color[0] = 180; */ 44 /* color[1] = 180; */ 45 /* color[2] = 180; */ 46 /* } */ 47 if (close) { 48 int c = 255-(int)(dist/(point_dist/sqrt(2)) * 255.0); 49 color[0] = c; 50 color[1] = c; 51 color[2] = c; 52 } 53 else { 54 color[0] = 0; 55 color[1] = 0; 56 color[2] = 0; 57 } 58 fwrite(color, 1, 3, fp); 59 close = 0; 60 } 61 fclose(fp); 62 } 63 64 int 65 grid_lookup(struct grid_info *info, int x, int y) { 66 if (x > info->cellw || y > info->cellw || x < 0 || y < 0) 67 return -2; 68 69 int grid_ix = x * info->cellw + y; 70 assert(grid_ix < info->cells && "blah"); 71 int ix = info->grid[grid_ix]; 72 73 return ix; 74 } 75 76 void 77 save_samples(struct point *samples, int seed, int n_samples) { 78 FILE *fp = fopen("samples.bin", "wb"); /* b - binary mode */ 79 fwrite(&seed, sizeof(seed), 1, fp); 80 fwrite(&n_samples, sizeof(n_samples), 1, fp); 81 fwrite(samples, sizeof(*samples), n_samples, fp); 82 fclose(fp); 83 } 84 85 struct point * 86 load_samples(int *seed, int *n_samples) { 87 FILE *fp = fopen("data/samples-200x200.bin", "rb"); 88 int res; 89 (void)res; 90 int localseed; 91 seed = seed ? seed : &localseed; 92 res = fread(seed, sizeof(*seed), 1, fp); 93 assert(res); 94 res = fread(n_samples, sizeof(*n_samples), 1, fp); 95 assert(res); 96 struct point *samples = malloc(*n_samples * sizeof(*samples)); 97 res = fread(samples, sizeof(*samples), *n_samples, fp); 98 assert(res); 99 fclose(fp); 100 return samples; 101 } 102 103 static int 104 add_sample(struct grid_info *info, struct point *candidate) { 105 int grid_x = floor(candidate->x / info->cell_size); 106 int grid_y = floor(candidate->y / info->cell_size); 107 108 int i = info->n_samples++; 109 int grid_ix = grid_ix = grid_x * info->cellw + grid_y; 110 int ai = info->n_active++; 111 112 assert(info->n_samples < info->cells); 113 assert(grid_ix < info->cells); 114 assert(ai < info->cells); 115 116 info->active[ai] = i; 117 info->samples[i].x = candidate->x; 118 info->samples[i].y = candidate->y; 119 assert(info->grid[grid_ix] == -1); 120 info->grid[grid_ix] = i; 121 return i; 122 } 123 124 static void 125 remove_active(struct grid_info *info, int i) { 126 if (info->n_active == 0) 127 return; 128 129 if (info->n_active == 1 && i == 0) { 130 info->n_active = 0; 131 return; 132 } 133 134 if (info->n_active == 1 && i != 0) 135 assert(!"trying to remove and index that doesn't exist"); 136 137 for (int j = i; j < info->n_active-1; ++j) { 138 assert(j >= 0); 139 assert(j+1 < info->cells); 140 info->active[j] = info->active[j+1]; 141 } 142 143 info->n_active--; 144 } 145 146 struct point * 147 poisson_disk_samples(const double point_dist, double size, 148 const int reject_limit, int *n_samples) { 149 struct point candidate; 150 struct grid_info info; 151 152 double cell_size = point_dist / sqrt(2); 153 154 int cellw = ceil(size/cell_size)+1; 155 int cells = cellw * cellw; 156 157 int *grid = malloc(cells * sizeof(*grid)); 158 int *active = malloc(cells * sizeof(*active)); 159 struct point *samples = malloc(cells * sizeof(*samples)); 160 161 info.n_samples = 0; 162 info.grid = grid; 163 info.samples = samples; 164 info.active = active; 165 info.n_active = 0; 166 info.cellw = cellw; 167 info.cells = cells; 168 info.size = size; 169 info.cell_size = cell_size; 170 171 for (int i = 0; i < cells; ++i) 172 grid[i] = -1; 173 174 candidate.x = rand_0to1() * size; 175 candidate.y = rand_0to1() * size; 176 177 add_sample(&info, &candidate); 178 179 assert(info.n_active == 1); 180 assert(info.n_samples == 1); 181 182 int nexti = -1; 183 while (info.n_active > 0) { 184 int rand_i = nexti != -1 ? nexti : rand() % info.n_active; 185 struct point *xi = &samples[rand_i]; 186 187 int all_rejected = 1; 188 189 for (int i = 0, tries = 0; i < reject_limit + tries; ++i) { 190 int sx = rand() % 2 ? 1 : -1; 191 int sy = rand() % 2 ? 1 : -1; 192 193 double rx = ((rand_0to1() * point_dist) + point_dist) * sx; 194 double ry = ((rand_0to1() * point_dist) + point_dist) * sy; 195 196 /* printf("xi %f xy %f rx %f ry %f\n", xi->x, xi->y, rx, ry); */ 197 198 candidate.x = xi->x + rx; 199 candidate.y = xi->y + ry; 200 201 int cx = floor(candidate.x / info.cell_size); 202 int cy = floor(candidate.y / info.cell_size); 203 204 if (candidate.x > size || candidate.x < 0 || candidate.y > size || candidate.y < 0) { 205 tries++; 206 continue; 207 } 208 209 int nearby[] = { 210 grid_lookup(&info, cx - 1, cy - 1), 211 grid_lookup(&info, cx + 0, cy - 1), 212 grid_lookup(&info, cx + 1, cy - 1), 213 214 grid_lookup(&info, cx - 1, cy + 0), 215 grid_lookup(&info, cx + 0, cy + 0), 216 grid_lookup(&info, cx + 1, cy + 0), 217 218 grid_lookup(&info, cx - 1, cy + 1), 219 grid_lookup(&info, cx + 0, cy + 1), 220 grid_lookup(&info, cx + 1, cy + 1), 221 }; 222 223 224 int found = 1; 225 226 if (grid_lookup(&info, cx + 0, cy + 0) != -1) { 227 found = 0; 228 } 229 else { 230 for (size_t j = 0; j < ARRAY_SIZE(nearby); ++j) { 231 int neari = nearby[j]; 232 233 if (neari < 0) 234 continue; 235 struct point *near = &samples[neari]; 236 237 assert(near); 238 239 double dx = near->x - candidate.x; 240 double dy = near->y - candidate.y; 241 double dist = sqrt(dx * dx + dy * dy); 242 if (dist < point_dist) { 243 found = 0; 244 break; 245 } 246 } 247 } 248 249 if (found) { 250 all_rejected = 0; 251 nexti = add_sample(&info, &candidate); 252 } 253 } // gen up to reject_limit samples 254 255 if (all_rejected) { 256 nexti = -1; 257 remove_active(&info, rand_i); 258 } 259 } 260 261 *n_samples = info.n_samples; 262 263 free(grid); 264 free(active); 265 return samples; 266 }