一、实验目的

1. 理解图像区域的轮廓提取，学会绘制图像中的轮廓。

2. 掌握python中测试算法时间的方法。

3. 分析不同轮廓提取算法的效果差异，掌握轮廓提取的优化方法（去噪、阈值调整、轮廓筛选等）

二、实验环境

• 操作系统：Windows 10
• 开发工具：PyCharm
• 依赖库：Python 3.8、OpenCV-Python 4.8.0、NumPy 1.24.3

三、实验原理

本次实验采用的轮廓提取算法原理如下：

1. Canny 边缘检测：通过高斯滤波去噪→计算梯度→非极大值抑制→双阈值边缘连接，实现连续、精准的边缘提取；
2. 梯度算子（Sobel/Prewitt/Roberts）：通过卷积核计算图像灰度的一阶导数，捕获灰度变化的边缘信息（Sobel 带权重、Prewitt 等权重、Roberts 为 2×2 对角线核）；
3. 拉普拉斯算法：计算图像灰度的二阶导数，检测灰度变化率的突变（对噪声敏感，需结合高斯滤波）；
4. 频域高通滤波：通过傅里叶变换将图像转换到频域，抑制低频（平滑区域）、保留高频（边缘），再逆变换回空间域得到轮廓。

四、实验内容

完整代码：

import cv2
import numpy as np
import time


# ---------------------- 读取并预处理图像 ----------------------
def load_image(path):
    img = cv2.imread(path)
    if img is None:
        raise ValueError("Image path error! Please check.")
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    gray = cv2.GaussianBlur(gray, (7, 7), 0)
    gray = cv2.medianBlur(gray, 5)
    return img, gray


# 替换为你的图像路径
img_original, gray_img = load_image("001.bmp")
h, w = gray_img.shape[:2]


# ---------------------- 通用轮廓提取 ----------------------
def get_contours(edge_img, original_img, color, draw_features=True):
    _, binary = cv2.threshold(edge_img, 40, 255, cv2.THRESH_BINARY)
    kernel = np.ones((5, 5), np.uint8)
    binary = cv2.morphologyEx(binary, cv2.MORPH_CLOSE, kernel, iterations=2)
    binary = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel, iterations=2)

    contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    valid_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > 200]

    result = original_img.copy()
    cv2.drawContours(result, valid_contours, -1, color, 2)

    if draw_features:
        for cnt in valid_contours:
            area = cv2.contourArea(cnt)
            if area > 500:
                M = cv2.moments(cnt)
                if M["m00"] != 0:
                    cX, cY = int(M["m10"] / M["m00"]), int(M["m01"] / M["m00"])
                    cv2.putText(result, f"Area: {int(area)}", (cX - 40, cY - 8),
                                cv2.FONT_HERSHEY_SIMPLEX, 0.3, (0, 0, 0), 1)
    return result


# ---------------------- 各算法实现 ----------------------
def canny_contour(gray, original):
    start = time.time()
    edges = cv2.Canny(gray, 60, 180)
    return get_contours(edges, original, (0, 255, 255)), time.time() - start


def sobel_contour(gray, original):
    start = time.time()
    sobel_x = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=5)
    sobel_y = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=5)
    sobel_abs = cv2.convertScaleAbs(cv2.addWeighted(sobel_x, 0.5, sobel_y, 0.5, 0))
    return get_contours(sobel_abs, original, (0, 255, 0)), time.time() - start


def prewitt_contour(gray, original):
    start = time.time()
    kernel_x = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], np.float32)
    kernel_y = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]], np.float32)
    prewitt_x = cv2.filter2D(gray, cv2.CV_64F, kernel_x)
    prewitt_y = cv2.filter2D(gray, cv2.CV_64F, kernel_y)
    prewitt_abs = cv2.convertScaleAbs(cv2.addWeighted(prewitt_x, 0.5, prewitt_y, 0.5, 0))
    return get_contours(prewitt_abs, original, (255, 0, 0)), time.time() - start


