skimage图像处理（全）

---

一、简介

Scikit-image 官方网站
使用建议：搜索想要了解的函数，会有非常详细的解释以及函数源码。

skimage（全称scikit-image）：是一个基于 scipy 和 numpy 构建的，专注于图像处理和计算机视觉的Python库。提供了一系列用于图像处理函数，包括图像滤波、边缘检测、图像分割、特征提取、图像配准等。

• scipy 是一个基于 numpy 构建的，具有更广泛通用的数学计算和科学计算的Python库。提供了优化、信号处理、统计、线性代数、数值积分、插值、稀疏矩阵等领域的函数。
• 优势：同时支持2D和3D图像处理（与单帧图像处理相比，极大降低耗时问题，但内存占比高）。相比于 OpenCV、Pillow、scipy 等专注于2D图像处理库，skimage 更加专业和全面。

二、安装

（1）pip 安装：pip install scikit-image
（2）conda 安装：conda install scikit-image

import skimage
print(skimage.__version__)  # 打印版本号

# 使用命令行，查看版本号
python -c "import skimage; print(skimage.__version__)"

三、模块简介：API reference

skimage包由许多的子模块组成，各个子模块提供不同的功能。

"""
主要子模块列表如下：
	（1）skimage.io               读取、保存和显示（图片或视频）
	（2）skimage.data             提供一些测试图片和样本数据
	（3）skimage.color            颜色空间转换	（如：rgb2gray、gray2rgb）
	（4）skimage.draw             图形绘制		（如：线条、矩形、圆等）
	（5）skimage.filters          图像过滤		（如：增强、降噪、边缘检测、纹理分析等）
	（6）skimage.transform        几何变换		（如：修改图像尺寸、旋转、裁剪、仿射变换等）
	（7）skimage.morphology       形态学操作		（如：腐蚀、膨胀、开闭运算等）
	（8）skimage.exposure         强度调整		（如：亮度调整、直方图均衡等）
	（9）skimage.feature          特征提取		（如：直方图、角点检测、模板匹配等）
	（10）skimage.graph			  图形结构		（如：最小生成树、连通性、骨架等）
	（11）skimage.measure         属性测量		（如：相似性或等高线等）
	（12）skimage.segmentation    图像分割		（将图像分成不同的部分或对象）
	（13）skimage.restoration     图像恢复和去噪	（如：去卷积、反卷积等）
	（14）skimage.util            通用函数		（如：类型转换等）
"""

四、项目实战

4.1、2D图像处理

python skimage 2D图像处理

4.1.1、打印图像属性

import skimage

img = skimage.io.imread('image.jpg')    # 读取图像
print("类型 =", type(img))
print("尺寸 =", img.shape)
print("高度 =", img.shape[0])
print("宽度 =", img.shape[1])
print("通道数 =", img.shape[2])
print("像素总个数 =", img.size)
print("最大像素值 =", img.max())
print("最小像素值 =", img.min())
print("像素平均值 =", img.mean())

4.1.2、读取 / 显示 / 保存图像：skimage.io.imread() + skimage.io.imshow() + skimage.io.imsave()

import skimage

img = skimage.io.imread('image.jpg')  # 读取图像
skimage.io.imshow(img)                # 显示图像
skimage.io.show()                     # 阻塞程序，以便图像可以正常显示。
skimage.io.imsave('image1.png', img)  # 保存图像

4.1.3、颜色空间转换：skimage.color.rgba2rgb + skimage.color.rgb2gray()

更多颜色空间转换函数，请看官网

import skimage

rgba = skimage.io.imread('image.jpg')                                        # 读取图像
rgb = skimage.color.rgba2rgb(rgba, background=(1, 1, 1), channel_axis=-1)    # rgba转rgb
gray = skimage.color.rgb2gray(rgb, channel_axis=-1)                          # rgb转gray
skimage.io.imshow(gray)     # 显示图像
skimage.io.show()           # 阻塞程序，以便图像可以正常显示。

4.2、3D图像处理

4.2.1、相位一致性：skimage.feature.hessian_matrix() + skimage.feature.hessian_matrix_eigvals()

