円を少し進化させてみる

バッファを埋めた塗りつぶしの円を、少しずつ変化させて目的の画像に近づけるという処理をやってみた。

int mutate_circle (circle *cir)  {
  cir->col = rand() & 0x00ffffff;
  cir->alpha = rand() % 100;
  cir->rad = rand() % 30 + 10;
  return 0;
}

変化させてみて、それがターゲット画像に近くなったかどうかを判定してどんどん変えていく。直径を変えてしまうと、不要な円がどんどん縮んで行く感じで、ちょっと反則というか面白くない感じもあるけど、とりあえず。

判定は、

double diff_buffers (const buf *buf1, const buf *buf2)
  :
  :
	res += sqrt((r1 - r2)*(r1 - r2) + (g1 - g2)*(g1 - g2) + (b1 - b2)*(b1 - b2));
  return res;
}

と色の距離で定義してみた。計算の重たそうなsqrtとか不要かも。

320×240ドットで単純な日の丸をターゲットに5000世代ほど進化させたら、ちょっとそれっぽくなった。

オリジナル
ターゲット
結果画像

知人が賞賛していたので、深く考えずにAmazonで買った「メタヒューリスティクスの数理」（久保幹雄／J.P.ペドロソ）という本をパラパラと眺めていたら、どうも今ぼくがやろうとしている問題は「組み合わせ最適化問題」という名前がついているジャンルの問題じゃないかという気がしている。組み合わせ数が膨大、それの評価関数も結構重たいとかなってくると、とても全組み合わせを試して最適解なんて計算できない。それで、色々とヒューリスティックな方法をやるわけだけど、それでも、現実問題で有効らしい方法論が過去数十年で蓄積しているという。そうした議論をまとめた本。とりあえず適当な出発点から近傍を探して小山を登り、登りきったらちょっと遠くへジャンプみたいな話とか、アリの群生にならって有向グラフに揮発性のフェロモンを落としつつ解の間を移動して最適解を探る戦略とか、計算機科学とかいうよりも生物学っぽくて面白い。

320×240ドット、200個の円、5000回ループで100秒ほどかかる。もうちょっと高速化したいので、pthreadで並列化とかしてみたい。並列化以前にやるべきこともあるような気もするけど。

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

# define NUMBER_OF_CIRCLES 200
# 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;

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

  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);
  }

  for (i = 0; i < height; i++) {
    for (j = 0; j < width; j++) {
      //      b->raw[i][j] = 0x00ff00ff; // purple-ish white
      b->raw[i][j] = BGCOLOR;
    }
  }
  return b;
}

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;
}

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, int x, int y, int rad, int color, int opc) {
  int i, j;
  int upper, lower, left, right;
  float back, fore;
  unsigned char r, g, b;
  unsigned char r2, g2, b2;

  if ((opc < 0) || (opc >100)) {
      return 1;
  }

  r2 = color & 0x000000ff;
  g2 = (color & 0x0000ff00) >>8;
  b2 = (color & 0x00ff0000) >> 16;

  fore = (float)opc / 100;
  back = (100 - (float)opc) / 100;

  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)  {
  cir->col = rand() & 0x00ffffff;
  cir->alpha = rand() % 100;
  cir->rad = rand() % 30 + 10;
  return 0;
}

double diff_buffers (const buf *buf1, const buf *buf2) {
  double res = 0;
  int i, j;
  unsigned char r1, g1, b1, r2, g2, b2;

  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));
    }
  }
  return res;
}

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

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

  srand(time(NULL));

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

  fill_circle(target, 160, 120, 50, 0x000000ff, 100);

  for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
    cir1[i].x = rand() % b1->width;
    cir1[i].y = rand() % b1->height;
    cir1[i].rad = rand() % 30 + 10;
    cir1[i].col = rand() & 0x00ffffff;
    cir1[i].alpha = 80;
    cir2[i] = cir1[i];
  }

  // initial pic

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

  // evolve pic

  for (j = 0; j < 3000; j++) {
    m = rand() % NUMBER_OF_CIRCLES;
    mutate_circle(&cir1[m]);

    clear_buf(b1);
    for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
      fill_circle(b1, cir1[i].x, cir1[i].y, cir1[i].rad, cir1[i].col, cir1[i].alpha);
    }
    
    if (diff_buffers(b1, target) < diff_buffers(b2, target)) {
      cir2[m] = cir1[m];
      clear_buf(b2);
      for (i = 0; i < NUMBER_OF_CIRCLES; i++) {
	fill_circle(b2, cir2[i].x, cir2[i].y, cir2[i].rad, cir2[i].col, cir2[i].alpha);
      }
    } else {
      cir1[m] = cir2[m];
    }
  }

  // result pic

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

  return 0;
}