def roberts_contour(gray, original):
    start = time.time()
    kernel1 = np.array([[1, 0], [0, -1]], np.float32)
    kernel2 = np.array([[0, 1], [-1, 0]], np.float32)
    roberts1 = cv2.filter2D(gray, cv2.CV_64F, kernel1)
    roberts2 = cv2.filter2D(gray, cv2.CV_64F, kernel2)
    roberts_abs = cv2.convertScaleAbs(cv2.addWeighted(roberts1, 0.5, roberts2, 0.5, 0))
    return get_contours(roberts_abs, original, (0, 0, 255)), time.time() - start


def laplacian_contour(gray, original):
    start = time.time()
    laplacian = cv2.Laplacian(gray, cv2.CV_64F, ksize=5)
    laplacian_abs = cv2.convertScaleAbs(laplacian)
    return get_contours(laplacian_abs, original, (255, 255, 0)), time.time() - start


def frequency_contour(gray, original):
    start = time.time()
    h, w = gray.shape
    f = np.fft.fft2(gray)
    f_shift = np.fft.fftshift(f)
    d0, n = 25, 2
    x, y = np.meshgrid(np.arange(w), np.arange(h))
    center_x, center_y = w // 2, h // 2
    distance = np.sqrt((x - center_x) ** 2 + (y - center_y) ** 2)
    butterworth_highpass = 1 / (1 + (d0 / (distance + 1e-6)) ** (2 * n))
    f_shift_filtered = f_shift * butterworth_highpass
    f_ishift = np.fft.ifftshift(f_shift_filtered)
    img_back = np.abs(np.fft.ifft2(f_ishift))
    img_back = cv2.normalize(img_back, None, 0, 255, cv2.NORM_MINMAX, cv2.CV_8U)
    return get_contours(img_back, original, (128, 0, 128)), time.time() - start


# ---------------------- 图像拼接 ----------------------
def stitch_images(images, row=2, col=3):
    min_h = min(img.shape[0] for img in images)
    min_w = min(img.shape[1] for img in images)
    resized = [cv2.resize(img, (min_w, min_h)) for img in images]
    rows = [cv2.hconcat(resized[i * col:(i + 1) * col]) for i in range(row)]
    combined = cv2.vconcat(rows)
    cv2.putText(combined, "Contour Extraction Comparison", (10, 30),
                cv2.FONT_HERSHEY_SIMPLEX, 1, (0, 0, 255), 2)
    return combined


# ---------------------- 主程序（优化退出逻辑） ----------------------
if __name__ == "__main__":
    try:
        # 运行算法
        res_canny, t_canny = canny_contour(gray_img, img_original)
        res_sobel, t_sobel = sobel_contour(gray_img, img_original)
        res_prewitt, t_prewitt = prewitt_contour(gray_img, img_original)
        res_roberts, t_roberts = roberts_contour(gray_img, img_original)
        res_laplacian, t_laplacian = laplacian_contour(gray_img, img_original)
        res_freq, t_freq = frequency_contour(gray_img, img_original)

        # 打印耗时
        print("【Algorithm Time Cost (s)】")
        print(f"Canny: {t_canny:.4f} | Sobel: {t_sobel:.4f} | Prewitt: {t_prewitt:.4f}")
        print(f"Roberts: {t_roberts:.4f} | Laplacian: {t_laplacian:.4f} | Frequency: {t_freq:.4f}")

        # 显示结果
        images = [img_original, res_canny, res_sobel, res_prewitt, res_roberts, res_laplacian]
        combined_img = stitch_images(images)
        cv2.imshow("Contour Comparison", combined_img)
        cv2.imshow("6. Frequency Filter", res_freq)

        # 等待用户按键（按任意键退出，更友好）
        cv2.waitKey(0)  # 0表示无限等待，直到有按键输入
        cv2.destroyAllWindows()  # 确保关闭所有窗口

    except Exception as e:
        print(f"Error: {e}")
        cv2.destroyAllWindows()  # 出错时也关闭窗口

