Numpy Python分形盒计数 – 分形维度

在本文中，我们将介绍Numpy Python中的分形盒计数方法，以及如何利用它来计算分形维度。分形维度是描述分形图形中复杂程度的一种度量。它通常是非整数，反映了图形的自相似性和自我重复性。

阅读更多：Numpy 教程

分形盒计数原理

分形盒计数方法是一种常用的计算分形维度的方法，它是基于一种简单的思想：用不同大小的盒子去覆盖分形图形。具体的步骤如下：

将分形图形覆盖上粗到细的等间距网格
在每个网格中判断是否覆盖了分形图形的部分
按照网格的大小来计算分形维度

下面是一个例子，我们先生成一个古典分形图形：科赫曲线，然后用分形盒计数法来计算它的分形维度。

import numpy as np
import matplotlib.pyplot as plt

def koch(n):
    if n == 0:
        X = [0, 1]
        Y = [0, 0]
    else:
        X, Y = koch(n-1)
        X1 = (2*X[0]+X[1])/3
        Y1 = (2*Y[0]+Y[1])/3
        X3 = (2*X[1]+X[0])/3
        Y3 = (2*Y[1]+Y[0])/3
        X2 = (X1+X3)/2 - np.sqrt(3)*(Y3-Y1)/6
        Y2 = (Y1+Y3)/2 - np.sqrt(3)*(X1-X3)/6
        X = [X[0], X1, X2, X3, X[1]]
        Y = [Y[0], Y1, Y2, Y3, Y[1]]
    return X, Y

X, Y = koch(4)

plt.plot(X, Y, 'k')
plt.xticks([])
plt.yticks([])
plt.show()

我们得到了一个科赫曲线的图形

我们定义一个函数：boxcount，用来计算分形盒计数：

def boxcount(x, y, boxsize=None):
    if boxsize == None:
        boxsize = (max(x)-min(x))/30

    xlow = np.min(np.floor(x/boxsize))*boxsize
    xhigh = np.max(np.ceil(x/boxsize))*boxsize
    ylow = np.min(np.floor(y/boxsize))*boxsize
    yhigh = np.max(np.ceil(y/boxsize))*boxsize

    xg, yg = np.meshgrid(np.arange(xlow, xhigh+boxsize, boxsize),
                         np.arange(ylow, yhigh+boxsize, boxsize))

    count = np.zeros_like(xg)

    for xi, yi in zip(x, y):
        xidx = int((xi-xlow)/boxsize)
        yidx = int((yi-ylow)/boxsize)
        count[yidx, xidx] = 1

    return xg, yg, count

这个函数接收分形图形的x, y坐标以及盒子的大小boxsize，然后返回由盒子组成的网格和对应的盒子中是否覆盖了图形的矩阵。

我们可以看到这里的盒子大小是按照max(x)-min(x)/30来计算的，默认大小为30。

然后我们调用boxcount函数来计算盒子的数目和盒子大小的对数，用来计算分形维度。

xg, yg, count = boxcount(X, Y, boxsize=None)

plt.imshow(count, cmap='gray_r')
plt.xticks([])
plt.yticks([])
plt.show()

boxsizes = np.logspace(-5, np.log10(np.max(xg-yg)), num=20, base=10.0)
boxcounts = []
for boxsize in boxsizes:
    xg, yg, count = boxcount(X, Y, boxsize=boxsize)
    boxcounts.append(np.sum(count > 0))

boxcounts = np.array(boxcounts)

plt.loglog(boxsizes, boxcounts, 'bo')
plt.xlabel('Box size (log scale)')
plt.ylabel('Number of boxes (log scale)')
plt.show()

plt.plot(np.log10(boxsizes), np.log10(boxcounts), 'bo')
plt.xlabel('Log(Box size)')
plt.ylabel('Log(Number of boxes)')
plt.show()

frac_dim, intercept = np.polyfit(np.log10(boxsizes), np.log10(boxcounts), 1)

print("Fractal dimension:", frac_dim)