简介：这套MATLAB工具资源专门面向需要掌握分形几何与盒维数计算的科研人员和工程师，提供1D、2D、3D三种空间下的分形盒维数计算实现，覆盖科赫曲线、曼德勃罗集、Sierpinski海绵等典型分形对象。包内共27个文件，以.m脚本、png结果图、gif演示动画为主，并有html说明文档与示例图片，整体压缩包约3.15MB，结构清晰，便于直接对照学习和二次开发。目前已有683人次浏览学习。借助box-counting核心算法和MATLAB实现，可完成对曲线、图像及三维结构的自相似性度量，并通过对数-对数图拟合得到维数D；脚本中包含了盒尺寸设置、计数统计、线性回归拟合等完整步骤，也提供了DLA、Apollonian垫片等演示范例。无论用于课程实验、论文复现还是工程中的复杂度分析，这套代码都能帮助读者快速上手，理解分形维数的计算原理与不同维度间的差异。

1. 分形盒维数：不是测长度，而是数盒子

一条不平滑的曲线，为什么用越小的尺子量，长度反而越长？分形盒维数（box-counting dimension）正是把这种“越测越细、数值越大”的现象变成一个稳定的斜率：用边长为 s 的盒子去覆盖目标物，统计非空盒子数 N(s)，当 s 不断减小时，N(s) 与 s 的负 D 次方成正比，这个 D 就是分形盒维数。它不是测量“占了多少面积”，而是测量“填充空间的能力”。我最近整理了一份能同时处理一维曲线、二维图像和三维体数据的分形盒维数计算代码，实际跑下来发现，真正影响结果准确度的不是核心公式，而是尺度范围怎么选、边界怎么界定、二值化怎么统一。这篇文章把实现细节和踩坑记录一起讲透，适合做信号特征提取、纹理识别、三维医学图像分析的人参考。

2. 一维信号：从曲线到网格计数的第一课

2.1 为什么一维最容易入门：线段的维度意义

所谓“一维信号”，比如时间序列、地震波形、股价曲线，本质上是一个一维数组 y[i]。单独看数组长度，它只有一维；但把它画在二维平面上，横轴是下标 i，纵轴是信号值 y，就变成了一条平面曲线。盒维数描述的就是这条曲线“霸占二维空间”的程度。

光滑的直线或正弦曲线，只要没有无限细的起伏，它的盒维数理论上等于 1；而白噪声、锯齿状的分形曲线，因为处处不光滑，填充空间的能力更强，盒维数就会介于 1 和 2 之间。这个区间很容易验证，所以我把一维作为理解盒维数的第一个实验对象。

盒子法的思想很朴素：用边长为 s 的正方形去铺满这条曲线所在的空间，统计有多少个正方形里至少包含一个曲线上的点。如果把 s 从大到小变化，非空盒子数 N(s) 会按幂律增长。取对数后：

log N(s) ≈ D · log(1/s) + const

斜率 D 就是盒维数。关键在于，这里的“点集”怎么从原始信号生成，以及 s 的取值区间怎么定。

2.2 实现步骤：归一化、网格划分、线性拟合

我常用的做法是把一维信号变成点集：横轴用采样点下标，纵轴用信号值，然后归一化到单位正方形里，再对不同尺度的盒子做网格统计。核心代码可以浓缩成下面这个函数：

import numpy as np def box_dimension_1d(signal, box_sizes=None, fit_ratio=0.6): # signal: 一维数组，长度 n x = np.arange(len(signal)) y = np.asarray(signal, dtype=float) # 归一化到单位正方形，避免横纵量纲不一致 x_norm = (x - x.min()) / (x.max() - x.min() + 1e-10) y_norm = (y - y.min()) / (y.max() - y.min() + 1e-10) # 盒子边长序列：默认从大到小等比递减 if box_sizes is None: min_size = max(2.0 / len(signal), 1e-3) box_sizes = np.exp(np.linspace(np.log(min_size), 0, 20))[::-1] counts = [] for s in box_sizes: # 用 floor 把坐标映射到网格索引 gx = np.floor(x_norm / s).astype(int) gy = np.floor(y_norm / s).astype(int) # 合并成点对去重，统计非空盒子数 pairs = set(zip(gx.tolist(), gy.tolist())) counts.append(len(pairs)) log_s = np.log(1.0 / np.asarray(box_sizes)) log_N = np.log(np.asarray(counts)) # 只取中间 60% 的点做线性拟合 n = len(log_s) start = int(n * (1 - fit_ratio) / 2) end = int(n * (1 + fit_ratio) / 2) slope, intercept = np.polyfit(log_s[start:end], log_N[start:end], 1) return slope

