什么是边缘检测 ?

边缘检测：边缘检测是图像处理中的一种技术，旨在识别和定位图像中亮度变化显著的区域。这些变化通常表现为强度、颜色或纹理的急剧变动，通常对应于物体的边界或轮廓。通过检测这些边缘信息，边缘检测能够帮助提取图像的结构特征，促进图像的分析和处理。

边缘检测的目的
（1）特征提取：帮助提取图像中重要的结构特征，便于进一步的分析和处理。
（2）图像分割：通过识别边缘来区分图像中的不同区域，实现图像的分割。
（3）增强图像质量：通过突出重要特征，提高图像的可读性和可分析性。

边缘检测的应用
（1）图像分析：用于对象检测、跟踪和识别。
（2）医学影像：在CT、MRI等医学图像中用于识别和分析结构。
（3）计算机视觉：为自动驾驶、视频监控等应用提供关键数据。
（4）图形设计：用于图像的自动轮廓提取和处理。

边缘检测的性能通常通过准确度、鲁棒性（对噪声的抵抗力）和计算效率来评估。选择合适的边缘检测算法需考虑具体应用的需求和图像特征。

一、边缘检测算子：Sobel算子、Scharr算子、Laplacian算子、Canny算子

OpenCV图像处理（全）

import cv2  # opencv    读取图像默认为BGR
import tifffileimage = tifffile.imread(r"F:\py\gray_image.tif")
################################################################################
sobel_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)  # 计算x方向的导数
sobel_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)  # 计算y方向的导数
sobel_magnitude = cv2.magnitude(sobel_x, sobel_y)  # 计算梯度幅值
sobel_magnitude_Abs = cv2.convertScaleAbs(sobel_magnitude)
################################################################################
scharr_x = cv2.Scharr(image, cv2.CV_64F, 1, 0)  # 计算x方向的导数
scharr_y = cv2.Scharr(image, cv2.CV_64F, 0, 1)  # 计算y方向的导数
scharr_magnitude = cv2.magnitude(scharr_x, scharr_y)  # 计算梯度幅值
scharr_magnitude_Abs = cv2.convertScaleAbs(scharr_magnitude)  # 计算绝对值
################################################################################
laplacian = cv2.Laplacian(image, cv2.CV_64F)
laplacian_Abs = cv2.convertScaleAbs(laplacian)
################################################################################
image_canny50_100 = cv2.Canny(image, 50, 100)
image_canny100_200 = cv2.Canny(image, 100, 200)
################################################################################
import napari
viewer = napari.Viewer()  # 创建napari视图
viewer.add_image(image, name="image")
viewer.add_image(sobel_magnitude_Abs, name="sobel_magnitude_Abs")
viewer.add_image(scharr_magnitude_Abs, name="scharr_magnitude_Abs")
viewer.add_image(laplacian_Abs, name="laplacian_Abs")
viewer.add_image(image_canny50_100, name="image_canny50_100")
viewer.add_image(image_canny100_200, name="image_canny100_200")
viewer.grid.enabled = True  # 启用网格视图模式
napari.run()  # 启动napari事件循环，使得窗口保持打开并可交互。

二、梯度计算 + 顶帽 + 黑帽 + 拉普拉斯金字塔

