边缘检测的100种方法
文章目录
- 什么是边缘检测 ?
- 一、边缘检测算子:Sobel算子、Scharr算子、Laplacian算子、Canny算子
- 二、梯度计算 + 顶帽 + 黑帽 + 拉普拉斯金字塔
- 三、相位一致性(Phase Congruency,PC)
- 3.1、底层代码(2D)
- 3.2、skimage.feature(2D / 3D)
什么是边缘检测 ?
边缘检测:边缘检测是图像处理中的一种技术,旨在识别和定位图像中亮度变化显著的区域。这些变化通常表现为强度、颜色或纹理的急剧变动,通常对应于物体的边界或轮廓。通过检测这些边缘信息,边缘检测能够帮助提取图像的结构特征,促进图像的分析和处理。
边缘检测的目的
(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)