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 r = random();
315 struct triangle *tri = &drawing->tri[r % drawing->num_tris];
317 r /= drawing->num_tris;
321 tri->coord[r/2].x = random() % master->width;
323 tri->coord[r/2].y = random() % master->height;
325 tri->color[r - 6] = random() % 256;
327 calc_multipliers(tri);
330 /* Breed two drawings together, and throw in a mutation. */
331 static struct drawing *breed_drawing(const void *ctx,
332 const struct drawing *a,
333 const struct drawing *b,
334 const struct image *master)
337 struct drawing *drawing;
338 unsigned int r = random(), randmax = RAND_MAX;
340 assert(a->num_tris == b->num_tris);
341 drawing = new_drawing(ctx, a->num_tris);
343 for (i = 0; i < a->num_tris; i++) {
347 drawing->tri[i] = a->tri[i];
351 drawing->tri[i] = b->tri[i];
361 mutate_drawing(drawing, master);
362 score_drawing(drawing, master);
366 /* This is our anti-thread. It does the time-consuming operation of
367 * breeding two drawings together and scoring the result. */
368 static void *breeder(struct at_pool *atp, const struct image *master)
370 const struct drawing *a, *b;
372 /* For simplicity, controller just hands us two pointers in
373 * separate writes. It could put them in the pool for us to
375 while ((a = at_read_parent(atp)) != NULL) {
376 struct drawing *child;
377 b = at_read_parent(atp);
379 child = breed_drawing(at_pool_ctx(atp), a, b, master);
380 at_tell_parent(atp, child);
382 /* Unused: we never exit. */
386 /* We keep a very rough count of how much work each athread has. This
387 * works fine since fairly all rendering takes about the same time.
389 * Better alternative would be to put all the pending work somewhere
390 * in the shared area and notify any idle thread. The threads would
391 * keep looking in that shared area until they can't see any more
392 * work, then they'd at_tell_parent() back. */
393 struct athread_work {
395 unsigned int pending;
398 /* It's assumed that there isn't more than num results pending. */
399 static unsigned gather_results(struct athread_work *athreads,
400 unsigned int num_threads,
401 struct drawing **drawing,
405 unsigned int i, maxfd = 0, done = 0;
406 struct timeval zero = { .tv_sec = 0, .tv_usec = 0 };
408 /* If it mattered, we could cache this fd mask and maxfd. */
409 for (i = 0; i < num_threads; i++) {
410 if (at_fd(athreads[i].at) > maxfd)
411 maxfd = at_fd(athreads[i].at);
417 for (i = 0; i < num_threads; i++)
418 FD_SET(at_fd(athreads[i].at), &set);
420 if (select(maxfd+1, &set, NULL, NULL, block ? NULL : &zero) < 0)
421 err(1, "Selecting on antithread results");
423 for (i = 0; i < num_threads; i++) {
424 if (!FD_ISSET(at_fd(athreads[i].at), &set))
426 *drawing = at_read(athreads[i].at);
428 err(1, "Error with thread %u", i);
431 athreads[i].pending--;
434 } while (block && num);
439 /* Hand work to an antithread to breed these two together. */
440 static void tell_some_breeder(struct athread_work *athreads,
441 unsigned int num_threads,
442 const struct drawing *a, const struct drawing *b)
444 unsigned int i, best = 0;
446 /* Find least loaded thread. */
447 for (i = 1; i < num_threads; i++) {
448 if (athreads[i].pending < athreads[best].pending)
452 at_tell(athreads[best].at, a);
453 at_tell(athreads[best].at, b);
454 athreads[best].pending++;
457 /* We seed initial triangles colours from the master image. */
458 static const unsigned char *initial_random_color(const struct image *master)
460 return master->buffer + (random() % (master->height * master->width))*3;
463 /* Create an initial random drawing. */
464 static struct drawing *random_drawing(const void *ctx,
465 const struct image *master,
466 unsigned int num_tris)
468 struct drawing *drawing = new_drawing(ctx, num_tris);
471 for (i = 0; i < drawing->num_tris; i++) {
473 struct triangle *tri = &drawing->tri[i];
474 for (c = 0; c < 3; c++) {
475 tri->coord[c].x = random() % master->width;
476 tri->coord[c].y = random() % master->height;
478 memcpy(tri->color, initial_random_color(master), 3);
479 tri->color[3] = (random() % 255) + 1;
480 calc_multipliers(tri);
482 score_drawing(drawing, master);
486 /* Read in a drawing from the saved state file. */
487 static struct drawing *read_drawing(const void *ctx, FILE *in,
488 const struct image *master)
490 struct drawing *drawing;
493 if (fscanf(in, "%u triangles:\n", &i) != 1)
494 errx(1, "Reading saved state");
495 drawing = new_drawing(ctx, i);
496 for (i = 0; i < drawing->num_tris; i++) {
497 unsigned int color[4];
498 if (fscanf(in, "%u,%u,%u,%u,%u,%u,%u,%u,%u,%u\n",
499 &drawing->tri[i].coord[0].x,
500 &drawing->tri[i].coord[0].y,
501 &drawing->tri[i].coord[1].x,
502 &drawing->tri[i].coord[1].y,
503 &drawing->tri[i].coord[2].x,
504 &drawing->tri[i].coord[2].y,
505 &color[0], &color[1], &color[2], &color[3]) != 10)
506 errx(1, "Reading saved state");
507 drawing->tri[i].color[0] = color[0];
508 drawing->tri[i].color[1] = color[1];
509 drawing->tri[i].color[2] = color[2];
510 drawing->tri[i].color[3] = color[3];
511 calc_multipliers(&drawing->tri[i]);
513 score_drawing(drawing, master);
517 /* Comparison function for sorting drawings best to worst. */
518 static int compare_drawing_scores(const void *_a, const void *_b)
520 struct drawing **a = (void *)_a, **b = (void *)_b;
522 return (*a)->score - (*b)->score;
525 /* Save one drawing to state file */
526 static void dump_drawings(struct drawing **drawing, const char *outname)
530 char *tmpout = talloc_asprintf(NULL, "%s.tmp", outname);
532 out = fopen(tmpout, "w");
534 err(1, "Opening %s", tmpout);
535 fprintf(out, "POPULATION_SIZE=%u\n", POPULATION_SIZE);
536 for (i = 0; i < POPULATION_SIZE; i++) {
537 fprintf(out, "%u triangles:\n", drawing[i]->num_tris);
538 for (j = 0; j < drawing[i]->num_tris; j++) {
539 fprintf(out, "%u,%u,%u,%u,%u,%u,%u,%u,%u,%u\n",
540 drawing[i]->tri[i].coord[0].x,
541 drawing[i]->tri[i].coord[0].y,
542 drawing[i]->tri[i].coord[1].x,
543 drawing[i]->tri[i].coord[1].y,
544 drawing[i]->tri[i].coord[2].x,
545 drawing[i]->tri[i].coord[2].y,
546 drawing[i]->tri[j].color[0],
547 drawing[i]->tri[j].color[1],
548 drawing[i]->tri[j].color[2],
549 drawing[i]->tri[j].color[3]);
552 if (fclose(out) != 0)
553 err(1, "Failure closing %s", tmpout);
555 if (rename(tmpout, outname) != 0)
556 err(1, "Renaming %s over %s", tmpout, outname);
560 /* Save state file. */
561 static void dump_state(struct drawing *drawing[POPULATION_SIZE],
562 const struct image *master,
564 const char *outstate,
567 char *out = talloc_asprintf(NULL, "%s.%08u.jpg", outpic, gen);
569 printf("Dumping gen %u to %s & %s\n", gen, out, outstate);
570 dump_drawings(drawing, outstate);
571 image = image_of_drawing(out, drawing[0], master);
572 write_jpeg_file(image, out, 80);
576 /* Biassed coinflip moves us towards top performers. I didn't spend
577 * too much time on it, but 1/32 seems to give decent results (see
578 * breeding-algos.gnumeric). */
579 static struct drawing *select_good_drawing(struct drawing *drawing[],
580 unsigned int population)
585 return select_good_drawing(drawing, population/2);
586 return select_good_drawing(drawing + population/2, population/2);
589 static void usage(void)
591 errx(1, "usage: <infile> <outfile> <statefile> <numtriangles> <numcpus> [<instatefile>]");
594 int main(int argc, char *argv[])
596 struct image *master;
597 unsigned int gen, since_prev_best, num_threads, i;
598 struct drawing *drawing[POPULATION_SIZE];
599 unsigned long prev_best, poolsize;
601 struct athread_work *athreads;
603 if (argc != 6 && argc != 7)
606 /* Room for triangles and master image, with some spare.
607 * We should really read master image header first, so we can be
608 * more precise than "about 3MB". ccan/alloc also needs some
609 * more work to be more efficient. */
610 poolsize = (POPULATION_SIZE + POPULATION_SIZE/4) * (sizeof(struct drawing) + atoi(argv[4]) * sizeof(struct triangle)) * 2 + 1000 * 1000 * 3;
611 atp = at_pool(poolsize);
613 err(1, "Creating pool of %lu bytes", poolsize);
615 /* Auto-free the pool and anything hanging off it (eg. threads). */
616 talloc_steal(talloc_autofree_context(), atp);
619 master = read_jpeg_file(at_pool_ctx(atp), argv[1]);
622 printf("Creating initial population");
624 for (i = 0; i < POPULATION_SIZE; i++) {
625 drawing[i] = random_drawing(at_pool_ctx(atp),
626 master, atoi(argv[4]));
634 state = fopen(argv[5], "r");
636 err(1, "Opening %s", argv[5]);
638 fgets(header, 100, state);
639 printf("Loading initial population from %s: %s", argv[5],
641 for (i = 0; i < POPULATION_SIZE; i++) {
642 drawing[i] = read_drawing(at_pool_ctx(atp),
649 num_threads = atoi(argv[5]);
653 /* Hang the threads off the pool (not *in* the pool). */
654 athreads = talloc_array(atp, struct athread_work, num_threads);
655 for (i = 0; i < num_threads; i++) {
656 athreads[i].pending = 0;
657 athreads[i].at = at_run(atp, breeder, master);
659 err(1, "Creating antithread %u", i);
663 /* Worse than theoretically worst case. */
664 prev_best = master->height * master->stride * 256;
666 for (gen = 0; since_prev_best < PLATEAU_GENS; gen++) {
667 unsigned int j, done = 0;
668 struct drawing *new[POPULATION_SIZE/4];
670 qsort(drawing, POPULATION_SIZE, sizeof(drawing[0]),
671 compare_drawing_scores);
673 printf("Best %lu, worst %lu\n",
674 drawing[0]->score, drawing[POPULATION_SIZE-1]->score);
676 /* Probability of being chosen to breed depends on
677 * rank. We breed over lowest 1/4 population. */
678 for (j = 0; j < POPULATION_SIZE / 4; j++) {
679 struct drawing *best1, *best2;
681 best1 = select_good_drawing(drawing, POPULATION_SIZE);
682 best2 = select_good_drawing(drawing, POPULATION_SIZE);
684 tell_some_breeder(athreads, num_threads, best1, best2);
686 /* We reap during loop, so return pipes don't fill.
687 * See "Better alternative" above. */
688 done += gather_results(athreads, num_threads,
689 new + done, j - done, false);
692 /* Collate final results. */
693 gather_results(athreads, num_threads, new+done, j-done, true);
695 /* Overwrite bottom 1/4 */
696 for (j = POPULATION_SIZE * 3 / 4; j < POPULATION_SIZE; j++) {
697 talloc_free(drawing[j]);
698 drawing[j] = new[j - POPULATION_SIZE * 3 / 4];
701 /* We dump on every 1% improvement in score. */
702 if (drawing[0]->score < prev_best * 0.99) {
704 dump_state(drawing, master, argv[2], argv[3], gen);
706 prev_best = drawing[0]->score;
717 /* Dump final state */
718 printf("No improvement over %lu for %u gens\n",
719 prev_best, since_prev_best);
720 dump_state(drawing, master, argv[2], argv[3], gen);