简介针对分形维数计算需求这份压缩包提供了一套基于Box Counting方法的MATLAB实现适合分形几何、图像分析与复杂系统研究方向的初学者及科研人员参考。压缩包仅2KB共包含3个m文件分别覆盖核心计数算法、Sierpinski三角形测试样本生成与分形维数封装调用整体结构紧凑便于单独抽取或整体运行。算法流程完整呈现了从初始网格覆盖、逐级细分、非空盒子统计到对数拟合估算维数的关键步骤读者结合Sierpinski三角形这一经典分形可直观校验不同网格粒度下的计数结果代码注释清晰、体量轻也适合逐行研读和参数调优。同时数据输入接口方便替换为其他图像或数值集可扩展到海岸线、云层结构等自然分形场景。目前已有191人浏览学习这款轻量脚本在教学演示与快速实验中具有不错的参考价值。1. 从 box_count.zip 说起盒计数维要量化的问题一个 zip 压缩包叫 box_count.zip解压后装好 numpy 就能算数这说的就是 box-counting 维数估计用边长越来越小的盒子覆盖图像前景数盒子数量拟合双对数斜率输出一个 0 到 2 之间的数 D。这个数对规则对象是整数直线是 1、平面是 2对分形是小数Koch 雪花约 1.26。它解决的问题是“这张图的结构复杂到什么程度”材料断面、道路网、血管、多孔介质都在用。下面从定义、实现、参数到校准给出一套能自己复现的最小方案读完可以直接拿任意二值图跑出维数并确认这个数是结构特征而不是拟合噪声。2. 盒计数维的定义N(ε)、ln N(ε) 与标度律2.1 从覆盖到极限Minkowski 维数怎么定义对平面点集 S用边长为 ε 的正方形盒子覆盖它记最少盒子数为 N(ε)。盒计数维也叫 Minkowski-Bouligand 维数定义为D lim(ε→0) ln N(ε) / ln(1/ε)实际代码里有两个近似一是“最少盒子数”很难严格证明通常用固定网格覆盖的数量代替二者在同一标度下只差常数倍取对数后常数项被斜率吸收不影响 D二是极限取不到图像分辨率就是 ε 的下限。所以工程问题变成在可观测的 ε 区间内N(ε) 是否呈现幂律幂指数是多少。用三个例子建立直觉。一条长度为 L 的水平线段N(ε) ≈ L/εD 1。一个边长为 L 的实心正方形N(ε) ≈ (L/ε)²D 2。Koch 雪花每迭代一次边长变为原来的 4/3覆盖它所需的盒子数按非整数幂增长D ln4/ln3 ≈ 1.2619。盒计数维的价值就在这里用连续的数刻画“介于线和面之间”的结构而不是靠肉眼说“它有点像分形”。2.2 幂律、对数坐标与斜率符号对自相似对象放大后的结构与整体同构覆盖数满足 N(ε) ∝ ε^(-D)。两边取对数得到ln N(ε) -D ln ε ln C在 ln ε 做自变量的坐标系里这是一条直线斜率是 -D。因此有两种等价写法对 (ln ε, ln N) 回归后取斜率相反数或对 (ln(1/ε), ln N) 回归直接读正斜率。我固定用后一种因为和定义式逐项对应下游写报告、对拍别人脚本时少一层换算。这里有个经常对不上的坑很多开源脚本输出负斜率另一些输出绝对值两者都叫“维数”。拿 A 脚本的结果去跟 B 脚本的报告比先确认坐标约定否则 1.26 和 -1.26 会让人白折腾半天。用实心方块可以快速校准这条约定下面是 512×512 方块在八个尺寸下的理想覆盖数import numpy as np sizes np.array([2, 4, 8, 16, 32, 64, 128, 256], dtypefloat) ns np.array([65536, 16384, 4096, 1024, 256, 64, 16, 4]) # (512/s)^2 x, y np.log(1.0 / sizes), np.log(ns) D, _ np.polyfit(x, y, 1) print(D) # 2.0000这就是实心方块的标定基准逻辑说明ns 是理想覆盖数也就是边长为 s 的盒子铺满 512×512 方块时需要的数量对数拟合出来的斜率必须严格等于 2。你的实现跑出 2.00 说明坐标约定和回归逻辑都对跑出别的数问题一定出在计数或拟合上。2.3 回归区间四个常见的拟合陷阱回归不是把所有点扔进去就完事四个最常见的偏差来源第一ε 1 的点必须去掉。N(1) 等于前景像素总数只反映面积不反映结构在对数坐标里它离其他点最远杠杆作用最大会明显拉低斜率。第二ε 接近图像边长时 N(ε) 跌到 1 附近进入饱和区这些点同样失真。第三点数太少三四个点做最小二乘单点噪声直接决定结果。第四图像不是正方形时超过短边一半的盒子覆盖的是窄条区域退化成近似一维计数。实际取点我一般这样定ε 从 2 开始按 2^k 递增取到 min(H,W)/2 为止。图像越小的图能用的档位越少常见参考如下图像短边推荐 size 序列有效点数1282, 4, 8, 16, 32, 6465122, 4, 8, 16, 32, 64, 128, 256810242 到 512共 9 档9注意max-size 并非越大越“精细”。超过 min(H,W)/2 后回归进入饱和区斜率向 0 偏移宁小勿大。3. 用 Python 实现 box-countingreshape 版核心函数与命令行3.1 zip 解压后的标准结构与最小运行命令box_count.zip 解压后最常见的布局是三件套box_count.py 一个自包含脚本requirements.txt 锁定 numpy 和 PillowREADME.md 写参数约定。先建环境跑一次最小命令python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy pillow python box_count.py sample.png --threshold 170 --csv result.csv命令里的 --threshold 是灰度二值化阈值灰度小于 170 的像素视为前景--csv 把每一档的 (size, count) 原始数据落盘。第一次跑一定要带 --csv因为只看最后输出的 D无法判断回归有没有踩到 2.3 说的饱和区。跑通之后再把 sample.png 换成自己的图。3.2 核心函数 box_count一次 reshape 统计全部盒子计数过程可以压缩成“裁剪 reshape any”三步不需要嵌套循环。完整可运行版本如下import numpy as np from PIL import Image def load_binary(path: str, threshold: int 128) - np.ndarray: 读图并二值化返回 bool 数组True 表示前景。 img Image.open(path).convert(L) arr np.array(img, dtypenp.uint8) return arr threshold def box_count(binary: np.ndarray, sizes) - list: 对每个边长 s返回 (s, 非空盒子数) 列表。 h, w binary.shape out [] for s in sizes: if s 2 or s h or s w: continue hc, wc h - h % s, w - w % s # 裁剪到 s 的整数倍 crop binary[:hc, :wc] # 高分成 hc//s 块每块内部 s 行宽同理 cells crop.reshape(hc // s, s, wc // s, s) occupied cells.any(axis(1, 3)) # 消掉盒内两个方向 out.append((s, int(occupied.sum()))) return out def fit_dimension(counts: list) - tuple: 对 (s, N) 做 log-log 回归返回 (D, intercept)。 sizes np.array([float(s) for s, _ in counts]) ns np.array([float(c) for _, c in counts]) valid (ns 1) (sizes 1) # 丢弃饱和点与 ε1 x np.log(1.0 / sizes[valid]) # 自变量 ln(1/ε) y np.log(ns[valid]) # 因变量 ln N(ε) slope, intercept np.polyfit(x, y, 1) return slope, intercept # slope 就是 D逻辑说明crop.reshape 把 (hc, wc) 数组看成四维张量四个轴依次是“行块号、盒内行、列块号、盒内列”用 any(axis(1, 3)) 一次性消掉两个盒内维度剩下的 (行块数, 列块数) 矩阵里 True 的数量就是 N(ε)。reshape 按行主序展开行块和列块的切分顺序正是先高后宽不需要额外转置。如果改成链式写法 any(axis1).any(axis2) 也等价但注意 3D 版本不能随便链式见第 6 章。参数说明sizes 通常给 2 的幂序列threshold 决定前景粗细曲线边缘的抗锯齿像素如果低于阈值会被算进前景D 系统性偏大valid 过滤只处理极端点过滤后只剩两三个点说明 max-size 设得太大。提示hc/wc 裁剪和 reshape 都只产生视图不复制底层数组内存开销接近原图本身4096×4096 的 bool 图约 16 MB可以放心跑。3.3 效率对比向量化 any 与三层 for 循环常见误用是写三层 for 循环对每个盒子遍历 s² 个像素判占用复杂度约等于 H·W×档位数但访存完全随机、缓存命中差reshape 版的 any 在连续内存块上做归约每个元素只读一次。同一台机器上的量级大致如下图像尺寸循环版耗时reshape 版耗时512×512约 2 秒约 15 毫秒1024×1024约 8 秒约 40 毫秒4096×4096约 2 分钟约 0.6 秒对单张图这点差距无所谓但批量处理几百张图时向量化直接决定任务能不能当晚跑完。这也是 box-counting 适合打成 zip 分发的原因核心逻辑就这么几十行依赖面窄换机器不会出现“性能玄学”。3.4 把 (size, count) 画出来调参前的第一件事D 只是一个回归斜率丢掉了很多信息。用 csv 落盘的数据画一张散点图很多问题一眼可见import csv import numpy as np import matplotlib.pyplot as plt rows list(csv.reader(open(result.csv)))[1:] # 跳过表头 sizes np.array([float(r[0]) for r in rows]) counts np.array([float(r[1]) for r in rows]) x np.log(1.0 / sizes) y np.log(counts) plt.plot(x, y, o-) plt.xlabel(ln(1/ε)) plt.ylabel(ln N(ε)) plt.show()实心方块的图应该是近似等间距的 8 个点最左端略微下弯是正常饱和最右端三点如果歪掉优先查阈值而不是拟合代码。matplotlib 只是检查工具建议不放进 box_count.py 的导入链服务器上没有图形环境时脚本不能因此报错。4. box size、边缘处理与阈值四个影响盒计数结果的参数4.1 box size 序列2 的幂还是线性递增理论上只要求 ε 趋于 0没规定怎么采样。实际影响在双对数坐标里2 的幂序列使 ln(1/ε) 等距分布回归时每个点的杠杆差不多线性递增序列会把大部分点挤在大盒子端小盒子端只剩一两个点拟合斜率偏向中段尺度对细小结构不敏感所以默认用 2^k。实操上我还会加一个缩比检查把同一张图缩到一半分辨率再算一次两次 D 的差超过 0.05说明最小特征尺寸和最小盒子已经可比需要换更高分辨率素材而不是调回归区间。显微镜图像尤其要养成这个习惯因为采样分辨率决定了可观测的最小尺度也决定了维数这个数算到哪一档截止。4.2 边缘处理裁剪、填充与输入规整图像边长不能被 s 整除时裁剪丢弃右/下边缘填充补背景。两种做法的偏差方向相反处理方式对 D 的偏差方向适用场景裁剪偏低丢边缘结构前景远离边缘计算要快补 0 填充偏高边界产生空盒子前景贴近边缘不能丢像素规整到 2 的幂几乎无偏采集阶段的最优解更稳的做法是让输入边长本身就是 2 的幂采集 ROI 时约定尺寸或处理前用Image.resize((512, 512))规整。注意缩放会引入抗锯齿必须先二值化再缩放缩放后若还想调阈值重新做一次 threshold。不要默认填充因为填充的偏差随盒子尺度漂移很难在回归阶段修正。4.3 二值化阈值与骨架化的边界阈值直接决定前景“粗细”。一条抗锯齿直线阈值取 100 会把半透明过渡像素算进前景等效线宽变大小盒子端 N(ε) 上升D 从 1.0 漂向 1.2。补救方法分两类灰度直方图双峰时用 Otsu 全局阈值skimage.filters.threshold_otsu 一行搞定对象本来就是线条结构道路网、血管、骨架时先 skeletonize 再计数否则线宽占多个像素会让 D 系统性偏高。骨架化也有代价细化算法会制造 2×2 方块和交叉点伪影这些局部结构在 ε 缩小时被放大。因此骨架化之后建议把最小 box size 从 2 提到 4让伪影落在回归区间外。判断标准始终是 5.3 的 log-log 图小盒子端上翘多半就是这类伪影。4.4 多网格偏移平均消除网格对齐偏差固定网格计数有个已知毛病45° 对角线这类结构在 s 较大时可能正好骑在网格线上非空盒子数周期性骤减log-log 图出现锯齿回归斜率随之摆动。补救方法是多网格法对同一个 s把网格原点平移几个偏移量分别计数取平均def box_count_offset(binary, s, offsetsNone): if offsets is None: offsets (0, s // 2) h, w binary.shape total 0 for dx in offsets: hc h - (h - dx) % s # 从 dx 起能整分的最大高度 wc w - (w - dx) % s # 宽度方向用同一偏移简单起见 crop binary[dx:dx hc, dx:dx wc] cells crop.reshape(hc // s, s, wc // s, s) total int(cells.any(axis(1, 3)).sum()) return s, total / len(offsets)逻辑说明偏移 dx 表示网格原点落在第 dx 行/列hc 取满足 hc ≡ dx (mod s) 的最大高度切片从 dx 开始取 hc 个元素末尾不残留。两次计数的平均值作为 N(ε) 参与回归对角线结构的锯齿被抹平。参数说明offsets 扩到 4 个四个象限后效果基本饱和每多一个偏移只是多一次 reshapeany总耗时线性增长对 4096×4096 也不敏感。如果回归点在某档位出现明显锯齿先单独验证这一档的偏移效果再决定是否全局启用。5. 用已知分形校准 box-countingCantor 集、Koch 曲线与 Sierpinski 地毯5.1 三个理论维数已知的生成函数验证实现的标准做法精确生成理论维数已知的结构跑同一套 box_count看结果偏差。三个生成函数只依赖 numpydef sierpinski_carpet(n5): n 阶 Sierpinski 地毯边长 3**n理论 D ln8/ln3 ≈ 1.8928。 carpet np.ones((1, 1), dtypebool) kernel np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtypebool) for _ in range(n): carpet np.kron(kernel, carpet) # 旧像素扩展成 3x3 块中心挖空 return carpet def cantor_1d(n6): n 阶 Cantor 集长度 3**n理论 D ln2/ln3 ≈ 0.6309。 line np.ones(3 ** n, dtypebool) step 3 ** n for _ in range(n): step // 3 for start in range(0, 3 ** n, 3 * step): line[start step:start 2 * step] False return line def koch_points(depth, a, b): 递归生成 Koch 曲线点列depth0 时就是端点。 a, b np.asarray(a, dtypefloat), np.asarray(b, dtypefloat) if depth 0: return [a, b] v (b - a) / 3.0 p1, p2 a v, a 2 * v apex p1 np.array([v[0]/2 - v[1]*np.sqrt(3)/2, v[0]*np.sqrt(3)/2 v[1]/2]) pts [] for pa, pb in [(a, p1), (p1, apex), (apex, p2), (p2, b)]: pts koch_points(depth - 1, pa, pb)[:-1] return pts [b]逻辑说明np.kron(kernel, carpet) 把 carpet 的每个元素替换成 kernel 标定过的 3×3 块kernel 为 1 的位置放整块 carpet为 0 的位置全 False迭代 n 次得到 3^n 边长中心孔洞逐级保留。Cantor 集在每段迭代中去掉中间 1/3start 步长按当前段长 3*step 递进。Koch 曲线用递归在每段上构造等边三角形凸包depth4 时点数是 4^41 257 个。Koch 点列画到图像上再二值化PIL 自带画线不需要额外依赖from PIL import Image, ImageDraw def koch_image(depth4, size729): pts koch_points(depth, (0.0, 0.0), (1.0, 0.0)) img Image.new(L, (size, size), 255) dr ImageDraw.Draw(img) prev None for x, y in pts: cur (int(x * (size - 4)) 2, int(y * (size - 4)) 2) if prev is not None: dr.line([prev, cur], fill0, width1) prev cur return np.array(img) 128注意画布的 y 方向只用到约 0.29 的比例上半部分留白是正常的空盒子不计入 N(ε)不影响计数但盒子尺寸超过曲线实际跨度后 count 会迅速降到 1所以 Koch 样本的 box size 上限取 128 而不是 256。5.2 校准结果与判定标准用 3.2 的函数跑上述样本典型的实测范围如下具体值随 numpy 版本和画线抗锯齿略有浮动样本理论 D回归区间典型实测范围单像素水平直线1.00002..5121.001.0145° 对角线1.00002..2560.981.00Cantor 集6 阶0.63092..1280.620.64Sierpinski 地毯5 阶1.89282..1281.881.90Koch 曲线4 阶1.26192..1281.251.27判定标准实心方块必须严格收敛到 2.00这是实现 bug 的试金石直线必须贴近 1.00两个标准分形偏差在 ±0.02 内算合格。超过这个范围先回第 4 章查边缘处理和阈值不要急着改回归区间。对角线这类结构对网格对齐敏感没做 4.4 的偏移平均时 D 会周期性波动这是方法特性不是 bug。5.3 快速目检log-log 图上的两种系统偏差给一个配套绘图函数把拟合直线和原始点叠在一起看def plot_fit(counts): import matplotlib.pyplot as plt sizes np.array([float(s) for s, _ in counts]) ns np.array([float(c) for _, c in counts]) ok ns 1 x, y np.log(1.0 / sizes[ok]), np.log(ns[ok]) slope, ic np.polyfit(x, y, 1) xs np.linspace(x.min(), x.max(), 50) plt.plot(x, y, o, labelfD{slope:.3f}) plt.plot(xs, slope * xs ic, --, colorgray) plt.xlabel(ln(1/ε)) plt.ylabel(ln N(ε)) plt.legend() plt.show()看这张图只需要两个判断。点列在直线下方系统性下弯是大盒子端进入饱和砍掉最大的几个尺寸点列在小盒子端上翘是抗锯齿或骨架伪影在抬升 N(ε)把最小尺寸从 2 提到 4或回去调阈值。两种情况都别靠删点硬凑根因在输入不在拟合。6. 升级到 3D 体数据六维 reshape 与开机自检6.1 从四维 reshape 到六维三维盒计数盒计数不限于二维图。CT 扫描、多孔介质、泡沫材料的二值体数据用同样的覆盖思想只是把盒子换成正方体D 落在 0 到 3 之间。实现只需把第 3 章的四维 reshape 扩成六维并且注意 any 的轴编号要一次给全def box_count_3d(binary3d: np.ndarray, sizes): 输入 (H, W, D) 的 bool 数组输出 (s, count) 列表。 h, w, d binary3d.shape out [] for s in sizes: hc, wc, dc h - h % s, w - w % s, d - d % s crop binary3d[:hc, :wc, :dc] cells crop.reshape(hc // s, s, wc // s, s, dc // s, s) occupied cells.any(axis(1, 3, 5)) # 一次消掉三个盒内方向 out.append((s, int(occupied.sum()))) return out六轴顺序是“行块号、盒内行、列块号、盒内列、深块号、盒内深”。这里不能写成交替链式 any(axis1).any(axis3).any(axis5)因为每次归约后轴编号会左移直接链式会把块维度误消掉。两个 3D 特有约束体素必须各向同性医学影像 z 轴层厚和面内分辨率不同时先重采样成立方体素边缘剩余超过 10% 时先规整到 2 的幂不要填充理由同 4.2。6.2 给 zip 分发的脚本加一个自检模式最后说交付习惯。box_count.zip 这类脚本到了别人机器上最常见的问题不是装不起来而是“装起来但不敢信结果”。我一般会在main里加 --self-test用 5.1 的生成函数做断言入口命令python box_count.py --self-test对应分支如下# argparse 里加 # ap.add_argument(--self-test, actionstore_true) if args.self_test: carpet sierpinski_carpet(5) counts box_count(carpet, [2 ** k for k in range(1, 8)]) d, _ fit_dimension(counts) assert abs(d - 1.8928) 0.03, fself-test failed: D{d:.4f} print(fself-test OK, D{d:.4f})自检用生成函数而不是素材图片因为生成只依赖 numpy解压后零素材直接验证环境。断言阈值放 ±0.03比 5.2 的判定标准略宽给不同 numpy 版本留余量跑过自检再处理用户自己的图D 的可信度就从“脚本给的数”变成“自己验证过的数”。本文还有配套的精品资源点击获取
Python实现盒计数维：从图像到分形维数的完整指南
简介针对分形维数计算需求这份压缩包提供了一套基于Box Counting方法的MATLAB实现适合分形几何、图像分析与复杂系统研究方向的初学者及科研人员参考。压缩包仅2KB共包含3个m文件分别覆盖核心计数算法、Sierpinski三角形测试样本生成与分形维数封装调用整体结构紧凑便于单独抽取或整体运行。算法流程完整呈现了从初始网格覆盖、逐级细分、非空盒子统计到对数拟合估算维数的关键步骤读者结合Sierpinski三角形这一经典分形可直观校验不同网格粒度下的计数结果代码注释清晰、体量轻也适合逐行研读和参数调优。同时数据输入接口方便替换为其他图像或数值集可扩展到海岸线、云层结构等自然分形场景。目前已有191人浏览学习这款轻量脚本在教学演示与快速实验中具有不错的参考价值。1. 从 box_count.zip 说起盒计数维要量化的问题一个 zip 压缩包叫 box_count.zip解压后装好 numpy 就能算数这说的就是 box-counting 维数估计用边长越来越小的盒子覆盖图像前景数盒子数量拟合双对数斜率输出一个 0 到 2 之间的数 D。这个数对规则对象是整数直线是 1、平面是 2对分形是小数Koch 雪花约 1.26。它解决的问题是“这张图的结构复杂到什么程度”材料断面、道路网、血管、多孔介质都在用。下面从定义、实现、参数到校准给出一套能自己复现的最小方案读完可以直接拿任意二值图跑出维数并确认这个数是结构特征而不是拟合噪声。2. 盒计数维的定义N(ε)、ln N(ε) 与标度律2.1 从覆盖到极限Minkowski 维数怎么定义对平面点集 S用边长为 ε 的正方形盒子覆盖它记最少盒子数为 N(ε)。盒计数维也叫 Minkowski-Bouligand 维数定义为D lim(ε→0) ln N(ε) / ln(1/ε)实际代码里有两个近似一是“最少盒子数”很难严格证明通常用固定网格覆盖的数量代替二者在同一标度下只差常数倍取对数后常数项被斜率吸收不影响 D二是极限取不到图像分辨率就是 ε 的下限。所以工程问题变成在可观测的 ε 区间内N(ε) 是否呈现幂律幂指数是多少。用三个例子建立直觉。一条长度为 L 的水平线段N(ε) ≈ L/εD 1。一个边长为 L 的实心正方形N(ε) ≈ (L/ε)²D 2。Koch 雪花每迭代一次边长变为原来的 4/3覆盖它所需的盒子数按非整数幂增长D ln4/ln3 ≈ 1.2619。盒计数维的价值就在这里用连续的数刻画“介于线和面之间”的结构而不是靠肉眼说“它有点像分形”。2.2 幂律、对数坐标与斜率符号对自相似对象放大后的结构与整体同构覆盖数满足 N(ε) ∝ ε^(-D)。两边取对数得到ln N(ε) -D ln ε ln C在 ln ε 做自变量的坐标系里这是一条直线斜率是 -D。因此有两种等价写法对 (ln ε, ln N) 回归后取斜率相反数或对 (ln(1/ε), ln N) 回归直接读正斜率。我固定用后一种因为和定义式逐项对应下游写报告、对拍别人脚本时少一层换算。这里有个经常对不上的坑很多开源脚本输出负斜率另一些输出绝对值两者都叫“维数”。拿 A 脚本的结果去跟 B 脚本的报告比先确认坐标约定否则 1.26 和 -1.26 会让人白折腾半天。用实心方块可以快速校准这条约定下面是 512×512 方块在八个尺寸下的理想覆盖数import numpy as np sizes np.array([2, 4, 8, 16, 32, 64, 128, 256], dtypefloat) ns np.array([65536, 16384, 4096, 1024, 256, 64, 16, 4]) # (512/s)^2 x, y np.log(1.0 / sizes), np.log(ns) D, _ np.polyfit(x, y, 1) print(D) # 2.0000这就是实心方块的标定基准逻辑说明ns 是理想覆盖数也就是边长为 s 的盒子铺满 512×512 方块时需要的数量对数拟合出来的斜率必须严格等于 2。你的实现跑出 2.00 说明坐标约定和回归逻辑都对跑出别的数问题一定出在计数或拟合上。2.3 回归区间四个常见的拟合陷阱回归不是把所有点扔进去就完事四个最常见的偏差来源第一ε 1 的点必须去掉。N(1) 等于前景像素总数只反映面积不反映结构在对数坐标里它离其他点最远杠杆作用最大会明显拉低斜率。第二ε 接近图像边长时 N(ε) 跌到 1 附近进入饱和区这些点同样失真。第三点数太少三四个点做最小二乘单点噪声直接决定结果。第四图像不是正方形时超过短边一半的盒子覆盖的是窄条区域退化成近似一维计数。实际取点我一般这样定ε 从 2 开始按 2^k 递增取到 min(H,W)/2 为止。图像越小的图能用的档位越少常见参考如下图像短边推荐 size 序列有效点数1282, 4, 8, 16, 32, 6465122, 4, 8, 16, 32, 64, 128, 256810242 到 512共 9 档9注意max-size 并非越大越“精细”。超过 min(H,W)/2 后回归进入饱和区斜率向 0 偏移宁小勿大。3. 用 Python 实现 box-countingreshape 版核心函数与命令行3.1 zip 解压后的标准结构与最小运行命令box_count.zip 解压后最常见的布局是三件套box_count.py 一个自包含脚本requirements.txt 锁定 numpy 和 PillowREADME.md 写参数约定。先建环境跑一次最小命令python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy pillow python box_count.py sample.png --threshold 170 --csv result.csv命令里的 --threshold 是灰度二值化阈值灰度小于 170 的像素视为前景--csv 把每一档的 (size, count) 原始数据落盘。第一次跑一定要带 --csv因为只看最后输出的 D无法判断回归有没有踩到 2.3 说的饱和区。跑通之后再把 sample.png 换成自己的图。3.2 核心函数 box_count一次 reshape 统计全部盒子计数过程可以压缩成“裁剪 reshape any”三步不需要嵌套循环。完整可运行版本如下import numpy as np from PIL import Image def load_binary(path: str, threshold: int 128) - np.ndarray: 读图并二值化返回 bool 数组True 表示前景。 img Image.open(path).convert(L) arr np.array(img, dtypenp.uint8) return arr threshold def box_count(binary: np.ndarray, sizes) - list: 对每个边长 s返回 (s, 非空盒子数) 列表。 h, w binary.shape out [] for s in sizes: if s 2 or s h or s w: continue hc, wc h - h % s, w - w % s # 裁剪到 s 的整数倍 crop binary[:hc, :wc] # 高分成 hc//s 块每块内部 s 行宽同理 cells crop.reshape(hc // s, s, wc // s, s) occupied cells.any(axis(1, 3)) # 消掉盒内两个方向 out.append((s, int(occupied.sum()))) return out def fit_dimension(counts: list) - tuple: 对 (s, N) 做 log-log 回归返回 (D, intercept)。 sizes np.array([float(s) for s, _ in counts]) ns np.array([float(c) for _, c in counts]) valid (ns 1) (sizes 1) # 丢弃饱和点与 ε1 x np.log(1.0 / sizes[valid]) # 自变量 ln(1/ε) y np.log(ns[valid]) # 因变量 ln N(ε) slope, intercept np.polyfit(x, y, 1) return slope, intercept # slope 就是 D逻辑说明crop.reshape 把 (hc, wc) 数组看成四维张量四个轴依次是“行块号、盒内行、列块号、盒内列”用 any(axis(1, 3)) 一次性消掉两个盒内维度剩下的 (行块数, 列块数) 矩阵里 True 的数量就是 N(ε)。reshape 按行主序展开行块和列块的切分顺序正是先高后宽不需要额外转置。如果改成链式写法 any(axis1).any(axis2) 也等价但注意 3D 版本不能随便链式见第 6 章。参数说明sizes 通常给 2 的幂序列threshold 决定前景粗细曲线边缘的抗锯齿像素如果低于阈值会被算进前景D 系统性偏大valid 过滤只处理极端点过滤后只剩两三个点说明 max-size 设得太大。提示hc/wc 裁剪和 reshape 都只产生视图不复制底层数组内存开销接近原图本身4096×4096 的 bool 图约 16 MB可以放心跑。3.3 效率对比向量化 any 与三层 for 循环常见误用是写三层 for 循环对每个盒子遍历 s² 个像素判占用复杂度约等于 H·W×档位数但访存完全随机、缓存命中差reshape 版的 any 在连续内存块上做归约每个元素只读一次。同一台机器上的量级大致如下图像尺寸循环版耗时reshape 版耗时512×512约 2 秒约 15 毫秒1024×1024约 8 秒约 40 毫秒4096×4096约 2 分钟约 0.6 秒对单张图这点差距无所谓但批量处理几百张图时向量化直接决定任务能不能当晚跑完。这也是 box-counting 适合打成 zip 分发的原因核心逻辑就这么几十行依赖面窄换机器不会出现“性能玄学”。3.4 把 (size, count) 画出来调参前的第一件事D 只是一个回归斜率丢掉了很多信息。用 csv 落盘的数据画一张散点图很多问题一眼可见import csv import numpy as np import matplotlib.pyplot as plt rows list(csv.reader(open(result.csv)))[1:] # 跳过表头 sizes np.array([float(r[0]) for r in rows]) counts np.array([float(r[1]) for r in rows]) x np.log(1.0 / sizes) y np.log(counts) plt.plot(x, y, o-) plt.xlabel(ln(1/ε)) plt.ylabel(ln N(ε)) plt.show()实心方块的图应该是近似等间距的 8 个点最左端略微下弯是正常饱和最右端三点如果歪掉优先查阈值而不是拟合代码。matplotlib 只是检查工具建议不放进 box_count.py 的导入链服务器上没有图形环境时脚本不能因此报错。4. box size、边缘处理与阈值四个影响盒计数结果的参数4.1 box size 序列2 的幂还是线性递增理论上只要求 ε 趋于 0没规定怎么采样。实际影响在双对数坐标里2 的幂序列使 ln(1/ε) 等距分布回归时每个点的杠杆差不多线性递增序列会把大部分点挤在大盒子端小盒子端只剩一两个点拟合斜率偏向中段尺度对细小结构不敏感所以默认用 2^k。实操上我还会加一个缩比检查把同一张图缩到一半分辨率再算一次两次 D 的差超过 0.05说明最小特征尺寸和最小盒子已经可比需要换更高分辨率素材而不是调回归区间。显微镜图像尤其要养成这个习惯因为采样分辨率决定了可观测的最小尺度也决定了维数这个数算到哪一档截止。4.2 边缘处理裁剪、填充与输入规整图像边长不能被 s 整除时裁剪丢弃右/下边缘填充补背景。两种做法的偏差方向相反处理方式对 D 的偏差方向适用场景裁剪偏低丢边缘结构前景远离边缘计算要快补 0 填充偏高边界产生空盒子前景贴近边缘不能丢像素规整到 2 的幂几乎无偏采集阶段的最优解更稳的做法是让输入边长本身就是 2 的幂采集 ROI 时约定尺寸或处理前用Image.resize((512, 512))规整。注意缩放会引入抗锯齿必须先二值化再缩放缩放后若还想调阈值重新做一次 threshold。不要默认填充因为填充的偏差随盒子尺度漂移很难在回归阶段修正。4.3 二值化阈值与骨架化的边界阈值直接决定前景“粗细”。一条抗锯齿直线阈值取 100 会把半透明过渡像素算进前景等效线宽变大小盒子端 N(ε) 上升D 从 1.0 漂向 1.2。补救方法分两类灰度直方图双峰时用 Otsu 全局阈值skimage.filters.threshold_otsu 一行搞定对象本来就是线条结构道路网、血管、骨架时先 skeletonize 再计数否则线宽占多个像素会让 D 系统性偏高。骨架化也有代价细化算法会制造 2×2 方块和交叉点伪影这些局部结构在 ε 缩小时被放大。因此骨架化之后建议把最小 box size 从 2 提到 4让伪影落在回归区间外。判断标准始终是 5.3 的 log-log 图小盒子端上翘多半就是这类伪影。4.4 多网格偏移平均消除网格对齐偏差固定网格计数有个已知毛病45° 对角线这类结构在 s 较大时可能正好骑在网格线上非空盒子数周期性骤减log-log 图出现锯齿回归斜率随之摆动。补救方法是多网格法对同一个 s把网格原点平移几个偏移量分别计数取平均def box_count_offset(binary, s, offsetsNone): if offsets is None: offsets (0, s // 2) h, w binary.shape total 0 for dx in offsets: hc h - (h - dx) % s # 从 dx 起能整分的最大高度 wc w - (w - dx) % s # 宽度方向用同一偏移简单起见 crop binary[dx:dx hc, dx:dx wc] cells crop.reshape(hc // s, s, wc // s, s) total int(cells.any(axis(1, 3)).sum()) return s, total / len(offsets)逻辑说明偏移 dx 表示网格原点落在第 dx 行/列hc 取满足 hc ≡ dx (mod s) 的最大高度切片从 dx 开始取 hc 个元素末尾不残留。两次计数的平均值作为 N(ε) 参与回归对角线结构的锯齿被抹平。参数说明offsets 扩到 4 个四个象限后效果基本饱和每多一个偏移只是多一次 reshapeany总耗时线性增长对 4096×4096 也不敏感。如果回归点在某档位出现明显锯齿先单独验证这一档的偏移效果再决定是否全局启用。5. 用已知分形校准 box-countingCantor 集、Koch 曲线与 Sierpinski 地毯5.1 三个理论维数已知的生成函数验证实现的标准做法精确生成理论维数已知的结构跑同一套 box_count看结果偏差。三个生成函数只依赖 numpydef sierpinski_carpet(n5): n 阶 Sierpinski 地毯边长 3**n理论 D ln8/ln3 ≈ 1.8928。 carpet np.ones((1, 1), dtypebool) kernel np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]], dtypebool) for _ in range(n): carpet np.kron(kernel, carpet) # 旧像素扩展成 3x3 块中心挖空 return carpet def cantor_1d(n6): n 阶 Cantor 集长度 3**n理论 D ln2/ln3 ≈ 0.6309。 line np.ones(3 ** n, dtypebool) step 3 ** n for _ in range(n): step // 3 for start in range(0, 3 ** n, 3 * step): line[start step:start 2 * step] False return line def koch_points(depth, a, b): 递归生成 Koch 曲线点列depth0 时就是端点。 a, b np.asarray(a, dtypefloat), np.asarray(b, dtypefloat) if depth 0: return [a, b] v (b - a) / 3.0 p1, p2 a v, a 2 * v apex p1 np.array([v[0]/2 - v[1]*np.sqrt(3)/2, v[0]*np.sqrt(3)/2 v[1]/2]) pts [] for pa, pb in [(a, p1), (p1, apex), (apex, p2), (p2, b)]: pts koch_points(depth - 1, pa, pb)[:-1] return pts [b]逻辑说明np.kron(kernel, carpet) 把 carpet 的每个元素替换成 kernel 标定过的 3×3 块kernel 为 1 的位置放整块 carpet为 0 的位置全 False迭代 n 次得到 3^n 边长中心孔洞逐级保留。Cantor 集在每段迭代中去掉中间 1/3start 步长按当前段长 3*step 递进。Koch 曲线用递归在每段上构造等边三角形凸包depth4 时点数是 4^41 257 个。Koch 点列画到图像上再二值化PIL 自带画线不需要额外依赖from PIL import Image, ImageDraw def koch_image(depth4, size729): pts koch_points(depth, (0.0, 0.0), (1.0, 0.0)) img Image.new(L, (size, size), 255) dr ImageDraw.Draw(img) prev None for x, y in pts: cur (int(x * (size - 4)) 2, int(y * (size - 4)) 2) if prev is not None: dr.line([prev, cur], fill0, width1) prev cur return np.array(img) 128注意画布的 y 方向只用到约 0.29 的比例上半部分留白是正常的空盒子不计入 N(ε)不影响计数但盒子尺寸超过曲线实际跨度后 count 会迅速降到 1所以 Koch 样本的 box size 上限取 128 而不是 256。5.2 校准结果与判定标准用 3.2 的函数跑上述样本典型的实测范围如下具体值随 numpy 版本和画线抗锯齿略有浮动样本理论 D回归区间典型实测范围单像素水平直线1.00002..5121.001.0145° 对角线1.00002..2560.981.00Cantor 集6 阶0.63092..1280.620.64Sierpinski 地毯5 阶1.89282..1281.881.90Koch 曲线4 阶1.26192..1281.251.27判定标准实心方块必须严格收敛到 2.00这是实现 bug 的试金石直线必须贴近 1.00两个标准分形偏差在 ±0.02 内算合格。超过这个范围先回第 4 章查边缘处理和阈值不要急着改回归区间。对角线这类结构对网格对齐敏感没做 4.4 的偏移平均时 D 会周期性波动这是方法特性不是 bug。5.3 快速目检log-log 图上的两种系统偏差给一个配套绘图函数把拟合直线和原始点叠在一起看def plot_fit(counts): import matplotlib.pyplot as plt sizes np.array([float(s) for s, _ in counts]) ns np.array([float(c) for _, c in counts]) ok ns 1 x, y np.log(1.0 / sizes[ok]), np.log(ns[ok]) slope, ic np.polyfit(x, y, 1) xs np.linspace(x.min(), x.max(), 50) plt.plot(x, y, o, labelfD{slope:.3f}) plt.plot(xs, slope * xs ic, --, colorgray) plt.xlabel(ln(1/ε)) plt.ylabel(ln N(ε)) plt.legend() plt.show()看这张图只需要两个判断。点列在直线下方系统性下弯是大盒子端进入饱和砍掉最大的几个尺寸点列在小盒子端上翘是抗锯齿或骨架伪影在抬升 N(ε)把最小尺寸从 2 提到 4或回去调阈值。两种情况都别靠删点硬凑根因在输入不在拟合。6. 升级到 3D 体数据六维 reshape 与开机自检6.1 从四维 reshape 到六维三维盒计数盒计数不限于二维图。CT 扫描、多孔介质、泡沫材料的二值体数据用同样的覆盖思想只是把盒子换成正方体D 落在 0 到 3 之间。实现只需把第 3 章的四维 reshape 扩成六维并且注意 any 的轴编号要一次给全def box_count_3d(binary3d: np.ndarray, sizes): 输入 (H, W, D) 的 bool 数组输出 (s, count) 列表。 h, w, d binary3d.shape out [] for s in sizes: hc, wc, dc h - h % s, w - w % s, d - d % s crop binary3d[:hc, :wc, :dc] cells crop.reshape(hc // s, s, wc // s, s, dc // s, s) occupied cells.any(axis(1, 3, 5)) # 一次消掉三个盒内方向 out.append((s, int(occupied.sum()))) return out六轴顺序是“行块号、盒内行、列块号、盒内列、深块号、盒内深”。这里不能写成交替链式 any(axis1).any(axis3).any(axis5)因为每次归约后轴编号会左移直接链式会把块维度误消掉。两个 3D 特有约束体素必须各向同性医学影像 z 轴层厚和面内分辨率不同时先重采样成立方体素边缘剩余超过 10% 时先规整到 2 的幂不要填充理由同 4.2。6.2 给 zip 分发的脚本加一个自检模式最后说交付习惯。box_count.zip 这类脚本到了别人机器上最常见的问题不是装不起来而是“装起来但不敢信结果”。我一般会在main里加 --self-test用 5.1 的生成函数做断言入口命令python box_count.py --self-test对应分支如下# argparse 里加 # ap.add_argument(--self-test, actionstore_true) if args.self_test: carpet sierpinski_carpet(5) counts box_count(carpet, [2 ** k for k in range(1, 8)]) d, _ fit_dimension(counts) assert abs(d - 1.8928) 0.03, fself-test failed: D{d:.4f} print(fself-test OK, D{d:.4f})自检用生成函数而不是素材图片因为生成只依赖 numpy解压后零素材直接验证环境。断言阈值放 ±0.03比 5.2 的判定标准略宽给不同 numpy 版本留余量跑过自检再处理用户自己的图D 的可信度就从“脚本给的数”变成“自己验证过的数”。本文还有配套的精品资源点击获取