逻辑分为四步：归一化、生成尺度序列、统计非空盒子、线性拟合。归一化这一步很关键，如果直接使用原始像素坐标或信号幅值，横纵轴单位不同，盒子就不是正方形，维数自然失去意义。

参数方面，box_sizes默认从min_size = 2 / len(signal)开始，到 1 为止，取 20 个对数等分点。min_size的含义是“至少大于两个采样点的间距”，因为更小的盒子会导致每个点各占一个格子，曲线内部结构被当作噪声处理。fit_ratio=0.6表示只保留 log-log 散点中间 60% 的部分做拟合，两端通常因为尺度上下限而偏离直线。

2.3 参数选择：尺度范围决定结果可信度

盒维数计算的“玄学”大部分集中在尺度范围的选择。同样一条曲线，如果把最小盒子边长设成 1/1000，拟合结果可能接近 1.8；如果设成 1/50，结果可能只有 1.2。原因在于，尺度下限低于曲线的采样分辨率时，曲线看起来像一堆孤立的点，维数会虚高；尺度上限太大时，整个曲线被装进很少的几个盒子，增长趋势饱和，维数又会偏低。

我一般的做法是：先画出 log-log 散点图，看中段有没有明显的线性区间。最小盒子边长的下限取平均点间距的 2～3 倍；上限取目标包围盒长边的 1/5 左右，因为再大了基本就是一两个盒子覆盖全图。如果信号长度很长，比如几万个点，可以多取一些尺度，比如 30～40 个，让拟合更稳。

还可以用一个粗糙的验证方法：对一段纯正弦信号，盒维数应当接近 1.0；对一段随机噪声，盒维数应当接近 2.0。如果代码算出来偏离这两个值太多，优先检查归一化和尺度序列，而不是怀疑算法。

3. 二维图像：二值化、网格计算与边界效应

3.1 灰度图到二值图：盒维数到底是在测谁

二维盒维数最常见的输入是图像，但很少有图像天生只有黑白两种颜色。如果直接拿灰度图做盒计数，你会很快发现阈值一变，维数跟着剧烈变化。我曾在某图像处理 Demo 里踩过这个坑：同一张分形图，阈值从 128 改成 100，维数从 1.3 跳到 1.7，后来才意识到，盒维数真正稳定描述的是“形状”或“目标集合”，而不是灰度值本身。

因此，计算二维图像盒维数前，必须先想清楚你要测的是前景形状，还是灰度表面纹理。如果测前景形状，做法是二值化；如果测纹理，应该用灰度表面法或差分盒维数，而不是简单地把灰度值当作是否非空的依据。这篇文章里聚焦二值图的盒维数。

二值化可以选固定阈值、Otsu 大津法或自适应阈值。对大多数图像分割场景，Otsu 是合理默认值；但如果目标像素占比极低，比如一张图里只有一个细长裂缝，Otsu 可能失效，需要手动设定阈值或基于局部背景调整。

3.2 二维网格计数：用 set 消除重复计数

输入是一张二值图，背景为 0，前景为 1。我们需要把所有前景像素的坐标提取出来，归一化到单位正方形，然后像一维情况那样做网格统计。注意一个问题：如果图像不是正方形，直接除以各自的宽高会把像素坐标拉伸成不同比例，网格会变成矩形，盒子不再是正方形，维数会被干扰。

我通常把长边归一化为 1，短边按相同比例缩放。这样盒子仍然是正方形，只是短边方向会有部分空白区域，但这并不会影响前景目标的覆盖计数，因为空区域本来就没有点。

代码实现如下：