（一）、Canny 边缘检测算法（代码实现详解）

核心原理映射

Canny 算法通过 “去噪→梯度计算→非极大值抑制→双阈值连接” 四步提取边缘，代码中通过cv2.Canny()函数直接实现核心逻辑，后续通过轮廓提取函数处理边缘图。

核心代码及解析

def canny_contour(gray, original):
    start = time.time()
    # 1. Canny边缘检测（核心步骤）
    edges = cv2.Canny(gray, threshold1=60, threshold2=180)  # 双阈值设置
    # 2. 从边缘图提取有效轮廓（通用函数处理）
    result = get_contours(edges, original, (0,255,255))  # 黄色轮廓
    return result, time.time()-start

关键代码解析：

1. cv2.Canny(gray, threshold1=60, threshold2=180)：

   • gray：输入灰度图（需提前去噪，代码中已通过GaussianBlur+medianBlur处理）；
   • threshold1（低阈值）：控制弱边缘连接，低于此值的像素被丢弃；
   • threshold2（高阈值）：控制强边缘，高于此值的像素被视为 “确定边缘”；
   • 函数内部自动完成 “高斯去噪→梯度计算→非极大值抑制→双阈值连接”，输出二值化边缘图（边缘为 255，背景为 0）。

2. 后续通过get_contours函数对边缘图进行形态学优化（闭合间隙、去除小噪点）和轮廓筛选（过滤面积 < 200 的小轮廓），最终得到清晰轮廓。

参数调优技巧：

• 若边缘断裂：降低threshold1（如从 60→40），让更多弱边缘被保留并连接；
• 若冗余边缘多：提高threshold2（如从 180→200），过滤更多噪声导致的弱边缘。

（二）、Sobel 算子（代码实现详解）

核心原理映射

Sobel 通过 x、y 方向的加权卷积核计算梯度，代码中用cv2.Sobel()分别计算两个方向梯度，再合并得到边缘强度图。

核心代码及解析

def sobel_contour(gray, original):
    start = time.time()
    # 1. 计算x、y方向梯度
    sobel_x = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=5)  # x方向梯度（垂直边缘）
    sobel_y = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=5)  # y方向梯度（水平边缘）
    # 2. 梯度合并与绝对值转换（避免负梯度丢失）
    sobel_abs = cv2.convertScaleAbs(cv2.addWeighted(sobel_x, 0.5, sobel_y, 0.5, 0))
    # 3. 提取轮廓
    result = get_contours(sobel_abs, original, (0,255,0))  # 绿色轮廓
    return result, time.time()-start

关键代码解析：

1. cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=5)：

   • cv2.CV_64F：输出数据类型（支持负梯度，因边缘可能有明暗变化导致梯度为负）；
   • dx=1, dy=0：计算 x 方向梯度（对垂直边缘敏感，如物体左右边缘）；
   • ksize=5：卷积核大小（3/5/7，越大边缘越平滑，抗噪性越强）。

2. cv2.convertScaleAbs(...)：将负梯度转为绝对值（边缘强度非负），并转换为uint8类型（便于后续处理）。

3. cv2.addWeighted(sobel_x, 0.5, sobel_y, 0.5, 0)：合并 x、y 方向梯度（各占 50% 权重），得到整体边缘强度。

参数调优技巧：

• 若边缘过粗：减小ksize（如从 5→3），保留更多细节；
• 若水平边缘不明显：提高 y 方向权重（如addWeighted(sobel_x, 0.3, sobel_y, 0.7, 0)）。

（三）、Prewitt 算子（代码实现详解）

核心原理映射

Prewitt 与 Sobel 类似，但用等权重卷积核计算梯度，代码中通过自定义核cv2.filter2D()实现卷积。

核心代码及解析