相位一致性（Phase Congruency, PC）：在信号处理或系统分析中，保持信号或系统的相位信息在不同频率下保持不变的性质。

• 定义：在频率域分析中，信号（或图像）可以表示为幅度谱和相位谱。
  • 频谱（频率谱）：表示信号在不同频率上的强度。
  • 幅度谱（Amplitude Spectrum）：表示信号在频率域上的振幅或幅度分布；
  • 相位谱（Phase Spectrum）：表示信号在频率域上的相位信息； —— 决定了信号或图像的位置和特征的空间分布
• 作用：在频域分析中，通常对信号进行傅里叶变换来获取其频谱，然后可以从频谱中提取相位谱和幅度谱。若不保持相位一致性，则频谱的相位信息被破坏，导致图像模糊、失真或者无法还原原始信号。
• 应用：图像配准、图像拼接和图像重建。
  基于相位一致性的边缘检测

相位一致性算法步骤：

• （1）归一化处理：将像素值归一化到0到1之间，以便后续计算的稳定性。
• （2）矩阵计算：对图像进行Hessian矩阵计算，获取每个像素点的Hessian矩阵。Hessian矩阵描述了图像的局部二阶导数信息，可以用于描述图像的曲率和边缘信息。
• （3）特征值计算：从Hessian矩阵中计算特征值和特征向量。特征值表示了Hessian矩阵在不同方向上的曲率，特征向量表示了对应于特征值的方向。在相位一致性算法中，我们主要关注特征值的大小。
• （4）相位一致性计算：通过对特征值的计算，得到每个像素点处的相位一致性值。相位一致性值反映了图像局部结构和纹理信息的一致性程度。边缘和纹理区域通常具有较高的相位一致性值，而平坦区域或噪声区域通常具有较低的相位一致性值。
• （5）后处理：根据相位一致性值进行一些后处理操作，如阈值化或非极大值抑制，以便进一步提取感兴趣的边缘和纹理特征。

import numpy as np
import tifffile
import napari
import time

from skimage.feature import hessian_matrix, hessian_matrix_eigvals, canny
from skimage.filters import median


def phase_congruency(image, k=2.0, f=2.0):
    image = median(image)

    edges = np.zeros_like(image)
    for i in range(image.shape[0]):
        edges[i] = canny(image[i], low_threshold=50, high_threshold=200)

    start_time = time.time()  # 记录开始时间
    #######################################################
    H = hessian_matrix(image, sigma=2.0, mode='reflect')
    eigvals = hessian_matrix_eigvals(H)
    #######################################################
    end_time = time.time()  # 记录结束时间
    runtime = end_time - start_time  # 计算运行时间
    print(runtime)

    lambda1 = eigvals[0, ...] * -1  # First eigenvalue
    lambda2 = eigvals[1, ...] * -1  # Second eigenvalue
    lambda3 = eigvals[2, ...] * -1  # Third eigenvalue

    phase_congruency = np.sqrt(lambda1 ** 2 + lambda2 ** 2 + lambda3 ** 2) / (np.abs(lambda1) + np.abs(lambda2) + np.abs(lambda3))
    # Enhance edges using Canny result for each slice
    pc_image = phase_congruency*0.5 + (edges)
    pc = pc_image ** k * f
    return pc


def main():
    image_path = r'SampleImage.tif'
    image = tifffile.imread(image_path)
    image = image.astype(np.float32)
    image /= np.max(image)
    pc_image = phase_congruency(image)

    viewer = napari.Viewer()
    viewer.add_image(pc_image, name='Phase Congruency', colormap='gray')
    viewer.add_image(image, name='raw Image', colormap='gray')
    napari.run()


if __name__ == '__main__':
    main()