import numpy as np def box_dimension_2d(binary_image, box_sizes=None, fit_ratio=0.6): # binary_image: shape (H, W)，前景为 1 H, W = binary_image.shape ys, xs = np.nonzero(binary_image) if len(xs) == 0: return np.nan # 按长边归一化，保持盒子为正方形 max_len = max(H, W) x_norm = xs / max_len y_norm = ys / max_len if box_sizes is None: # 最小盒子边长取 2 个像素的长度 min_size = 2.0 / max_len box_sizes = np.exp(np.linspace(np.log(min_size), 0, 20))[::-1] counts = [] for s in box_sizes: gx = np.floor(x_norm / s).astype(int) gy = np.floor(y_norm / s).astype(int) counts.append(len(set(zip(gx.tolist(), gy.tolist())))) log_s = np.log(1.0 / np.asarray(box_sizes)) log_N = np.log(np.asarray(counts)) n = len(log_s) start = int(n * (1 - fit_ratio) / 2) end = int(n * (1 + fit_ratio) / 2) slope, _ = np.polyfit(log_s[start:end], log_N[start:end], 1) return slope

这段代码的核心与一维版本完全一致，区别在于坐标来源不同。np.nonzero返回的是二维数组中所有非零元素的行列坐标，ys是行号，xs是列号，对应图像坐标的 y、x。归一化时除以max_len，而不是分别除以 H 和 W，是为了保证两个轴使用同一个缩放系数。

如果你有大量图像要批量计算，不建议在这段循环里频繁调用set(zip(...))，因为 Python 的 set 在处理几十万点时会偏慢。更快的做法是用np.unique(coords, axis=0)或np.ravel_multi_index把二维坐标映射成一维 ID 后去重。后面三维部分会展示更省内存的写法。

3.3 边界效应：小尺度下的像素化误差

二维盒维数的第一个大坑是“点到度”的边界效应。当盒子边长 s 小于等于图像的一个像素时，每个前景像素都独占一个盒子，N(s) 不再随 s 减小而增长，log-log 曲线左端会明显向下弯曲，拟合出的斜率偏低。

要判断小尺度是否越过分辨率极限，可以打印每个 s 对应的 counts：如果出现连续几个 s 的非空盒子数几乎不变，说明已经进入像素化区间。我通常把最小尺度设成至少 2 个像素的长度，并且在进行线性拟合时只取中间 50%～70% 的尺度点。

另一个容易忽略的问题是图像边缘。如果前景目标刚好横跨网格线，它可能被同时分进两个盒子，造成非空盒计数偏大。但这种边界误差在每个尺度下都会出现，且方向一致，对斜率的影响相对有限。真正需要警惕的是网格原点偏移：同样的目标在原点位置平移几个像素，计数结果可能不同。解决方法是把坐标统一减掉一个随机偏移，或者多次随机小偏移计算后取平均。

4. 三维体数据：体系计数与内存优化

4.1 三维数组的读取和二值化

三维体数据常见于医学 CT、工业扫描、材料微观结构，通常是一个三维数组 shape 为 (D, H, W)，每个体系（voxel）是灰度值。和二维一样，盒维数描述的是“哪些体素属于目标”，所以必须先二值化。

举个例子，模拟项目 X 里有一段血管扫描数据，目标亮度明显高于背景，直接设定一个固定阈值vol > 200就能得到 mask。但有些数据对比度不高，我会先做中值滤波去掉噪点，再用 Otsu 阈值。注意：三维 Otsu 的实现较复杂，最简单的做法是先把三维数组展平，对全部灰度值做 Otsu，再还原成三维 mask。这种全局阈值对大多数体数据够用了。

三维盒维数的计算逻辑与二维完全一致，只是从二维网格变成了三维网格：把每个非零体素的坐标归一化到单位立方体，然后用不同边长的立方体盒子去覆盖，统计非空盒子数。

4.2 三维网格计数实现：用 unique 的行进行去重

直接照搬二维的set(zip(...))在三维时容易卡成蜗牛，因为体素数量可能达到百万量级。更合理的做法是提前把非零体素的坐标变成一个二维数组，每一行是一个体素的 (x, y, z)，然后用np.unique对行去重。下面是参考实现：

