知远漫谈头像
关注
Python NumPy - 实战 图像的读取与像素值处理封面图

Python NumPy - 实战 图像的读取与像素值处理

在这里插入图片描述

👋 大家好,欢迎来到我的技术博客!
📚 在这里,我会分享学习笔记、实战经验与技术思考,力求用简单的方式讲清楚复杂的问题。
🎯 本文将围绕NumPy这个话题展开,希望能为你带来一些启发或实用的参考。
🌱 无论你是刚入门的新手,还是正在进阶的开发者,希望你都能有所收获!


Python NumPy实战:图像的读取与像素值处理 📷

在数字图像处理的世界中,Python和NumPy的组合为我们提供了强大的工具集来操作和分析图像数据。从简单的像素级操作到复杂的图像变换,掌握这些技能对于计算机视觉、机器学习和数据分析领域的从业者来说至关重要。

🔍 什么是数字图像?

数字图像是由像素组成的二维矩阵,每个像素包含颜色信息。根据颜色深度的不同,图像可以分为:

  • 灰度图像:每个像素只有一个亮度值,通常范围是0-255
  • 彩色图像:通常使用RGB色彩空间,每个像素有红、绿、蓝三个通道值
  • 二值图像:每个像素只有0或1两个值
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# 创建一个简单的示例图像
def create_sample_image():
    """创建一个示例图像用于演示"""
    # 创建一个100x100的随机灰度图像
    gray_img = np.random.randint(0, 256, (100, 100), dtype=np.uint8)
    
    # 创建一个彩色图像
    color_img = np.random.randint(0, 256, (100, 100, 3), dtype=np.uint8)
    
    return gray_img, color_img

# 展示图像的基本属性
gray_sample, color_sample = create_sample_image()
print(f"灰度图像形状: {gray_sample.shape}")
print(f"彩色图像形状: {color_sample.shape}")
print(f"灰度图像数据类型: {gray_sample.dtype}")
print(f"彩色图像数据类型: {color_sample.dtype}")

📚 NumPy基础回顾

在深入图像处理之前,让我们快速回顾一些NumPy的核心概念,这些将在后续的图像操作中频繁使用。

数组创建与基本操作

# 创建不同类型的数组
zeros_array = np.zeros((3, 4))  # 全零数组
ones_array = np.ones((2, 3, 4))  # 全一数组
random_array = np.random.rand(5, 5)  # 随机数组
identity_matrix = np.eye(4)  # 单位矩阵

print("全零数组:")
print(zeros_array)
print("\n单位矩阵:")
print(identity_matrix)

# 数组索引和切片
sample_array = np.arange(24).reshape(4, 6)
print(f"\n原始数组:\n{sample_array}")

# 基本索引
print(f"第一行: {sample_array[0]}")
print(f"第三列: {sample_array[:, 2]}")

# 切片操作
print(f"前两行前三列:\n{sample_array[:2, :3]}")

数学运算和统计函数

# 数学运算
arr1 = np.array([[1, 2], [3, 4]])
arr2 = np.array([[5, 6], [7, 8]])

print("数组加法:")
print(arr1 + arr2)

print("数组乘法(元素级):")
print(arr1 * arr2)

print("矩阵乘法:")
print(np.dot(arr1, arr2))

# 统计函数
data = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
print(f"均值: {np.mean(data)}")
print(f"标准差: {np.std(data)}")
print(f"最大值: {np.max(data)}")
print(f"最小值: {np.min(data)}")

🖼️ 图像读取与显示

现在让我们开始真正的图像处理之旅。首先需要了解如何读取和显示图像。

import cv2
from PIL import Image
import matplotlib.pyplot as plt

def read_image_with_pil(image_path):
    """使用PIL读取图像"""
    try:
        img = Image.open(image_path)
        return np.array(img)
    except Exception as e:
        print(f"读取图像时出错: {e}")
        return None

def read_image_with_opencv(image_path):
    """使用OpenCV读取图像"""
    try:
        img = cv2.imread(image_path)
        if img is not None:
            # OpenCV默认BGR格式,转换为RGB
            img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        return img
    except Exception as e:
        print(f"读取图像时出错: {e}")
        return None

# 创建测试图像
test_img = np.random.randint(0, 256, (200, 200, 3), dtype=np.uint8)
plt.figure(figsize=(8, 6))
plt.imshow(test_img)
plt.title("测试图像")
plt.axis('off')
plt.show()

图像文件

选择库

PIL/Pillow

OpenCV

matplotlib

读取为numpy数组

图像处理操作

保存/显示结果

🎨 像素值基础操作

图像处理的核心是对像素值的操作。让我们探索一些基本的像素处理技术。

灰度化处理

将彩色图像转换为灰度图像是图像处理中的常见操作:

def rgb_to_grayscale(rgb_image):
    """
    将RGB图像转换为灰度图像
    使用标准权重: R*0.299 + G*0.587 + B*0.114
    """
    if len(rgb_image.shape) == 3:
        r, g, b = rgb_image[:,:,0], rgb_image[:,:,1], rgb_image[:,:,2]
        gray = 0.299 * r + 0.587 * g + 0.114 * b
        return gray.astype(np.uint8)
    else:
        return rgb_image

# 创建示例彩色图像
sample_rgb = np.random.randint(0, 256, (100, 100, 3), dtype=np.uint8)
gray_version = rgb_to_grayscale(sample_rgb)

print(f"原图像形状: {sample_rgb.shape}")
print(f"灰度图像形状: {gray_version.shape}")

# 显示对比
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].imshow(sample_rgb)
axes[0].set_title("原彩色图像")
axes[0].axis('off')

axes[1].imshow(gray_version, cmap='gray')
axes[1].set_title("灰度图像")
axes[1].axis('off')
plt.show()

亮度调整

调整图像亮度是最基本的图像增强操作之一:

def adjust_brightness(image, factor):
    """
    调整图像亮度
    factor > 1: 变亮
    factor < 1: 变暗
    factor = 1: 不变
    """
    # 确保结果在0-255范围内
    adjusted = np.clip(image * factor, 0, 255)
    return adjusted.astype(np.uint8)