"""##################################################################################################
函数作用：用于计算图像局部二阶导数的方法，得到图像的Hessian矩阵。Hessian矩阵可用于检测图像中的特定结构，如边缘、角点和斑点等。
函数说明：H = skimage.feature.hessian_matrix(image, sigma=1.0, order='rc')
输入参数：
        （1）image：输入的灰度图像。该图像必须是二维的（单通道）。
        （2）sigma：用于控制高斯滤波的标准差，影响图像平滑效果。默认值为1.0。
                - 较小的sigma值会使高斯滤波器的响应更加锐利，可以提取细节特征，但同时也会增加噪声的敏感性。
                - 较大的sigma值会使高斯滤波器的响应更加平滑，可以提取更大尺度的特征，但可能会导致细节丢失。
        （3）order：表示输出的Hessian矩阵的排列顺序。
                'rc'：按照rc(row-major, 行主)顺序排列，即按行优先排列。
                'xy'：按照xy顺序排列，即按照x和y方向的导数排列。
输出参数：
        H：图像的Hessian矩阵。对于二维图像，输出是一个形状为(2, 2, M, N)的数组，其中M和N是输入图像的高度和宽度。
##################################################################################################
函数作用：用于计算图像的Hessian矩阵的特征值（Eigenvalues），特征值是Hessian矩阵的特征向量对应的值，可以用于描述图像局部结构的特性。
函数说明：Eigenvalues = skimage.feature.hessian_matrix_eigvals(H)
输入参数：   H：图像的Hessian矩阵。———— 对于二维图像，H是一个形状为(2, 2, M, N)的数组，其中M和N是输入图像的高度和宽度。
输出参数：   Hessian矩阵的特征值数组。———— 对于每个像素，输出是一个形状为(2, M, N)的数组，其中M和N是输入图像的高度和宽度。
                                        Eigenvalues[0]对应第一个特征值，Eigenvalues[1]对应第二个特征值。
##################################################################################################"""

4.2.2、连通区域：skimage.filters.threshold_otsu() + skimage.measure.label() + skimage.measure.regionprops()

from skimage import filters, measure
import napari
import tifffile
import numpy as np

# （1）读取图像
image_stack = tifffile.imread(r'F:\py\561result-1-part.tif')  # [通道, 高度, 宽度]
gray_image = image_stack.copy()  # 复制图像
################################################################
# （2）获取连通区域
threshold = filters.threshold_otsu(gray_image)  # 大津阈值法自动获取灰度阈值
binary_image = gray_image > threshold  # 二值化图像
labeled_image, return_num = measure.label(binary_image, return_num=True)  # 标记连通区域
regions = measure.regionprops(labeled_image)  # 获取连通区域的属性
print("区域数量 =", return_num)
################################################################
# （3）筛选指定形状 + 绘图（原图+标记图像）
max_size = 12  # 筛选形状的边长最大值
min_size = 3   # 筛选形状的边长最小值
num_size = 0   # 符合要求的形状数量
for region in regions:
    # （3.1）获取区域属性
    bbox = region.bbox  # 获取细胞的边界框
    mask = region.image  # 获取区域的像素值掩码

    # （3.2）获取区域边界框（最小行、最大行、最小列、最大列、最小切片、最大切片）———— 每个值都对应原图中的位置（region）
    min_slice_box, min_row_box, min_col_box, max_slice_box, max_row_box, max_col_box = bbox
    cur_cell_width_box = max_col_box - min_col_box
    cur_cell_height_box = max_row_box - min_row_box
    cur_cell_depth_box = max_slice_box - min_slice_box
    cur_cell_size_box = [cur_cell_width_box, cur_cell_height_box, cur_cell_depth_box]

    # （3.3）提取区域数据
    region_data1 = image_stack[min_slice_box: max_slice_box, min_row_box: max_row_box, min_col_box: max_col_box]
    region_data2 = np.where(mask, region_data1, 0)  # 将区域内的像素值提取出来，只保留mask中为True的像素位置的值

    # （3.4）筛选符合要求的形状
    if max(cur_cell_size_box) <= max_size and min(cur_cell_size_box) >= min_size:
        num_size += 1
        gray_image[min_slice_box:max_slice_box + 1, min_row_box:max_row_box + 1, min_col_box:max_col_box + 1] = 65535  # 8-bit:255    16-bit:65535
        labeled_image[min_slice_box:max_slice_box + 1, min_row_box:max_row_box + 1, min_col_box:max_col_box + 1] = 65535

