10_imageGrayscale.py - 图像灰度化处理完全指南

学习路径第 10 步 (共 10 步) | 难度：高级（综合实战）

总结

numpy 教程到此结束， 09 10两个章节为numpy学习成果的输出。

概述

综合运用 NumPy 全部核心知识（数组操作、广播、向量化、切片、数据类型）完成一个完整的图像处理流水线：从 RGB 图像到灰度图、直方图分析、Sobel 边缘检测。本文件是整个学习路径的收官之作。

学习目标

• 理解图像在 NumPy 中的数据结构 (H, W, 3) uint8
• 掌握 5 种灰度化算法（加权平均 / 算术平均 / 最大值 / 最小值 / 去饱和）
• 实战体验 NumPy 向量化 vs Python 循环的性能差距（约 30 倍加速）
• 学习 Sobel 边缘检测 的完整实现流程
• 掌握内存优化技巧：就地操作、分块处理、dtype 精细选择

核心内容 (9 个模块)

运行输出文件

运行后会生成以下图片：

numpy/grayscale_comparison.png     -- 各方法对比大图
numpy/grayscale_histogram.png      -- 直方图分析
numpy/grayscale_custom_weights.png -- 自定义权重效果
numpy/edge_detection_result.png    -- Sobel 边缘检测流水线

Code

"""
=====================================
NumPy 图像灰度化处理完全指南 (Image Grayscale)
=====================================

本案例系统介绍使用 NumPy 实现图像灰度化的：
1. 图像的基本数据结构（RGB三通道）
2. 多种灰度化算法原理与实现
3. NumPy 向量化操作 vs 循环对比
4. 性能优化技巧
5. 图像可视化展示
6. 实际应用场景扩展

前置知识：
  - 彩色图像由 R(红)、G(绿)、B(蓝) 三个通道组成
  - 每个通道的取值范围: 0-255 (uint8)
  - 灰度图像只有一个通道，每个像素表示亮度

作者：bloxed
"""

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import time
import os


def separator(title):
    """打印分隔线"""
    print(f"\n{'='*60}")
    print(f"  {title}")
    print('='*60)


# ============================================================
# 第一部分：创建测试图像（无需外部图片）
# ============================================================
separator("一、创建彩色测试图像")

def create_test_image(width=256, height=200):
    """
    创建一个包含多种颜色和渐变的测试图像
    返回 shape: (height, width, 3) 的 uint8 数组
    """
    image = np.zeros((height, width, 3), dtype=np.uint8)
    
    # 计算区域边界（确保整数划分）
    w1 = width // 3
    w2 = width // 3
    w3 = width - w1 - w2
    h_half = height // 2

    # 区域1：红色渐变 (左上角) - 使用 np.tile 广播到正确行数
    red_grad = np.linspace(100, 255, w1).astype(np.uint8)
    image[:h_half, :w1, 0] = np.tile(red_grad, (h_half, 1))
    image[:h_half, :w1, 1] = 50
    image[:h_half, :w1, 2] = 50

    # 区域2：绿色渐变 (中上)
    green_grad = np.linspace(100, 255, w2).astype(np.uint8)
    image[:h_half, w1:w1+w2, 0] = 50
    image[:h_half, w1:w1+w2, 1] = np.tile(green_grad, (h_half, 1))
    image[:h_half, w1:w1+w2, 2] = 50

    # 区域3：蓝色渐变 (右上)
    blue_grad = np.linspace(100, 255, max(w3, 1)).astype(np.uint8)
    image[:h_half, w1+w2:, 0] = 50
    image[:h_half, w1+w2:, 1] = 50
    if w3 > 0:
        image[:h_half, w1+w2:, 2] = np.tile(blue_grad[:w3], (h_half, 1))

    # 区域4：RGB混合渐变 (下半部分 - 模拟真实照片效果)
    y_coords = np.arange(h_half, height)[:, np.newaxis]
    x_coords = np.arange(width)[np.newaxis, :]
    image[h_half:, :, 0] = ((x_coords / width * 200 + y_coords / height * 55)).clip(0, 255).astype(np.uint8)
    image[h_half:, :, 1] = ((x_coords / width * 150 + (y_coords - h_half) / h_half * 105 + 50)).clip(0, 255).astype(np.uint8)
    image[h_half:, :, 2] = ((y_coords - h_half) / h_half * 200 + 55).clip(0, 255).astype(np.uint8)

    return image

