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.0、测试图像：skimage.data

from skimage import data

# 2D图像
camera = data.camera()
coins = data.coins()
brick = data.brick()
chessboard_GRAY = data.checkerboard()

# 3D图像
cat = data.cat()
chelsea = data.chelsea()
coffee = data.coffee()

astronaut = data.astronaut()
rocket = data.rocket()
retina = data.retina()

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()

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

• JPEG (.jpg)
  • 通常是三通道（RGB）：主要使用红色、绿色和蓝色（RGB）三个通道来存储颜色信息。
  • 不支持透明度：JPEG不具备Alpha通道，因此不支持透明背景。
  • YCbCr色彩空间：在压缩过程中，JPEG图像通常使用YCbCr色彩空间，将图像分为亮度（Y）和两个色度分量（Cb、Cr），便于更高效的压缩。
  • 单通道图像（灰度图像）：JPEG格式可以存储单通道灰度图像，适用于只需黑白信息的情况。

• PNG (.png)
  • 通常是三通道（RGB）或四通道（RGBA）：可以存储RGB三通道或包含Alpha通道的四通道（RGBA）图像。
  • 支持透明度通道：与JPEG不同，PNG支持透明图像，因此常用于需要透明背景的图像场景。
  • 双通道图像（灰度图像与Alpha通道）：支持将灰度图像与Alpha通道结合，形成两通道（GA）图像。
  • 单通道图像（灰度图像）：PNG可以是单通道灰度图像，适用于黑白图像。

import skimage
import numpy as np

# （1）读取四通道图像RGBA=[245, 612, 4]
rgba = skimage.io.imread(r"F:\py\color_image.png")

# （2）若为三通道图像RGB =[245, 612, 3]，则手动添加透明度（Alpha通道）
if rgba.shape[-1] == 3:
    alpha = 1  # 透明度=0.5
    # 创建 Alpha 通道 A=[245, 612, 1]
    alpha_channel = np.full((rgba.shape[0], rgba.shape[1], 1), fill_value=int(255*alpha), dtype=np.uint8)
    rgba = np.concatenate([rgba, alpha_channel], axis=-1)  # 将 Alpha 通道与 RGB 图像合并，生成 RGBA 图像

# （3）颜色空间转换
rgb = skimage.color.rgba2rgb(rgba)  # rgba转rgb
gray = skimage.color.rgb2gray(rgb)  # rgb转gray

skimage.io.imshow(rgb)              # 显示图像
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()

import matplotlib.pyplot as plt
from skimage import data, filters, measure, color

# （1）加载示例图像（硬币图像）
image = data.coins()

# （2）图像处理
thresh_value = filters.threshold_otsu(image)    # 使用Otsu方法计算全局阈值
binary_image = image > thresh_value             # 生成二值图像（大于阈值的像素为前景，小于阈值的为背景）
label_image = measure.label(binary_image)       # 对二值图像进行标签标记
regions = measure.regionprops(label_image)      # 使用regionprops函数提取标记区域的属性
print(f"thresh_value = {thresh_value}")

# （3）输出每个区域的属性
for region in regions:
    print(f"Region {region.label}:")
    print(f" - Area: {region.area}")
    print(f" - Bounding Box: {region.bbox}")
    print(f" - Centroid: {region.centroid}")

# （4）可视化结果
plt.subplot(131), plt.imshow(image, 'gray'), plt.title('image')
plt.subplot(132), plt.imshow(label_image, 'gray'), plt.title('label_image')
plt.subplot(133), plt.imshow(color.label2rgb(label_image, image=image, bg_label=0)), plt.title('rgb_image')
plt.show()
# color.label2rgb(label_image, image=image, bg_label=0): 用于将标记图像转换为伪彩色图像
# （将每个标记区域以不同颜色显示，并且可以叠加在原始灰度图像上进行对比。）


"""###########################################################################
函数介绍：大津阈值法 - 根据图像的直方图自适应获取最佳的二值化阈值（前景+背景）
函数说明：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×2 个像素的线性加权平均计算插值结果。
                        优点：计算速度较快，生成的图像比最近邻更平滑；
                        缺点：细节表现一般，可能会出现轻微的模糊。
                2: 双二次插值（Bi-quadratic）
                        特点：使用 3×3 邻域内的像素通过二次多项式拟合进行插值。
                        优点：在保持较高平滑度的同时能改善细节；
                        缺点：计算量比双线性略大，效果不如双三次常见。
                3: 双三次插值（Bi-cubic）
                        特点：采用 4×4 邻域像素，通过三次多项式插值计算结果。
                        优点：能较好地平衡图像平滑度与细节保留，是较为常用且效果优良的插值方法；
                        缺点：计算量较双线性大，处理速度稍慢。
                4: 双四次插值（Bi-quartic）
                        特点：使用较大邻域（例如 5×5）并采用四次多项式进行插值。
                        优点：能获得更高的平滑度和连续性；
                        缺点：计算复杂度进一步增加，处理速度较慢，且在某些情况下可能过于平滑。
                5: 双五次插值（Bi-quintic）
                        特点：使用更大邻域（例如 6x6）并采用五次多项式进行插值。
                        优点：在理论上能够更准确地重构连续图像；
                        缺点：计算量大，处理速度慢，而且过高阶的插值有时可能引入不必要的过渡效果。
		（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阶数越高，精度越高，但耗时越长。
##############################################################################################################"""