import cv2                          # opencv    读取图像默认为BGR
import numpy as np
import tifffileimage = tifffile.imread(r"F:\py\gray_image.tif")
################################################################################
kernel = np.ones((10, 10), np.uint8)  # 初始化卷积核（np.ones: 生成一个数值全为1的3x3数组）
image_grad = cv2.morphologyEx(image, cv2.MORPH_GRADIENT, kernel)   # 梯度计算
image_top  = cv2.morphologyEx(image, cv2.MORPH_TOPHAT, kernel)     # 顶帽运算
image_black = cv2.morphologyEx(image, cv2.MORPH_BLACKHAT, kernel)  # 黑帽运算image_up = cv2.pyrUp(image)            # 高斯金字塔：图像上采样（放大一倍）
image_up_down = cv2.pyrDown(image_up)  # 拉普拉斯金字塔：先上采样，再下采样
image_up_down_resize = cv2.resize(image_up_down, (image.shape[1], image.shape[0]))  # 确保尺度一致
image_up_down_laplacian = cv2.subtract(image, image_up_down_resize)  # 拉普拉斯金字塔：当前层减去上采样image_down = cv2.pyrDown(image)        # 高斯金字塔：图像下采样（缩小一倍）
image_down_up = cv2.pyrUp(image_down)  # 拉普拉斯金字塔：先下采样，再上采样
image_down_up_resize = cv2.resize(image_down_up, (image.shape[1], image.shape[0]))  # 确保尺度一致
image_down_up_laplacian = cv2.subtract(image, image_down_up_resize)  # 拉普拉斯金字塔：当前层减去上采样
################################################################################
import napari
viewer = napari.Viewer()  # 创建napari视图
viewer.add_image(image, name="image")
viewer.add_image(image_grad, name="image_grad")
viewer.add_image(image_top, name="image_top")
viewer.add_image(image_black, name="image_black")
viewer.add_image(image_up_down_laplacian, name="image_up_down_laplacian")
viewer.add_image(image_down_up_laplacian, name="image_down_up_laplacian")
viewer.grid.enabled = True  # 启用网格视图模式
napari.run()  # 启动napari事件循环，使得窗口保持打开并可交互。

三、相位一致性（Phase Congruency，PC）

相位一致性（Phase Congruency），用于检测图像中的结构信息，如边缘、角点、纹理等。

3.1、底层代码（2D）

Phase Based Feature Detection and Phase Congruency