def prewitt_contour(gray, original):
    start = time.time()
    # 1. 定义Prewitt卷积核（x方向检测垂直边缘，y方向检测水平边缘）
    kernel_x = np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], np.float32)
    kernel_y = np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]], np.float32)
    # 2. 卷积计算梯度
    prewitt_x = cv2.filter2D(gray, cv2.CV_64F, kernel_x)  # x方向梯度
    prewitt_y = cv2.filter2D(gray, cv2.CV_64F, kernel_y)  # y方向梯度
    # 3. 梯度合并与绝对值转换
    prewitt_abs = cv2.convertScaleAbs(cv2.addWeighted(prewitt_x, 0.5, prewitt_y, 0.5, 0))
    # 4. 提取轮廓
    result = get_contours(prewitt_abs, original, (255,0,0))  # 蓝色轮廓
    return result, time.time()-start

关键代码解析：

1. 自定义卷积核：

   • kernel_x：x 方向核（左右像素差，检测垂直边缘），所有非中心像素权重均为 ±1（等权重，区别于 Sobel 的加权）；
   • kernel_y：y 方向核（上下像素差，检测水平边缘），权重同样为 ±1。

2. cv2.filter2D(gray, cv2.CV_64F, kernel_x)：用自定义核对图像进行卷积，得到 x 方向梯度（原理同 Sobel，但核权重不同）。

与 Sobel 的代码差异：

• Sobel 用cv2.Sobel()内置函数，Prewitt 需自定义核并用cv2.filter2D()实现；
• Sobel 核权重随距离中心变化（如中心邻域权重为 ±2），Prewitt 核权重均等（±1），因此 Prewitt 抗噪性稍弱（代码中需加强去噪）。

（四）、Roberts 算子（代码实现详解）

核心原理映射

Roberts 用 2×2 核计算对角线方向梯度，代码中通过自定义 2×2 核卷积实现，对陡峭边缘敏感。

核心代码及解析

def roberts_contour(gray, original):
    start = time.time()
    # 1. 定义Roberts卷积核（检测对角线边缘）
    kernel1 = np.array([[1, 0], [0, -1]], np.float32)  # 45°方向
    kernel2 = np.array([[0, 1], [-1, 0]], np.float32)  # 135°方向
    # 2. 卷积计算梯度
    roberts1 = cv2.filter2D(gray, cv2.CV_64F, kernel1)
    roberts2 = cv2.filter2D(gray, cv2.CV_64F, kernel2)
    # 3. 梯度合并与绝对值转换
    roberts_abs = cv2.convertScaleAbs(cv2.addWeighted(roberts1, 0.5, roberts2, 0.5, 0))
    # 4. 提取轮廓
    result = get_contours(roberts_abs, original, (0,0,255))  # 红色轮廓
    return result, time.time()-start

关键代码解析：

1. 2×2 卷积核：

   • kernel1：检测 45° 方向边缘（左上角 - 右下角像素差）；
   • kernel2：检测 135° 方向边缘（右上角 - 左下角像素差）；核尺寸小（2×2），计算速度快，但抗噪性差（代码中需用medianBlur加强去噪）。

2. 由于核尺寸小，梯度对细节敏感，代码中get_contours需提高min_area（如 200）过滤噪声导致的微小轮廓。

（五）、拉普拉斯算法（代码实现详解）

核心原理映射

拉普拉斯通过二阶导数检测边缘（零交叉点），代码中用cv2.Laplacian()计算二阶导数，因对噪声敏感，需配合强去噪。

核心代码及解析

def laplacian_contour(gray, original):
    start = time.time()
    # 1. 计算拉普拉斯二阶导数（已提前通过高斯滤波去噪）
    laplacian = cv2.Laplacian(gray, cv2.CV_64F, ksize=5)  # ksize=5增强平滑
    # 2. 绝对值转换（二阶导数有正负）
    laplacian_abs = cv2.convertScaleAbs(laplacian)
    # 3. 提取轮廓
    result = get_contours(laplacian_abs, original, (255,255,0))  # 青蓝色轮廓
    return result, time.time()-start

关键代码解析：