def add_brightness_offset(image, offset):
    """
    通过增加偏移量调整亮度
    offset > 0: 变亮
    offset < 0: 变暗
    """
    adjusted = np.clip(image.astype(np.int16) + offset, 0, 255)
    return adjusted.astype(np.uint8)

# 创建测试图像
test_image = np.random.randint(50, 200, (100, 100, 3), dtype=np.uint8)

# 应用不同的亮度调整
brighter = adjust_brightness(test_image, 1.5)
darker = adjust_brightness(test_image, 0.7)
offset_bright = add_brightness_offset(test_image, 50)
offset_dark = add_brightness_offset(test_image, -30)

# 显示结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
images = [test_image, brighter, darker, test_image, offset_bright, offset_dark]
titles = ['原图', '亮度×1.5', '亮度×0.7', '原图', '亮度+50', '亮度-30']

for i, (img, title) in enumerate(zip(images, titles)):
    row, col = i // 3, i % 3
    axes[row, col].imshow(img)
    axes[row, col].set_title(title)
    axes[row, col].axis('off')

plt.tight_layout()
plt.show()

对比度调整

对比度调整可以增强图像的视觉效果:

def adjust_contrast(image, factor):
    """
    调整图像对比度
    factor > 1: 增强对比度
    factor < 1: 降低对比度
    factor = 1: 不变
    """
    # 计算平均亮度
    mean_value = np.mean(image)
    
    # 围绕平均值调整对比度
    adjusted = mean_value + factor * (image - mean_value)
    
    # 确保结果在有效范围内
    adjusted = np.clip(adjusted, 0, 255)
    return adjusted.astype(np.uint8)

# 测试对比度调整
original = np.random.randint(50, 200, (100, 100, 3), dtype=np.uint8)
high_contrast = adjust_contrast(original, 1.5)
low_contrast = adjust_contrast(original, 0.5)

# 显示结果
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(original)
axes[0].set_title("原图")
axes[0].axis('off')

axes[1].imshow(high_contrast)
axes[1].set_title("高对比度 (factor=1.5)")
axes[1].axis('off')

axes[2].imshow(low_contrast)
axes[2].set_title("低对比度 (factor=0.5)")
axes[2].axis('off')

plt.show()

🔧 高级像素操作技巧

掌握了基础操作后,让我们探索一些更高级的像素处理技术。

直方图均衡化

直方图均衡化是一种常用的图像增强技术,可以改善图像的对比度:

def histogram_equalization(image):
    """
    对灰度图像进行直方图均衡化
    """
    if len(image.shape) == 3:
        image = rgb_to_grayscale(image)
    
    # 计算直方图
    hist, bins = np.histogram(image.flatten(), 256, [0, 256])
    
    # 计算累积分布函数
    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max() / cdf.max()
    
    # 构建查找表
    cdf_m = np.ma.masked_equal(cdf, 0)
    cdf_m = (cdf_m - cdf_m.min()) * 255 / (cdf_m.max() - cdf_m.min())
    cdf = np.ma.filled(cdf_m, 0).astype('uint8')
    
    # 应用查找表
    equalized = cdf[image]
    
    return equalized

# 创建测试图像
test_gray = np.random.randint(0, 256, (100, 100), dtype=np.uint8)

# 应用直方图均衡化
equalized_img = histogram_equalization(test_gray)

# 显示直方图对比
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 原始图像和均衡化图像
axes[0, 0].imshow(test_gray, cmap='gray')
axes[0, 0].set_title("原始图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(equalized_img, cmap='gray')
axes[0, 1].set_title("均衡化后图像")
axes[0, 1].axis('off')

# 直方图
axes[1, 0].hist(test_gray.flatten(), bins=256, range=[0, 256], alpha=0.7)
axes[1, 0].set_title("原始图像直方图")
axes[1, 0].set_xlabel("像素值")
axes[1, 0].set_ylabel("频次")

axes[1, 1].hist(equalized_img.flatten(), bins=256, range=[0, 256], alpha=0.7)
axes[1, 1].set_title("均衡化后直方图")
axes[1, 1].set_xlabel("像素值")
axes[1, 1].set_ylabel("频次")

plt.tight_layout()
plt.show()

阈值分割

阈值分割是将图像转换为二值图像的重要技术:

def simple_threshold(image, threshold):
    """
    简单阈值分割
    """
    binary = np.where(image >= threshold, 255, 0)
    return binary.astype(np.uint8)

def otsu_threshold(image):
    """
    Otsu自动阈值算法
    """
    if len(image.shape) == 3:
        image = rgb_to_grayscale(image)
    
    # 计算直方图
    hist, _ = np.histogram(image.flatten(), bins=256, range=[0, 256])
    
    # 归一化直方图
    hist_norm = hist / hist.sum()
    
    # 计算累积直方图
    Q = hist_norm.cumsum()
    
    # 初始化参数
    bins = np.arange(256)
    
    fn_min = np.inf
    thresh = -1
    
    for i in range(1, 256):
        p1, p2 = np.hsplit(hist_norm, [i])  # 概率
        q1, q2 = Q[i], Q[255] - Q[i]  # 类别和
        
        if q1 < 1.e-6 or q2 < 1.e-6:
            continue
            
        b1, b2 = np.hsplit(bins, [i])  # 权重
        
        # 寻找均值和方差
        m1, m2 = np.sum(p1 * b1) / q1, np.sum(p2 * b2) / q2
        v1, v2 = np.sum(((b1 - m1) ** 2) * p1) / q1, np.sum(((b2 - m2) ** 2) * p2) / q2
        
        # 计算类内方差
        fn = v1 * q1 + v2 * q2
        if fn < fn_min:
            fn_min = fn
            thresh = i
    
    return simple_threshold(image, thresh)

# 创建测试图像
gradient_img = np.tile(np.linspace(0, 255, 100), (100, 1)).astype(np.uint8)

# 应用不同的阈值方法
simple_thresh = simple_threshold(gradient_img, 128)
otsu_thresh = otsu_threshold(gradient_img)

# 显示结果
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].imshow(gradient_img, cmap='gray')
axes[0].set_title("原始梯度图像")
axes[0].axis('off')

