分形维数法概念

分形维数法是一种用于描述和分析复杂形状和结构的方法。这些形状和结构可能看起来很不规则或混乱，但实际上，它们在某种程度上有一定的规律性。分形维数帮助我们量化这种规律性和复杂性。

具体展开来说，分形就是几何形状，它的某些部分与整体相似，即使放大了看，部分的形状依然和整体类似。例如，雪花、海岸线和树木的分枝都是一种分形。
传统的维数有一维（线）、二维（平面）和三维（立体）。分形维数是一种非整数维数，用来描述那些复杂度介于这些传统维数之间的形状。例如，一个分形的维数可能是1.5或2.7，这表示它比一条线复杂，但还没有达到一个平面的复杂程度。

分形维数可以用来描述自然界中的各种复杂形态，如山脉的轮廓、云的形状、植物的生长模式等；在图像处理、信号分析和材料科学等领域，分形维数可以分析和理解复杂系统和数据。

分形维数有多种计算方法，其中一种常见的方法是盒子计数法，这种方法通过在不同尺度下计算覆盖分形所需的盒子数量来估算维数。

Python程序实现

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, color, filters
from matplotlib.widgets import RectangleSelector
from osgeo import gdal


def differential_box_counting(img, min_box_size=5, num_boxes=20):
    assert len(img.shape) == 2
    img = (img - img.min()) / (img.max() - img.min())
    N = img.shape[0]
    sizes = np.logspace(np.log2(min_box_size), np.log2(N), num=num_boxes, base=2).astype(int)
    counts = []

    for size in sizes:
        reduced_image = measure_block_counting(img, size)
        counts.append(reduced_image)

    counts = np.array(counts)
    coeffs = np.polyfit(np.log(sizes), np.log(counts), 1)
    return -coeffs[0], sizes, counts


def measure_block_counting(Z, box_size):
    block_count = 0
    for i in range(0, Z.shape[0], box_size):
        for j in range(0, Z.shape[1], box_size):
            block = Z[i:i + box_size, j:j + box_size]
            if np.any(block):
                block_count += 1
    return block_count


def read_tif_image(image_path):
    dataset = gdal.Open(image_path)
    band = dataset.GetRasterBand(1)
    image = band.ReadAsArray()
    return image


class SelectRegion:
    def __init__(self, ax, img):
        self.ax = ax
        self.img = img
        self.rect_selector = RectangleSelector(ax, self.onselect, interactive=True)
        self.selected_region = None

    def onselect(self, eclick, erelease):
        x1, y1 = int(eclick.xdata), int(eclick.ydata)
        x2, y2 = int(erelease.xdata), int(erelease.ydata)
        self.selected_region = self.img[y1:y2, x1:x2]
        self.rect_selector.set_active(False)
        plt.close()


def main():
    img_path = 'img path'  # 替换为你的TIF文件路径
    img = read_tif_image(img_path)
    if img.ndim == 3 and img.shape[2] == 4:
        img = img[:, :, :3]
    if img.ndim == 3:
        img = color.rgb2gray(img)

    # 保存原图用于显示
    original_img = img.copy()

    # 边缘提取
    img = filters.sobel(img)

    fig, ax = plt.subplots(1, 1, figsize=(8, 8))
    ax.imshow(original_img, cmap='gray')  # 显示原图
    select_region = SelectRegion(ax, original_img)

    plt.show()

    if select_region.selected_region is not None:
        selected_img = select_region.selected_region

        # 使用预处理后的图像进行分形维数计算
        selected_processed_img = filters.sobel(selected_img)

        fd, sizes, counts = differential_box_counting(selected_processed_img, min_box_size=5, num_boxes=20)

        print(f'Fractal Dimension: {fd}')

        fig, axs = plt.subplots(1, 2, figsize=(12, 6))
        axs[0].imshow(selected_img, cmap='gray')  # 显示选定的原图区域
        axs[0].set_title(f'Selected Region\nFractal Dimension: {fd}')

        axs[1].plot(np.log(sizes), np.log(counts), 'o-', label='Data')
        axs[1].plot(np.log(sizes), np.polyval(np.polyfit(np.log(sizes), np.log(counts), 1), np.log(sizes)), 'r--', label=f'Fit: D={fd:.2f}')
        axs[1].set_xlabel('log(Box Size)')
        axs[1].set_ylabel('log(Count)')
        axs[1].legend()
        axs[1].set_title('Fractal Dimension Fitting')

        plt.tight_layout()
        plt.show()


if __name__ == '__main__':
    main()

补充路径之后，程序的具体流程如下：
打开一张图片，锚定一个区域：

可以自动得出区域的分形维数和维数（斜率）随盒子尺寸变化的拟合情况：

xlabel是盒子的尺寸，也就是2的n次方。拿其中一个点说，盒子大小为2的4次方（16）个像素点时，对应的要覆盖这个分形结构的盒子数（二值化之后包含非0像素的盒子数）。（因为我的程序是打开一张图片截取感兴趣的区域这样的思路，所以横坐标用这样的方式来设置，也有一些程序横坐标是分数值，这样是针对真个图片来说的，比如四分之一，表示盒子的边长为整个图片的四分之一）。
具体的分形过程可视化两张过程图看看吧，直观一些。



本文只介绍了盒计数分形法，比较常用。
具体原理和其他方法的介绍可以多看看相关文章。

吕泽昆.图像分形维数的计盒方法改进[J].电子技术与软件工程,2020,(11):144-145.
张志，董福安，伍友利.二维灰度图像的分形维数计算[J].计算机应用,2005,(12): 2853-2854
…