C#でマンデルブロ集合を描く

2019年に発行された『(新版)独習C#』の初版第4刷を購入し再び読んで結構忘れてしまっていたのを復習し、フラクタルについては以下の本で学び、マンデルブロ集合を描くコードを書いたので記録に残しておこう。

フラクタル紀行』芹沢浩(著)、森北出版株式会社
   ISBN 4-627-01560-7
  1993年出版

マンデルブロ集合の発散しない領域の中で周期的に変動する領域が色々存在するということを本書で知った。単に発散するかしないかだけではなかった。

f:id:willwealth:20201025221738j:plain
マンデルブロ集合には周期の異なる領域がある

周期Nの領域の値をジュリア集合を描くときの写像の定数項として使用すると周期Nのジュリア集合が描かれるというのが興味深かった。例えば赤の領域の値を定数項として選ぶと赤のジュリア集合が描かれます。うまく説明できないので興味のある方は上記の本をお読みください。コブとコブが接触している脇の箇所が面白いというアドバイスもあり、周期Nが大きくなるカオス的な領域は色々と複素数的な意匠を楽しめます。
さて、この本の中では私の知らない Turbo C のコードになっていたがC言語なのでロジックはC#に写しやすかったが、フォームアプリケーションは厄介なので、Visual Studioで自動生成されたからコードを流用した。
適当に作ったフォームアプリケーションから、

  • Form1.cs
  • Form1.Designer.cs
  • Program.cs

の3つのファイルからコピペしてまず骨格部分を作り、描画のコードを追加した。
Pythonだと標準で複素数が使えるが、Numericsを使えばC#でも使えるということを知らなかったので、今回はNumericsを使わず複素数クラス(割り算なし)を作成した。

using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading.Tasks;
using System.Windows.Forms;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Drawing.Imaging;

namespace WindowsFormsApplication1
{

    //
    //  Complex Number
    //
    public class Complex
    {
        public double X;
        public double Y;

        public Complex(double x, double y)
        {
            this.X = x;
            this.Y = y;
        }

        public double AbsSqr()
        {
            return (this.X * this.X + this.Y * this.Y);
        }

        public static Complex operator +(Complex a, Complex b)
        {
            return new Complex(a.X + b.X, a.Y + b.Y);
        }
        public static Complex operator -(Complex a, Complex b)
        {
            return new Complex(a.X - b.X, a.Y - b.Y);
        }

        public static Complex operator *(Complex a, Complex b)
        {
            return new Complex(a.X * b.X - a.Y * b.Y, a.X * b.Y + a.Y * b.X);
        }
    }

    partial class Form1
    {
        /// <summary>
        /// 必要なデザイナー変数です。
        /// </summary>
        private System.ComponentModel.IContainer components = null;

        /// <summary>
        /// 使用中のリソースをすべてクリーンアップします。
        /// </summary>
        /// <param name="disposing">マネージ リソースが破棄される場合 true、破棄されない場合は false です。</param>
        protected override void Dispose(bool disposing)
        {
            if (disposing && (components != null))
            {
                components.Dispose();
            }
            base.Dispose(disposing);
        }

        #region Windows フォーム デザイナーで生成されたコード

        /// <summary>
        /// デザイナー サポートに必要なメソッドです。このメソッドの内容を
        /// コード エディターで変更しないでください。
        /// </summary>
        private void InitializeComponent()
        {
            this.components = new System.ComponentModel.Container();
            this.AutoScaleMode = System.Windows.Forms.AutoScaleMode.Font;
            this.Text = "Sample Program";
        }

        #endregion
    }

    public partial class Form1 : Form
    {
        Image image;
        int width  = 800;
        int height = 500;

        public Complex Pixel2Complex(int i, int j, Complex c, double unit)
        {
            double a, b;
            a = c.X + unit * i;
            b = c.Y - unit * j;
            return new Complex(a, b);
        }