# （4）保存结果
print("shape_num =", num_size)
labeled_image = labeled_image.astype(gray_image.dtype)  # 将数据类型转换为指定数据类型.如:uint8
tifffile.imwrite("gray_image.tif", gray_image)  # 保存图像为TIFF
tifffile.imwrite("labeled_image.tif", labeled_image)  # 保存图像为TIFF

# （5）在napari中显示图像
viewer = napari.Viewer()  # 创建napari视图
viewer.add_image(image_stack, name="image_stack")  # 添加图像（原始图像）
viewer.add_image(gray_image, name="gray_image")  # 添加图像（在原始图像中，显示标记图像）
viewer.add_image(labeled_image, name="labeled_image")  # 添加图像（在二值化图像中，显示标记图像）
napari.run()  # 显示napari图形界面

"""###########################################################################
函数介绍：大津阈值法 - 根据图像的直方图自适应获取最佳的二值化阈值（前景+背景）
函数说明：threshold = skimage.filters.threshold_otsu(image, nbins=256)
输入参数：
        image：待阈值处理的输入图像。可以是灰度图像或多通道图像。灰度图像将直接使用，多通道图像会首先转换为灰度图像。
        nbins：（可选参数）用于指定图像直方图的 bin 数量，默认为 256。增加 bin 数量可能会更准确地计算阈值，但会增加计算开销。
输出参数：
        threshold：计算得到的 Otsu 阈值，用于对输入图像进行二值化处理。
###########################################################################
函数介绍：连通区域标记 - 对二值化图像的每个连通区域进行唯一标记，并用不同的整数值表示。
函数说明：labeled_image, num_labels = skimage.measure.label(image, connectivity=1, return_num=True)
输入参数：
        image：输入的二值化图像。其像素值只有 True 或 False，表示图像的前景和背景。
        connectivity：（可选参数）用于指定连通区域的连接方式。默认为 1，表示使用 4 连通（4-connectivity），即上下左右四个方向相邻的像素被认为是连通的。可以设置为 2，表示使用 8 连通（8-connectivity），即包括对角线方向的像素。
        return_num：（可选参数）控制是否返回连通区域的数量。默认为 False，不返回。如果设置为 True，将返回连通区域的数量。
输出参数：
        labeled_image：连通区域标记后的图像
        num_labels：如果设置了 return_num=True，将返回连通区域的数量。
###########################################################################
函数介绍：属性计算 - 获取连通区域的属性：面积 + 周长 + 重心 + 边界框。
函数说明：props = skimage.measure.regionprops(label_image, intensity_image=None, cache=True)
输入参数：
        label_image：连通区域标记后的图像，通常是通过 skimage.measure.label() 函数得到的标记图像。
        intensity_image：（可选参数）输入的原始图像（灰度图像或彩色图像），用于计算连通区域的强度相关属性。默认为 None，表示不计算强度属性。
        cache：（可选参数）控制是否缓存标记图像。默认为 True，表示对标记图像进行缓存，以便后续的处理。如果设置为 False，则每次调用 regionprops 函数时都会重新计算标记图像。
输出参数：
        regions：一个包含连通区域属性的列表。每个连通区域都对应一个字典，其中包含了该区域的各种属性。

获取属性：
        for region in regions:
			label 	= region.label		区域的唯一标签或标识符。用于标识不同的连通区域。
			image 	= region.image		区域的像素值掩码，通常是一个与原始图像相同大小的布尔数组，True表示该区域非零值。
			bbox 		= region.bbox		表示区域边界框的坐标范围(min_slice, min_row, min_col, max_slice, max_row, max_col)
			coords 		= region.coords		区域内非零像素的坐标列表。
			centroid 	= region.centroid	区域的中心点坐标(centroid_slice, centroid_row, centroid_col)

			perimeter 	= region.perimeter			区域的周长。
			area 		= region.area				区域的面积，表示区域内非零像素的数量。
			convex_area	= region.convex_area		区域的凸包面积。
			equivalent	= region.equivalent_diameter与区域相同面积的圆的直径。

			orientation = region.orientation		区域的主要轴的方向。
			major 		= region.major_axis_length	区域的主要轴的长度。
			minor 		= region.minor_axis_length	区域的次要轴的长度。		
			moments 	= region.moments			区域的矩属性，用于计算区域的矩，包括一阶矩、二阶矩、中心矩、归一化矩等。
			eccentricity = region.eccentricity		区域的离心率，表示区域的形状。
			solidity 	= region.solidity			区域的实心度，表示区域内部像素占整个凸包的比例。
###########################################################################"""

