利用ogrid和广播机制写的低通滤波器

ogrid和numpy广播机制

向量化低通滤波详解

我将详细解释使用向量化操作实现低通滤波的代码，并介绍其中使用的函数和原理。

向量化低通滤波代码

1
2
3

y, x = np.ogrid[:image_f.shape[0], :image_f.shape[1]]
mask = np.sqrt((x - center[1])**2 + (y - center[0])**2) <= 50
image_f_lp = image_f * mask

这三行代码完全替代了原代码中的双重循环，更高效地实现了低通滤波。下面我将逐行解释：

1. 创建坐标网格

1	`y, x = np.ogrid[:image_f.shape[0], :image_f.shape[1]]`

功能说明：

这行代码创建了两个数组，分别表示图像中每个点的行坐标(y)和列坐标(x)。

函数详解：

np.ogrid：创建开放网格索引数组
- 返回两个数组：行坐标数组和列坐标数组
- 行坐标数组形状为 (H, 1)，其中 H 是图像高度
- 列坐标数组形状为 (1, W)，其中 W 是图像宽度
- 示例：对于 3×3 图像：
  1
  2
  3
  y, x = np.ogrid[0:3, 0:3] # y = [[0], [1], [2]] # x = [[0, 1, 2]]

为什么使用 ogrid 而不是 mgrid？

np.mgrid 会创建完整的坐标网格，返回形状为 (H, W) 的数组
np.ogrid 返回稀疏表示，更节省内存
通过 NumPy 的广播机制，稀疏网格可以高效计算

2. 创建距离掩模

1	`mask = np.sqrt((x - center[1])2 + (y - center[0])2) <= 50`

功能说明：

这行代码计算每个点到中心点的欧几里得距离，并创建布尔掩模，距离小于等于50的点为True，否则为False。

数学原理：

欧几里得距离公式：$d = \sqrt{(x - cx)^2 + (y - cy)^2}$
- x：列坐标
- y：行坐标
- center[1]：中心点的列索引
- center[0]：中心点的行索引

广播机制详解：

(x - center[1])：列坐标减去中心列坐标
- 形状 (1, W) 的数组减去标量，结果形状 (1, W)
(y - center[0])：行坐标减去中心行坐标
- 形状 (H, 1) 的数组减去标量，结果形状 (H, 1)

当计算平方和时，NumPy 自动广播到 (H, W) 形状

# 示例：2×2 图像，中心在 (0.5, 0.5)
x_diff = x - 0.5  # [[-0.5, 0.5]]
y_diff = y - 0.5  # [[-0.5], [0.5]]
squared_sum = x_diff**2 + y_diff**2
# 广播计算：
# [[(-0.5)**2 + (-0.5)**2, (-0.5)**2 + (0.5)**2],
#  [(0.5)**2 + (-0.5)**2, (0.5)**2 + (0.5)**2]]
# = [[0.25+0.25, 0.25+0.25], [0.25+0.25, 0.25+0.25]]
# = [[0.5, 0.5], [0.5, 0.5]]

布尔掩模：

<= 50：将距离与阈值比较，创建布尔数组
结果形状与图像相同 (H, W)
True 表示保留该频率分量，False 表示滤除

3. 应用掩模

1	`image_f_lp = image_f * mask`

功能说明：

将频率域图像与布尔掩模相乘，实现低通滤波。

工作原理：

在频率域，中心区域对应低频分量
掩模中：
- 中心区域（距离≤50）为True（相当于1）
- 边缘区域（距离>50）为False（相当于0）
乘法操作：
- 保留中心低频分量
- 滤除边缘高频分量

为什么高效？

向量化操作利用CPU/GPU的并行计算能力
避免Python循环开销
使用NumPy底层优化（C语言实现）

可视化理解

原始图像与频谱：

[原始图像]         [频谱幅度]
  ┌───────┐        ┌───────┐
  │       │        │   ●   │ ← 中心低频
  │   ●   │        │ /???\ │
  │       │        │|     || ← 高频成分
  └───────┘        └───────┘