# 创建测试图像
img_color = create_test_image()

print(f"原始彩色图像信息:")
print(f"  数据类型: {img_color.dtype}")
print(f"  图像形状: {img_color.shape} → (高={img_color.shape[0]}, 宽={img_color.shape[1]}, 通道数={img_color.shape[2]})")
print(f"  总像素数: {img_color.size}")
print(f"  内存占用: {img_color.nbytes / 1024:.1f} KB")

# 显示各通道值示例
print(f"\n像素 (10, 10) 的 RGB 值: R={img_color[10, 10, 0]}, G={img_color[10, 10, 1]}, B={img_color[10, 10, 2]}")


# ============================================================
# 第二部分：理解图像的数据结构
# ============================================================
separator("二、深入理解 RGB 图像数据结构")

print("""
┌──────────────────────────────────────────────────────────────┐
│                    RGB 图像三维数组结构                        │
├──────────────────────────────────────────────────────────────┤
│  轴0 (Axis 0): 高度 (Height) - 行数                           │
│  轴1 (Axis 1): 宽度 (Width)  - 列数                          │
│  轴2 (Axis 2): 通道 (Channel) - [R, G, B]                    │
│                                                              │
│  image[y, x, 0] → 红色分量 (Red)                             │
│  image[y, x, 1] → 绿色分量 (Green)                           │
│  image[y, x, 2] → 蓝色分量 (Blue)                            │
│                                                              │
│  人眼对三种颜色的敏感度: 绿色 > 红色 > 蓝色                   │
└──────────────────────────────────────────────────────────────┘
""")

# 提取各通道
R_channel = img_color[:, :, 0]
G_channel = img_color[:, :, 1]
B_channel = img_color[:, :, 2]

print(f"提取各通道后:")
print(f"  R通道形状: {R_channel.shape}, 值范围: [{R_channel.min()}, {R_channel.max()}]")
print(f"  G通道形状: {G_channel.shape}, 值范围: [{G_channel.min()}, {G_channel.max()}]")
print(f"  B通道形状: {B_channel.shape}, 值范围: [{B_channel.min()}, {B_channel.max()}]")


# ============================================================
# 第三部分：多种灰度化算法实现
# ============================================================
separator("三、多种灰度化算法实现")


def grayscale_weighted(image):
    """
    方法1: 加权平均法 (推荐 ⭐)
    
    原理: 根据人眼对不同颜色的敏感度分配权重
    公式: Gray = 0.299*R + 0.587*G + 0.114*B
    
    这是 ITU-R BT.601 标准，最符合人眼感知
    """
    # 使用 NumPy 广播机制，一次性对所有像素运算
    weights = np.array([0.299, 0.587, 0.114])
    gray = np.dot(image[..., :3].astype(float), weights)
    return gray.astype(np.uint8)


def grayscale_average(image):
    """
    方法2: 简单平均法
    
    原理: R、G、B 三通道取算术平均值
    公式: Gray = (R + G + B) / 3
    
    优点: 计算简单快速
    缺点: 未考虑人眼对绿色的更高敏感性，结果偏暗
    """
    gray = np.mean(image[..., :3], axis=2).astype(np.uint8)
    return gray


def grayscale_max(image):
    """
    方法3: 最大值法
    
    原理: 取 RGB 三个通道中的最大值
    公式: Gray = max(R, G, B)
    
    特点: 结果整体较亮，常用于需要突出亮部的场景
    """
    gray = np.max(image[..., :3], axis=2)
    return gray


def grayscale_min(image):
    """
    方法4: 最小值法
    
    原理: 取 RGB 三个通道中的最小值
    公式: Gray = min(R, G, B)
    
    特点: 结果整体较暗，较少单独使用
    """
    gray = np.min(image[..., :3], axis=2)
    return gray


