分形造型（一）

一、引言

欧氏几何方法

以规则几何图形为研究对象，物体形状由方程来描述

欧氏几何主要描述工具

点，直线，平滑的曲线、平面

欧氏几何不能模拟自然景物

山脉、云、树、草、火、雾、波浪等

为解决复杂图形生成问题，分形造型应运而生

二、分形的历史

分形（Fractal）一词，是曼德布罗特创造出来的，其原意具有不规则、支离破碎等意义

分形几何学是一门以非规则几何形态为研究对象的几何学

分形理论创始人是美籍法国数学家B.B.曼德布罗特。他1924年11月20日生于波兰华沙。1985年获巴纳德（Barnard）奖章。1986年获富兰克林（Franklin）奖章。1967年发表于美国《科学》杂志上的“英国的海岸线有多长”的划时代论文，是他的分形思想萌芽的重要标志。1975年冬他创造了fractal一词。著作《分形：形状、机遇和维数》法文版于同年出版。

1967年B.B.Mandelbrot提出了“英国的海岸线有多长？”的问题，以1km为单位测量海岸线，得到的近似长度将短于1km的迂回曲折都忽略掉了，若以1m为单位测量，则能测出被忽略掉的迂回曲折，长度将变大，测量单位进一步变小，测得的长度将愈来愈大，这些愈来愈大的长度将趋近于一个确定值，这个极限值就是海岸线的长度。海岸线的长度是不确定的，长度依赖于测量单位。

三、分形几何与传统几何相比

1. 从整体上看，分形几何图形是处处不规则的。例如，海岸线和山川形状，从远距离观察，其形状是极不规则的。

2. 在不同尺度上，图形的规则性又是相同的。上述的海岸线和山川形状，从近距离观察，其局部形状又和整体形态相似，它们从整体到局部，都是自相似的。

四、分数维D定义

N相当于把原图复制出若干份

如果某图形是由把原图缩小为\(\frac 1 S\)的相似的N个图形所组成，有：

\[D=\frac {lnN}{lnS} \]

的关系成立，则指数D称为相似性维数，D可以是整数，也可以是分数。

五、经典分形

Koch曲线

import turtle as t


def koch(size, n):
    if n == 0:
        t.fd(size)
    else:
        for angle in [0, 60, -120, 60]:
            t.left(angle)
            koch(size / 3, n - 1)


if __name__ == '__main__':
    t.setup(1000, 1000)
    t.pen(speed=0, pendown=False, pencolor='blue')

    a, n = 400, 4
    t.goto(-a / 2, a / 2 / pow(3, 0.5))
    t.pd()
    for i in range(3):
        koch(a, n)
        t.right(120)
    t.ht()
    t.done()

原来的一条边变成4条边，所以N=4

现在的一个小三角的边长是原来的\(\frac 1 3\)，所以S=\(\frac 1 3\)

所以\(D=\frac {ln4}{ln3}=1.2618\)，D介于1和2之间

Cantor集

Cantor集是由不断去掉线段中间的三分之一开区间而得出。

#!/apps/bin/python
from tkinter import *


class Cantor(Frame):
    limit = 1

    def __init__(self, master=None):
        Frame.__init__(self, master)
        self.draw = Canvas(self, width=800, height=500)
        self.grid()
        self.createWidgets()

    def createWidgets(self):
        # self.quitButton = Button(self, text="QUIT", command=self.quit)
        # self.quitButton.pack(side=BOTTOM, fill=BOTH)
        self.draw.pack(side=LEFT)
        self.drawCanvas(100, 100, 700, 100)

    def drawCanvas(self, ax, ay, bx, by):
        self.draw.create_line(ax, ay, bx, by)
        if (bx - ax) > self.limit:
            cx = ax + (bx - ax) / 3
            cy = ay + 50
            dx = bx - (bx - ax) / 3 
            dy = by + 50
            ay = ay + 50 
            by = by + 50 
            self.drawCanvas(ax, ay, cx, cy)
            self.drawCanvas(dx, dy, bx, by)


app = Cantor()
app.master.title("Cantor (recursive)")
app.mainloop()

一条线变成2条线，所以D=2，长度是原来的\(\frac 1 3\)，所以S=\(\frac 1 3\)

所以

\[D=\frac {ln2}{ln3}=0.6309 \]