掩模创建：

[坐标网格]         [距离计算]        [布尔掩模]
  x: 0 1 2         d: √0 √1 √2      mask: T T F
  y: 0              √1 √2 √3               T T F
     1              √2 √3 √4               F F F

滤波效果：

[原始频谱]   ×   [掩模]   =   [滤波后频谱]
  a b c         T T F         a b 0
  d e f         T T F         d e 0
  g h i         F F F         0 0 0

完整优化代码

import numpy as np
from scipy.fft import fft2, fftshift, ifft2, ifftshift
import matplotlib.pyplot as plt
from PIL import Image

# 读取图像
with Image.open("house.bmp") as im:
    image = np.array(im)
    if image.ndim == 3:  # 如果是彩色图，转为灰度
        image = np.mean(image, axis=2).astype(np.uint8)

# 计算FFT并中心化
image_f = fftshift(fft2(image))

# 计算中心点
center = [s // 2 for s in image_f.shape]

# 向量化低通滤波
rows, cols = image_f.shape
y, x = np.ogrid[:rows, :cols]  # 创建坐标网格
distance = np.sqrt((x - center[1])**2 + (y - center[0])**2)  # 计算距离
mask = distance <= 50  # 创建圆形掩模
image_f_lp = image_f * mask  # 应用掩模

# 反变换
image_lp = np.real(ifft2(ifftshift(image_f_lp)))  # 注意使用ifftshift移回

# 归一化并保存
image_lp_normalized = (image_lp - np.min(image_lp)) / (np.max(image_lp) - np.min(image_lp)) * 255
image_lp_uint8 = image_lp_normalized.astype(np.uint8)
Image.fromarray(image_lp_uint8).save("output.png")

# 可视化
fig, axs = plt.subplots(2, 2, figsize=(10, 8))

# 原始频谱（对数幅度）
axs[0, 0].imshow(np.log(1 + np.abs(image_f)), cmap='viridis')
axs[0, 0].set_title('Original Spectrum')
axs[0, 0].axis('off')

# 滤波后频谱
axs[1, 0].imshow(np.log(1 + np.abs(image_f_lp)), cmap='viridis')
axs[1, 0].set_title('Low-Pass Spectrum')
axs[1, 0].axis('off')

# 原始图像
axs[0, 1].imshow(image, cmap='gray')
axs[0, 1].set_title('Original Image')
axs[0, 1].axis('off')

# 滤波后图像
axs[1, 1].imshow(image_lp_uint8, cmap='gray')
axs[1, 1].set_title('Filtered Image')
axs[1, 1].axis('off')

plt.tight_layout()
plt.savefig("filter_comparison.png", dpi=150)
plt.show()

广播机制详解

NumPy广播机制详解

NumPy的广播机制（Broadcasting）是NumPy中最强大且最重要的特性之一，它允许不同形状的数组进行算术运算。广播机制的核心思想是：自动扩展较小的数组，使其形状与较大的数组匹配，而无需实际复制数据。

为什么需要广播机制？

想象一下，你想将一个1000×1000的矩阵与一个1×1000的行向量相加。如果没有广播机制，你需要：

将行向量复制1000次以匹配矩阵的行数
然后执行加法

广播机制自动完成这个扩展过程，无需显式复制数据，既节省内存又提高计算效率。

广播规则

广播遵循严格的规则来决定如何操作数组。规则只有两条：

维度对齐：从最右边的维度开始向左比较两个数组的形状
维度兼容：每个维度必须满足以下条件之一：
- 两个维度大小相等
- 其中一个维度大小为1
- 其中一个数组在该维度上不存在（维度数为0）

如果所有维度都兼容，则可以广播；否则会引发错误。

广播示例

示例1：标量与数组

import numpy as np

a = np.array([[1, 2, 3], 
              [4, 5, 6]])  # 形状 (2, 3)
b = 10                      # 标量 (形状 ())

c = a + b  # 广播：将b扩展为 [[10, 10, 10], [10, 10, 10]]
print(c)
"""
输出:
[[11 12 13]
 [14 15 16]]
"""

示例2：向量与矩阵

a = np.array([[1, 2, 3], 
              [4, 5, 6]])  # 形状 (2, 3)
b = np.array([10, 20, 30])  # 形状 (3,)

c = a + b  # 广播：将b扩展为 [[10, 20, 30], [10, 20, 30]]
print(c)
"""
输出:
[[11 22 33]
 [14 25 36]]
"""

示例3：不同维度

a = np.array([[[1], [2], [3]]])  # 形状 (1, 3, 1)
b = np.array([10, 20, 30])       # 形状 (3,)

c = a + b  # 广播步骤：
           # 1. 对齐维度：a(1,3,1) 和 b(3) -> b变为(1,3)
           # 2. 扩展：a扩展为(1,3,3), b扩展为(1,3,3)
print(c)
"""
输出:
[[[11 21 31]
  [12 22 32]
  [13 23 33]]]
"""

广播在低通滤波中的应用

回到我们的低通滤波代码：

y, x = np.ogrid[:rows, :cols]  # y形状(rows,1), x形状(1,cols)
center = [cy, cx]              # cy和cx是标量

# 计算距离
distance = np.sqrt((x - cx)**2 + (y - cy)**2)

这里发生了多层广播：

第一层广播：`x - cx`

x 形状: (1, cols)
cx 形状: () 标量
结果: cx被广播为 (1, cols) 的形状

第二层广播：`y - cy`

y 形状: (rows, 1)
cy 形状: () 标量
结果: cy被广播为 (rows, 1) 的形状

第三层广播：`(x - cx)2 + (y - cy)2`

(x - cx)**2 形状: (1, cols)
(y - cy)**2 形状: (rows, 1)
加法操作:
- 行维度: 1 vs rows -> 1被广播为rows
- 列维度: cols vs 1 -> 1被广播为cols
最终形状: (rows, cols)

广播机制的优势

内存效率：不创建实际的数据副本
代码简洁：避免显式循环
计算高效：利用NumPy底层优化（C实现）
可读性强：数学表达式更直观

广播机制的限制

不是所有形状都能广播：

1
2
3

a = np.ones((3, 4))
b = np.ones((2, 3))
c = a + b  # 报错：ValueError: operands could not be broadcast together with shapes (3,4) (2,3)

原因：第一个维度3 vs 2（都不为1且不相等）

广播机制可视化理解

考虑一个简单例子：

a = np.array([[0, 0, 0],    # 形状 (3,3)
              [10, 10, 10],
              [20, 20, 20]])
              
b = np.array([1, 2, 3])     # 形状 (3,)

c = a + b

广播过程：

a: 3×3        b: 3         →  b广播为3×3
[[0, 0, 0],    [1,2,3]      [[1,2,3],
 [10,10,10],   → 扩展为      [1,2,3],
 [20,20,20]]                [1,2,3]]
 
结果：
[[0+1, 0+2, 0+3],   = [[1,2,3],
 [10+1,10+2,10+3],     [11,12,13],
 [20+1,20+2,20+3]]     [21,22,23]]

在低通滤波中的具体应用

在低通滤波代码中：

# 创建坐标网格
y, x = np.ogrid[:rows, :cols]  # y: (rows,1), x: (1,cols)

# 计算每个点到中心的距离
distance = np.sqrt((x - center[1])**2 + (y - center[0])**2)

假设：

图像大小：3×3
中心点：(1,1)

计算过程：

x坐标: [[0,1,2]]   (1×3)
y坐标: [[0],[1],[2]] (3×1)

x - cx = [[0-1, 1-1, 2-1]] = [[-1,0,1]]
y - cy = [[0-1],[1-1],[2-1]] = [[-1],[0],[1]]

(x - cx)**2 = [[1,0,1]]
(y - cy)**2 = [[1],[0],[1]]

相加（广播）：
[[1,0,1]   +   [[1,1,1]     =  [[1+1, 0+1, 1+1],
           →     [0,0,0]        [1+0, 0+0, 1+0],
 [1,0,1]]        [1,1,1]]       [1+1, 0+1, 1+1]]
 
结果：
[[2,1,2],
 [1,0,1],
 [2,1,2]]
 
距离 = sqrt([[2,1,2],
             [1,0,1],
             [2,1,2]])

这就是广播机制如何高效计算每个点到中心的距离，而无需显式循环！

总结

NumPy广播机制是：

一种自动扩展数组以匹配形状的智能机制
基于严格的维度对齐和兼容规则
内存高效（无数据复制）
代码简洁且计算高效
广泛应用于科学计算、图像处理等领域

理解广播机制是掌握NumPy和高效科学计算的关键。它允许我们用简洁的数学表达式代替复杂的循环，同时保持高性能。

ogrid详解

坐标顺序的深度解析：为什么是 `y, x` 而不是 `x, y`

这个问题触及了NumPy数组索引和图像坐标系统的核心差异，也是许多初学者容易混淆的地方。让我们深入探讨一下。

核心原因：图像坐标系统 vs 笛卡尔坐标系统

1. 图像坐标系统（矩阵表示法）

在图像处理和NumPy中：

第一个维度（行索引）表示垂直方向（Y轴）
第二个维度（列索引）表示水平方向（X轴）

(0,0)-----> X (列索引)
  |
  |
  V
Y (行索引)

2. 笛卡尔坐标系统（数学表示法）

在传统数学中：

X轴表示水平方向
Y轴表示垂直方向

Y
^
|
|
(0,0)-----> X

`np.ogrid` 的工作原理

np.ogrid 返回两个数组：

第一个数组：行坐标（Y方向）
第二个数组：列坐标（X方向）

1	`y, x = np.ogrid[:height, :width]`

这里：

y 是垂直坐标（行索引）
x 是水平坐标（列索引）

为什么不能交换顺序？

错误方式：

1 2	`# 错误！交换了坐标顺序 x, y = np.ogrid[:image_f.shape[0], :image_f.shape[1]]`

结果：

x 变成了行坐标（Y方向）
y 变成了列坐标（X方向）

这会导致：

坐标含义颠倒
距离计算错误
滤波效果完全错误

具体示例分析

假设我们有一个3×3的图像：

正确方式：

y, x = np.ogrid[0:3, 0:3]

print("Y (行坐标):")
print(y)  # [[0], [1], [2]]

print("X (列坐标):")
print(x)  # [[0, 1, 2]]

错误方式：

x, y = np.ogrid[0:3, 0:3]  # 错误顺序

print("X (被错误赋值为行坐标):")
print(x)  # [[0], [1], [2]] 但变量名是x，实际是Y值

print("Y (被错误赋值为列坐标):")
print(y)  # [[0, 1, 2]] 但变量名是y，实际是X值

在距离计算中的应用

计算点(1,1)到其他点的距离：

正确计算：

center_y, center_x = 1, 1  # 中心点坐标

# 正确使用y,x
distance = np.sqrt((x - center_x)**2 + (y - center_y)**2)

错误计算：

# 错误顺序赋值
center_x, center_y = 1, 1  # 中心点坐标

# 错误使用x,y（实际含义颠倒）
distance = np.sqrt((x - center_x)**2 + (y - center_y)**2)

可视化对比

正确顺序的结果：

距离矩阵：
[[1.41, 1.00, 1.41],
 [1.00, 0.00, 1.00],
 [1.41, 1.00, 1.41]]

错误顺序的结果：

距离矩阵：
[[0.00, 1.00, 2.00],
 [1.00, 1.41, 2.24],
 [2.00, 2.24, 2.83]]

为什么NumPy采用这种设计？

数组内存布局：
- NumPy数组按行优先存储
- 第一维度（行）在内存中连续存储
图像处理惯例：
- OpenCV：image[y, x]
- Matplotlib：imshow 第一个维度是垂直方向
- PIL：Image 对象使用 (width, height) 但转换为数组后是 (height, width)
数学矩阵表示：
- 矩阵元素 $a_{ij}$ 中，i是行索引，j是列索引
- 对应图像中，i是垂直位置，j是水平位置

如何避免混淆？

最佳实践：

# 明确命名变量
height, width = image_f.shape
vertical_coords, horizontal_coords = np.ogrid[:height, :width]

# 或使用更清晰的命名
rows, cols = image_f.shape
y_coords, x_coords = np.ogrid[:rows, :cols]

# 计算距离
distance = np.sqrt((x_coords - center_x)**2 + (y_coords - center_y)**2)

完整正确代码：

# 计算中心点（注意顺序）
center_y = image_f.shape[0] // 2  # 行中心（垂直）
center_x = image_f.shape[1] // 2  # 列中心（水平）

# 创建坐标网格
y, x = np.ogrid[:image_f.shape[0], :image_f.shape[1]]

# 计算距离
distance = np.sqrt((x - center_x)**2 + (y - center_y)**2)

# 创建掩模
mask = distance <= 50

总结

顺序重要性：
- y, x = np.ogrid[...] 不是随意选择的
- 它反映了图像处理中的坐标惯例
根本原因：
- 图像数组是 [行, 列] 结构
- 行对应垂直方向（Y）
- 列对应水平方向（X）
错误后果：
- 交换顺序会导致坐标系统混乱
- 所有后续计算（如距离）都会错误
- 滤波效果完全不符合预期
记忆技巧：
- 记住 “Y 在前，X 在后”
- 对应 “行在前，列在后”
- 或联想 “高度（Y）在前，宽度（X）在后”

正确理解和使用坐标顺序是图像处理的基础，也是避免难以调试错误的关键。NumPy的这种设计虽然初看起来反直觉（与传统笛卡尔坐标不同），但符合数组内存布局和图像处理领域的通用惯例。

python > numpy

利用ogrid和广播机制写的低通滤波器

https://mingzaitown.github.io/2025/10/29/python/numpy/利用ogrid和广播机制写的低通滤波器/

作者

MingZai

发布于

2025年10月30日

许可协议

变量作用域上一篇

7.2 下一篇

利用ogrid和广播机制写的低通滤波器

ogrid和numpy广播机制

向量化低通滤波详解

向量化低通滤波代码

1. 创建坐标网格

功能说明：

函数详解：

为什么使用 ogrid 而不是 mgrid？

2. 创建距离掩模

功能说明：

数学原理：

广播机制详解：

布尔掩模：

3. 应用掩模

功能说明：

工作原理：

为什么高效？

可视化理解

原始图像与频谱：

掩模创建：

滤波效果：

完整优化代码

广播机制详解

NumPy广播机制详解

为什么需要广播机制？

广播规则

广播示例

示例1：标量与数组

示例2：向量与矩阵

示例3：不同维度

广播在低通滤波中的应用

第一层广播：x - cx

第二层广播：y - cy

第三层广播：(x - cx)**2 + (y - cy)**2

广播机制的优势

广播机制的限制

广播机制可视化理解

在低通滤波中的具体应用

总结

ogrid详解

坐标顺序的深度解析：为什么是 y, x 而不是 x, y

核心原因：图像坐标系统 vs 笛卡尔坐标系统

1. 图像坐标系统（矩阵表示法）

2. 笛卡尔坐标系统（数学表示法）

np.ogrid 的工作原理

为什么不能交换顺序？

错误方式：

结果：

具体示例分析

正确方式：

错误方式：

在距离计算中的应用

正确计算：

错误计算：

可视化对比

正确顺序的结果：

错误顺序的结果：

为什么NumPy采用这种设计？

如何避免混淆？

最佳实践：

完整正确代码：

总结

第一层广播：`x - cx`

第二层广播：`y - cy`

第三层广播：`(x - cx)2 + (y - cy)2`

坐标顺序的深度解析：为什么是 `y, x` 而不是 `x, y`

`np.ogrid` 的工作原理