import sys
sys.dont_write_bytecode = True  # Python 解释器在运行脚本时不生成 .pyc 文件
#####################################################################
import math
import numpy as np
import cv2  # Faster Fourier transforms than NumPy and Scikit-Image
import tifffiledef phase_congruency(input_image, num_scales, num_angles):"""函数说明：计算图像的相位一致性输入参数：input_image: 输入图像，应为二维的NumPy数组。num_scales: 环形带通滤波器的尺度数量。如：4num_angles: 环形带通滤波器的角度数量。如：6输出参数：M（归一化的相位一致性矩阵，在0~1范围之间）"""############################################################################################################# （1）初始化参数和变量min_wave_length = 3   # 最小波长，用于滤波器设计scale_multiplier = 2.1  # 滤波器比例因子gaussian_std = 0.55  # 频率域高斯滤波器的标准差noise_threshold_param = 2.0  # 噪声阈值参数sigmoid_cutoff = 0.5  # Sigmoidal加权函数的截断值sigmoid_slope = 10  # Sigmoidal加权函数的斜率noise_estimation_method = -1  # 噪声估计方法（-1表示使用median方法）epsilon = .0001  # 防止分母为零的小值############################################################################################################# （2）构建环形带通滤波器：由多个尺度的Log-Gabor滤波器和角度高斯函数组合而成，在不同的频率和方向上对图像进行滤波。f_transform = cv2.dft(np.float32(input_image), flags=cv2.DFT_COMPLEX_OUTPUT)  # 离散傅里叶变换，得到频域复数图像nrows, ncols = input_image.shape  # 获取输入图像的行数和列数# 初始化一些中间结果矩阵和累加器zero_matrix = np.zeros((nrows, ncols))complex_response = np.zeros((nrows, ncols, num_scales, num_angles), dtype=complex)phase_congruency_matrix = np.zeros((nrows, ncols, num_angles))cov_x2 = np.zeros((nrows, ncols))cov_y2 = np.zeros((nrows, ncols))cov_xy = np.zeros((nrows, ncols))energy_vector = np.zeros((nrows, ncols, 3))pc_sum = np.zeros((nrows, ncols))# 矩阵的半径cy = math.floor(nrows / 2)  # 图像中心的y坐标cx = math.floor(ncols / 2)  # 图像中心的x坐标y_coords, x_coords = np.mgrid[0:nrows, 0:ncols]y_coords = (y_coords - cy) / nrows  # 归一化y坐标，以中心为原点x_coords = (x_coords - cx) / ncols  # 归一化x坐标，以中心为原点radius_matrix = np.sqrt(x_coords ** 2 + y_coords ** 2)  # 构建图像中每个像素点到中心的距离矩阵radius_matrix[cy, cx] = 1  # 中心点的半径设为1，避免分母为零theta_matrix = np.arctan2(-y_coords, x_coords)  # 构建每个像素点到中心的角度矩阵（范围：-pi 到 pi）sin_theta = np.sin(theta_matrix)  # sin(theta)矩阵cos_theta = np.cos(theta_matrix)  # cos(theta)矩阵# 初始化一组环形带通滤波器annular_bandpass_filters = np.empty((nrows, ncols, num_scales))# 设置滤波器参数d_theta_on_sigma = 1.3filter_orientations = np.arange(start=0, stop=math.pi - math.pi / num_angles, step=math.pi / num_angles)# 角度高斯函数的标准偏差，用于在频率面构建滤波器theta_sigma = math.pi / num_angles / d_theta_on_sigma# 初始化滤波器和奇偶波的变量bandpass_filters = np.empty((nrows, ncols, num_scales, num_angles))even_wavelets = np.empty((nrows, ncols, num_scales, num_angles))odd_wavelets = np.empty((nrows, ncols, num_scales, num_angles))# 实现对数Gabor传递函数的方法filter_order = 15  # filter 'sharpness'lowcut = .45norm_radius = radius_matrix / (abs(x_coords).max() * 2)lowpass_butterworth = 1.0 / (1.0 + (norm_radius / lowcut) ** (2 * filter_order))for s in np.arange(num_scales):wavelength = min_wave_length * scale_multiplier ** s  # 计算每个尺度的波长center_frequency = 1.0 / wavelength  # 计算每个尺度的中心频率log_gabor_filter = np.exp((-(np.log(radius_matrix / center_frequency)) ** 2) / (2 * math.log(gaussian_std) ** 2))  # 计算环形带通滤波器annular_bandpass_filters[:, :, s] = log_gabor_filter * lowpass_butterworth  # 应用低通滤波器annular_bandpass_filters[cy, cx, s] = 0  # 将中心点设为零，去除直流分量############################################################################################################# （3）循环处理不同方向的滤波器：对于每个角度方向，将图像与滤波器进行卷积，并计算响应的幅度和相位。然后，根据幅度和相位计算能量。for o in np.arange(num_angles):angle = o * math.pi / num_angles  # 计算当前方向的滤波角度ds = sin_theta * math.cos(angle) - cos_theta * math.sin(angle)  # 计算sin差异dc = cos_theta * math.cos(angle) + sin_theta * math.sin(angle)  # 计算cos差异dtheta = np.abs(np.arctan2(ds, dc))  # 计算绝对角距离dtheta = np.minimum(dtheta * num_angles / 2, math.pi)  # 归一化角距离到0到pi之间spread = (np.cos(dtheta) + 1) / 2  # 计算角度分布权重，归一化到0到1之间# 初始化累加器矩阵sum_even_response = np.zeros((nrows, ncols))sum_odd_response = np.zeros((nrows, ncols))sum_amplitude_response = np.zeros((nrows, ncols))energy = np.zeros((nrows, ncols))max_amplitude = []# 循环处理不同尺度的滤波器for s in np.arange(num_scales):# 使用角度分布权重乘以径向和角度滤波器得到滤波器current_filter = annular_bandpass_filters[:, :, s] * spread# 对滤波器进行傅里叶变换critical_filter_shift = np.fft.ifftshift(current_filter)critical_filter_shift_cv = np.empty((nrows, ncols, 2))for ip in range(2):critical_filter_shift_cv[:, :, ip] = critical_filter_shift# 使用滤波器对输入图像进行卷积，得到奇偶响应结果matrix_eo = cv2.idft(critical_filter_shift_cv * f_transform)complex_response[:, :, s, o] = matrix_eo[:, :, 1] + 1j * matrix_eo[:, :, 0]# 计算奇偶响应结果的振幅amplitude_response = cv2.magnitude(matrix_eo[:, :, 0], matrix_eo[:, :, 1])# 累加振幅响应结果sum_amplitude_response += amplitude_responsesum_even_response += matrix_eo[:, :, 1]sum_odd_response += matrix_eo[:, :, 0]# 在最小尺度下，通过振幅响应结果估计噪声特性if s == 0:tau = np.median(sum_amplitude_response) / math.sqrt(math.log(4))  # 使用中值法估计噪声统计信息max_amplitude = amplitude_responseelse:max_amplitude = np.maximum(max_amplitude, amplitude_response)  # 记录跨尺度的最大分量的振幅。这是确定频率扩展加权的必要信息。# 累加奇偶响应的能量energy_vector[:, :, 0] += sum_even_responseenergy_vector[:, :, 1] += math.cos(angle) * sum_odd_responseenergy_vector[:, :, 2] += math.sin(angle) * sum_odd_response# 计算能量X_energy = np.sqrt(sum_even_response ** 2 + sum_odd_response ** 2) + epsilonmean_even_response = sum_even_response / X_energymean_odd_response = sum_odd_response / X_energyfor s in np.arange(num_scales):E = complex_response[:, :, s, o].real  # 提取实部O = complex_response[:, :, s, o].imag  # 提取虚部energy += E * mean_even_response + O * mean_odd_response - np.abs(E * mean_odd_response - O * mean_even_response)############################################################################################################# （4）去除噪声：估计噪声水平，并根据噪声阈值去除一部分低能量信号（噪声和弱特征）。total_tau = tau * (1 - (1 / scale_multiplier) ** num_scales) / (1 - (1 / scale_multiplier))est_noise_mean = total_tau * math.sqrt(math.pi / 2)est_noise_sigma = total_tau * math.sqrt((4 - math.pi) / 2)threshold = est_noise_mean + noise_threshold_param * est_noise_sigma  # 根据噪声均值和标准差估计噪声energy = np.maximum(energy - threshold, 0)  # 去除噪声# 根据振幅响应计算能量width = (sum_amplitude_response / (max_amplitude + epsilon) - 1) / (num_scales - 1)weight = 1.0 / (1 + np.exp((sigmoid_cutoff - width) * sigmoid_slope))  # Sigmoidal加权函数phase_congruency_matrix[:, :, o] = weight * energy / sum_amplitude_response  # 归一化相位一致性pc_sum += phase_congruency_matrix[:, :, o]  # 累加到最终相位一致性# 计算协方差cov_x = phase_congruency_matrix[:, :, o] * math.cos(angle)cov_y = phase_congruency_matrix[:, :, o] * math.sin(angle)cov_x2 += cov_x ** 2cov_y2 += cov_y ** 2cov_xy += cov_x * cov_y############################################################################################################# （5）计算相位一致性：根据不同方向的滤波器响应，计算相位一致性的协方差矩阵，然后求得最大矩和最小矩。cov_x2 /= (num_angles / 2)cov_y2 /= (num_angles / 2)cov_xy = 4 * cov_xy / num_anglesdenominator = np.sqrt(cov_xy ** 2 + (cov_x2 - cov_y2) ** 2) + epsilonmax_energy_matrix = (cov_y2 + cov_x2 + denominator) / 2  # 归一化最大值min_energy_matrix = (cov_y2 + cov_x2 - denominator) / 2  # 归一化最小值############################################################################################################# （6）计算图像特征（方向和特征相位）：根据滤波器响应的能量分量，计算主方向和特征相位。orientation_matrix = np.arctan2(energy_vector[:, :, 2], energy_vector[:, :, 1])  # 方向矩阵orientation_matrix[orientation_matrix < 0] += math.pi  # 确保方向值为正orientation_matrix = np.round(orientation_matrix * 180 / math.pi)  # 方向转为度odd_vector_magnitude = np.sqrt(energy_vector[:, :, 1] ** 2 + energy_vector[:, :, 2] ** 2)  # 奇偶向量的模feature_type = np.arctan2(energy_vector[:, :, 0], odd_vector_magnitude)  # 特征类型（也可用于方向性分析）return max_energy_matrix  # 归一化后的相位一致性矩阵if __name__ == '__main__':image = tifffile.imread(r"F:\py\gray_image.tif")image_pc = phase_congruency(image, num_scales=4, num_angles=6)import napariviewer = napari.Viewer()  # 创建napari视图viewer.add_image(image, name="image")viewer.add_image(image_pc, name="phase_congruency")viewer.grid.enabled = True  # 启用网格视图模式napari.run()  # 启动napari事件循环，使得窗口保持打开并可交互。