Sierpinski垫片

将一个等边三角形四等分，得到四个小等边三角形，去掉中间一个，保留它的三条边，将剩下的三个小等边三角形在分别进行四等分，并分别去掉中间一个，保留它们的边。

"""
    功能：绘制谢尔宾斯基三角形
    环境：python3.7
    日期：2019/1/14 21:49
    作者：指尖魔法师
    版本：1.0
"""
import turtle as t


def sanjiaoxing(san):
    """
    传入三个点坐标，绘制三角形
    """
    t.penup()
    t.goto(san[0])
    t.pendown()
    t.goto(san[1])
    t.goto(san[2])
    t.goto(san[0])


def get_mid(a, b):
    """
    计算返回2个点的中间点坐标
    """
    x = (a[0] + b[0]) / 2
    y = (a[1] + b[1]) / 2
    return [x, y]


def draw_san(size, i):
    """
    绘制谢尔宾斯基三角形函数
    :param size: 三个点坐标列表
    :param i: 递归次数
    """
    # 绘制三角形
    sanjiaoxing(size)
    if i > 0:
        # 绘制左边小三角形
        size2 = [size[0], get_mid(size[0], size[1]), get_mid(size[0], size[2])]
        draw_san(size2, i - 1)

        # 绘制上边的小三角形
        size3 = [get_mid(size[0], size[2]), get_mid(size[1], size[2]), size[2]]
        draw_san(size3, i - 1)

        # 绘制右边的小三角形
        size4 = [get_mid(size[0], size[1]), size[1], get_mid(size[1], size[2])]
        draw_san(size4, i - 1)


def main():
    """
    主函数
    """
    # 打印图形标题
    t.penup()
    t.left(90)
    t.forward(350)
    t.pendown()
    t.write("谢尔宾斯基三角形", False, align="center", font=("宋体", 20, "normal"))
    t.speed(5)

    # 初始三角形坐标
    points = [[-200, 0], [200, 0], [0, 300]]
    # 递归5次
    count = 5
    # 调用绘制谢尔宾斯基三角形函数
    draw_san(points, count)

    t.exitonclick()


if __name__ == '__main__':
    main()

\[D=\frac {ln3} {ln2}=1.5850 \]

Sierpinski海绵

import javax.swing.*;
import java.awt.*;
import java.util.stream.IntStream;

public class Sponge extends JFrame implements Runnable {
    public static void main(String[] args) {
        Sponge mf = new Sponge();
        mf.initUI();
    }

    Point[] pn = new Point[20];

    public void paint(Graphics g) {
        super.paint(g);
        //正方形的顶点坐标
        int x = 100, y = 500;
        int d = 120, dx = 81, dy = 60;

        Point p0 = new Point(x, y);
        drawBottomLevel(g, p0, d, dx, dy);
        Point p1 = new Point(x, y - d);
        drawMidLevel(g, p1, d, dx, dy);
        Point p2 = new Point(x, y - 2 * d);
        drawTopLevel(g, p2, d, dx, dy);

        IntStream.range(0, pn.length).forEach(i -> {
            Point ptem = pn[i];
            if (ptem != null) {
                Point ptem0 = new Point(pn[i].x, pn[i].y + 2 * d / 3);
                drawBottomLevel(g, ptem0, d / 3, dx / 3, dy / 3);

                Point ptem1 = new Point(pn[i].x, pn[i].y + d / 3);
                drawMidLevel(g, ptem1, d / 3, dx / 3, dy / 3);

                Point ptem2 = new Point(pn[i].x, pn[i].y);
                drawTopLevel(g, ptem2, d / 3, dx / 3, dy / 3);
            }
        });
    }

    //根据起点p0得到另外6个点
    public Point[] getPointByP0(Point p0, int d, int dx, int dy) {

        Point p1 = new Point(p0.x + dx, p0.y - dy);
        Point p2 = new Point(p0.x + dx + d, p0.y - dy);
        Point p3 = new Point(p0.x + d, p0.y);
        Point p4 = new Point(p0.x + d, p0.y + d);
        Point p5 = new Point(p0.x, p0.y + d);
        Point p6 = new Point(p0.x + dx + d, p0.y + d - dy);

        Point[] ps = new Point[7];
        ps[0] = p0;
        ps[1] = p1;
        ps[2] = p2;
        ps[3] = p3;
        ps[4] = p4;
        ps[5] = p5;
        ps[6] = p6;

        return ps;
    }