        void DrawFractal()
        {
            Complex leftUpPosition = new Complex(-2.35,1.25);
            double pixelUnit = 0.005;
            int iteration = 200;

            Complex[] z = new Complex[iteration + 1];
            Complex zc;
            Brush [] palet =
                {
                    new SolidBrush(Color.Black),                   // 0
                    new SolidBrush(Color.Blue),                    // 1
                    new SolidBrush(Color.FromArgb(140, 255, 255)), // 2
                    new SolidBrush(Color.FromArgb(255, 170, 170)), // 3
                    new SolidBrush(Color.Red),                     // 4
                    new SolidBrush(Color.Purple),                  // 5
                    new SolidBrush(Color.Yellow),                  // 6
                    new SolidBrush(Color.YellowGreen),             // 7
                    new SolidBrush(Color.FromArgb(255, 0, 170)),   // 8
                    new SolidBrush(Color.SeaGreen),                // 9
                    new SolidBrush(Color.FromArgb(255, 200, 0)),   // 10 
                    new SolidBrush(Color.Silver),                  // 11 
                    new SolidBrush(Color.FromArgb(255, 160, 0)),   // 12
                    new SolidBrush(Color.SpringGreen),             // 13 
                    new SolidBrush(Color.FromArgb(255, 200, 255)), // 14
                    new SolidBrush(Color.FromArgb(0, 170, 170)),   // 15
                    new SolidBrush(Color.FromArgb(0, 250, 0)),     // 16
                    new SolidBrush(Color.FromArgb(0, 200, 0)),     // 17
                    new SolidBrush(Color.FromArgb(0, 150, 0)),     // 18
                    new SolidBrush(Color.FromArgb(0, 100, 0)),     // 19
                    new SolidBrush(Color.FromArgb(0, 50, 0)),      // 20
                    new SolidBrush(Color.White)                    // 21
                };

            int MaxCycle = palet.Length - 1;

            Graphics g = Graphics.FromImage(image);

            int imax = image.Width;
            int jmax = image.Height;
            double n_progress = imax / 100.0;
            double s;
            int i;
            int j;
            int k;

            for (i = 0; i < imax; i++)
            {
                for (j = 0; j < jmax; j++)
                {
                    z[0] = Pixel2Complex(i, j, leftUpPosition, pixelUnit);
                    zc = z[0];

                    for (k = 0; k < iteration; k++)
                    {
                        z[k + 1] = z[k]*z[k] + zc;

                        s = z[k].AbsSqr();
                        if ( s > 4.0 ) {
                            g.FillRectangle(palet[0], i, j, 1, 1);
                            break;
                        }
                    }

                    if (k == iteration)
                    {
                        int n;
                        for (n = 1; n < MaxCycle; n++)
                        {
                            s = ( z[iteration] - z[iteration - n]).AbsSqr();
                            if ( s < pixelUnit )
                            {
                                g.FillRectangle(palet[n], i, j, 1, 1);
                                break;
                            }
                        }
                        if (n == MaxCycle)
                        {
                            g.FillRectangle(palet[MaxCycle], i, j, 1, 1);
                        }
                    }
                }
                // Disp Progress
                Console.Write("{0,7:d0}%", (int)((double)i/n_progress) + 1);
                Console.SetCursorPosition(0, Console.CursorTop); // カーソル位置を初期化
            }
            Console.Write("\n");
        }

        public Form1()
        {
            InitializeComponent();

            image = new Bitmap(width,height);
            DrawFractal();
        }

        protected override void OnPaint(PaintEventArgs e)
        {
            base.OnPaint(e);
            e.Graphics.DrawImage(image, 0, 0);
        }
    }

    static class Program
    {
        /// <summary>
        /// アプリケーションのメイン エントリ ポイントです。
        /// </summary>
        [STAThread]
        static void Main()
        {
            Application.EnableVisualStyles();
            Application.SetCompatibleTextRenderingDefault(false);
            Application.Run(new Form1());
        }
    }
}