3.2、skimage.feature（2D / 3D）

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

import numpy as np
import tifffile
import napariimport skimage.feature
import skimage.filtersdef calculate_phase_congruency(input_image, canny_thresholds=(50, 200), sigma=2.0, phase_power=2.0, amplification_factor=2.0):"""计算图像的相位一致性特征，并与 Canny 边缘检测结果合成。:param input_image: 输入图像，支持2D或3D。:param phase_power: 相位幂，用于控制相位一致性的非线性调节。:param amplification_factor: 放大因子，用于相位一致性特征增强。:param sigma: Hessian 矩阵的平滑参数。:param canny_thresholds: Canny 边缘检测的低阈值和高阈值。:return: 相位一致性和 Canny 边缘检测合成图。"""# （1）中值滤波以减少噪声filtered_image = skimage.filters.median(input_image)# （2）Canny 边缘检（2D / 3D）if filtered_image.ndim == 2:edges = skimage.feature.canny(filtered_image,low_threshold=canny_thresholds[0],high_threshold=canny_thresholds[1])elif filtered_image.ndim == 3:edges = np.zeros_like(filtered_image)for slice_index in range(filtered_image.shape[0]):edges[slice_index] = skimage.feature.canny(filtered_image[slice_index],low_threshold=canny_thresholds[0],high_threshold=canny_thresholds[1])else:raise ValueError("Only 2D or 3D images are supported.")# （3）计算 Hessian 矩阵和特征值hessian_matrices = skimage.feature.hessian_matrix(filtered_image, sigma=sigma, mode='reflect')eigenvalues = skimage.feature.hessian_matrix_eigvals(hessian_matrices)# （4）计算相位一致性特征图（2D / 3D）if filtered_image.ndim == 2:lambda1, lambda2 = eigenvaluesphase_congruency_map = np.sqrt(lambda1 ** 2 + lambda2 ** 2) / \(np.abs(lambda1) + np.abs(lambda2) + 1e-8)elif filtered_image.ndim == 3:lambda1, lambda2, lambda3 = eigenvalues * -1phase_congruency_map = np.sqrt(lambda1 ** 2 + lambda2 ** 2 + lambda3 ** 2) / \(np.abs(lambda1) + np.abs(lambda2) + np.abs(lambda3) + 1e-8)# （5）合成 Canny 边缘检测与相位一致性特征combined_phase_congruency = (phase_congruency_map * 0.5 + edges) ** phase_power * amplification_factorimport napariviewer = napari.Viewer()  # 创建napari视图viewer.add_image(image, name="image")viewer.add_image(edges, name="edges")viewer.add_image(phase_congruency_map, name="phase_congruency_map")viewer.add_image(combined_phase_congruency, name="combined_phase_congruency")viewer.grid.enabled = True  # 启用网格视图模式napari.run()  # 启动napari事件循环，使得窗口保持打开并可交互。return combined_phase_congruencyif __name__ == '__main__':image = tifffile.imread(r"F:\py\gray_image.tif")image_pc = calculate_phase_congruency(image)