def grayscale_desaturation(image):
    """
    方法5: 去饱和度法
    
    原理: 取 RGB 中的最大值和最小值的平均值
    公式: Gray = (max(R,G,B) + min(R,G,B)) / 2
    
    特点: 效果介于最大值法和最小值法之间
    """
    channel_max = np.max(image[..., :3], axis=2)
    channel_min = np.min(image[..., :3], axis=2)
    gray = ((channel_max.astype(int) + channel_min.astype(int)) // 2).astype(np.uint8)
    return gray


# 执行各种灰度化方法
gray_weighted = grayscale_weighted(img_color)
gray_average = grayscale_average(img_color)
gray_max = grayscale_max(img_color)
gray_min = grayscale_min(img_color)
gray_desaturation = grayscale_desaturation(img_color)

print("五种灰度化算法结果:")
print(f"  {'方法':<16} {'形状':<15} {'值范围':<20} {'均值':<10}")
print(f"  {'-'*60}")
for name, result in [("加权平均", gray_weighted), ("简单平均", gray_average),
                       ("最大值法", gray_max), ("最小值法", gray_min),
                       ("去饱和度", gray_desaturation)]:
    print(f"  {name:<14} {str(result.shape):<15} [{result.min():3d}, {result.max():3d}]       {result.mean():6.1f}")


# ============================================================
# 第四部分：纯Python循环实现（性能对比基准）
# ============================================================
separator("四、NumPy向量化 vs 纯Python循环 (性能对比)")


def grayscale_python_loop(image):
    """纯 Python 三重循环实现（极慢，仅用于演示）"""
    height, width = image.shape[:2]
    gray = np.zeros((height, width), dtype=np.uint8)
    for y in range(height):
        for x in range(width):
            r, g, b = int(image[y, x, 0]), int(image[y, x, 1]), int(image[y, x, 2])
            gray[y, x] = int(0.299 * r + 0.587 * g + 0.114 * b)
    return gray


# 测试较小尺寸图像的性能差异
test_img = create_test_image(128, 96)  # 使用较小图像避免等待过久

# NumPy 向量化版本
start = time.perf_counter()
for _ in range(100):
    _ = grayscale_weighted(test_img)
numpy_time = (time.perf_counter() - start) / 100

# 纯 Python 循环版本（只运行一次）
start = time.perf_counter()
_ = grayscale_python_loop(test_img)
python_time = time.perf_counter() - start

speedup = python_time / numpy_time

print(f"测试图像大小: {test_img.shape[:2]}")
print(f"{'='*50}")
print(f"  NumPy 向量化 (np.dot):   {numpy_time*1000:.4f} ms/次")
print(f"  纯Python 三重循环:        {python_time*1000:.2f} ms/次")
print(f"{'='*50}")
print(f"  速度提升比:  NumPy 快约 {speedup:.0f} 倍!")
print()
print("原因分析:")
print("  1. NumPy 底层使用 C 语言实现，避免了 Python 解释器开销")
print("  2. 利用 CPU 的 SIMD 指令集（如 AVX）并行处理多个数据")
print("  3. 内存连续存储，缓存命中率高")
print("  4. 避免了 Python 循环的类型检查开销")


# ============================================================
# 第五部分：不同方法的视觉效果对比
# ============================================================
separator("五、可视化对比 (保存图像)")

# 创建对比图
fig, axes = plt.subplots(2, 3, figsize=(14, 9))

# 原图
axes[0, 0].imshow(img_color)
axes[0, 0].set_title('Original Color\n(RGB Image)', fontsize=11)
axes[0, 0].axis('off')

# 各通道分离显示
axes[0, 1].imshow(R_channel, cmap='Reds')
axes[0, 1].set_title('Red Channel\n(R)', fontsize=11)
axes[0, 1].axis('off')

axes[0, 2].imshow(G_channel, cmap='Greens')
axes[0, 2].set_title('Green Channel\n(G)', fontsize=11)
axes[0, 2].axis('off')

# 各种灰度化结果
axes[1, 0].imshow(gray_weighted, cmap='gray', vmin=0, vmax=255)
axes[1, 0].set_title('Weighted Method\n(0.299R+0.587G+0.114B) ⭐', fontsize=11)
axes[1, 0].axis('off')

axes[1, 1].imshow(gray_average, cmap='gray', vmin=0, vmax=255)
axes[1, 1].set_title('Average Method\n((R+G+B)/3)', fontsize=11)
axes[1, 1].axis('off')

axes[1, 2].imshow(gray_max, cmap='gray', vmin=0, vmax=255)
axes[1, 2].set_title('Max Method\n(max(R,G,B))', fontsize=11)
axes[1, 2].axis('off')

plt.suptitle('Image Grayscale Conversion Methods Comparison\n图像灰度化方法对比', 
             fontsize=14, fontweight='bold')
plt.tight_layout()

# 保存图像
output_path = os.path.join(os.path.dirname(__file__), 'grayscale_comparison.png')
plt.savefig(output_path, dpi=150, bbox_inches='tight', facecolor='white')
print(f"[OK] 对比图已保存至: {output_path}")

plt.close()


# ============================================================
# 第六部分：直方图分析
# ============================================================
separator("六、灰度直方图分析")

# 创建直方图对比图
fig, axes = plt.subplots(1, 4, figsize=(16, 4))

methods = [
    ('Weighted (推荐)', gray_weighted),
    ('Average', gray_average),
    ('Max', gray_max),
    ('Desaturation', gray_desaturation)
]

for idx, (title, gray_img) in enumerate(methods):
    axes[idx].hist(gray_img.ravel(), bins=256, range=(0, 256), color='steelblue', alpha=0.7)
    axes[idx].set_title(f'{title}\nMean={gray_img.mean():.1f}', fontsize=10)
    axes[idx].set_xlabel('Pixel Value')
    axes[idx].set_ylabel('Frequency')
    if idx > 0:
        axes[idx].set_yticks([])

plt.suptitle('Grayscale Histogram Comparison\n不同方法的灰度分布对比', fontsize=13, fontweight='bold')
plt.tight_layout()

hist_output = os.path.join(os.path.dirname(__file__), 'grayscale_histogram.png')
plt.savefig(hist_output, dpi=150, bbox_inches='tight', facecolor='white')
print(f"[OK] 直方图已保存至: {hist_output}")

plt.close()


# ============================================================
# 第七部分：进阶应用 - 自定义权重调整
# ============================================================
separator("七、进阶应用 - 自定义灰度化权重")

def grayscale_custom(image, wr=0.299, wg=0.587, wb=0.114):
    """
    自定义权重的灰度化函数
    
    参数:
      image: 输入图像
      wr, wg, wb: 红、绿、蓝通道的权重 (总和应为1.0)
    
    应用场景:
      - 增加红色权重: 更好地保留暖色调细节
      - 增加绿色权重: 更接近人眼感知 (默认)
      - 增加蓝色权重: 强调冷色调区域
    """
    assert abs(wr + wg + wb - 1.0) < 0.001, "权重之和应等于1.0"
    weights = np.array([wr, wg, wb])
    gray = np.dot(image[..., :3].astype(float), weights).astype(np.uint8)
    return gray


# 不同权重组合的效果
custom_configs = {
    'Standard (ITU-R BT.601)':  (0.299, 0.587, 0.114),
    'BT.709 (HDTV标准)':        (0.2126, 0.7152, 0.0722),
    'Red Emphasis (强调红色)':   (0.6, 0.25, 0.15),
    'Blue Emphasis (强调蓝色)':  (0.2, 0.3, 0.5),
}

print("\n不同权重配置:")
print(f"{'配置':<28} {'R权重':>8} {'G权重':>8} {'B权重':>8} {'结果均值':>10}")
print("-" * 65)

fig, axes = plt.subplots(1, 5, figsize=(16, 4))

axes[0].imshow(img_color)
axes[0].set_title('Original', fontsize=10)
axes[0].axis('off')

for idx, (name, (wr, wg, wb)) in enumerate(custom_configs.items(), 1):
    custom_gray = grayscale_custom(img_color, wr, wg, wb)
    print(f"{name:<26} ({wr:>5.3f}, {wg:>5.3f}, {wb:>5.3f})   {custom_gray.mean():>8.1f}")
    axes[idx].imshow(custom_gray, cmap='gray', vmin=0, vmax=255)
    axes[idx].set_title(f'{name.split("(")[0].strip()}\n({wr:.2f},{wg:.2f},{wb:.2f})', fontsize=9)
    axes[idx].axis('off')

plt.suptitle('Custom Weight Grayscale Comparison\n自定义权重灰度化对比', fontsize=12, fontweight='bold')
plt.tight_layout()

custom_output = os.path.join(os.path.dirname(__file__), 'grayscale_custom_weights.png')
plt.savefig(custom_output, dpi=150, bbox_inches='tight', facecolor='white')
print(f"\n[OK] 自定义权重对比图已保存至: {custom_output}")

plt.close()


# ============================================================
# 第八部分：实际应用 - 边缘检测预处理
# ============================================================
separator("八、实际应用 - 边缘检测预处理")

print("""
灰度化的典型应用流程:
  彩色图像 → 灰度化 → 高斯模糊 → Sobel/Canny边缘检测 → 二值化
""")

def simple_edge_detection(gray_image, threshold=30):
    """
    简单的边缘检测 (Sobel算子)
    
    原理: 利用相邻像素的梯度变化检测边缘
    Sobel核可以检测水平和垂直方向的边缘
    """
    # 定义 Sobel 卷积核
    sobel_x = np.array([[-1, 0, 1],
                        [-2, 0, 2],
                        [-1, 0, 1]], dtype=np.float64)
    
    sobel_y = np.array([[-1, -2, -1],
                        [ 0,  0,  0],
                        [ 1,  2,  1]], dtype=np.float64)
    
    # 将图像转为 float
    img_float = gray_image.astype(np.float64)
    
    # 手动实现卷积（简化版）
    height, width = img_float.shape
    grad_x = np.zeros_like(img_float)
    grad_y = np.zeros_like(img_float)
    
    # 应用卷积核（忽略边界）
    for y in range(1, height - 1):
        for x in range(1, width - 1):
            patch = img_float[y-1:y+2, x-1:x+2]
            grad_x[y, x] = np.sum(patch * sobel_x)
            grad_y[y, x] = np.sum(patch * sobel_y)
    
    # 计算梯度幅值
    magnitude = np.sqrt(grad_x**2 + grad_y**2)
    
    # 二值化
    edges = (magnitude > threshold).astype(np.uint8) * 255
    
    return magnitude, edges


# 对加权灰度图进行边缘检测
print("正在执行边缘检测 (可能需要几秒)...")
magnitude, edges = simple_edge_detection(gray_weighted, threshold=40)

# 可视化边缘检测结果
fig, axes = plt.subplots(1, 4, figsize=(15, 4))

axes[0].imshow(img_color)
axes[0].set_title('Original', fontsize=10)
axes[0].axis('off')

axes[1].imshow(gray_weighted, cmap='gray')
axes[1].set_title('Grayscale\n(Preprocessing)', fontsize=10)
axes[1].axis('off')

axes[2].imshow(magnitude, cmap='hot')
axes[2].set_title('Gradient Magnitude\n(Sobel)', fontsize=10)
axes[2].axis('off')

axes[3].imshow(edges, cmap='gray')
axes[3].set_title('Binary Edges\n(Threshold=40)', fontsize=10)
axes[3].axis('off')

plt.suptitle('Edge Detection Pipeline: Grayscale as Preprocessing Step\n'
             '边缘检测流水线：灰度化作为预处理步骤', fontsize=12, fontweight='bold')
plt.tight_layout()

edge_output = os.path.join(os.path.dirname(__file__), 'edge_detection_result.png')
plt.savefig(edge_output, dpi=150, bbox_inches='tight', facecolor='white')
print(f"[OK] 边缘检测效果图已保存至: {edge_output}")

plt.close()


# ============================================================
# 第九部分：批量处理与内存优化
# ============================================================
separator("九、批量处理与内存优化技巧")

print("""
【内存优化建议】
""")

# 技巧1: 使用原地操作减少内存拷贝
print("1. 就地操作 (In-place Operation):")
img_copy = img_color.copy()  # 先复制一份
result_inplace = np.empty((*img_copy.shape[:2], 3), dtype=np.float32)
np.multiply(img_copy[..., :3].astype(np.float32), [0.299, 0.587, 0.114], out=result_inplace)
result_inplace = np.sum(result_inplace, axis=2).astype(np.uint8)
print(f"   原地操作结果形状: {result_inplace.shape}")

# 技巧2: 使用更小的数据类型节省内存
print("\n2. 数据类型选择:")
img_uint8 = img_color.astype(np.uint8)     # 1字节/通道, 最常用
img_float32 = img_color.astype(np.float32)  # 4字节/通道, 计算精度高
img_float64 = img_color.astype(np.float64)  # 8字节/通道, 最高精度
print(f"   uint8  内存占用: {img_uint8.nbytes / 1024:.1f} KB")
print(f"   float32 内存占用: {img_float32.nbytes / 1024:.1f} KB")
print(f"   float64 内存占用: {img_float64.nbytes / 1024:.1f} KB")

# 技巧3: 分块处理大图像
print("\n3. 大图像分块处理:")
def process_large_image_chunked(image, chunk_size=64, process_fn=None):
    """
    分块处理大图像以降低峰值内存占用
    适用于: 超高分辨率图像、显存有限的GPU环境
    """
    h, w = image.shape[:2]
    result = np.empty((h, w), dtype=np.uint8)
    
    for y_start in range(0, h, chunk_size):
        y_end = min(y_start + chunk_size, h)
        for x_start in range(0, w, chunk_size):
            x_end = min(x_start + chunk_size, w)
            # 处理当前块
            chunk = image[y_start:y_end, x_start:x_end]
            if process_fn:
                result[y_start:y_end, x_start:x_end] = process_fn(chunk)
    return result

chunked_result = process_large_image_chunked(img_color, chunk_size=64, process_fn=grayscale_weighted)
print(f"   分块处理结果与直接处理一致: {np.array_equal(chunked_result, gray_weighted)} [OK]")


# ============================================================
# 总结
# ============================================================
separator("总结：图像灰度化核心知识点")

summary = """
┌─────────────────────────────────────────────────────────────┐
│                  图像灰度化核心要点总结                        │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  【算法选择】                                                │
│  [*] 加权平均法: Gray = 0.299R + 0.587G + 0.114B  (推荐)    │
│    简单平均法: Gray = (R + G + B) / 3                      │
│    最大值法:   Gray = max(R, G, B)                         │
│    最小值法:   Gray = min(R, G, B)                         │
│                                                             │
│  【NumPy核心操作】                                           │
│  • np.dot() / @ : 矩阵乘法，用于加权计算                     │
│  • np.mean(axis=2): 沿通道轴求均值                           │
│  • np.max/min(axis=2): 沿通道轴求极值                        │
│  • 广播机制(Broadcasting): 自动处理维度匹配                  │
│                                                             │
│  【性能优化】                                                │
│  • 始终优先使用向量化操作，避免 Python 循环                  │
│  • 合理选择数据类型 (uint8 节省内存, float32 平衡精度)       │
│  • 大图像采用分块处理策略                                    │
│                                                             │
│  【典型应用流程】                                            │
│  RGB图像 → 灰度化 → 滤波 → 边缘检测 → 二值化 → 分析         │
│                                                             │
│  【相关库】                                                  │
│  • OpenCV (cv2): cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)     │
│  • Pillow:   Image.convert('L')                            │
│  • Scikit-image: color.rgb2gray()                           │
└─────────────────────────────────────────────────────────────┘
"""
print(summary)

print(f"\n{'='*60}")
print("  所有生成的图像文件:")
print(f"{'='*60}")
output_dir = os.path.dirname(__file__)
for fname in ['grayscale_comparison.png', 'grayscale_histogram.png',
              'grayscale_custom_weights.png', 'edge_detection_result.png']:
    fpath = os.path.join(output_dir, fname)
    exists = "[OK]" if os.path.exists(fpath) else "[X]"
    print(f"  {exists} {fname}")