import numpy as np def box_dimension_3d(vol, box_sizes=None, fit_ratio=0.6): # vol: shape (D, H, W) 的三维数组，非零为目标 zz, yy, xx = np.nonzero(vol) if len(xx) == 0: return np.nan max_len = max(vol.shape) coords = np.stack([ xx / max_len, yy / max_len, zz / max_len ], axis=1) # shape (n_points, 3) if box_sizes is None: # 最小盒子边长取 2 个体素的长度 min_size = 2.0 / max_len box_sizes = np.exp(np.linspace(np.log(min_size), 0, 20))[::-1] counts = [] for s in box_sizes: # 计算每个体素落在哪个盒子 grid = np.floor(coords / s).astype(np.int64) # 对盒子索引行去重，统计非空立方体数量 unique_grid = np.unique(grid, axis=0) counts.append(len(unique_grid)) log_s = np.log(1.0 / np.asarray(box_sizes)) log_N = np.log(np.asarray(counts)) n = len(log_s) start = int(n * (1 - fit_ratio) / 2) end = int(n * (1 + fit_ratio) / 2) slope, _ = np.polyfit(log_s[start:end], log_N[start:end], 1) return slope

这段代码里，np.unique(grid, axis=0)会按行去重，返回所有出现过的盒子索引。相比set(zip(...))，它在 numpy 底层用排序实现，内存占用可控，速度也快很多。

注意coords / s得到的是浮点数，floor后转 int 作为盒子索引。如果s很小，盒子数量会非常多，unique排序的时间会上升。一个折中是把盒子边长序列的步长设得稍微大一点，比如每个尺度是前一个的 0.7 倍，而不是 0.5 倍，这样总尺度数减少，拟合依然有效。

4.3 内存优化：用一维索引代替三维行去重

当非零体素数量超过几百万时，即使np.unique排序也可能产生较大临时内存。一个更省内存的方法是先把三维坐标映射到一维索引，再进行去重。具体做法是：把盒子网格索引的三个数字分别乘以合理的步长后相加。

例如，设定盒子最大索引范围n_x, n_y, n_z，则：

# 假设 gx, gy, gz 都是 int64 数组 linear_idx = gx * n_y * n_z + gy * n_z + gz unique_count = len(np.unique(linear_idx))

这种做法的好处是不需要保存一个 m×3 的数组，只需要保存一个 m 长度的一维数组，去重时内存开销小很多。注意n_y * n_z要足够大，避免不同 (gx, gy, gz) 映射到同一个线性索引，也就是乘法步长要大于坐标范围的最大值。

另外，如果体数据非常稠密，比如目标占比超过一半，用np.nonzero本身就会返回大量坐标数组。这时可以考虑分块计算：把体数据切成几个子块，每个子块分别计算盒子覆盖数，最终合并时要注意去重，因为边界上的盒子可能跨子块。简单做法是把子块重叠一个盒子边长，再做去重。这部分属于优化细节，大多数场景下做好坐标归一化，已经够用了。

5. 避坑与常见问题：为什么我算出来的维数总是偏大或偏小

5.1 现象：log-log 曲线两端弯曲，拟合斜率不稳定

画出来的散点图不是一条直线，左右两端都往下弯，只有中间一段看起来才像直线。

原因：小尺度过小，每个目标点占据独立盒子，N(s) 不再增长；大尺度过大，所有点集中到很少几个盒子，同样破坏幂律关系。二者都会让拟合斜率失真。

解决：不要用全部尺度做线性拟合。我习惯先把尺度序列取 30 个左右，通过fit_ratio=0.6只保留中间 60% 的点参与拟合；如果散点中间段也不直，就把尺度上下限缩窄，重新计算。一个自查方法是打印每个尺度的N(s)，查看有没有连续几个尺度数值不变，如果有，说明已经越过分辨率极限。

5.2 现象：同一幅图像旋转 30 度后，维数变化超过 0.1

同一张二值图，旋转后再算盒维数，两次结果差异明显。

原因：盒子网格是轴对齐的，目标图案与网格的相位和方向关系会影响非空盒子数。尤其对于直线边缘、规则纹理，这种影响更大。