4.2.3、调整图像（缩放+插值）：skimage.transform.resize()

若order=0（最近邻插值），则输出图像中的每个像素都是原始图像中最近邻像素的值，没有使用插值。

import numpy as np
from skimage import transform
import time

original_image = np.random.random((100, 200, 300))  # 随机生成示例数据
target_size = (200, 400, 600)  # 设置目标长宽高
###############################################
# order阶数越高，精度越高，但耗时越长。
###############################################
for ii in range(6):
    t = time.time()
    # 调整图像: 缩放+插值
    resized_image = transform.resize(original_image, target_size, mode='constant', cval=0, order=ii)
    print("总共耗时:", time.time() - t)
print("原始图像 =", original_image.shape)
print("缩放图像 =", resized_image.shape)
"""
原始图像 = (100, 200, 300)
缩放图像 = (200, 400, 600)

总共耗时: 0.3540208339691162
总共耗时: 2.0999197959899902
总共耗时: 4.99407434463501
总共耗时: 10.33707880973816
总共耗时: 19.263009786605835
总共耗时: 32.0386061668396
"""

"""##############################################################################################################
函数说明：skimage.transform.resize(image, output_shape, mode='reflect', cval=0, order=1, clip=True,
                                preserve_range=False, anti_aliasing=None, anti_aliasing_sigma=None)
输入参数：
		（1）image：待调整大小的输入图像。
                2D图像 = (height, width)；
                3D图像 = (depth, height, width)。
		（2）output_shape：目标输出图像的大小。
                2D图像 = (new_height, new_width)；
                3D图像 = (new_depth, new_height, new_width)。				
		（3）mode（可选）：插值方法。
				'constant': 	在图像边界外使用常数值填充。
				'edge': 		在图像边界外使用图像边缘像素值填充。
				'symmetric': 	在图像边界外使用图像的对称像素值填充。
				'reflect': 		在图像边界外使用图像的镜像像素值填充。
				'wrap': 		在图像边界外使用图像的环绕像素值填充。
		（4）cval（可选）：当mode='constant'时使用，指定填充的常数值（默认0）。
		（5）order（可选）：插值阶数（默认1）。阶数越高，精度越高，但耗时越长。
                0: 最近邻插值（Nearest-neighbor）
                1: 双线性插值（Bi-linear）
                2: 双二次插值（Bi-quadratic）
                3: 双三次插值（Bi-cubic）
                4: 双四次插值（Bi-quartic）
                5: 双五次插值（Bi-quintic）
		（6）clip（可选）：	            	指定是否将输出图像剪裁到指定输出范围内（默认True）。如：将输出图像限制在0到1之间。
		（7）preserve_range（可选）：		指定是否保持图像数据类型的范围。True表示输出图像保持与原始图像的数据类型。False表示输出图像为浮点数类型（float64）。
		（8）anti_aliasing（可选）：			指定是否启用抗锯齿处理。True表示将应用抗锯齿过滤器来减少下采样导致的混叠效应。
		（9）anti_aliasing_sigma（可选）：	指定抗锯齿过滤器的标准差 [默认0.01 * max(new_shape)]。如果anti_aliasing设置为True，则该参数用于指定抗锯齿过滤器的标准差。
输出参数：
		返回调整后的图像。
##############################################################################################################
请注意：
（1）skimage.transform.resize函数通常用于简单的插值操作。对于更复杂的3D图像处理任务，可能需要使用其他库，如SimpleITK、ITK或PyTorch等。
（2）order阶数越高，精度越高，但耗时越长。
##############################################################################################################"""