    int count = 0;//递归次数

    public void draw(Graphics g, Point p0, Point p1, Point p2, Point p3, Point p4, Point p5, Point p6) {
        if (count < 20) {
            pn[count] = p0;
        }
        count++;

        Graphics2D gD = (Graphics2D) g;

        //填充面的颜色
        Polygon pon1 = new Polygon();
        //填充第一个可见面
        pon1.addPoint(p0.x, p0.y);
        pon1.addPoint(p1.x, p1.y);
        pon1.addPoint(p2.x, p2.y);
        pon1.addPoint(p3.x, p3.y);
        gD.setColor(new Color(0xFFC0CB));
        gD.fillPolygon(pon1);
        //填充第二个可见面
        Polygon pon2 = new Polygon();
        pon2.addPoint(p0.x, p0.y);
        pon2.addPoint(p3.x, p3.y);
        pon2.addPoint(p4.x, p4.y);
        pon2.addPoint(p5.x, p5.y);
        gD.setColor(new Color(0x39c5bb));
        gD.fillPolygon(pon2);
        //填充第三个可见面
        Polygon pon3 = new Polygon();
        pon3.addPoint(p2.x, p2.y);
        pon3.addPoint(p3.x, p3.y);
        pon3.addPoint(p4.x, p4.y);
        pon3.addPoint(p6.x, p6.y);
        gD.setColor(new Color(0xFFA500));
        gD.fillPolygon(pon3);
    }