解决：对图像做多种随机偏移或多次小角度旋转，计算多个维数后取平均。旋转后需要重新二值化和插值，建议使用抗锯齿插值并固定阈值。如果对计算效率要求高，也可以只对坐标施加一个随机平移扰动，重复 5～10 次后取均值。

5.3 现象：灰度图换一个阈值，维数从 1.3 变成 1.7

同一张灰度图，手动调阈值后维数大幅波动。

原因：盒维数统计的是“非空盒子”，如果目标分割不明确，阈值变化直接改变了前景面积分布，维数自然跟着变。灰度图本身并不是二值集合，不能直接拿灰度值做盒计数。

解决：明确你要测的是形状还是灰度表面。测形状就必须先分割成稳定的二值图，阈值选择应基于先验背景或 Otsu；测表面纹理可以改用灰度差分盒维数，把灰度值当作高度场，再对高度场做网格覆盖。切忌一边测形状一边用不同阈值反复试。

5.4 现象：三维 512×512×512 数据计算时内存爆掉

运行三维代码时，程序直接内存溢出或卡到无响应。

原因：np.nonzero返回三个坐标数组，每个元素都是 int64，如果目标体素占一半以上，光坐标数组就有几百 MB。后续np.unique还会产生排序副本，内存压力更大。

解决：先做稀疏化，只保留目标体素坐标。如果目标占比实在太高，就要分块计算。另一个技巧是先用vol = vol.astype(np.bool_)压缩布尔存储，再用上面提到的线性索引去重法减少临时数组维度。同时把盒子边长下界从 1 个体系提高到 2～3 个体系，能显著减少计数点数。

5.5 现象：两个开源脚本算同一个数据集，结果差 0.15

同一个二值图，A脚本算出来 1.4，B脚本算出来 1.55，看起来都差不多，但结果就是不一致。

原因：不同实现对坐标归一化的边界定义不一样。有的把坐标范围放在 [0,1]，有的放在 [0,1)，有的用floor有的用ceil，还有的网格原点从左上角像素中心开始。这些微小差异在小尺度下会被放大。

解决：统一约定：坐标除以包围盒长度后，范围控制在 [0,1) 左闭右开；使用floor作为取整方式；盒子边长从最小两个像素开始。我写代码时会在归一化前对坐标做coord = coord - coord.min()，再除以(coord.max() - coord.min()) * 1.0000001，确保最大值落在下一个盒子边界时不会因为浮点精度误入新盒子。

6. 让盒维数算得更稳：自适应尺度选择与已知分形验证

当你已经能稳定计算一维、二维、三维盒维数，接下来真正提升可信度的是两件事：尺度自适应和结果校验。

尺度自适应，是指不要固定 20 个等比尺度，而是根据数据点密度动态决定最小盒子边长。我一般用包围盒对角线的 1/500 作为最小尺度下限，再根据拟合残差把明显偏离直线的点剔除。具体做法是先做一次初始拟合，计算每个散点的残差，去掉残差超过 2 倍标准差的点，再重新拟合一遍。这样得到的斜率更稳定，尤其是三维数据，不会因为少数极端尺度破坏整体趋势。

下表是我平时用的三种尺度策略对比：

验证方法比优化算法更值得投入。我最常做的校验不是对照论文数据，而是直接用已知维数的图形测试代码：一条直线段的盒维数应当接近 1.0，一个实心正方形接近 2.0，一个实心立方体接近 3.0。更严格一点的测试是用 Sierpinski 地毯标准图案，理论维数约 1.8928。如果亲手生成的 Sierpinski 图案算出来在 1.85～1.93 区间，基本可以放心。

代码里可以这样加一个快速自检：

def validate_dimension_calc(func, data, expected, tol=0.05): result = func(data) print(f"computed={result:.4f}, expected={expected:.4f}") assert abs(result - expected) < tol, "维数校验未通过"

从那以后，我每次拿到新的体数据或图像集，都会先跑一遍直线、正方形、立方体三个基准测试，确认代码没有在归一化、边界处理或拟合尺度上出问题，再上真实数据。这一步看似浪费时间，实际上省掉了我后面大量返工。希望这篇笔记里的实现和避坑经验能让你少走一圈弯路，计算分形盒维数时心里更有底。

本文还有配套的精品资源，点击获取