axes[1].imshow(simple_thresh, cmap='gray')
axes[1].set_title("简单阈值分割 (threshold=128)")
axes[1].axis('off')

axes[2].imshow(otsu_thresh, cmap='gray')
axes[2].set_title("Otsu自动阈值分割")
axes[2].axis('off')

plt.show()

🌈 颜色空间转换

不同的颜色空间适用于不同的应用场景。让我们看看如何在各种颜色空间之间转换。

def rgb_to_hsv(rgb):
    """
    RGB转HSV颜色空间
    """
    rgb = rgb.astype('float') / 255.0
    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
    
    max_val = np.maximum(np.maximum(r, g), b)
    min_val = np.minimum(np.minimum(r, g), b)
    diff = max_val - min_val
    
    # 计算色调(Hue)
    h = np.zeros_like(max_val)
    mask = diff != 0
    
    # 红色主导
    red_mask = (max_val == r) & mask
    h[red_mask] = (60 * ((g[red_mask] - b[red_mask]) / diff[red_mask]) + 360) % 360
    
    # 绿色主导
    green_mask = (max_val == g) & mask
    h[green_mask] = (60 * ((b[green_mask] - r[green_mask]) / diff[green_mask]) + 120) % 360
    
    # 蓝色主导
    blue_mask = (max_val == b) & mask
    h[blue_mask] = (60 * ((r[blue_mask] - g[blue_mask]) / diff[blue_mask]) + 240) % 360
    
    # 计算饱和度(Saturation)
    s = np.zeros_like(max_val)
    s[max_val != 0] = diff[max_val != 0] / max_val[max_val != 0]
    
    # 计算明度(Value)
    v = max_val
    
    hsv = np.stack([h/360, s, v], axis=2)
    return (hsv * 255).astype(np.uint8)

def hsv_to_rgb(hsv):
    """
    HSV转RGB颜色空间
    """
    hsv = hsv.astype('float')
    h, s, v = hsv[:,:,0]/255*360, hsv[:,:,1]/255, hsv[:,:,2]/255
    
    c = v * s
    x = c * (1 - np.abs((h/60) % 2 - 1))
    m = v - c
    
    rgb_prime = np.zeros((h.shape[0], h.shape[1], 3))
    
    mask1 = (h >= 0) & (h < 60)
    rgb_prime[mask1] = np.stack([c[mask1], x[mask1], np.zeros_like(c[mask1])], axis=1)
    
    mask2 = (h >= 60) & (h < 120)
    rgb_prime[mask2] = np.stack([x[mask2], c[mask2], np.zeros_like(c[mask2])], axis=1)
    
    mask3 = (h >= 120) & (h < 180)
    rgb_prime[mask3] = np.stack([np.zeros_like(c[mask3]), c[mask3], x[mask3]], axis=1)
    
    mask4 = (h >= 180) & (h < 240)
    rgb_prime[mask4] = np.stack([np.zeros_like(c[mask4]), x[mask4], c[mask4]], axis=1)
    
    mask5 = (h >= 240) & (h < 300)
    rgb_prime[mask5] = np.stack([x[mask5], np.zeros_like(c[mask5]), c[mask5]], axis=1)
    
    mask6 = (h >= 300) & (h < 360)
    rgb_prime[mask6] = np.stack([c[mask6], np.zeros_like(c[mask6]), x[mask6]], axis=1)
    
    rgb = (rgb_prime + m[:,:,np.newaxis]) * 255
    return np.clip(rgb, 0, 255).astype(np.uint8)

# 创建测试图像
test_rgb = np.random.randint(0, 256, (100, 100, 3), dtype=np.uint8)

# 转换到HSV再转回RGB
hsv_img = rgb_to_hsv(test_rgb)
recovered_rgb = hsv_to_rgb(hsv_img)

# 分离HSV通道
h_channel = hsv_img[:,:,0]
s_channel = hsv_img[:,:,1]
v_channel = hsv_img[:,:,2]

