200個の円で任意の写真を近似

100個とか200個の円で320×240ドットの領域を埋めて、それらのサイズや色、位置、透明度をランダムに変えていくことでターゲットの画像に近づけるという処理。既存画像（ppm形式）をロードする関数を書いたので、実際にどんなもんかと思って数千世代とか1万世代ぐらいやってみたけど、どうもイマイチな感じ。目を細めれば、確かにガンバってる様子が分かるという程度。たぶん元画像が精細すぎるというか、比率の問題が大きそうな気がする。使うべきなのは円じゃなくて三角か四角なのかもしれない。

結果を眺めていると、ある程度以上の世代を経ると、どうも行き詰まっている。ということは、世代が進むに連れてだんだん変異がアクセプトされる率が下がるということで、これは直感にも合う。

実際、どのぐらい進化というか変化が起こるのかを数えてみたのが以下のグラフ。1度に変異させるパラメータを制限したほうが、より変化が起こりやすいんじゃないかという直感も確かめられた。一番進化ペースが落ちないのは直径を固定した場合で、何もかも変化させると進化のペースが落ちてしまい、10世代で何の変化もないというようなことが起こるようになる。

結構処理が重たくて数千世代で2、3分とかそういうオーダー。もうちょい速ければ、結果画像を見ながら、いろいろアプローチを試す気にもなるのになと思って一部の処理をマルチスレッド化してみた。まず簡単そうなところで、2つの画像の距離を測定するdiff_buffers関数を2つ同時に走らせる。

