マンデルブロ集合

突然マンデルブロー集合を図示したくなったので暇つぶしに．
マンデルブロー集合とは

という漸化式において，複素数Znがn→∞で発散しないような複素数cのことである．
マンデルブロー集合はガウス平面上でフラクタルとして描画され，相似形の模様が様々なスケールで現れる．これを自己相似性という．

フラクタルは非常に興味深い．
単純な処理を繰り返すことで，自然界に見られる複雑さを表現することができる．
一見複雑に見える模様も，実は自己相似性をもったシンプルな模様の繰り返しであることが多い．
簡単に思いつくところでは木の枝，葉脈，雪の結晶など，フラクタルとして有名なものは，ロマネスコと呼ばれるブロッコリーの仲間や，オウムガイの貝殻の断面図などが挙げられる．

以下マンデルブロー集合の描画の話．
あるc=x+iyに対して十分な回数だけ漸化式を計算し，Znがある閾値以下となった場合，cはマンデルブロー集合に含まれるとした．このcをガウス平面上に黒点としてプロットする．
画像のサイズを300x300とし，x，yを0から300まで計算すると以下のような結果を得た．

左隅に黒い点が見える．
計算式を間違ったのかと思い，コードを見直したがおかしい箇所は見当たらない．
よく考えると，左上を(0.0)とした第1象限を描画しているのであった．
x，yの範囲を-150から150に変更した．

中央に横になったTが見える．
これはスケールの問題であろう．
x,yを-1.5から1.5までにして再度計算したものが下図となる．
今度は美しいマンデルブロー集合を見ることが出来た．

あんなに単純な漸化式から，無限に拡大可能な複雑な模様を生み出すことができるということに感動する．
発散の速度によって色を変化させるという処理を行えば，さらに美しい模様を得ることもできる．

以下，全く洗練されていないソースコード．

#include <opencv2/highgui/highgui.hpp>
#include <cmath> 
#include <iostream>

using namespace std;

const int WINDOW_SIZE = 300;

struct complex1{
  double x;
  double y;
};

struct mycolor{
  int R;
  int G;
  int B;
};

int main ( int argc, char* argv[]){

  struct mycolor color;

  double roi = (double)(atoi(argv[1]));
  double scale = (double)(atof(argv[2]));

  int cnt_1 = 0;
  int cnt_2=0;
  struct complex1 z_n1;
  struct complex1 z_n2;
  struct complex1 c;
  z_n1.x = 0.0; z_n1.y = 0.0;
  z_n2.x = 0.0; z_n2.y = 0.0;
  IplImage* img = cvCreateImage(cvSize(300, 300), IPL_DEPTH_8U, 3);

  int inf_flag = 0.0;
  
  for( c.x=roi; c.x < roi+300.0; c.x += 1.0 ){
    for( c.y=-150.0 ; c.y < 150.0 ; c.y += 1.0 ){
      cnt_1 = 0;
      inf_flag = 0;
      z_n1.x = 0.0; z_n1.y = 0.0;
      z_n2.x = 0.0; z_n2.y = 0.0;

      while(cnt_1 < 2000){
	z_n2.x = z_n1.x * z_n1.x - z_n1.y * z_n1.y + c.x*scale;
	z_n2.y = 2 * z_n1.x * z_n1.y + c.y*scale;
	z_n1 = z_n2;
    
	cnt_1++;

	 if(sqrt(z_n1.x * z_n1.x + z_n1.y * z_n1.y) > 200000.0 ){
	   cout << "hassan" << endl;
	   inf_flag = 1;
	   break;
	 }
      }
      cout << cnt_1 << endl;

      if(inf_flag == 1){
	color.R = (int)(255/cnt_1*1);
	color.G = (int)(255/cnt_1*3);
	color.B = (int)(255/cnt_1*7);
	
      
	if(color.R >255){
	  color.R = 255;
	}if(color.G >255){
	  color.G = 255;
	}if(color.B >255){
	  color.B = 255; 
	}
      }else{
	color.R = 0;
	color.G = 0;
	color.B = 0;
      }

      uchar* ptr =(uchar*)(img->imageData + (int)(c.y + 150.0) * img->widthStep);
     
      ptr[3*(int)(c.x-roi)+0] = color.B;	
      ptr[3*(int)(c.x-roi)+1] = color.G;
      ptr[3*(int)(c.x-roi)+2] = color.R;
      
      if(cnt_2%100 == 0){
	cvNamedWindow ("img", CV_WINDOW_AUTOSIZE);
	cvShowImage ("img", img);
	cvWaitKey(10);
      }
      cnt_2++;
    }
  } 
  
  cvNamedWindow ("img", CV_WINDOW_AUTOSIZE);
  cvShowImage ("img", img);

  cvWaitKey(0);
  cvSaveImage("mndl.png", img);
  cvDestroyWindow("img");

  return 0;
}