    public void drawBottomLevel(Graphics g, Point p0, int d, int dx, int dy) {
        Point[] ps1 = getPointByP0(p0, d, dx, dy);
        Point[] ps2 = getPointByP0(ps1[3], d, dx, dy);
        Point[] ps3 = getPointByP0(ps2[3], d, dx, dy);
        Point[] ps4 = getPointByP0(ps3[1], d, dx, dy);
        Point[] ps5 = getPointByP0(ps4[1], d, dx, dy);
        Point[] ps6 = getPointByP0(ps1[1], d, dx, dy);
        Point[] ps7 = getPointByP0(ps6[1], d, dx, dy);
        Point[] ps8 = getPointByP0(ps7[3], d, dx, dy);

        draw(g, ps7[0], ps7[1], ps7[2], ps7[3], ps7[4], ps7[5], ps7[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps8[0], ps8[1], ps8[2], ps8[3], ps8[4], ps8[5], ps8[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps5[0], ps5[1], ps5[2], ps5[3], ps5[4], ps5[5], ps5[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps4[0], ps4[1], ps4[2], ps4[3], ps4[4], ps4[5], ps4[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps6[0], ps6[1], ps6[2], ps6[3], ps6[4], ps6[5], ps6[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps1[0], ps1[1], ps1[2], ps1[3], ps1[4], ps1[5], ps1[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps2[0], ps2[1], ps2[2], ps2[3], ps2[4], ps2[5], ps2[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
        draw(g, ps3[0], ps3[1], ps3[2], ps3[3], ps3[4], ps3[5], ps3[6]);
        try {
            Thread.sleep(100);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
    }


    public void drawMidLevel(Graphics g, Point p0, int d, int dx, int dy) {
        Point[] ps1 = getPointByP0(p0, d, dx, dy);
        Point[] ps2 = getPointByP0(ps1[3], d, dx, dy);
        Point[] ps3 = getPointByP0(ps2[3], d, dx, dy);
        Point[] ps4 = getPointByP0(ps3[1], d, dx, dy);
        Point[] ps5 = getPointByP0(ps4[1], d, dx, dy);
        Point[] ps6 = getPointByP0(ps1[1], d, dx, dy);
        Point[] ps7 = getPointByP0(ps6[1], d, dx, dy);
        Point[] ps8 = getPointByP0(ps7[3], d, dx, dy);
        
        draw(g, ps7[0], ps7[1], ps7[2], ps7[3], ps7[4], ps7[5], ps7[6]);
        draw(g, ps5[0], ps5[1], ps5[2], ps5[3], ps5[4], ps5[5], ps5[6]);
        draw(g, ps1[0], ps1[1], ps1[2], ps1[3], ps1[4], ps1[5], ps1[6]);
        draw(g, ps3[0], ps3[1], ps3[2], ps3[3], ps3[4], ps3[5], ps3[6]);
    }


    public void drawTopLevel(Graphics g, Point p0, int d, int dx, int dy) {
        Point[] ps1 = getPointByP0(p0, d, dx, dy);
        Point[] ps2 = getPointByP0(ps1[3], d, dx, dy);
        Point[] ps3 = getPointByP0(ps2[3], d, dx, dy);
        Point[] ps4 = getPointByP0(ps3[1], d, dx, dy);
        Point[] ps5 = getPointByP0(ps4[1], d, dx, dy);
        Point[] ps6 = getPointByP0(ps1[1], d, dx, dy);
        Point[] ps7 = getPointByP0(ps6[1], d, dx, dy);
        Point[] ps8 = getPointByP0(ps7[3], d, dx, dy);

        draw(g, ps7[0], ps7[1], ps7[2], ps7[3], ps7[4], ps7[5], ps7[6]);
        draw(g, ps8[0], ps8[1], ps8[2], ps8[3], ps8[4], ps8[5], ps8[6]);
        draw(g, ps5[0], ps5[1], ps5[2], ps5[3], ps5[4], ps5[5], ps5[6]);
        draw(g, ps4[0], ps4[1], ps4[2], ps4[3], ps4[4], ps4[5], ps4[6]);
        draw(g, ps6[0], ps6[1], ps6[2], ps6[3], ps6[4], ps6[5], ps6[6]);
        draw(g, ps1[0], ps1[1], ps1[2], ps1[3], ps1[4], ps1[5], ps1[6]);
        draw(g, ps2[0], ps2[1], ps2[2], ps2[3], ps2[4], ps2[5], ps2[6]);
        draw(g, ps3[0], ps3[1], ps3[2], ps3[3], ps3[4], ps3[5], ps3[6]);
    }


    public void initUI() {
        this.setTitle("门格海绵");
        this.setSize(1000, 800);
        this.setDefaultCloseOperation(3);
        this.setResizable(false);
        this.setLocationRelativeTo(null);
        this.setVisible(true);
    }

    @Override
    public void run() {
    }
}

Julia集

Python实现

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
import matplotlib.backend_bases as mbb
import time


class MandelBrot():
    def __init__(self, x0, x1, y0, y1, n):
        self.oriAxis = np.array([x0, x1, y0, y1])  # 初始坐标
        self.axis = self.oriAxis
        self.nx, self.ny, self.nMax = n, n, n  # x,y方向的网格划分个数
        self.nIter = 100  # 迭代次数
        self.n0 = 0  # 预迭代次数
        self.z = genZ(self.oriAxis, self.nx, self.ny)
        self.DrawMandelbrot()

    def DrawMandelbrot(self):
        mBrot = getJulia(self.z, self.z, self.nIter)

        self.fig, ax = plt.subplots()
        plt.imshow(mBrot, cmap=cm.jet, extent=self.axis)
        plt.gca().set_axis_off()

        self.fig.canvas.mpl_disconnect(self.fig.canvas.manager.key_press_handler_id)
        self.fig.canvas.mpl_connect('button_press_event', self.OnMouse)
        self.fig.canvas.mpl_connect('button_release_event', self.OnRelease)
        self.fig.canvas.mpl_connect('scroll_event', self.OnScroll)

        plt.show()

    def DrawJulia(self, c0):

        z = genZ([-2, 2, -2, 2], 800, 800)
        julia = getJulia(z, c0, self.nIter)

        jFig, jAx = plt.subplots()
        plt.cla()
        plt.imshow(julia, cmap=cm.jet, extent=self.axis)
        plt.gca().set_axis_off()
        plt.show()
        jFig.canvas.draw_idle()

    # 滚轮缩放
    def OnScroll(self, evt):
        x0, y0 = evt.xdata, evt.ydata

        if evt.button == "up":
            self.axis = (self.axis + [x0, x0, y0, y0]) / 2
        elif evt.button == 'down':
            self.axis = 2 * self.axis - [x0, x0, y0, y0]

        z = genZ(self.axis, self.nx, self.ny)
        mBrot = getJulia(z, z, self.nIter)
        plt.cla()
        plt.imshow(mBrot, cmap=cm.jet, extent=self.axis)
        plt.gca().set_axis_off()

        mBrot[mBrot < 1] == self.n0 + self.nIter
        self.n0 = int(np.min(mBrot))
        self.fig.canvas.draw_idle()
        pass

    def OnMouse(self, evt):
        self.xStart = evt.xdata
        self.yStart = evt.ydata
        self.fig.canvas.draw_idle()

    def OnRelease(self, evt):
        x0, y0, x1, y1 = self.xStart, self.yStart, evt.xdata, evt.ydata
        if evt.button == mbb.MouseButton.LEFT:
            self.DrawJulia(x1 + y1 * 1j)  # 如果释放的是左键，那么就绘制Julia集并返回
            return
        # 右键拖动，可以对Mandelbrot集进行真实的放大
        self.axis = np.array([min(x0, x1), max(x0, x1),
                              min(y0, y1), max(y0, y1)])

        nxny = self.axis[[1, 3]] - self.axis[[0, 2]]
        self.nx, self.ny = (nxny / max(nxny) * self.nMax).astype(int)
        z = genZ(self.axis, self.nx, self.ny)

        n = 100  # n为迭代次数
        mBrot = getJulia(z, z, n)
        plt.cla()
        plt.imshow(mBrot, cmap=cm.jet, extent=self.axis)
        plt.gca().set_axis_off()

        mBrot[mBrot < 1] == self.n0 + n
        self.n0 = int(np.min(mBrot))
        self.fig.canvas.draw_idle()


def genZ(axis, nx, ny):
    x0, x1, y0, y1 = axis
    x = np.linspace(x0, x1, nx)
    y = np.linspace(y0, y1, ny)
    real, img = np.meshgrid(x, y)
    z = real + img * 1j
    return z


def getJulia(z, c, n, n0=0, m=2):
    t = time.time()
    c = np.zeros_like(z) + c
    out = abs(z)
    for _ in range(n0):
        z = z * z + c
    for i in range(n0, n0 + n):
        absz = abs(z)
        z[absz > m] = 0
        c[absz > m] = 0
        out[absz > m] = i
        z = z * z + c
    print("time:", time.time() - t)
    return out


if __name__ == "__main__":
    x, y = 0, 0
    brot = MandelBrot(-2, 1, -1.5, 1.5, 1000)

介绍

在复平面上，水平的轴线代表实数，垂直的轴线代表虚数。每个Julia集合（有无限多个点）都决定一个常数C，它是一个复数。在复平面上任意取一个点，其值是复数Z。将其代入下面方程中进行反复迭代运算：

\[z_{n+1}=z_{n}^2+C \]

Julia集的生成过程

假设显示器的分辨率是\(a*b\)点，可显示的颜色为\(K+1\)种，分别以数字\(0,1,2,\cdots,K\)表示，且以0表示黑色。

从\(Z_i\)到\(Z_{i+1}\)的迭代过程就是

\[x_{i+1}=x_i^2-y_i^2+p\\y_{i+1}=2x_iy_i+q \]

步骤0：

选定参数\(x_{min}=y_{min}=-1.5, x_{max}=y_{max}=1.5\)，又取M=100

令

\[\Delta x= \frac {x_{max}-x_{min}}{a-1}\\\Delta y=\frac {y_{max}-y_{min}}{b-1} \]

对所有的点\((n_x,n_y),n_x=0,1,\cdots,a-1\)及\(n_y=0,1,\cdots,b-1\)完成下面的循环

步骤1：

令

\[x_0=x_{min}+n_x\cdot \Delta x\\y_0=y_{min}+n_y\cdot \Delta y\\k=0 \]

步骤2：

由迭代过程从\((x_i,y_i)\)算出\((x_{i+1},y_{i+1})\)，并计算\(k:=k+1\)

步骤3：

计算\(r=x_i^2+y_i^2\)。

如果\(r>M\)，则选择颜色\(k\)，转至步骤4；

如果\(k=K\)，则选择颜色\(0\)，转至步骤4；

如果\(r≤M\)，且\(k<K\)，则转至步骤2。

步骤4：

对点\((n_x,n_y)\)显示颜色\(k\)并转至下一点，再从头作步骤1