简介这套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_sizesNone, fit_ratio0.6): # signal: 一维数组长度 n x np.arange(len(signal)) y np.asarray(signal, dtypefloat) # 归一化到单位正方形避免横纵量纲不一致 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_ratio0.6表示只保留 log-log 散点中间 60% 的部分做拟合两端通常因为尺度上下限而偏离直线。2.3 参数选择尺度范围决定结果可信度盒维数计算的“玄学”大部分集中在尺度范围的选择。同样一条曲线如果把最小盒子边长设成 1/1000拟合结果可能接近 1.8如果设成 1/50结果可能只有 1.2。原因在于尺度下限低于曲线的采样分辨率时曲线看起来像一堆孤立的点维数会虚高尺度上限太大时整个曲线被装进很少的几个盒子增长趋势饱和维数又会偏低。我一般的做法是先画出 log-log 散点图看中段有没有明显的线性区间。最小盒子边长的下限取平均点间距的 23 倍上限取目标包围盒长边的 1/5 左右因为再大了基本就是一两个盒子覆盖全图。如果信号长度很长比如几万个点可以多取一些尺度比如 3040 个让拟合更稳。还可以用一个粗糙的验证方法对一段纯正弦信号盒维数应当接近 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_sizesNone, fit_ratio0.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, axis0)或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_sizesNone, fit_ratio0.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 ], axis1) # 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, axis0) 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, axis0)会按行去重返回所有出现过的盒子索引。相比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_ratio0.6只保留中间 60% 的点参与拟合如果散点中间段也不直就把尺度上下限缩窄重新计算。一个自查方法是打印每个尺度的N(s)查看有没有连续几个尺度数值不变如果有说明已经越过分辨率极限。5.2 现象同一幅图像旋转 30 度后维数变化超过 0.1同一张二值图旋转后再算盒维数两次结果差异明显。原因盒子网格是轴对齐的目标图案与网格的相位和方向关系会影响非空盒子数。尤其对于直线边缘、规则纹理这种影响更大。解决对图像做多种随机偏移或多次小角度旋转计算多个维数后取平均。旋转后需要重新二值化和插值建议使用抗锯齿插值并固定阈值。如果对计算效率要求高也可以只对坐标施加一个随机平移扰动重复 510 次后取均值。5.3 现象灰度图换一个阈值维数从 1.3 变成 1.7同一张灰度图手动调阈值后维数大幅波动。原因盒维数统计的是“非空盒子”如果目标分割不明确阈值变化直接改变了前景面积分布维数自然跟着变。灰度图本身并不是二值集合不能直接拿灰度值做盒计数。解决明确你要测的是形状还是灰度表面。测形状就必须先分割成稳定的二值图阈值选择应基于先验背景或 Otsu测表面纹理可以改用灰度差分盒维数把灰度值当作高度场再对高度场做网格覆盖。切忌一边测形状一边用不同阈值反复试。5.4 现象三维 512×512×512 数据计算时内存爆掉运行三维代码时程序直接内存溢出或卡到无响应。原因np.nonzero返回三个坐标数组每个元素都是 int64如果目标体素占一半以上光坐标数组就有几百 MB。后续np.unique还会产生排序副本内存压力更大。解决先做稀疏化只保留目标体素坐标。如果目标占比实在太高就要分块计算。另一个技巧是先用vol vol.astype(np.bool_)压缩布尔存储再用上面提到的线性索引去重法减少临时数组维度。同时把盒子边长下界从 1 个体系提高到 23 个体系能显著减少计数点数。5.5 现象两个开源脚本算同一个数据集结果差 0.15同一个二值图A脚本算出来 1.4B脚本算出来 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 倍标准差的点再重新拟合一遍。这样得到的斜率更稳定尤其是三维数据不会因为少数极端尺度破坏整体趋势。下表是我平时用的三种尺度策略对比策略适用场景优点注意点固定等比尺度快速批量计算代码简单结果可复现需要人工判断尺度范围自适应尺度数据点分布不均匀自动避开分辨率极限要额外计算点间距统计多次随机偏移平均目标方向性较强降低网格相位干扰计算量大需重复 510 次验证方法比优化算法更值得投入。我最常做的校验不是对照论文数据而是直接用已知维数的图形测试代码一条直线段的盒维数应当接近 1.0一个实心正方形接近 2.0一个实心立方体接近 3.0。更严格一点的测试是用 Sierpinski 地毯标准图案理论维数约 1.8928。如果亲手生成的 Sierpinski 图案算出来在 1.851.93 区间基本可以放心。代码里可以这样加一个快速自检def validate_dimension_calc(func, data, expected, tol0.05): result func(data) print(fcomputed{result:.4f}, expected{expected:.4f}) assert abs(result - expected) tol, 维数校验未通过从那以后我每次拿到新的体数据或图像集都会先跑一遍直线、正方形、立方体三个基准测试确认代码没有在归一化、边界处理或拟合尺度上出问题再上真实数据。这一步看似浪费时间实际上省掉了我后面大量返工。希望这篇笔记里的实现和避坑经验能让你少走一圈弯路计算分形盒维数时心里更有底。本文还有配套的精品资源点击获取