1. cv2.Laplacian(gray, cv2.CV_64F, ksize=5)：

   • ksize=5：拉普拉斯核大小（必须为正奇数），越大平滑效果越强（抑制噪声）；
   • 函数计算二阶导数（∇²f = ∂²f/∂x² + ∂²f/∂y²），边缘处会出现正负交替的 “零交叉点”。

2. 代码中提前通过GaussianBlur(gray, (7,7), 0)加强去噪，否则laplacian_abs会包含大量噪声边缘。

参数调优技巧：

• 若噪声边缘多：增大ksize（如从 5→7）或增强高斯滤波核（如(9,9)）；
• 若边缘过弱：降低get_contours中的二值化阈值（如从 40→30）。

（六）、频域高通滤波法（代码实现详解）

核心原理映射

通过傅里叶变换将图像转为频域，用高通滤波器保留高频（边缘），再逆变换回空间域，代码中需手动实现傅里叶变换、滤波、逆变换三步。

核心代码及解析

def frequency_contour(gray, original):
    start = time.time()
    h, w = gray.shape
    # 1. 傅里叶变换+中心化（低频移到中心）
    f = np.fft.fft2(gray)  # 二维傅里叶变换
    f_shift = np.fft.fftshift(f)  # 低频中心化
    # 2. 构建巴特沃斯高通滤波器
    d0 = 25  # 截止频率（越小保留高频越多）
    n = 2    # 阶数（越大边缘越锐利）
    x, y = np.meshgrid(np.arange(w), np.arange(h))  # 生成坐标网格
    center_x, center_y = w//2, h//2  # 频域中心
    distance = np.sqrt((x-center_x)**2 + (y-center_y)** 2)  # 各点到中心的距离
    butterworth_highpass = 1 / (1 + (d0/(distance+1e-6))**(2*n))  # 高通滤波器
    # 3. 频域滤波（相乘）
    f_shift_filtered = f_shift * butterworth_highpass
    # 4. 逆傅里叶变换（转回空间域）
    f_ishift = np.fft.ifftshift(f_shift_filtered)  # 逆中心化
    img_back = np.abs(np.fft.ifft2(f_ishift))  # 逆变换+取幅度
    img_back = cv2.normalize(img_back, None, 0, 255, cv2.NORM_MINMAX, cv2.CV_8U)  # 归一化
    # 5. 提取轮廓
    result = get_contours(img_back, original, (128,0,128))  # 紫色轮廓
    return result, time.time()-start

关键代码解析：

1. 傅里叶变换步骤：

   • np.fft.fft2(gray)：将空间域图像转为频域（复数矩阵，包含幅度和相位）；
   • np.fft.fftshift(f)：将低频分量移到中心（原始频域低频在四角），便于设计滤波器。

2. 巴特沃斯高通滤波器：

   • distance：频域中每个点到中心的距离（低频在中心，距离近；高频在边缘，距离远）；
   • 滤波器公式确保 “距离越远（高频），保留越强”，d0越小，更多高频被保留（边缘更密集）。

3. 逆变换与归一化：

   • np.fft.ifft2(f_ishift)：将滤波后的频域图像转回空间域；
   • cv2.normalize(...)：将频域转换后的灰度值缩放到 0-255（便于后续二值化）。

参数调优技巧：

• 若边缘模糊：减小d0（如从 25→15），保留更多高频；
• 若噪声多：增大d0（如从 25→35），抑制高频噪声；
• 若边缘过锐：减小n（如从 2→1），让滤波器过渡更平滑。

五、实验总结

算法特性：Canny 算法综合效果最优（边缘连续、效率高），频域高通滤波抗噪性强但耗时久，梯度算子适合不同复杂度场景；
代码实现：OpenCV 内置函数（如cv2.Canny()、cv2.Sobel()）简化了空间域算法实现，频域算法需手动处理傅里叶变换，复杂度较高；
调优关键：预处理（去噪）、阈值调整、轮廓筛选是提升效果的核心，需根据图像特点（噪声、边缘强度）动态调整参数。