# 显示结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].imshow(test_rgb)
axes[0, 0].set_title("原RGB图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(recovered_rgb)
axes[0, 1].set_title("HSV转换后恢复的RGB")
axes[0, 1].axis('off')

axes[0, 2].imshow(np.abs(test_rgb.astype(int) - recovered_rgb.astype(int)), cmap='hot')
axes[0, 2].set_title("转换误差")
axes[0, 2].axis('off')

axes[1, 0].imshow(h_channel, cmap='hsv')
axes[1, 0].set_title("Hue通道")
axes[1, 0].axis('off')

axes[1, 1].imshow(s_channel, cmap='gray')
axes[1, 1].set_title("Saturation通道")
axes[1, 1].axis('off')

axes[1, 2].imshow(v_channel, cmap='gray')
axes[1, 2].set_title("Value通道")
axes[1, 2].axis('off')

plt.tight_layout()
plt.show()

🧮 几何变换

几何变换包括平移、旋转、缩放等操作,在图像处理中非常重要。

def translate_image(image, tx, ty):
    """
    平移图像
    tx: x方向平移量
    ty: y方向平移量
    """
    rows, cols = image.shape[:2]
    
    # 创建输出图像
    if len(image.shape) == 3:
        translated = np.zeros_like(image)
    else:
        translated = np.zeros_like(image)
    
    # 计算新坐标
    for i in range(rows):
        for j in range(cols):
            new_i = i - ty
            new_j = j - tx
            
            # 检查边界
            if 0 <= new_i < rows and 0 <= new_j < cols:
                if len(image.shape) == 3:
                    translated[i, j, :] = image[new_i, new_j, :]
                else:
                    translated[i, j] = image[new_i, new_j]
    
    return translated

def rotate_image(image, angle):
    """
    旋转变换
    angle: 旋转角度(度)
    """
    # 转换为弧度
    angle_rad = np.radians(angle)
    cos_angle = np.cos(angle_rad)
    sin_angle = np.sin(angle_rad)
    
    rows, cols = image.shape[:2]
    
    # 计算中心点
    center_x, center_y = cols // 2, rows // 2
    
    # 创建输出图像
    if len(image.shape) == 3:
        rotated = np.zeros_like(image)
    else:
        rotated = np.zeros_like(image)
    
    # 执行旋转变换
    for i in range(rows):
        for j in range(cols):
            # 相对于中心点的坐标
            x = j - center_x
            y = i - center_y
            
            # 逆向变换计算原坐标
            src_x = x * cos_angle - y * sin_angle + center_x
            src_y = x * sin_angle + y * cos_angle + center_y
            
            # 检查是否在原图像范围内
            if 0 <= src_x < cols and 0 <= src_y < rows:
                src_x_int = int(src_x)
                src_y_int = int(src_y)
                
                if len(image.shape) == 3:
                    rotated[i, j, :] = image[src_y_int, src_x_int, :]
                else:
                    rotated[i, j] = image[src_y_int, src_x_int]
    
    return rotated

def scale_image(image, scale_x, scale_y):
    """
    缩放图像
    scale_x: x方向缩放因子
    scale_y: y方向缩放因子
    """
    rows, cols = image.shape[:2]
    new_rows = int(rows * scale_y)
    new_cols = int(cols * scale_x)
    
    # 创建输出图像
    if len(image.shape) == 3:
        scaled = np.zeros((new_rows, new_cols, image.shape[2]), dtype=image.dtype)
    else:
        scaled = np.zeros((new_rows, new_cols), dtype=image.dtype)
    
    # 执行缩放变换
    for i in range(new_rows):
        for j in range(new_cols):
            # 计算对应原图像的位置
            src_i = int(i / scale_y)
            src_j = int(j / scale_x)
            
            # 边界检查
            if 0 <= src_i < rows and 0 <= src_j < cols:
                if len(image.shape) == 3:
                    scaled[i, j, :] = image[src_i, src_j, :]
                else:
                    scaled[i, j] = image[src_i, src_j]
    
    return scaled

# 创建测试图像
test_pattern = np.zeros((100, 100, 3), dtype=np.uint8)
test_pattern[20:80, 20:80] = [255, 0, 0]  # 红色方块
test_pattern[30:70, 30:70] = [0, 255, 0]  # 绿色方块

# 应用几何变换
translated_img = translate_image(test_pattern, 10, 15)
rotated_img = rotate_image(test_pattern, 30)
scaled_img = scale_image(test_pattern, 1.5, 1.5)

# 显示结果
fig, axes = plt.subplots(2, 2, figsize=(12, 12))

axes[0, 0].imshow(test_pattern)
axes[0, 0].set_title("原图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(translated_img)
axes[0, 1].set_title("平移变换 (tx=10, ty=15)")
axes[0, 1].axis('off')

axes[1, 0].imshow(rotated_img)
axes[1, 0].set_title("旋转变换 (30°)")
axes[1, 0].axis('off')

axes[1, 1].imshow(scaled_img)
axes[1, 1].set_title("缩放变换 (1.5x)")
axes[1, 1].axis('off')

plt.tight_layout()
plt.show()

输入图像

几何变换

平移

旋转

缩放

仿射变换

输出图像

📊 图像滤波与卷积

滤波是图像处理中的核心概念,用于去噪、边缘检测等任务。

基础卷积操作

def convolve_2d(image, kernel):
    """
    2D卷积操作
    """
    # 获取图像和核的尺寸
    img_rows, img_cols = image.shape
    k_rows, k_cols = kernel.shape
    
    # 计算填充大小
    pad_rows = k_rows // 2
    pad_cols = k_cols // 2
    
    # 对图像进行填充
    padded_img = np.pad(image, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode='constant')
    
    # 创建输出图像
    output = np.zeros_like(image)
    
    # 执行卷积
    for i in range(img_rows):
        for j in range(img_cols):
            # 提取局部区域
            region = padded_img[i:i+k_rows, j:j+k_cols]
            # 计算卷积
            output[i, j] = np.sum(region * kernel)
    
    return output

# 定义常用滤波器核
blur_kernel = np.ones((3, 3)) / 9  # 平滑滤波器
sharpen_kernel = np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]])  # 锐化滤波器
edge_kernel = np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]])  # 边缘检测滤波器

# 创建测试图像
test_image = np.random.randint(0, 256, (50, 50), dtype=np.uint8)

# 应用不同的滤波器
blurred = convolve_2d(test_image, blur_kernel)
sharpened = convolve_2d(test_image, sharpen_kernel)
edged = convolve_2d(test_image, edge_kernel)

# 确保结果在有效范围内
blurred = np.clip(blurred, 0, 255).astype(np.uint8)
sharpened = np.clip(sharpened, 0, 255).astype(np.uint8)
edged = np.clip(edged, 0, 255).astype(np.uint8)

# 显示结果
fig, axes = plt.subplots(2, 2, figsize=(12, 12))

