RGB图像类梯度轮廓检测的Python高效实现优化请求

需求概述 #

将RGB图像视为向量场$f$的离散化（定义域为像素二维网格，值域为RGB三维向量），需完成以下计算：

1. 计算标量场$g$：每个像素$x$的$g(x)$是其邻域像素$y$与$x$的RGB向量差的最大欧几里得范数。
2. 计算二值场$h$：当$g(x) \geq t$时$h(x,t)=1$，否则为0；通过调整阈值$t$，$h$可逼近图像轮廓。

现有实现依赖双重循环逐像素处理，速度远慢于库函数卷积操作，且结果不够理想，需基于Numpy向量化技术优化。

优化思路 #

1. 用数组广播替代双重循环：利用Numpy的广播机制，一次性处理所有像素的邻域计算，避免逐像素遍历的开销。
2. 统一边界处理：通过np.pad给图像添加1层边界（复制边缘像素），避免逐像素判断邻域范围，简化计算逻辑。
3. 向量化计算欧几里得范数：直接对偏移后的图像与原图像做差，计算平方和后开根号，无需循环计算每个向量的点积。
4. 替换np.vectorize：用布尔索引直接生成二值图，比np.vectorize效率更高。

优化后代码 #

import numpy as np
import matplotlib.pyplot as plt
import os

directories = ['images']
max_norm = np.sqrt(3 * 255**2)  # 预计算最大可能范数

# 定义3x3邻域的偏移量（排除中心像素，因为差为0不影响最大值）
offsets = [(-1, -1), (-1, 0), (-1, 1),
           (0, -1),          (0, 1),
           (1, -1),  (1, 0), (1, 1)]

for directory in directories:
    for filename in os.scandir(directory):
        if filename.is_file():
            print(filename.path)
            image = plt.imread(filename.path).astype(np.float64)
            plt.imshow(image / 255)  # 显示原图像（plt.imread读入的是0-255的uint8）
            plt.show()
            
            height, width = image.shape[:2]
            # 给图像加1层边界，处理边缘像素的邻域
            padded_image = np.pad(image, pad_width=((1,1), (1,1), (0,0)), mode='edge')
            
            # 初始化最大范数数组
            max_diff_norm = np.zeros((height, width), dtype=np.float64)
            
            # 遍历所有邻域偏移，计算每个偏移下的差的范数，更新最大值
            for dr, dc in offsets:
                # 取偏移后的图像区域，与原图像对齐
                shifted = padded_image[1+dr : 1+dr+height, 1+dc : 1+dc+width]
                # 计算RGB向量差的平方和，开根号得到范数
                diff_norm = np.sqrt(np.sum((shifted - image)**2, axis=2))
                # 更新当前像素的最大范数
                max_diff_norm = np.maximum(max_diff_norm, diff_norm)
            
            # 归一化到0-1范围
            gradient_image = max_diff_norm / max_norm
            
            # 计算阈值并生成二值图（用布尔索引替代vectorize）
            mu = gradient_image.mean()
            gradient_image = (gradient_image >= mu).astype(np.uint8)
            
            plt.imshow(gradient_image, cmap='gray')
            plt.show()

关键优化点说明 #

• 邻域偏移处理：通过对图像做不同偏移，一次性计算所有像素对应邻域位置的向量差，避免逐像素提取邻域。
• 边界填充：使用mode='edge'的padding，保证边缘像素的邻域是合理的边缘扩展，无需额外判断边界范围。
• 向量化范数计算：np.sum((shifted - image)**2, axis=2)直接对所有像素的RGB分量计算平方和，再开根号，效率远高于循环计算。
• 布尔索引生成二值图：(gradient_image >= mu).astype(np.uint8)直接将满足条件的像素设为1，否则为0，比np.vectorize快数倍。

内容的提问来源于stack exchange，提问作者Alberto