はじめに
Wikipedia によると
複素力学系における充填ジュリア集合（じゅうてんジュリアしゅうごう、英: filled Julia set または 英: filled-in Julia set）は、多項式の複素函数を繰り返し適用したときに無限に発散しない複素数の集合である。
とのこと. 本稿では数学的な話には深入りせず, 画像や動画の作成を試みる.
計算方法
を複素数とする. の更新を以下で定義する.
と書くと, 更新式は以下のように書いても同じである.
, の範囲の複素平面 を のグリッドに区切る.
グリッド上の各点 に対応する複素数を と書くと
である.
更新式を使えば から を計算するのは容易である.更新する度に を計算し, となれば発散したとみなして点 を黒で塗りつぶし, 発散しない場合は最大で まで計算する. 最後まで発散しなければ点 を白で塗りつぶす.
以上の手順によって の各点が黒から白で着色され, 充填ジュリア集合の画像が得られる.
プログラム
具体的なプログラムを以下に示す.
#include <stdio.h>
#include <math.h>
#define RE_MAX 1.5
#define RE_MIN -1.5
#define IM_MAX 1.5
#define IM_MIN -1.5
#define MAX_ITER 50
#define BLACK 0x000000
#define WHITE 0xffffff
#define XMAX 256
#define YMAX 256
long gr_screen[YMAX][XMAX];
void putbytes(FILE *f, int n, unsigned long x)
{
while (--n >= 0) {
fputc(x & 255, f); x >>= 8;
}
}
void gr_dot(int x, int y, long color)
{
if (x >= 0 && x < XMAX && y >= 0 && y < YMAX)
gr_screen[y][x] = color;
}
void gr_clear(long color)
{
int x, y;
for (x = 0; x < XMAX; x++)
for (y = 0; y < YMAX; y++)
gr_dot(x, y, color);
}
void gr_BMP(const char *filename)
{
int x, y;
FILE *f = fopen(filename, "wb");
fputs("BM", f); /* ファイルタイプ */
putbytes(f, 4, XMAX * YMAX * 4 + 54); /* ファイルサイズ */
putbytes(f, 4, 0); /* 予約領域 */
putbytes(f, 4, 54); /* ファイル先端から画像までのオフセット */
putbytes(f, 4, 40); /* 情報ヘッダサイズ */
putbytes(f, 4, XMAX); /* 画像の幅 */
putbytes(f, 4, YMAX); /* 画像の高さ */
putbytes(f, 2, 1); /* プレーン数（1） */
putbytes(f, 2, 32); /* 色ビット数 */
putbytes(f, 4, 0); /* 圧縮形式（無圧縮） */
putbytes(f, 4, XMAX * YMAX * 4); /* 画像データサイズ */
putbytes(f, 4, 3780); /* 水平解像度（dot/m） */
putbytes(f, 4, 3780); /* 垂直解像度（dot/m） */
putbytes(f, 4, 0); /* 格納パレット数 */
putbytes(f, 4, 0); /* 重要色数 */
for (y = 0; y < YMAX; y++)
for (x = 0; x < XMAX; x++)
putbytes(f, 4, gr_screen[y][x]); /* 画像データ */
fclose(f);
}
int main()
{
const char* fname = "test.bmp";
int w,h,n,flag;
double x,y,new_x,new_y;
double re_c = -0.06;
double im_c = 0.68;
double alpha = 100000000000000000000.0;
for(w=0;w<XMAX;w++){
for(h=0;h<YMAX;h++){
x = RE_MIN+(RE_MAX-RE_MIN)*(double)w/XMAX;
y = IM_MIN+(IM_MAX-IM_MIN)*(double)h/YMAX;
flag = 0;
for(n=0;n<MAX_ITER;n++){
new_x = x*x-y*y+re_c;
new_y = 2*x*y+im_c;
if(sqrt(new_x*new_x+new_y*new_y)>alpha){
flag = 1;
break;
}
x = new_x;
y = new_y;
}
if(flag){
gr_dot(w,h,BLACK);
}else{
gr_dot(w,h,WHITE);
}
}
}
gr_BMP(fname);
return 0;
}
画像はBMP形式とし, 画像作成のために 文献 のコードを使用した.
また, プログラムの作成にあたって以下を参考とした.
- https://ja.wikipedia.org/wiki/%E5%85%85%E5%A1%AB%E3%82%B8%E3%83%A5%E3%83%AA%E3%82%A2%E9%9B%86%E5%90%88
- https://watlab-blog.com/2020/06/28/julia-set/
実行結果
のとき, 以下の結果が得られた.
プログラムを修正する
発散する点をすべて黒で塗りつぶすのではなく, 発散するまでに必要とした反復回数 をHSL色空間の色相(hue)に対応づける.
具体的には, 以下のように修正を行う.
/* 変更前 */
if(flag){
gr_dot(w,h,BLACK);
}else{
gr_dot(w,h,WHITE);
}
/* 変更後 */
if(flag){
hue = 180.0*n/MAX_ITER;
gr_dot(w, h, hsl2rgb(hue,saturation,lightness));
}else{
gr_dot(w,h,BLACK);
ここで hsl2rgb() はHSL色空間の変数 hue, saturation, lightness を RGBに変換する関数である. 適切に色の変換を行うため, 以下のコードを追加する.
/*
* h : hue(angle) 0.0~360.0
* s : saturation 0.0~1.0
* l : lightness 0.0~1.0
*/
int hsl2rgb(double h, double s, double l)
{
double a = s*(1-fabs(2*l-1))/2;
double vmax = l + a;
double vmin = l - a;
double r,g,b;
int R,G,B;
double diff = vmax-vmin;
if(h<60){
r = vmax;
g = diff*h/60+vmin;
b = vmin;
}else if(h<120){
r = (120-h)/60*diff + vmin;
g = vmax;
b = vmin;
}else if(h<180){
r = vmin;
g = vmax;
b = (h-120)/60*diff + vmin;
}else if(h<240){
r = vmin;
g = (240-h)/60*diff + vmin;
b = vmax;
}else if(h<300){
r = (h-240)/60*diff + vmin;
g = vmin;
b = vmax;
}else if(h<360){
r = vmax;
g = vmin;
b = (360-h)/60*diff + vmin;
}else{
r = vmax;
g = vmax;
b = vmax;
}
R = (int)(r*255);
G = (int)(g*255);
B = (int)(b*255);
return (R<<16) + (G<<8) +B;
}
なお, 関数 hsl2rgb の作成にあたって以下を参考とした.
- https://www.peko-step.com/tool/hslrgb.html
- https://note.kiriukun.com/entry/20181206-rgb-and-hsl-conversion
実行結果
修正したプログラムを実行すると以下の結果が得られた.
発散するまでに必要とした反復回数が少ないほど色相が小さい(=赤色に近い)ため, 発散しやすい領域が画像からはっきり読み取れる.
連番画像と動画
を固定し, を から まで動かす. このとき,
とすれば, の値に対応した充填ジュリア集合の画像が得られる.
を変化させたときのアニメーションを作成するため, 以下のプログラムを作成した. プログラムを実行すると images ディレクトリ配下に連番画像が作成される.
#include <stdio.h>
#include <math.h>
#define RE_MAX 1.5
#define RE_MIN -1.5
#define IM_MAX 1.5
#define IM_MIN -1.5
#define MAX_ITER 50
#define FRAMES 480
#define BLACK 0x000000
#define WHITE 0xffffff
#define XMAX 1080
#define YMAX 1920
long gr_screen[YMAX][XMAX];
void putbytes(FILE *f, int n, unsigned long x)
{
while (--n >= 0) {
fputc(x & 255, f); x >>= 8;
}
}
void gr_dot(int x, int y, long color)
{
if (x >= 0 && x < XMAX && y >= 0 && y < YMAX)
gr_screen[y][x] = color;
}
void gr_clear(long color)
{
int x, y;
for (x = 0; x < XMAX; x++)
for (y = 0; y < YMAX; y++)
gr_dot(x, y, color);
}
void gr_BMP(const char *filename)
{
int x, y;
FILE *f = fopen(filename, "wb");
fputs("BM", f); /* ファイルタイプ */
putbytes(f, 4, XMAX * YMAX * 4 + 54); /* ファイルサイズ */
putbytes(f, 4, 0); /* 予約領域 */
putbytes(f, 4, 54); /* ファイル先端から画像までのオフセット */
putbytes(f, 4, 40); /* 情報ヘッダサイズ */
putbytes(f, 4, XMAX); /* 画像の幅 */
putbytes(f, 4, YMAX); /* 画像の高さ */
putbytes(f, 2, 1); /* プレーン数（1） */
putbytes(f, 2, 32); /* 色ビット数 */
putbytes(f, 4, 0); /* 圧縮形式（無圧縮） */
putbytes(f, 4, XMAX * YMAX * 4); /* 画像データサイズ */
putbytes(f, 4, 3780); /* 水平解像度（dot/m） */
putbytes(f, 4, 3780); /* 垂直解像度（dot/m） */
putbytes(f, 4, 0); /* 格納パレット数 */
putbytes(f, 4, 0); /* 重要色数 */
for (y = 0; y < YMAX; y++)
for (x = 0; x < XMAX; x++)
putbytes(f, 4, gr_screen[y][x]); /* 画像データ */
fclose(f);
}
/*
* h : hue(angle) 0.0~360.0
* s : saturation 0.0~1.0
* l : lightness 0.0~1.0
*/
int hsl2rgb(double h, double s, double l)
{
double a = s*(1-fabs(2*l-1))/2;
double vmax = l + a;
double vmin = l - a;
double r,g,b;
int R,G,B;
double diff = vmax-vmin;
if(h<60){
r = vmax;
g = diff*h/60+vmin;
b = vmin;
}else if(h<120){
r = (120-h)/60*diff + vmin;
g = vmax;
b = vmin;
}else if(h<180){
r = vmin;
g = vmax;
b = (h-120)/60*diff + vmin;
}else if(h<240){
r = vmin;
g = (240-h)/60*diff + vmin;
b = vmax;
}else if(h<300){
r = (h-240)/60*diff + vmin;
g = vmin;
b = vmax;
}else if(h<360){
r = vmax;
g = vmin;
b = (360-h)/60*diff + vmin;
}else{
r = vmax;
g = vmax;
b = vmax;
}
R = (int)(r*255);
G = (int)(g*255);
B = (int)(b*255);
return (R<<16) + (G<<8) +B;
}
int main()
{
char fname[100];
int w,h,n,m,flag;
double x,y,new_x,new_y;
double re_c,im_c;
double abs_c = 0.7;
double alpha = 100000000000000000000.0;
double theta;
double hue;
double saturation = 1.0;
double lightness = 0.5;
for(m=0;m<FRAMES;m++){
sprintf(fname, "./images/img%03d.bmp", m);
theta = 2*M_PI*m/FRAMES;
re_c = abs_c*cos(theta);
im_c = abs_c*sin(theta);
for(w=0;w<XMAX;w++){
for(h=0;h<YMAX;h++){
x = RE_MIN+(RE_MAX-RE_MIN)*(double)w/XMAX;
y = IM_MIN+(IM_MAX-IM_MIN)*(double)h/YMAX;
flag = 0;
for(n=0;n<MAX_ITER;n++){
new_x = x*x-y*y+re_c;
new_y = 2*x*y+im_c;
if(sqrt(new_x*new_x+new_y*new_y)>alpha){
flag = 1;
break;
}
x = new_x;
y = new_y;
}
if(flag){
hue = 180.0*n/MAX_ITER;
gr_dot(w, h, hsl2rgb(hue,saturation,lightness));
}else{
gr_dot(w,h,BLACK);
}
}
}
gr_BMP(fname);
printf("save %s\n",fname);
gr_clear(BLACK);
}
return 0;
}
連番画像を作成したのち, 以下のコマンドによって動画を作成した.
$ ffmpeg -r 30 -i ./images/img%03d.bmp -vcodec libx264 -pix_fmt yuv420p -r 60 out.mp4
作成した動画は以下の通りである.
参考文献
[1] 奥村晴彦, [改定新版]C言語による標準アルゴリズム事典, 技術評論社.