axes[0, 0].imshow(test_image, cmap='gray')
axes[0, 0].set_title("原图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(blurred, cmap='gray')
axes[0, 1].set_title("模糊滤波")
axes[0, 1].axis('off')

axes[1, 0].imshow(sharpened, cmap='gray')
axes[1, 0].set_title("锐化滤波")
axes[1, 0].axis('off')

axes[1, 1].imshow(edged, cmap='gray')
axes[1, 1].set_title("边缘检测")
axes[1, 1].axis('off')

plt.tight_layout()
plt.show()

高斯滤波

高斯滤波是一种重要的平滑滤波器,能够有效去除噪声:

def gaussian_kernel(size, sigma):
    """
    生成高斯核
    size: 核大小(必须为奇数)
    sigma: 标准差
    """
    ax = np.arange(-size // 2 + 1., size // 2 + 1.)
    xx, yy = np.meshgrid(ax, ax)
    kernel = np.exp(-(xx**2 + yy**2) / (2 * sigma**2))
    return kernel / np.sum(kernel)

def gaussian_filter(image, kernel_size, sigma):
    """
    高斯滤波
    """
    kernel = gaussian_kernel(kernel_size, sigma)
    return convolve_2d(image, kernel)

# 生成带有噪声的测试图像
clean_image = np.zeros((100, 100), dtype=np.uint8)
clean_image[30:70, 30:70] = 255  # 白色方块

# 添加噪声
noisy_image = clean_image + np.random.normal(0, 20, clean_image.shape)
noisy_image = np.clip(noisy_image, 0, 255).astype(np.uint8)

# 应用不同参数的高斯滤波
gaussian_filtered1 = gaussian_filter(noisy_image, 5, 1.0)
gaussian_filtered2 = gaussian_filter(noisy_image, 9, 2.0)

# 显示结果
fig, axes = plt.subplots(2, 2, figsize=(12, 12))

axes[0, 0].imshow(clean_image, cmap='gray')
axes[0, 0].set_title("干净图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(noisy_image, cmap='gray')
axes[0, 1].set_title("含噪声图像")
axes[0, 1].axis('off')

axes[1, 0].imshow(gaussian_filtered1, cmap='gray')
axes[1, 0].set_title("高斯滤波 (5x5, σ=1.0)")
axes[1, 0].axis('off')

axes[1, 1].imshow(gaussian_filtered2, cmap='gray')
axes[1, 1].set_title("高斯滤波 (9x9, σ=2.0)")
axes[1, 1].axis('off')

plt.tight_layout()
plt.show()

# 显示高斯核
kernel_5_1 = gaussian_kernel(5, 1.0)
kernel_9_2 = gaussian_kernel(9, 2.0)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
im1 = axes[0].imshow(kernel_5_1, cmap='hot')
axes[0].set_title("5x5 高斯核 (σ=1.0)")
plt.colorbar(im1, ax=axes[0])

im2 = axes[1].imshow(kernel_9_2, cmap='hot')
axes[1].set_title("9x9 高斯核 (σ=2.0)")
plt.colorbar(im2, ax=axes[1])

plt.show()

🔍 边缘检测技术

边缘检测是计算机视觉中的重要技术,用于识别图像中的物体边界。

Sobel算子

Sobel算子是经典的边缘检测算子:

def sobel_edge_detection(image):
    """
    使用Sobel算子进行边缘检测
    """
    # Sobel算子核
    sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    
    # 计算x和y方向的梯度
    grad_x = convolve_2d(image, sobel_x)
    grad_y = convolve_2d(image, sobel_y)
    
    # 计算梯度幅值
    magnitude = np.sqrt(grad_x**2 + grad_y**2)
    
    # 计算梯度方向
    direction = np.arctan2(grad_y, grad_x)
    
    return magnitude, direction, grad_x, grad_y

# 创建测试图像
test_shape = np.zeros((100, 100), dtype=np.uint8)
test_shape[20:80, 20:80] = 255  # 正方形
test_shape[30:70, 30:70] = 0    # 内部正方形

# 应用Sobel边缘检测
magnitude, direction, grad_x, grad_y = sobel_edge_detection(test_shape)

# 显示结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].imshow(test_shape, cmap='gray')
axes[0, 0].set_title("原图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(np.abs(grad_x), cmap='gray')
axes[0, 1].set_title("X方向梯度")
axes[0, 1].axis('off')

axes[0, 2].imshow(np.abs(grad_y), cmap='gray')
axes[0, 2].set_title("Y方向梯度")
axes[0, 2].axis('off')

axes[1, 0].imshow(magnitude, cmap='gray')
axes[1, 0].set_title("梯度幅值")
axes[1, 0].axis('off')

axes[1, 1].imshow(direction, cmap='hsv')
axes[1, 1].set_title("梯度方向")
axes[1, 1].axis('off')

axes[1, 2].imshow(np.clip(magnitude, 0, 255), cmap='gray')
axes[1, 2].set_title("边缘检测结果")
axes[1, 2].axis('off')

plt.tight_layout()
plt.show()

Canny边缘检测

Canny边缘检测是更为先进的边缘检测算法:

def canny_edge_detection(image, low_threshold=50, high_threshold=150):
    """
    简化的Canny边缘检测实现
    """
    # 1. 高斯滤波降噪
    smoothed = gaussian_filter(image, 5, 1.0)
    
    # 2. 计算梯度
    sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    
    grad_x = convolve_2d(smoothed, sobel_x)
    grad_y = convolve_2d(smoothed, sobel_y)
    
    magnitude = np.sqrt(grad_x**2 + grad_y**2)
    direction = np.arctan2(grad_y, grad_x)
    
    # 3. 非极大值抑制
    suppressed = non_maximum_suppression(magnitude, direction)
    
    # 4. 双阈值检测和边缘连接
    edges = double_threshold(suppressed, low_threshold, high_threshold)
    
    return edges

def non_maximum_suppression(magnitude, direction):
    """
    非极大值抑制
    """
    rows, cols = magnitude.shape
    suppressed = np.zeros_like(magnitude)
    
    # 将角度量化为四个方向
    angle = np.rad2deg(direction) % 180
    
    for i in range(1, rows-1):
        for j in range(1, cols-1):
            # 根据梯度方向比较相邻像素
            if (0 <= angle[i,j] < 22.5) or (157.5 <= angle[i,j] <= 180):
                neighbors = [magnitude[i, j-1], magnitude[i, j+1]]
            elif 22.5 <= angle[i,j] < 67.5:
                neighbors = [magnitude[i-1, j+1], magnitude[i+1, j-1]]
            elif 67.5 <= angle[i,j] < 112.5:
                neighbors = [magnitude[i-1, j], magnitude[i+1, j]]
            else:  # 112.5 <= angle[i,j] < 157.5
                neighbors = [magnitude[i-1, j-1], magnitude[i+1, j+1]]
            
            # 如果当前像素是局部最大值,则保留
            if magnitude[i,j] >= max(neighbors):
                suppressed[i,j] = magnitude[i,j]
    
    return suppressed

def double_threshold(image, low_thresh, high_thresh):
    """
    双阈值检测
    """
    strong_edges = (image >= high_thresh)
    weak_edges = (image >= low_thresh) & (image < high_thresh)
    
    # 简单的边缘连接:如果弱边缘与强边缘相连则保留
    edges = np.zeros_like(image, dtype=bool)
    edges[strong_edges] = True
    
    # 检查弱边缘是否与强边缘相连
    for i in range(1, image.shape[0]-1):
        for j in range(1, image.shape[1]-1):
            if weak_edges[i,j]:
                # 检查8邻域是否有强边缘
                if np.any(strong_edges[i-1:i+2, j-1:j+2]):
                    edges[i,j] = True
    
    return edges.astype(np.uint8) * 255

# 创建复杂测试图像
complex_shape = np.zeros((100, 100), dtype=np.uint8)
# 圆形
y, x = np.ogrid[:100, :100]
mask = (x-50)**2 + (y-50)**2 <= 25**2
complex_shape[mask] = 255

# 应用Canny边缘检测
canny_edges = canny_edge_detection(complex_shape, 30, 100)

# 显示结果
fig, axes = plt.subplots(1, 2, figsize=(12, 6))

axes[0].imshow(complex_shape, cmap='gray')
axes[0].set_title("原图像")
axes[0].axis('off')

axes[1].imshow(canny_edges, cmap='gray')
axes[1].set_title("Canny边缘检测")
axes[1].axis('off')

plt.show()

🎯 形态学操作

形态学操作主要用于二值图像处理,包括腐蚀、膨胀等操作。

def erosion(image, kernel):
    """
    腐蚀操作
    """
    rows, cols = image.shape
    k_rows, k_cols = kernel.shape
    
    pad_rows = k_rows // 2
    pad_cols = k_cols // 2
    
    padded_img = np.pad(image, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode='constant')
    result = np.zeros_like(image)
    
    for i in range(rows):
        for j in range(cols):
            region = padded_img[i:i+k_rows, j:j+k_cols]
            # 如果结构元素为1的所有位置在原图像中都为1,则结果为1
            if np.all((kernel == 1) <= (region == 255)):
                result[i, j] = 255
    
    return result

def dilation(image, kernel):
    """
    膨胀操作
    """
    rows, cols = image.shape
    k_rows, k_cols = kernel.shape
    
    pad_rows = k_rows // 2
    pad_cols = k_cols // 2
    
    padded_img = np.pad(image, ((pad_rows, pad_rows), (pad_cols, pad_cols)), mode='constant')
    result = np.zeros_like(image)
    
    for i in range(rows):
        for j in range(cols):
            region = padded_img[i:i+k_rows, j:j+k_cols]
            # 如果结构元素为1的任意位置在原图像中为1,则结果为1
            if np.any((kernel == 1) & (region == 255)):
                result[i, j] = 255
    
    return result

def opening(image, kernel):
    """
    开运算:先腐蚀后膨胀
    """
    eroded = erosion(image, kernel)
    opened = dilation(eroded, kernel)
    return opened

def closing(image, kernel):
    """
    闭运算:先膨胀后腐蚀
    """
    dilated = dilation(image, kernel)
    closed = erosion(dilated, kernel)
    return closed

# 创建测试二值图像
binary_test = np.zeros((100, 100), dtype=np.uint8)
binary_test[30:70, 30:70] = 255  # 白色方块
binary_test[40:60, 40:60] = 0    # 黑色内部

# 添加一些噪声点
binary_test[10, 10] = 255
binary_test[80, 80] = 255

# 定义结构元素
structuring_element = np.ones((3, 3), dtype=np.uint8)

# 应用形态学操作
eroded_img = erosion(binary_test, structuring_element)
dilated_img = dilation(binary_test, structuring_element)
opened_img = opening(binary_test, structuring_element)
closed_img = closing(binary_test, structuring_element)

# 显示结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].imshow(binary_test, cmap='gray')
axes[0, 0].set_title("原二值图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(eroded_img, cmap='gray')
axes[0, 1].set_title("腐蚀操作")
axes[0, 1].axis('off')

axes[0, 2].imshow(dilated_img, cmap='gray')
axes[0, 2].set_title("膨胀操作")
axes[0, 2].axis('off')

axes[1, 0].imshow(opened_img, cmap='gray')
axes[1, 0].set_title("开运算")
axes[1, 0].axis('off')

axes[1, 1].imshow(closed_img, cmap='gray')
axes[1, 1].set_title("闭运算")
axes[1, 1].axis('off')

# 显示结构元素
axes[1, 2].imshow(structuring_element, cmap='gray')
axes[1, 2].set_title("结构元素")
axes[1, 2].axis('off')

plt.tight_layout()
plt.show()

📈 性能优化技巧

在处理大型图像时,性能优化变得至关重要。以下是一些实用的优化技巧:

向量化操作

import time

def pixel_by_pixel_operation(image):
    """
    逐像素操作(较慢)
    """
    rows, cols = image.shape
    result = np.zeros_like(image)
    
    for i in range(rows):
        for j in range(cols):
            result[i, j] = image[i, j] * 2
    
    return result

def vectorized_operation(image):
    """
    向量化操作(较快)
    """
    return image * 2

# 性能比较
large_image = np.random.randint(0, 256, (1000, 1000), dtype=np.uint8)

# 测量逐像素操作时间
start_time = time.time()
result1 = pixel_by_pixel_operation(large_image)
time1 = time.time() - start_time

# 测量向量化操作时间
start_time = time.time()
result2 = vectorized_operation(large_image)
time2 = time.time() - start_time

print(f"逐像素操作耗时: {time1:.4f} 秒")
print(f"向量化操作耗时: {time2:.4f} 秒")
print(f"加速比: {time1/time2:.2f}x")

# 验证结果一致性
print(f"结果一致: {np.array_equal(result1, result2)}")

内存管理

def memory_efficient_processing(image):
    """
    内存高效的图像处理
    """
    # 使用就地操作减少内存分配
    image_copy = image.copy()
    
    # 就地修改而不是创建新数组
    image_copy *= 1.2  # 亮度调整
    np.clip(image_copy, 0, 255, out=image_copy)  # 就地裁剪
    
    return image_copy.astype(np.uint8)

def memory_intensive_processing(image):
    """
    内存密集型处理
    """
    # 创建多个中间数组
    temp1 = image * 1.2
    temp2 = np.clip(temp1, 0, 255)
    result = temp2.astype(np.uint8)
    
    return result

# 内存使用监控示例
def monitor_memory_usage():
    """
    监控内存使用情况的示例
    """
    import psutil
    import os
    
    process = psutil.Process(os.getpid())
    memory_before = process.memory_info().rss / 1024 / 1024  # MB
    
    large_img = np.random.randint(0, 256, (2000, 2000), dtype=np.uint8)
    
    memory_after_allocation = process.memory_info().rss / 1024 / 1024
    
    # 处理图像
    processed = memory_efficient_processing(large_img)
    
    memory_after_processing = process.memory_info().rss / 1024 / 1024
    
    print(f"初始内存: {memory_before:.2f} MB")
    print(f"分配图像后: {memory_after_allocation:.2f} MB")
    print(f"处理后: {memory_after_processing:.2f} MB")
    
    del large_img, processed

# monitor_memory_usage()  # 取消注释以运行内存监控

🛠️ 实际应用案例

让我们通过几个实际应用案例来展示这些技术的实用性。

图像去噪系统

class ImageDenoisingSystem:
    """
    图像去噪系统
    """
    
    def __init__(self):
        self.noise_level = 0.1
    
    def add_noise(self, image, noise_type='gaussian'):
        """
        向图像添加噪声
        """
        if noise_type == 'gaussian':
            noise = np.random.normal(0, self.noise_level * 255, image.shape)
            noisy_image = image + noise
            return np.clip(noisy_image, 0, 255).astype(np.uint8)
        
        elif noise_type == 'salt_pepper':
            noisy_image = image.copy()
            # 添加椒盐噪声
            num_salt = np.ceil(self.noise_level * image.size * 0.5)
            num_pepper = np.ceil(self.noise_level * image.size * 0.5)
            
            # 添加盐噪声
            coords = [np.random.randint(0, i-1, int(num_salt)) 
                     for i in image.shape]
            noisy_image[tuple(coords)] = 255
            
            # 添加胡椒噪声
            coords = [np.random.randint(0, i-1, int(num_pepper)) 
                     for i in image.shape]
            noisy_image[tuple(coords)] = 0
            
            return noisy_image
    
    def median_filter(self, image, kernel_size=3):
        """
        中值滤波去噪
        """
        from scipy.ndimage import median_filter
        if len(image.shape) == 3:
            # 对每个通道分别处理
            filtered = np.zeros_like(image)
            for i in range(image.shape[2]):
                filtered[:,:,i] = median_filter(image[:,:,i], size=kernel_size)
            return filtered
        else:
            return median_filter(image, size=kernel_size)
    
    def bilateral_filter(self, image, d=9, sigma_color=75, sigma_space=75):
        """
        双边滤波(简化版本)
        """
        # 这里使用简化的双边滤波实现
        # 实际应用中建议使用OpenCV的cv2.bilateralFilter
        return gaussian_filter(image, 5, 1.0)  # 简化替代

# 使用去噪系统
denoiser = ImageDenoisingSystem()

# 创建测试图像
test_clean = np.zeros((100, 100, 3), dtype=np.uint8)
test_clean[20:80, 20:80] = [255, 128, 64]  # 棕色方块

# 添加噪声
noisy_gaussian = denoiser.add_noise(test_clean, 'gaussian')
noisy_salt_pepper = denoiser.add_noise(test_clean, 'salt_pepper')

# 应用去噪
denoised_gaussian = denoiser.median_filter(noisy_gaussian)
denoised_salt_pepper = denoiser.median_filter(noisy_salt_pepper)

# 显示结果
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

axes[0, 0].imshow(test_clean)
axes[0, 0].set_title("干净图像")
axes[0, 0].axis('off')

axes[0, 1].imshow(noisy_gaussian)
axes[0, 1].set_title("高斯噪声图像")
axes[0, 1].axis('off')

axes[0, 2].imshow(denoised_gaussian)
axes[0, 2].set_title("去噪后图像(高斯噪声)")
axes[0, 2].axis('off')

axes[1, 0].imshow(test_clean)
axes[1, 0].set_title("干净图像")
axes[1, 0].axis('off')

axes[1, 1].imshow(noisy_salt_pepper)
axes[1, 1].set_title("椒盐噪声图像")
axes[1, 1].axis('off')

axes[1, 2].imshow(denoised_salt_pepper)
axes[1, 2].set_title("去噪后图像(椒盐噪声)")
axes[1, 2].axis('off')

plt.tight_layout()
plt.show()

图像特征提取器

class ImageFeatureExtractor:
    """
    图像特征提取器
    """
    
    def __init__(self):
        pass
    
    def extract_color_histogram(self, image, bins=32):
        """
        提取颜色直方图特征
        """
        if len(image.shape) == 3:
            histograms = []
            for i in range(image.shape[2]):
                hist, _ = np.histogram(image[:,:,i], bins=bins, range=[0, 256])
                histograms.append(hist)
            return np.concatenate(histograms)
        else:
            hist, _ = np.histogram(image, bins=bins, range=[0, 256])
            return hist
    
    def extract_texture_features(self, image):
        """
        提取纹理特征(基于灰度共生矩阵的简化版)
        """
        if len(image.shape) == 3:
            image = rgb_to_grayscale(image)
        
        # 计算基本纹理统计量
        features = {
            'mean': np.mean(image),
            'std': np.std(image),
            'contrast': np.mean((image[:-1, :-1] - image[1:, 1:])**2),
            'energy': np.sum((image/255.0)**2),
            'entropy': -np.sum((image/255.0) * np.log2(image/255.0 + 1e-10))
        }
        
        return np.array(list(features.values()))
    
    def extract_shape_features(self, image):
        """
        提取形状特征
        """
        if len(image.shape) == 3:
            image = rgb_to_grayscale(image)
        
        # 简单的形状特征
        binary = simple_threshold(image, 128)
        area = np.sum(binary > 0)
        perimeter = np.sum(np.abs(np.diff(binary, axis=0))) + np.sum(np.abs(np.diff(binary, axis=1)))
        
        # 紧凑度
        compactness = perimeter**2 / (4 * np.pi * area) if area > 0 else 0
        
        return np.array([area, perimeter, compactness])
    
    def extract_all_features(self, image):
        """
        提取所有特征
        """
        color_features = self.extract_color_histogram(image)
        texture_features = self.extract_texture_features(image)
        shape_features = self.extract_shape_features(image)
        
        return np.concatenate([color_features, texture_features, shape_features])

# 使用特征提取器
extractor = ImageFeatureExtractor()

# 创建测试图像
feature_test = np.random.randint(0, 256, (100, 100, 3), dtype=np.uint8)
feature_test[30:70, 30:70] = [255, 0, 0]  # 红色方块

# 提取特征
all_features = extractor.extract_all_features(feature_test)
color_hist = extractor.extract_color_histogram(feature_test)
texture_feats = extractor.extract_texture_features(feature_test)
shape_feats = extractor.extract_shape_features(feature_test)

print("特征提取结果:")
print(f"总特征数: {len(all_features)}")
print(f"颜色直方图特征数: {len(color_hist)}")
print(f"纹理特征数: {len(texture_feats)}")
print(f"形状特征数: {len(shape_feats)}")

print("\n纹理特征值:")
feature_names = ['Mean', 'Std', 'Contrast', 'Energy', 'Entropy']
for name, value in zip(feature_names, texture_feats):
    print(f"{name}: {value:.4f}")

print("\n形状特征值:")
shape_names = ['Area', 'Perimeter', 'Compactness']
for name, value in zip(shape_names, shape_feats):
    print(f"{name}: {value:.4f}")

原始图像

特征提取

颜色特征

纹理特征

形状特征

颜色直方图

纹理统计量

几何属性

机器学习模型

分类/识别结果

📚 学习资源推荐

为了进一步深入学习图像处理,以下是一些优质的学习资源:

在线教程和文档

  1. NumPy官方文档 - NumPy的权威参考资料,包含了所有函数的详细说明和使用示例。

  2. SciPy Lecture Notes - 免费的科学Python教程,涵盖了NumPy、SciPy、Matplotlib等库的使用。

  3. OpenCV-Python Tutorials - OpenCV的Python接口教程,适合学习更高级的计算机视觉技术。

书籍推荐

# 推荐书籍列表
recommended_books = [
    {
        "title": "Python计算机视觉编程",
        "author": "Jan Erik Solem",
        "description": "全面介绍使用Python进行计算机视觉开发"
    },
    {
        "title": "数字图像处理",
        "author": "Rafael C. Gonzalez",
        "description": "数字图像处理的经典教材"
    },
    {
        "title": "Learning OpenCV 4",
        "author": "Adrian Kaehler",
        "description": "OpenCV 4的实用指南"
    }
]

print("📚 推荐学习书籍:")
for book in recommended_books:
    print(f"- {book['title']} by {book['author']}")
    print(f"  {book['description']}\n")

实践项目建议

  1. 图像分类器:使用提取的特征训练简单的分类器
  2. 图像拼接:实现全景图像的自动拼接
  3. 目标检测:基于传统方法的目标检测系统
  4. 图像增强应用:开发完整的图像编辑工具

🔚 总结

通过本文的学习,我们深入了解了使用Python和NumPy进行图像处理的各种技术。从基础的像素操作到高级的特征提取,每一步都展示了NumPy在图像处理中的强大能力。

关键要点回顾:

✅ 基础操作:掌握了图像读取、显示和基本像素操作
✅ 颜色处理:学会了颜色空间转换和色彩调整技术
✅ 几何变换:实现了平移、旋转、缩放等几何操作
✅ 滤波技术:理解了卷积操作和各种滤波器的应用
✅ 边缘检测:掌握了Sobel和Canny等经典边缘检测算法
✅ 形态学操作:学会了腐蚀、膨胀等二值图像处理技术
✅ 性能优化:了解了向量化操作和内存管理的重要性

这些技能不仅在学术研究中有价值,在工业应用中也同样重要。无论是医学图像分析、卫星遥感处理还是工业质量检测,掌握这些图像处理技术都能为你打开新的可能性大门。

记住,实践是最好的老师。建议你尝试自己动手实现更多有趣的图像处理项目,不断深化对这些概念的理解。随着经验的积累,你会发现图像处理是一个既有趣又实用的领域!🌟


🙌 感谢你读到这里!
🔍 技术之路没有捷径,但每一次阅读、思考和实践,都在悄悄拉近你与目标的距离。
💡 如果本文对你有帮助,不妨 👍 点赞、📌 收藏、📤 分享 给更多需要的朋友!
💬 欢迎在评论区留下你的想法、疑问或建议,我会一一回复,我们一起交流、共同成长 🌿
🔔 关注我,不错过下一篇干货!我们下期再见!✨

转载自 CSDN-专业IT技术社区

原文链接:https://blog.csdn.net/qq_41187124/article/details/157775064

文章来源转载

评论

赞0

评论列表

微信小程序
QQ小程序

关于作者

点赞数:0
关注数:0
粉丝:0
文章:0
关注标签:0
加入于:--