if (diff_buffers(b1, target) < diff_buffers(b2, target)) {

という処理が倍速に近くなったはず。だけど計測してみたら、全体への影響は小さくて、せいぜい5、6％の高速化にとどまった。まあ当然か。

一番重たそうな円を描くプロセスも並列化してみたい。x方向とy方向の2重ループになっているのを並列化するために、y方向を1ステップずつじゃなくてNステップずつの処理にして、ループ開始位置をスレッドごとにずらして関数への引数で渡せばいいんじゃないかという予想……、で、とりあえず円オブジェクトをmainのスタックじゃなくて、ヒープに取るように変更してみた……、ところでハマった。構造体のコピーって、ポインタ経由ではできないっぽい。

typedef struct {
  int hoge;
  int hoge;
} circle;

circle c1[100], c2[100];
:
:
c2[i] = c1[i];

とやっていたものが、

typedef struct {
  int hoge;
  int hoge;
} circle

circle *c1[100], *c2[100];
&c2[i] = &c1[i];
||

とはできないらしい。

>|c|
memcpy(c2[i], c1[i], sizeof(circle));

とメモリをそのままコピーしてみた。うーん。

そもそも、「memcpy(c2[i], c1[i], sizeof(circle))」の意味で、最初は「c2[i] = c1[i]」とポインタの値だけをコピーしてたりしてしまっていた。当然エラーもでないし、一見動いてるように見えるしで、何が壊れたのか一瞬ワケが分からなかった。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <pthread.h>
#include <unistd.h>

# define NUMBER_OF_CIRCLES 400
# define GENERATION 6000
# define BGCOLOR 0x00ffffff
# define WIDTH 320
# define HEIGHT 240

typedef struct {
  int height;
  int width;
  int **raw;
} buf;

typedef struct {
  int x;
  int y;
  int rad;
  int col;
  int alpha;
} circle;

typedef struct {
  buf *buf1;
  buf *buf2;
  double res;
} diff_args;

int clear_buf (buf *b) {
  int i, j;
  for (i = 0; i < b->height; i++) {
    for (j = 0; j < b->width; j++) {
      b->raw[i][j] = BGCOLOR;
    }
  }
  return 0;
}

buf *init_buf (int width, int height) {
  buf *b;
  int i;

  b = malloc(sizeof(buf));
  b->height = height;
  b->width =width;

  b->raw = malloc(sizeof(int*) * height); 
  b->raw[0] = malloc(sizeof(int) * width * height); 
  for (i = 1; i < height; i ++) {
    b->raw[i] = b->raw[0] + (width * i);
  }

  clear_buf(b);
  return b;
}

buf *load_buf (char *filename) {
  FILE *fp;
  char s[255], size[255];
  int i, j;
  int width, height;
  buf *b;

  if ((fp = fopen(filename, "r")) == NULL) {
      printf("Can't open the file\n");
      exit(1);
  }

  // ppm header

  fgets(s,      255, fp); // P6
  fgets(size,   255, fp); // 320 240
  fgets(s,      255, fp); // 255

  width = atoi(size);
  height = atoi(strstr(size, " "));

  b = init_buf(width, height);

  for (i = 0; i < b->height; i++) {
    for (j = 0; j < b->width; j++) {
      fread(&b->raw[i][j], sizeof(char), 3, fp);
    }
  }
  fclose(fp);
  return b;
}

int save_buf (buf *b, char *filename) {
  FILE *fp;
  int i, j;

  if ((fp = fopen(filename, "w")) == NULL) {
      printf("Can't open the file\n");
      exit(1);
  }

  fprintf(fp, "P6\n%d %d\n255\n", b->width, b->height);

  for (i = 0; i < b->height; i++) {
    for (j = 0; j < b->width; j++) {
      fwrite(&b->raw[i][j], sizeof(char), 3, fp);
    }
  }
  fclose(fp);
  return 0;
}

int fill_circle (buf *buf, const circle *cir) {
  int i, j;
  int x, y, rad;
  int upper, lower, left, right;
  float back, fore;
  unsigned char r, g, b;
  unsigned char r2, g2, b2;

  if ((cir->alpha < 0) || (cir->alpha >100)) {
      return 1;
  }

  r2 = cir->col & 0x000000ff;
  g2 = (cir->col & 0x0000ff00) >>8;
  b2 = (cir->col & 0x00ff0000) >> 16;

  fore = (float)cir->alpha / 100;
  back = (100 - (float)cir->alpha) / 100;

  x = cir->x; y = cir->y; rad = cir->rad;

  upper = (y - rad) > 0 ? (y - rad) : 0;
  lower = (y + rad) > buf->height ? buf->height : (y + rad);
  left  = (x - rad) > 0 ? (x - rad) : 0;
  right = (x + rad) > buf->width ? buf->width : (x + rad);

  for (i = upper; i < lower; i++) {
    for (j = left; j < right; j++) {
      if ((rad*rad - (y-i)*(y-i) - (x-j)*(x-j)) > 0) { 
	r = buf->raw[i][j] & 0x000000ff;
	g = (buf->raw[i][j] & 0x0000ff00) >> 8;
	b = (buf->raw[i][j] & 0x00ff0000) >> 16;
	r = (r2 * fore) + (r * back);
	g = (g2 * fore) + (g * back);
	b = (b2 * fore) + (b * back);
	buf->raw[i][j] = (b << 16) + (g << 8) + r;
      }
    }
  }
  return 0;
}

int mutate_circle (circle *cir, int width, int height) {
  cir->col = rand() & 0x00ffffff;
  cir->alpha = rand() % 100;
  cir->rad = rand() % 40 + 5;
  //  cir->x = rand() % width;
  //  cir->y = rand() % height;
  return 0;
}

double p_diff_buffers (void *p) {
  diff_args *args;
  buf *buf1;
  buf *buf2;
  double res = 0;
  int i, j;
  unsigned char r1, g1, b1, r2, g2, b2;

  args = (diff_args *)p;
  buf1 = args->buf1;
  buf2 = args->buf2;

  if ((buf1->width != buf2->width) ||
      (buf1->height != buf2->height)) {
    printf("buffer sizes don't match\n");
    exit(1);
  }

  for (i = 0; i < buf1->height; i++) {
    for (j = 0; j < buf1->width; j++) {
	r1 = buf1->raw[i][j] & 0x000000ff;
	g1 = (buf1->raw[i][j] & 0x0000ff00) >> 8;
	b1 = (buf1->raw[i][j] & 0x00ff0000) >> 16;

	r2 = buf2->raw[i][j] & 0x000000ff;
	g2 = (buf2->raw[i][j] & 0x0000ff00) >> 8;
	b2 = (buf2->raw[i][j] & 0x00ff0000) >> 16;

	res += sqrt((r1 - r2)*(r1 - r2) + (g1 - g2)*(g1 - g2) + (b1 - b2)*(b1 - b2));
    }
  }
  args->res = res;
  pthread_exit(NULL);
}

int main (int argc,  char *argv[]) {
  buf *b1, *b2, *target;
  int i, j, m, e=0;

  char *filename = argv[1];
  circle *cir1[NUMBER_OF_CIRCLES], *cir2[NUMBER_OF_CIRCLES];

  pthread_t thread1, thread2;
  diff_args *args1, *args2;

  args1 = malloc(sizeof(diff_args));
  args2 = malloc(sizeof(diff_args));
  
  srand(time(NULL));

  target = load_buf("photo.ppm");

  b1 = init_buf(WIDTH, HEIGHT);
  b2 = init_buf(WIDTH, HEIGHT);

  for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
    cir1[i] = malloc(sizeof(circle));
    cir2[i] = malloc(sizeof(circle));
    cir1[i]->x = rand() % b1->width;
    cir1[i]->y = rand() % b1->height;
    cir1[i]->rad = rand() % 40 + 5;
    cir1[i]->col = rand() & 0x00ffffff;
    cir1[i]->alpha = 80;
    memcpy(cir2[i], cir1[i], sizeof(circle));
  }

  // initial pic

  for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
    fill_circle(b1, cir1[i]);
    fill_circle(b2, cir2[i]);
  }
  save_buf(b1, "1.ppm");
  save_buf(target, "target.ppm");

  // evolve pic

  for (j = 0; j < GENERATION; j++) {
    m = rand() % NUMBER_OF_CIRCLES;
    mutate_circle(cir1[m], b1->width, b1->height);

    clear_buf(b1);
    for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
      fill_circle(b1, cir1[i]);
    }
    
    args1->buf1 = b1; args1->buf2 = target;
    args2->buf1 = b2; args2->buf2 = target;
    pthread_create(&thread1, NULL, p_diff_buffers, (void *)args1);
    pthread_create(&thread2, NULL, p_diff_buffers, (void *)args2);
    pthread_join(thread1, NULL);
    pthread_join(thread2, NULL);
    //    printf("%f,%f\n", args1->res, args2->res);

    if (args1->res < args2->res) {
      memcpy(cir2[m], cir1[m], sizeof(circle));
      clear_buf(b2);
      for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
	fill_circle(b2, cir2[i]);
      }
      e++;
    } else{
      memcpy(cir1[m], cir2[m], sizeof(circle));
    }
    if ((j % 10) == 0) {
      printf ("%d,   %d\n", j, e);
    }
  }

  // result pic

  for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
    fill_circle(b1, cir1[i]);
  }
  save_buf(b1, "2.ppm");

  return 0;
}