1 /* Antithread example to approximate picture with triangles using
2 * genetic algorithm. Example for a 2 cpu system:
3 * ./arabella arabella.jpg out out.save 2
7 #include <ccan/talloc/talloc.h>
14 #include <sys/select.h>
15 #include <ccan/str/str.h>
16 #include <ccan/antithread/antithread.h>
17 #include <sys/types.h>
20 // Define this to run 100 times without dumping state
23 /* How many drawings in entire population. */
24 #define POPULATION_SIZE 1000
26 /* How many generations without 1% improvement before we terminate. */
27 #define PLATEAU_GENS 200
29 /* An image buffer to render into. */
31 unsigned height, width;
34 unsigned char *buffer;
37 /* A drawing is a (fixed) number of triangles. */
40 unsigned int num_tris;
44 /* Here's a triangle. */
51 unsigned char color[4];
56 /* Get pointer into image at a specific location. */
57 static unsigned char *image_at(struct image *image,
58 unsigned int x, unsigned int y)
60 return image->buffer + y * image->stride + x * 3;
63 /* Blend a dot into this location. */
64 static void add_dot(unsigned char *buf,
65 unsigned char mult, const uint16_t add[])
68 /* /256 isn't quite right, but it's much faster than /255 */
69 for (c = 0; c < 3; c++)
70 buf[c] = (buf[c] * mult + add[c]) / 256;
73 /* Code based on example taken from:
74 * http://www.devmaster.net/forums/showthread.php?t=1094 */
75 static void add_flat_triangle(struct image *image,
76 int x0, int y0, int x1, int y1, int x2, int y2,
77 unsigned char mult, const uint16_t add[])
82 // compute slopes for the two triangle legs
83 float dx0 = (float)(x2 - x0) / (y2 - y0);
84 float dx1 = (float)(x2 - x1) / (y2 - y1);
88 float lx = (float) x0, rx = (float) x1;
92 buf = image_at(image, 0, y0);
95 buf = image_at(image, 0, y2);
99 for (i=0; i < yRange; ++i) {
100 for (j=(int)(lx); j<(int)((rx) + 1.0f); ++j)
101 add_dot(buf + 3 * j, mult, add);
105 buf += image->stride;
109 static void swap(int *a, int *b)
116 static void paint_triangle(struct image *image, const struct triangle *tri)
118 int i0 = 0, i1 = 1, i2 = 2;
119 int x0, y0, x1, y1, x2, y2;
121 /* Could do this on triangle creation. */
122 // sort verts by height
123 if (tri->coord[i0].y > tri->coord[i1].y) swap(&i0, &i1);
124 if (tri->coord[i0].y > tri->coord[i2].y) swap(&i0, &i2);
125 if (tri->coord[i1].y > tri->coord[i2].y) swap(&i1, &i2);
127 x0 = tri->coord[i0].x, y0 = tri->coord[i0].y;
128 x1 = tri->coord[i1].x, y1 = tri->coord[i1].y;
129 x2 = tri->coord[i2].x, y2 = tri->coord[i2].y;
131 // test for easy cases, else split trinagle in two and render both halfs
133 if (x1 > x2) swap(&x1, &x2);
134 add_flat_triangle(image, x1, y1, x2, y2, x0, y0,
135 tri->mult, tri->add);
136 } else if (y0 == y1) {
137 if (x0 > x1) swap(&x0, &x1);
138 add_flat_triangle(image, x0, y0, x1, y1, x2, y2,
139 tri->mult, tri->add);
141 // compute x pos of the vert that builds the splitting line with x1
142 int tmp_x = x0 + (int)(0.5f + (float)(y1-y0) * (float)(x2-x0) / (float)(y2-y0));
144 if (x1 > tmp_x) swap(&x1, &tmp_x);
145 add_flat_triangle(image, x1, y1, tmp_x, y1, x0, y0,
146 tri->mult, tri->add);
147 add_flat_triangle(image, x1, y1, tmp_x, y1, x2, y2,
148 tri->mult, tri->add);
152 /* Create a new image, allocated off context. */
153 static struct image *new_image(const void *ctx,
154 unsigned width, unsigned height, unsigned stride)
156 struct image *image = talloc(ctx, struct image);
158 image->width = width;
159 image->height = height;
160 image->stride = stride;
161 image->buffer = talloc_zero_array(image, unsigned char,
166 /* Taken from JPEG example code. Quality is 1 to 100. */
167 static void write_jpeg_file(const struct image *image,
168 const char *filename, int quality)
170 struct jpeg_compress_struct cinfo;
171 struct jpeg_error_mgr jerr;
173 JSAMPROW row_pointer[1];
176 cinfo.err = jpeg_std_error(&jerr);
177 jpeg_create_compress(&cinfo);
179 if ((outfile = fopen(filename, "wb")) == NULL)
180 err(1, "can't open %s", filename);
182 jpeg_stdio_dest(&cinfo, outfile);
184 cinfo.image_width = image->width;
185 cinfo.image_height = image->height;
186 cinfo.input_components = 3;
187 cinfo.in_color_space = JCS_RGB;
188 jpeg_set_defaults(&cinfo);
189 jpeg_set_quality(&cinfo, quality, TRUE);
191 jpeg_start_compress(&cinfo, TRUE);
192 row_stride = image->width * 3;
194 while (cinfo.next_scanline < cinfo.image_height) {
195 row_pointer[0] = image->buffer + cinfo.next_scanline*row_stride;
196 (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
199 jpeg_finish_compress(&cinfo);
202 jpeg_destroy_compress(&cinfo);
205 static struct image *read_jpeg_file(const void *ctx, const char *filename)
207 struct jpeg_decompress_struct cinfo;
209 struct jpeg_error_mgr jerr;
213 if ((infile = fopen(filename, "rb")) == NULL)
214 err(1, "can't open %s", filename);
216 cinfo.err = jpeg_std_error(&jerr);
218 jpeg_create_decompress(&cinfo);
220 jpeg_stdio_src(&cinfo, infile);
222 jpeg_read_header(&cinfo, TRUE);
224 jpeg_start_decompress(&cinfo);
225 row_stride = cinfo.output_width * cinfo.output_components;
227 image = new_image(ctx,
228 cinfo.output_width, cinfo.output_height, row_stride);
229 while (cinfo.output_scanline < cinfo.output_height) {
230 JSAMPROW row = &image->buffer[cinfo.output_scanline*row_stride];
231 jpeg_read_scanlines(&cinfo, &row, 1);
234 (void) jpeg_finish_decompress(&cinfo);
235 jpeg_destroy_decompress(&cinfo);
240 /* Higher means closer to perfect match. We assume images same size. */
241 static unsigned long compare_images(const struct image *a,
242 const struct image *b)
244 unsigned long result = 0;
247 /* Huge images won't work here. We'd need to get cleverer. */
248 assert(a->height * a->stride < ULONG_MAX / 8);
250 for (i = 0; i < a->height * a->stride; i++) {
251 if (a->buffer[i] > b->buffer[i])
252 result += a->buffer[i] - b->buffer[i];
254 result += b->buffer[i] - a->buffer[i];
259 /* Precalculate the alpha adds and multiplies for this color/alpha combo. */
260 static void calc_multipliers(struct triangle *tri)
262 /* Multiply by 255 - alpha. */
263 tri->mult = (255 - tri->color[3]);
264 /* Add alpha * color */
265 tri->add[0] = (tri->color[0] * tri->color[3]);
266 tri->add[1] = (tri->color[1] * tri->color[3]);
267 tri->add[2] = (tri->color[2] * tri->color[3]);
270 /* Render the image of this drawing, and return it. */
271 static struct image *image_of_drawing(const void *ctx,
272 const struct drawing *drawing,
273 const struct image *master)
278 image = new_image(ctx, master->width, master->height, master->stride);
280 for (i = 0; i < drawing->num_tris; i++)
281 paint_triangle(image, &drawing->tri[i]);
285 /* Render the image and compare with the master. */
286 static void score_drawing(struct drawing *drawing,
287 const struct image *master)
291 /* We don't allocate it off the drawing, since we don't need
292 * it inside the shared area. */
293 image = image_of_drawing(NULL, drawing, master);
294 drawing->score = compare_images(image, master);
298 /* Create a new drawing allocated off context (which is the antithread
299 * pool context, so this is all allocated inside the pool so the
300 * antithreads can access it).
302 static struct drawing *new_drawing(const void *ctx, unsigned int num_tris)
304 struct drawing *drawing = talloc_zero(ctx, struct drawing);
305 drawing->num_tris = num_tris;
306 drawing->tri = talloc_array(drawing, struct triangle, num_tris);
310 /* Make a random change to the drawing: frob one triangle. */
311 static void mutate_drawing(struct drawing *drawing,
312 const struct image *master)
314 unsigned int i, r = random();
315 struct triangle *tri = &drawing->tri[r % drawing->num_tris];
317 r /= drawing->num_tris;
320 /* Move one corner in x or y dir. */
322 tri->coord[r/2].x = random() % master->width;
324 tri->coord[r/2].y = random() % master->height;
326 /* Change one aspect of color. */
327 tri->color[r - 6] = random() % 256;
328 } else if (r == 10) {
329 /* Completely move a triangle. */
330 for (i = 0; i < 3; i++) {
331 tri->coord[i].x = random() % master->width;
332 tri->coord[i].y = random() % master->height;
335 /* Completely change a triangle's colour. */
336 for (i = 0; i < 4; i++)
337 tri->color[i] = random() % 256;
339 calc_multipliers(tri);
342 /* Breed two drawings together, and throw in a mutation. */
343 static struct drawing *breed_drawing(const void *ctx,
344 const struct drawing *a,
345 const struct drawing *b,
346 const struct image *master)
349 struct drawing *drawing;
350 unsigned int r = random(), randmax = RAND_MAX;
352 assert(a->num_tris == b->num_tris);
353 drawing = new_drawing(ctx, a->num_tris);
355 for (i = 0; i < a->num_tris; i++) {
359 drawing->tri[i] = a->tri[i];
363 drawing->tri[i] = b->tri[i];
373 mutate_drawing(drawing, master);
374 score_drawing(drawing, master);
378 /* This is our anti-thread. It does the time-consuming operation of
379 * breeding two drawings together and scoring the result. */
380 static void *breeder(struct at_pool *atp, const struct image *master)
382 const struct drawing *a, *b;
384 /* For simplicity, controller just hands us two pointers in
385 * separate writes. It could put them in the pool for us to
387 while ((a = at_read_parent(atp)) != NULL) {
388 struct drawing *child;
389 b = at_read_parent(atp);
391 child = breed_drawing(at_pool_ctx(atp), a, b, master);
392 at_tell_parent(atp, child);
394 /* Unused: we never exit. */
398 /* We keep a very rough count of how much work each athread has. This
399 * works fine since fairly all rendering takes about the same time.
401 * Better alternative would be to put all the pending work somewhere
402 * in the shared area and notify any idle thread. The threads would
403 * keep looking in that shared area until they can't see any more
404 * work, then they'd at_tell_parent() back. */
405 struct athread_work {
407 unsigned int pending;
410 /* It's assumed that there isn't more than num results pending. */
411 static unsigned gather_results(struct athread_work *athreads,
412 unsigned int num_threads,
413 struct drawing **drawing,
417 unsigned int i, maxfd = 0, done = 0;
418 struct timeval zero = { .tv_sec = 0, .tv_usec = 0 };
420 /* If it mattered, we could cache this fd mask and maxfd. */
421 for (i = 0; i < num_threads; i++) {
422 if (at_fd(athreads[i].at) > maxfd)
423 maxfd = at_fd(athreads[i].at);
429 for (i = 0; i < num_threads; i++)
430 FD_SET(at_fd(athreads[i].at), &set);
432 if (select(maxfd+1, &set, NULL, NULL, block ? NULL : &zero) < 0)
433 err(1, "Selecting on antithread results");
435 for (i = 0; i < num_threads; i++) {
436 if (!FD_ISSET(at_fd(athreads[i].at), &set))
438 *drawing = at_read(athreads[i].at);
440 err(1, "Error with thread %u", i);
443 athreads[i].pending--;
446 } while (block && num);
451 /* Hand work to an antithread to breed these two together. */
452 static void tell_some_breeder(struct athread_work *athreads,
453 unsigned int num_threads,
454 const struct drawing *a, const struct drawing *b)
456 unsigned int i, best = 0;
458 /* Find least loaded thread. */
459 for (i = 1; i < num_threads; i++) {
460 if (athreads[i].pending < athreads[best].pending)
464 at_tell(athreads[best].at, a);
465 at_tell(athreads[best].at, b);
466 athreads[best].pending++;
469 /* We seed initial triangles colours from the master image. */
470 static const unsigned char *initial_random_color(const struct image *master)
472 return master->buffer + (random() % (master->height * master->width))*3;
475 /* Create an initial random drawing. */
476 static struct drawing *random_drawing(const void *ctx,
477 const struct image *master,
478 unsigned int num_tris)
480 struct drawing *drawing = new_drawing(ctx, num_tris);
483 for (i = 0; i < drawing->num_tris; i++) {
485 struct triangle *tri = &drawing->tri[i];
486 for (c = 0; c < 3; c++) {
487 tri->coord[c].x = random() % master->width;
488 tri->coord[c].y = random() % master->height;
490 memcpy(tri->color, initial_random_color(master), 3);
491 tri->color[3] = (random() % 255) + 1;
492 calc_multipliers(tri);
494 score_drawing(drawing, master);
498 /* Read in a drawing from the saved state file. */
499 static struct drawing *read_drawing(const void *ctx, FILE *in,
500 const struct image *master)
502 struct drawing *drawing;
505 if (fscanf(in, "%u triangles:\n", &i) != 1)
506 errx(1, "Reading saved state");
507 drawing = new_drawing(ctx, i);
508 for (i = 0; i < drawing->num_tris; i++) {
509 unsigned int color[4];
510 if (fscanf(in, "%u,%u,%u,%u,%u,%u,%u,%u,%u,%u\n",
511 &drawing->tri[i].coord[0].x,
512 &drawing->tri[i].coord[0].y,
513 &drawing->tri[i].coord[1].x,
514 &drawing->tri[i].coord[1].y,
515 &drawing->tri[i].coord[2].x,
516 &drawing->tri[i].coord[2].y,
517 &color[0], &color[1], &color[2], &color[3]) != 10)
518 errx(1, "Reading saved state");
519 drawing->tri[i].color[0] = color[0];
520 drawing->tri[i].color[1] = color[1];
521 drawing->tri[i].color[2] = color[2];
522 drawing->tri[i].color[3] = color[3];
523 calc_multipliers(&drawing->tri[i]);
525 score_drawing(drawing, master);
529 /* Comparison function for sorting drawings best to worst. */
530 static int compare_drawing_scores(const void *_a, const void *_b)
532 struct drawing **a = (void *)_a, **b = (void *)_b;
534 return (*a)->score - (*b)->score;
537 /* Save one drawing to state file */
538 static void dump_drawings(struct drawing **drawing, const char *outname)
542 char *tmpout = talloc_asprintf(NULL, "%s.tmp", outname);
544 out = fopen(tmpout, "w");
546 err(1, "Opening %s", tmpout);
547 fprintf(out, "POPULATION_SIZE=%u\n", POPULATION_SIZE);
548 for (i = 0; i < POPULATION_SIZE; i++) {
549 fprintf(out, "%u triangles:\n", drawing[i]->num_tris);
550 for (j = 0; j < drawing[i]->num_tris; j++) {
551 fprintf(out, "%u,%u,%u,%u,%u,%u,%u,%u,%u,%u\n",
552 drawing[i]->tri[j].coord[0].x,
553 drawing[i]->tri[j].coord[0].y,
554 drawing[i]->tri[j].coord[1].x,
555 drawing[i]->tri[j].coord[1].y,
556 drawing[i]->tri[j].coord[2].x,
557 drawing[i]->tri[j].coord[2].y,
558 drawing[i]->tri[j].color[0],
559 drawing[i]->tri[j].color[1],
560 drawing[i]->tri[j].color[2],
561 drawing[i]->tri[j].color[3]);
564 if (fclose(out) != 0)
565 err(1, "Failure closing %s", tmpout);
567 if (rename(tmpout, outname) != 0)
568 err(1, "Renaming %s over %s", tmpout, outname);
572 /* Save state file. */
573 static void dump_state(struct drawing *drawing[POPULATION_SIZE],
574 const struct image *master,
576 const char *outstate,
579 char *out = talloc_asprintf(NULL, "%s.%08u.jpg", outpic, gen);
581 printf("Dumping gen %u to %s & %s\n", gen, out, outstate);
582 dump_drawings(drawing, outstate);
583 image = image_of_drawing(out, drawing[0], master);
584 write_jpeg_file(image, out, 80);
588 /* Biassed coinflip moves us towards top performers. I didn't spend
589 * too much time on it, but 1/32 seems to give decent results (see
590 * breeding-algos.gnumeric). */
591 static struct drawing *select_good_drawing(struct drawing *drawing[],
592 unsigned int population)
597 return select_good_drawing(drawing, population/2);
598 return select_good_drawing(drawing + population/2, population/2);
601 static void usage(void)
603 errx(1, "usage: <infile> <outfile> <statefile> <numtriangles> <numcpus> [<instatefile>]");
606 int main(int argc, char *argv[])
608 struct image *master;
609 unsigned int gen, since_prev_best, num_threads, i;
610 struct drawing *drawing[POPULATION_SIZE];
611 unsigned long prev_best, poolsize;
613 struct athread_work *athreads;
615 if (argc != 6 && argc != 7)
618 /* Room for triangles and master image, with some spare.
619 * We should really read master image header first, so we can be
620 * more precise than "about 3MB". ccan/alloc also needs some
621 * more work to be more efficient. */
622 poolsize = (POPULATION_SIZE + POPULATION_SIZE/4) * (sizeof(struct drawing) + atoi(argv[4]) * sizeof(struct triangle)) * 2 + 1000 * 1000 * 3;
623 atp = at_pool(poolsize);
625 err(1, "Creating pool of %lu bytes", poolsize);
627 /* Auto-free the pool and anything hanging off it (eg. threads). */
628 talloc_steal(talloc_autofree_context(), atp);
631 master = read_jpeg_file(at_pool_ctx(atp), argv[1]);
634 printf("Creating initial population");
636 for (i = 0; i < POPULATION_SIZE; i++) {
637 drawing[i] = random_drawing(at_pool_ctx(atp),
638 master, atoi(argv[4]));
646 state = fopen(argv[6], "r");
648 err(1, "Opening %s", argv[6]);
650 fgets(header, 100, state);
651 printf("Loading initial population from %s: %s", argv[6],
653 for (i = 0; i < POPULATION_SIZE; i++) {
654 drawing[i] = read_drawing(at_pool_ctx(atp),
661 num_threads = atoi(argv[5]);
665 /* Hang the threads off the pool (not *in* the pool). */
666 athreads = talloc_array(atp, struct athread_work, num_threads);
667 for (i = 0; i < num_threads; i++) {
668 athreads[i].pending = 0;
669 athreads[i].at = at_run(atp, breeder, master);
671 err(1, "Creating antithread %u", i);
675 /* Worse than theoretically worst case. */
676 prev_best = master->height * master->stride * 256;
678 for (gen = 0; since_prev_best < PLATEAU_GENS; gen++) {
679 unsigned int j, done = 0;
680 struct drawing *new[POPULATION_SIZE/4];
682 qsort(drawing, POPULATION_SIZE, sizeof(drawing[0]),
683 compare_drawing_scores);
685 printf("Best %lu, worst %lu\n",
686 drawing[0]->score, drawing[POPULATION_SIZE-1]->score);
688 /* Probability of being chosen to breed depends on
689 * rank. We breed over lowest 1/4 population. */
690 for (j = 0; j < POPULATION_SIZE / 4; j++) {
691 struct drawing *best1, *best2;
693 best1 = select_good_drawing(drawing, POPULATION_SIZE);
694 best2 = select_good_drawing(drawing, POPULATION_SIZE);
696 tell_some_breeder(athreads, num_threads, best1, best2);
698 /* We reap during loop, so return pipes don't fill.
699 * See "Better alternative" above. */
700 done += gather_results(athreads, num_threads,
701 new + done, j - done, false);
704 /* Collate final results. */
705 gather_results(athreads, num_threads, new+done, j-done, true);
707 /* Overwrite bottom 1/4 */
708 for (j = POPULATION_SIZE * 3 / 4; j < POPULATION_SIZE; j++) {
709 talloc_free(drawing[j]);
710 drawing[j] = new[j - POPULATION_SIZE * 3 / 4];
713 /* We dump on every 1% improvement in score. */
714 if (drawing[0]->score < prev_best * 0.99) {
716 dump_state(drawing, master, argv[2], argv[3], gen);
718 prev_best = drawing[0]->score;
729 /* Dump final state */
730 printf("No improvement over %lu for %u gens\n",
731 prev_best, since_prev_best);
732 dump_state(drawing, master, argv[2], argv[3], gen);