Python 模块 NumPy 影像处理 AI协作
Python 模块 NumPy 影像处理 AI协作
模块——代码组织的艺术
AI 时代的思考：你真正需要掌握什么？
AI 能在几秒内写出一个完美的排序函数、一个正则表达式解析器、甚至一整段业务逻辑。
但有一件事 AI 经常犯错——代码该怎么组织。
当你让 AI 帮你写一个稍微复杂点的项目时，它经常：
- 把所有东西塞进一个文件
- 随意命名模块，边界不清
- 循环导入（circular import）导致运行报错
- 忘记
__init__.py，或者把不该暴露的东西全暴露出来 - 依赖管理一塌糊涂，
requirements.txt里写死了几十个不必要的包
代码组织能力，是 AI 时代工程师的核心竞争力之一。 AI 写函数，你设计架构。
模块怎么划分、依赖怎么管理、包怎么发布——这章就聊这些。
模块导入
import 的本质
在 C++ 里，#include 做的是文本替换——把头文件内容复制进来。
Python 的 import 不一样，它做的是执行并缓存：
- 找到目标模块文件
- 执行那个文件（是的，真的运行一遍）
- 把结果缓存在
sys.modules里 - 把模块对象绑定到当前命名空间
C++ 对比：
| C++ | Python |
|---|---|
#include <vector> |
import math |
| 编译时展开 | 运行时执行 |
| 每个编译单元独立 | 全局缓存，只执行一次 |
| namespace 隔离 | module 对象隔离 |
from...import：精准导入
from math import sqrt, pi # 只导入需要的名字
from os.path import join, exists
from typing import List, Dict, Optional # 类型注解常用
用这种方式，就可以直接写 sqrt(4) 而不是 math.sqrt(4)。
什么时候用哪种？
import math— 标准库、常用库，保留命名空间前缀更清晰from x import y— 某个名字用得很频繁，或者明确知道不会冲突import numpy as np— 社区惯例，遵守就好
别这样做：
这类似于 C++ 里 using namespace std; 写在头文件里——不是不行，但不推荐。
import 的搜索路径
Python 按以下顺序搜索模块：
- 内置模块（
sys、os这类） - 当前目录（或脚本所在目录）
PYTHONPATH环境变量指定的路径- 标准库路径
- site-packages（pip 安装的包在这里）
安装模块
pip 基础回顾
第1章已经介绍过 pip install，这里深入讲依赖管理。
pip install requests # 安装最新版
pip install requests==2.31.0 # 安装指定版本
pip install "requests>=2.28,<3.0" # 版本范围
pip install -U requests # 升级到最新版
pip uninstall requests # 卸载
pip list # 列出已安装的包
pip show requests # 查看某个包的详情
requirements.txt：依赖清单
这是 Python 项目的标配。作用就像 CMake 的 find_package + vcpkg 的 vcpkg.json——告诉别人（和 CI/CD）需要装哪些依赖。
生成 requirements.txt：
安装 requirements.txt 中的依赖：
requirements.txt 的内容长这样：
版本锁定的哲学
pip freeze 会输出所有包的精确版本，包括间接依赖（依赖的依赖）。这样做的好处是完全可复现，坏处是文件很长、版本升级麻烦。
实践中常见两种策略：
| 策略 | 写法 | 适合场景 |
|---|---|---|
| 精确锁定 | requests==2.31.0 |
生产部署、团队协作 |
| 宽松范围 | requests>=2.28 |
开源库、灵活升级 |
进阶工具：pip-tools
pip install pip-tools
# 写 requirements.in（只写直接依赖，不写版本或写范围）
# 生成精确锁定的 requirements.txt
pip-compile requirements.in
这样直接依赖和锁定文件分开管理，更专业。
自定义模块
最简单的模块
Python 里，一个 .py 文件就是一个模块。没有头文件，没有声明文件，就是这么简单。
# main.py
import utils
print(utils.greet("C++ programmer")) # Hello, C++ programmer!
print(utils.PI) # 3.14159
C++ 对比：
C++ 需要 utils.h（声明）+ utils.cpp（实现），Python 只需要 utils.py。
没有 header guard，没有 forward declaration，模块天然不会被重复执行（缓存机制）。
包（Package）：模块的集合
当项目变大，单个文件不够用了，就需要包——一个包含 __init__.py 的目录。
myproject/
main.py
mylib/
__init__.py ← 有这个文件，mylib/ 就是一个包
utils.py
math_helpers.py
io/
__init__.py
file_reader.py
file_writer.py
# 在 main.py 中使用
import mylib.utils
from mylib.math_helpers import calculate
from mylib.io.file_reader import read_csv
__init__.py 的作用
__init__.py 可以是空文件（只是标记这是个包），也可以包含代码：
# mylib/__init__.py
# 1. 控制 `from mylib import *` 时暴露哪些名字
__all__ = ["utils", "math_helpers"]
# 2. 简化导入路径（让用户不用知道内部结构）
from mylib.utils import greet
from mylib.math_helpers import calculate
# 3. 包级别的初始化代码
print("mylib loaded") # 一般不这样做，举例说明可以有代码
有了第2种用法，用户可以直接：
设计原则： __init__.py 是你包的公共接口。把常用的东西提升到这里，隐藏内部结构细节。这和 C++ 的公共头文件思想一致。
__all__：显式声明公共 API
# utils.py
__all__ = ["greet", "PI"] # 只有这两个是"公开的"
def greet(name):
return f"Hello, {name}!"
PI = 3.14159
def _internal_helper(): # 下划线开头，约定为"私有"
pass
这类似 C++ 的 public: / private:，但 Python 是约定而非强制。下划线开头的名字不会被 from x import * 导入，但你仍然可以手动导入——Python 信任程序员。
以模块方式运行代码
这是 Python 里非常常见的模式：
# utils.py
def greet(name):
return f"Hello, {name}!"
if __name__ == "__main__":
# 直接运行这个文件时，__name__ 是 "__main__"
# 被其他模块 import 时，__name__ 是 "utils"
print(greet("World"))
print("This only runs when executed directly")
实际用途：
- 测试/演示代码 — 可以直接跑文件看效果
- 命令行入口 — 既是库，也是工具
- 防止副作用 — 确保 import 时不会意外执行代码
模块引用
绝对导入 vs 相对导入
# file_reader.py 里导入 utils.py
# 绝对导入（推荐）
from mylib import utils
from mylib.utils import greet
# 相对导入（.表示当前包，..表示上级包）
from .. import utils # 上级包的 utils
from ..utils import greet # 上级包的 utils 里的 greet
from . import file_writer # 同级的 file_writer
相对导入的场景： 包内部互相引用时，相对导入更灵活——重命名包时不需要改所有内部引用。但绝对导入更易读，优先用绝对导入。
循环导入：最常见的坑
循环导入是 Python 新手（和 AI）最常遇到的问题：
运行时会报 ImportError: cannot import name 'func_a' from partially initialized module 'a'。
解决方法：
- 重构代码 — 最根本，把共同依赖提取到第三个模块
- 延迟导入 — 在函数内部 import，而不是模块顶层
- 只导入模块，不导入名字 —
import a而不是from a import func_a
C++ 类比： 循环 #include 有 header guard 保护，但循环依赖的设计问题是一样的——说明模块边界划分有问题。
扩展Python
Python 扩展的几种方式
有时候 Python 不够快，或者需要调用已有的 C/C++ 代码。Python 提供了多种扩展方式：
1. ctypes — 调用动态库（最简单）
import ctypes
# 加载动态库
lib = ctypes.CDLL("./mylib.so") # Linux/Mac
# lib = ctypes.CDLL("./mylib.dll") # Windows
# 调用函数
lib.add.argtypes = [ctypes.c_int, ctypes.c_int]
lib.add.restype = ctypes.c_int
result = lib.add(3, 4)
print(result) # 7
2. cffi — 更现代的 C 扩展接口
from cffi import FFI
ffi = FFI()
ffi.cdef("int add(int a, int b);")
lib = ffi.dlopen("./mylib.so")
print(lib.add(3, 4))
3. Cython — Python 转 C 的编译器
# mymodule.pyx (Cython 文件)
def fast_sum(list numbers):
cdef double total = 0
for n in numbers:
total += n
return total
编译后性能接近纯 C，语法是 Python 超集。
4. pybind11 — C++ 扩展的现代方案
// mymodule.cpp
#include <pybind11/pybind11.h>
namespace py = pybind11;
int add(int a, int b) { return a + b; }
PYBIND11_MODULE(mymodule, m) {
m.def("add", &add, "A function that adds two numbers");
}
pybind11 是 C++ 程序员最熟悉的方式，可以直接暴露 C++ 类、函数、枚举到 Python，性能无损。
虚拟环境
为什么需要虚拟环境？
想象一下：
- 项目 A 需要
numpy 1.20 - 项目 B 需要
numpy 1.24 - 你的系统只能安装一个版本
这就是虚拟环境要解决的问题。每个项目一个隔离的 Python 环境，互不干扰。
C++ 类比：
| C++ | Python |
|---|---|
| vcpkg / Conan | pip + venv |
| CMake find_package | import |
| build 目录隔离 | venv 目录隔离 |
| 每个项目独立编译依赖 | 每个项目独立安装依赖 |
venv 基础用法
# 创建虚拟环境
python -m venv .venv # .venv 是约定的目录名
# 激活虚拟环境
source .venv/bin/activate # Linux/Mac
.venv\Scripts\activate # Windows
# 此后所有 pip install 都只影响这个环境
pip install requests numpy
# 退出虚拟环境
deactivate
激活后，命令行提示符前会多一个 (.venv) 标记：
项目标准工作流
# 新项目开始
python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt
# 开发过程中安装新包
pip install new-package
pip freeze > requirements.txt # 更新依赖清单
# 提交代码（不要提交 .venv 目录！）
echo ".venv/" >> .gitignore
git add requirements.txt
git commit -m "Update dependencies"
进阶：pyenv + virtualenv
pyenv— 管理多个 Python 版本（类似 nvm for Node.js）virtualenv— 比内置 venv 功能更多
更现代的选择：uv
uv 是用 Rust 写的 Python 包管理器，速度快 10-100 倍：
# 安装 uv
curl -LsSf https://astral.sh/uv/install.sh | sh
# 创建虚拟环境并安装依赖
uv venv
uv pip install -r requirements.txt
# 甚至可以直接管理 Python 版本
uv python install 3.12
打包发布
现代 Python 打包：pyproject.toml
Python 打包方式经历了多次演变。现代做法是用 pyproject.toml（PEP 517/518 标准），取代老式的 setup.py。
项目结构：
mypackage/
pyproject.toml ← 项目元数据和构建配置
README.md
src/
mypackage/
__init__.py
utils.py
core.py
tests/
test_utils.py
pyproject.toml 示例：
[build-system]
requires = ["setuptools>=68", "wheel"]
build-backend = "setuptools.backends.legacy:build"
[project]
name = "mypackage"
version = "0.1.0"
description = "A sample Python package"
readme = "README.md"
license = {text = "MIT"}
requires-python = ">=3.9"
dependencies = [
"requests>=2.28",
"numpy>=1.24",
]
[project.optional-dependencies]
dev = [
"pytest>=7.0",
"black",
"mypy",
]
[project.scripts]
mypackage-cli = "mypackage.cli:main" # 命令行工具入口点
构建和安装
# 安装构建工具
pip install build
# 构建 wheel 和 sdist
python -m build
# 生成：
# dist/mypackage-0.1.0-py3-none-any.whl ← 二进制发行包
# dist/mypackage-0.1.0.tar.gz ← 源码发行包
# 本地安装（开发模式，修改立即生效）
pip install -e .
wheel 类比 C++ 的预编译库——用户安装时不需要编译。对于纯 Python 包，wheel 就是个打包好的 zip。对于含 C 扩展的包，wheel 包含了编译后的 .so/.pyd 文件。
发布到 PyPI
# 安装发布工具
pip install twine
# 发布到测试 PyPI（先在这里试试）
twine upload --repository testpypi dist/*
# 发布到正式 PyPI（需要先注册账号）
twine upload dist/*
# 别人就可以安装你的包了
pip install mypackage
PyPI 地址：https://pypi.org
发布前检查清单：
- [ ]
README.md写清楚了 - [ ] 版本号更新了
- [ ] 在 test PyPI 测试过安装
- [ ] LICENSE 文件存在
- [ ] 敏感信息（密钥、密码）没有打包进去
打包成可执行文件
有时候需要给不懂 Python 的用户分发工具，可以打包成单个可执行文件：
PyInstaller（最流行）：
pip install pyinstaller
# 打包成单个文件
pyinstaller --onefile main.py
# 生成：dist/main（Linux/Mac）或 dist/main.exe（Windows）
zipapp（Python 自带，轻量）：
# 把整个包打包成 .pyz 文件（需要目标机器有 Python）
python -m zipapp mypackage -m "mypackage.main:main" -o myapp.pyz
# 运行
python myapp.pyz
C++ 类比： PyInstaller 类似静态链接——把 Python 解释器和所有依赖都打包进去，文件大但完全独立。zipapp 类似动态链接——需要目标机器有 Python 运行时。
代码组织的哲学
最后，回到开头的问题：模块该怎么划分？
几个实用原则：
1. 单一职责
每个模块只做一件事。utils.py 里什么都放是反模式。
2. 按变化频率划分
经常变化的和稳定的分开。核心逻辑和 UI 逻辑分开。
3. 公共接口最小化__init__.py 里只暴露用户需要的，内部实现细节不要泄露出去。
4. 避免深层嵌套mylib.utils.helpers.string.format.advanced 这种路径说明结构有问题。
思考题：
- 你有一个项目，包含数据库操作、HTTP API、命令行界面三个部分。你会怎么划分模块结构？
- 如果两个模块都需要一个辅助函数，该放在哪里？
__init__.py里应该有多少代码？什么情况下应该保持为空文件？
小结
本章要点回顾：
| 概念 | Python 方式 | C++ 类比 |
|---|---|---|
| 模块 | .py 文件 |
.h + .cpp |
| 包 | 含 __init__.py 的目录 |
namespace + 目录 |
| 导入 | import x / from x import y |
#include |
| 依赖管理 | pip + requirements.txt | vcpkg / Conan |
| 环境隔离 | venv | 独立 build 目录 |
| 发布 | wheel + PyPI | 库文件 + 包管理器 |
最重要的三件事：
- 虚拟环境是标配——每个项目都用，别装全局
requirements.txt要提交——团队协作的基础- 模块边界要想清楚——这是 AI 替代不了你的地方
下一章我们聊 Python 的面向对象——你会发现它和 C++ 差异很大，但又有一些熟悉的感觉。
科学计算 — NumPy 与向量化思维
AI 时代，你真正需要掌握什么？
问 ChatGPT「如何用 NumPy 做矩阵乘法」，它会立刻给你正确答案。问它「如何对一组图像做归一化」，它也能写出漂亮的代码。
所以，NumPy 的 API 你不需要死记硬背。
但是，当 AI 生成了一段科学计算代码，你能判断它是否正确、是否高效吗？
- 这段代码用了 Python 循环还是向量化操作？性能差 100 倍你能发现吗？
- Broadcasting 的结果是你想要的 shape 吗？
- 这里用 C order 还是 Fortran order 更快？
np.dot在这里是矩阵乘法还是点积？
这才是你需要内化的东西：向量化思维、广播规则、内存模型。这些是 AI 写不出来的判断力。
为什么需要 NumPy？
Python list 的局限
C++ 程序员都懂数组：连续内存、缓存友好、SIMD 加速。Python 的 list 不是这样的：
# Python list 存的是"对象引用"，不是数据本身
nums = [1.0, 2.0, 3.0, 4.0]
# 每个元素是一个 Python float 对象，在堆上随机分布
# 遍历需要：解引用 → 检查类型 → 取值 → 下一步
// C++ vector 存的是连续的数据
std::vector<double> nums = {1.0, 2.0, 3.0, 4.0};
// 内存：[1.0][2.0][3.0][4.0] — 连续，SIMD 友好
这个差异导致什么结果？
import time
# 对 100 万个数求平方
n = 1_000_000
data = list(range(n))
# 方式 1：Python list comprehension
start = time.perf_counter()
result = [x * x for x in data]
t1 = time.perf_counter() - start
# 方式 2：NumPy
import numpy as np
arr = np.arange(n)
start = time.perf_counter()
result = arr * arr
t2 = time.perf_counter() - start
print(f"Python list: {t1*1000:.1f} ms")
print(f"NumPy: {t2*1000:.1f} ms")
print(f"加速比: {t1/t2:.0f}x")
# 典型输出：
# Python list: 85.3 ms
# NumPy: 1.2 ms
# 加速比: 71x
为什么 NumPy 这么快？三个原因： 1. 连续内存 — ndarray 就是 C 数组，缓存命中率高 2. SIMD 指令 — 底层调用 AVX/SSE，一次处理 4-8 个 double 3. 跳过解释器 — 运算在 C 层完成，Python 只负责调度
NumPy vs C++ 数组的关系
你可以把 np.ndarray 理解为 C++ 的 std::vector<T> 的 Python 包装，加上大量线性代数操作。
// C++ 视角
double* data = new double[n]; // 原始数组
// NumPy ndarray ≈ 这个指针 + shape + stride + dtype 的元数据
NumPy 基础
ndarray — 核心数据结构
import numpy as np
# 从 list 创建
a = np.array([1, 2, 3, 4, 5])
print(a.dtype) # int64
print(a.shape) # (5,)
print(a.ndim) # 1
# 二维数组（矩阵）
m = np.array([[1, 2, 3],
[4, 5, 6]])
print(m.shape) # (2, 3) — 2行3列
print(m.ndim) # 2
print(m.size) # 6 — 元素总数
dtype — 类型是你的责任
C++ 中类型是编译期确定的。NumPy 中你需要主动管理：
# 默认推断
a = np.array([1, 2, 3]) # int64
b = np.array([1.0, 2.0, 3.0]) # float64
c = np.array([1, 2, 3], dtype=np.float32) # 显式指定
# 常用 dtype
np.float32 # 单精度浮点，等同 C++ float
np.float64 # 双精度浮点，等同 C++ double（默认）
np.int32 # 32位整数
np.int64 # 64位整数
np.complex128 # 复数
# 类型转换
arr = np.array([1, 2, 3], dtype=np.int32)
arr_f = arr.astype(np.float64)
# ⚠️ 陷阱：整数溢出不报错
small = np.array([127], dtype=np.int8)
print(small + 1) # [128] → 溢出为 -128！
创建数组的常用方式
# 等差数列（类比 C++ for 循环填充）
np.arange(10) # [0, 1, 2, ..., 9]
np.arange(0, 1, 0.1) # [0.0, 0.1, ..., 0.9]
np.linspace(0, 1, 11) # [0.0, 0.1, ..., 1.0]，包含终点，共11个
# 特殊矩阵
np.zeros((3, 4)) # 3×4 全零矩阵
np.ones((2, 3)) # 2×3 全一矩阵
np.eye(4) # 4×4 单位矩阵
np.full((3, 3), 7.0) # 3×3 全为 7.0
# 随机数
np.random.rand(3, 4) # [0, 1) 均匀分布
np.random.randn(100) # 标准正态分布
np.random.randint(0, 10, size=(3, 3)) # 整数随机
# 从已有 shape 创建
a = np.array([[1, 2], [3, 4]])
np.zeros_like(a) # 与 a 同 shape 的零数组
np.ones_like(a)
向量化思维 — 本章核心
这是 C++ 程序员最需要转变的思维方式。
循环思维 vs 向量化思维
问题：对数组每个元素应用 f(x) = 2x² + 3x - 1
// C++ 思维：一个个处理
std::vector<double> y(n);
for (int i = 0; i < n; i++) {
y[i] = 2 * x[i]*x[i] + 3 * x[i] - 1;
}
# ❌ 直接翻译 C++ —— 可运行，但慢
x = np.linspace(0, 10, 1_000_000)
y = np.zeros_like(x)
for i in range(len(x)):
y[i] = 2 * x[i]**2 + 3 * x[i] - 1
# ✅ 向量化 —— 把操作作用于整个数组
y = 2 * x**2 + 3 * x - 1 # 就这一行！
向量化版本不只是更短，它快了 50-100 倍。
为什么向量化快？
Python 循环的每次迭代：
Python 解释器 dispatch
→ 检查 x[i] 的类型
→ 调用 Python float 的 __mul__
→ 创建临时 Python 对象
→ 垃圾回收
→ 重复 100 万次
NumPy 向量化：
Python 调度一次
→ C 层循环 + SIMD（AVX2 可以一次处理 4 个 double）
→ 结果直接写入预分配内存
向量化的本质：消灭 Python 层循环
# 场景：找出数组中所有大于均值的元素的平方和
data = np.random.randn(1_000_000)
# ❌ Python 循环
result = 0
mean = sum(data) / len(data)
for x in data:
if x > mean:
result += x * x
# ✅ 向量化（三行搞定，快 100x）
mask = data > data.mean() # 布尔数组
result = (data[mask] ** 2).sum()
通用函数 ufunc
NumPy 的数学函数都是向量化的（称为 ufunc）：
x = np.linspace(0, 2 * np.pi, 1000)
# ✅ 向量化数学函数
y = np.sin(x) # 不是 math.sin！
z = np.exp(-x) * np.cos(2 * x)
w = np.sqrt(x**2 + 1)
# ❌ 陷阱：用了 math 模块
import math
# math.sin(x) # 这会报错，math 函数不接受数组
# 应该用 np.sin(x)
索引和切片
基础索引
a = np.array([10, 20, 30, 40, 50])
# 与 Python list 类似
print(a[0]) # 10
print(a[-1]) # 50
print(a[1:4]) # [20, 30, 40]
print(a[::2]) # [10, 30, 50]
# 二维数组
m = np.arange(12).reshape(3, 4)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]]
print(m[1, 2]) # 6 — 第1行第2列（0-indexed）
print(m[1]) # [4, 5, 6, 7] — 整行
print(m[:, 2]) # [2, 6, 10] — 整列
print(m[0:2, 1:3]) # 子矩阵
切片是视图，不是拷贝！
这是 NumPy 与 Python list 的关键区别：
a = np.array([1, 2, 3, 4, 5])
# Python list：切片是拷贝
lst = [1, 2, 3, 4, 5]
sliced = lst[1:4]
sliced[0] = 99
print(lst) # [1, 2, 3, 4, 5] — 原 list 不变
# NumPy：切片是视图（view）
arr = np.array([1, 2, 3, 4, 5])
view = arr[1:4]
view[0] = 99
print(arr) # [ 1 99 3 4 5] — ⚠️ 原数组被修改了！
# 需要拷贝时，用 .copy()
copy = arr[1:4].copy()
copy[0] = 0
print(arr) # 不变
这是 NumPy 的设计哲学：避免不必要的内存拷贝，让你显式控制。C++ 程序员应该很熟悉这个理念——类似引用 vs 值语义。
花式索引（Fancy Indexing）
a = np.array([10, 20, 30, 40, 50])
# 用整数数组索引
idx = [0, 2, 4]
print(a[idx]) # [10, 30, 50]
# 二维索引
m = np.arange(12).reshape(3, 4)
rows = [0, 2]
cols = [1, 3]
print(m[rows, cols]) # [m[0,1], m[2,3]] = [1, 11]
# ⚠️ 花式索引返回拷贝，不是视图
fancy = a[[0, 2, 4]]
fancy[0] = 99
print(a) # 不变！
布尔索引 — 最强大的索引方式
data = np.array([3, -1, 4, -1, 5, -9, 2, 6])
# 生成布尔掩码
mask = data > 0
print(mask) # [ True False True False True False True True]
# 用掩码过滤
print(data[mask]) # [3, 4, 5, 2, 6]
# 简写（常用）
print(data[data > 0]) # 正数
print(data[data % 2 == 0]) # 偶数
# 组合条件（用 & | ~，不能用 and or not！）
print(data[(data > 0) & (data < 5)]) # 0 到 5 之间
# 配合赋值：替换所有负数为 0
data[data < 0] = 0
print(data) # [3, 0, 4, 0, 5, 0, 2, 6]
广播机制（Broadcasting）
Broadcasting 是 NumPy 最强大也最容易出错的特性。
直觉理解
a = np.array([1, 2, 3])
b = 10
# 标量与数组：b 被"广播"到与 a 相同 shape
print(a + b) # [11, 12, 13]
# 等价于 a + np.array([10, 10, 10])
# 二维 + 一维
m = np.array([[1, 2, 3],
[4, 5, 6]]) # shape (2, 3)
v = np.array([10, 20, 30]) # shape (3,)
print(m + v)
# [[11, 22, 33],
# [14, 25, 36]]
# v 被广播到每一行
Broadcasting 规则（重要！）
两个数组做运算时，从右边对齐维度，逐维比较： - 相同 → OK - 其中一个为 1 → 那个维度被广播（拉伸） - 其他 → 报错
shape (2, 3) + shape (3,)
↓ 右对齐
(2, 3) + (1, 3) ← (3,) 被当作 (1, 3)
↓ 第0维：2 vs 1，1被广播为2
(2, 3) + (2, 3) ← OK!
# 更复杂的例子
a = np.ones((3, 1, 4)) # shape (3, 1, 4)
b = np.ones((5, 4)) # shape (5, 4)
# 右对齐：(3, 1, 4) + (1, 5, 4)
# 广播后：(3, 5, 4)
print((a + b).shape) # (3, 5, 4)
实用场景：列归一化
# 每列减去列均值（常见的数据预处理操作）
data = np.array([[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0]]) # shape (3, 3)
col_mean = data.mean(axis=0) # shape (3,) — 每列的均值
normalized = data - col_mean # shape (3, 3) - shape (3,) → Broadcasting!
print(col_mean) # [4. 5. 6.]
print(normalized)
# [[-3. -3. -3.]
# [ 0. 0. 0.]
# [ 3. 3. 3.]]
常见陷阱
a = np.array([1, 2, 3]) # shape (3,)
b = np.array([[1], [2]]) # shape (2, 1)
# 这会成功，但结果是 (2, 3)，你想要的吗？
print((a + b).shape) # (2, 3)
# [[2, 3, 4],
# [3, 4, 5]]
# 如果你想要逐元素相加，shape 必须完全匹配或是 (1,) 或标量
常用变形操作
reshape — 改变形状不改变数据
a = np.arange(12) # [0, 1, 2, ..., 11]
# reshape 返回视图（尽可能）
m = a.reshape(3, 4) # 3×4 矩阵
print(m)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]]
# -1 让 NumPy 自动计算
a.reshape(2, -1) # (2, 6)，-1 自动推断为 6
a.reshape(-1) # 展平为一维，等同 a.ravel()
# 增加维度（用于 Broadcasting）
v = np.array([1, 2, 3]) # shape (3,)
col = v.reshape(-1, 1) # shape (3, 1) — 列向量
row = v.reshape(1, -1) # shape (1, 3) — 行向量
# 等价简写：
col = v[:, np.newaxis]
row = v[np.newaxis, :]
transpose — 转置
m = np.arange(6).reshape(2, 3)
print(m.shape) # (2, 3)
print(m.T.shape) # (3, 2) — 转置
# 多维数组
t = np.ones((2, 3, 4))
print(t.transpose(2, 0, 1).shape) # (4, 2, 3)
concatenate / stack — 拼接
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6], [7, 8]])
# 按行拼（竖向）
print(np.concatenate([a, b], axis=0))
# [[1, 2], [3, 4], [5, 6], [7, 8]]
# 按列拼（横向）
print(np.concatenate([a, b], axis=1))
# [[1, 2, 5, 6], [3, 4, 7, 8]]
# vstack / hstack（常用简写）
np.vstack([a, b]) # 等同 axis=0
np.hstack([a, b]) # 等同 axis=1
# stack：新增一个维度再拼
np.stack([a, b], axis=0).shape # (2, 2, 2)
np.stack([a, b], axis=2).shape # (2, 2, 2)
split — 分割
a = np.arange(12).reshape(3, 4)
# 按列分割成 2 份
left, right = np.hsplit(a, 2) # 每份 (3, 2)
# 按行分割
top, bottom = np.vsplit(a, [1]) # 分割点在第1行
# top: shape (1, 4), bottom: shape (2, 4)
线性代数
矩阵乘法
A = np.random.randn(3, 4)
B = np.random.randn(4, 5)
v = np.array([1.0, 2.0, 3.0])
w = np.array([4.0, 5.0, 6.0])
# Python 3.5+ 有 @ 运算符
C = A @ B # 矩阵乘法，shape (3, 5)
d = v @ w # 向量点积 = 32.0
# 等价写法
C = np.dot(A, B)
d = np.dot(v, w)
# ⚠️ 重要区分：
# np.dot(A, B) — 矩阵乘法（当 A, B 是二维时）
# A * B — 逐元素乘法（element-wise）！
print((A * np.ones_like(A)).shape) # (3, 4)，不是矩阵乘法
np.linalg — 线性代数工具箱
from numpy import linalg
A = np.array([[2.0, 1.0],
[1.0, 3.0]])
# 行列式
print(linalg.det(A)) # 5.0
# 逆矩阵
A_inv = linalg.inv(A)
print(A @ A_inv) # 单位矩阵（近似）
# 特征值和特征向量
eigenvalues, eigenvectors = linalg.eig(A)
print(eigenvalues) # [1.38, 3.61]（近似）
# 解线性方程组 Ax = b
b = np.array([3.0, 5.0])
x = linalg.solve(A, b) # 比 inv(A) @ b 更数值稳定
print(A @ x) # 验证：应等于 b
# SVD 分解
U, s, Vh = linalg.svd(A) # A = U @ diag(s) @ Vh
# 矩阵范数
print(linalg.norm(A)) # Frobenius 范数
print(linalg.norm(A, ord=2)) # 谱范数
与 C++ Eigen/BLAS 的对比
| 操作 | C++ Eigen | NumPy |
|---|---|---|
| 矩阵乘法 | A * B |
A @ B |
| 逐元素乘 | A.cwiseProduct(B) |
A * B |
| 转置 | A.transpose() |
A.T |
| 求逆 | A.inverse() |
linalg.inv(A) |
| 特征值 | EigenSolver |
linalg.eig(A) |
NumPy 底层调用 BLAS/LAPACK，大矩阵运算性能接近 C++ Eigen。
SciPy 简介
SciPy 是在 NumPy 基础上构建的科学计算库，相当于 C++ 里 Boost 之于标准库的关系。
常用子模块
# 优化（scipy.optimize）
from scipy import optimize
def f(x):
return (x - 2)**2 + 1
result = optimize.minimize_scalar(f)
print(result.x) # ≈ 2.0（最优解）
# 非线性方程求根
def equation(x):
return x**3 - x - 2
root = optimize.brentq(equation, 1, 2)
print(root) # ≈ 1.5214
# 信号处理（scipy.signal）
from scipy import signal
# 设计低通滤波器（类比 C++ 里用 DSP 库）
b, a = signal.butter(4, 0.1, btype='low')
filtered = signal.filtfilt(b, a, noisy_signal)
# 统计（scipy.stats）
from scipy import stats
data = np.random.randn(1000)
# t 检验
t_stat, p_value = stats.ttest_1samp(data, popmean=0)
print(f"p-value: {p_value:.4f}")
# 正态分布 PDF
x = np.linspace(-4, 4, 100)
pdf = stats.norm.pdf(x, loc=0, scale=1)
什么时候用 SciPy？
- 数值优化：梯度下降、最小化、方程求根
- 信号处理：滤波器设计、FFT、卷积
- 统计检验：t 检验、卡方检验、相关性分析
- 插值：
scipy.interpolate - 积分：
scipy.integrate（数值积分、ODE 求解器） - 稀疏矩阵：
scipy.sparse
性能提示
NumPy 快的时候
# ✅ 大数组的向量化运算
result = np.exp(large_array) * np.sin(another_array)
# ✅ 矩阵运算（调用 BLAS）
C = A @ B # 对大矩阵，接近 Eigen 的速度
# ✅ 聚合操作
total = array.sum()
max_val = array.max(axis=0)
NumPy 慢的时候
# ❌ 小数组（Python 函数调用开销 > 计算本身）
small = np.array([1.0, 2.0, 3.0])
# 对 3 个元素用 NumPy 不比 Python 快
# ❌ 需要条件逻辑的复杂操作
for i in range(n):
if condition[i]:
result[i] = f(arr[i])
else:
result[i] = g(arr[i])
# 改写为：
result = np.where(condition, f(arr), g(arr)) # ✅
# ❌ 递推/依赖上一步的结果
# 如 Fibonacci、状态机 → 无法向量化，考虑 Numba/Cython
内存布局：C order vs Fortran order
# C order（行优先，默认）：同行元素在内存中连续
# Fortran order（列优先）：同列元素在内存中连续
c_arr = np.array([[1, 2, 3], [4, 5, 6]], order='C')
f_arr = np.array([[1, 2, 3], [4, 5, 6]], order='F')
# 检查内存布局
print(c_arr.flags['C_CONTIGUOUS']) # True
print(f_arr.flags['F_CONTIGUOUS']) # True
# 按行遍历 C order 快（连续内存访问）
# 按列遍历 Fortran order 快
# 实践：如果你大量做 row-wise 操作，用 C order（默认即可）
# 如果大量做 column-wise 操作，考虑 Fortran order 或先 transpose
// 类比 C++ 二维数组：
// row-major（C++默认）: arr[row][col] → arr[row * cols + col]
// column-major（Fortran）: arr[row][col] → arr[col * rows + row]
避免中间数组
# ❌ 创建多个临时数组
a = np.ones(1_000_000)
b = np.ones(1_000_000)
c = np.ones(1_000_000)
result = a * b + c # 先创建 a*b 的临时数组，再加 c
# ✅ 用 out 参数（高级用法）
tmp = np.empty_like(a)
np.multiply(a, b, out=tmp)
np.add(tmp, c, out=tmp)
# ✅ 对于复杂表达式，考虑 numexpr
# import numexpr as ne
# result = ne.evaluate('a * b + c') # 单次遍历，无临时数组
快速诊断工具
# 检查是否连续内存（影响性能）
arr = np.ones((100, 100))
print(arr.flags)
# 查看内存占用
print(arr.nbytes) # 字节数
# 用 %timeit（Jupyter）或 timeit 模块
import timeit
t = timeit.timeit(lambda: arr.sum(axis=0), number=1000)
print(f"{t/1000*1000:.3f} ms per call")
小结：C++ 程序员的 NumPy 备忘录
| C++ 概念 | NumPy 对应 |
|---|---|
double arr[100] |
np.zeros(100) |
std::vector<double> |
np.array([...]) |
| 指针运算/引用 | 切片（视图） |
memcpy |
.copy() |
| SIMD 循环 | 向量化运算 |
| Eigen/BLAS | np.linalg / @ |
| 行优先内存 | C order（默认） |
| 列优先内存 | Fortran order |
核心心法：
- 消灭 Python 循环 — 凡是能用向量化替代的循环，都替代
- 切片是视图 — 需要独立副本时显式
.copy() - Broadcasting 从右对齐 — 脑子里想好 shape 再写代码
- dtype 要匹配 — float32 vs float64 会影响精度和速度
- AI 会写代码，但你要能判断它是否正确高效
下一步
第6章将介绍 Pandas — 当你的数据不只是数字矩阵，而是带标签的表格时，NumPy 的下一层抽象。
Python 应用实战 — 医学影像处理
🤖 AI 时代的思考：真正的编程能力是什么？
AI 可以帮你写代码。这是事实，而且这个能力还在快速增强。
但有一件事 AI 替代不了：把一个模糊的现实问题，转化为可执行的代码方案。
这个转化过程需要你： - 理解问题的本质（不是表面需求） - 把问题分解成可处理的子问题 - 知道用什么工具、用什么算法 - 判断结果是否正确
本章通过医学影像处理这个真实领域，带你走一遍完整的思路。不只是教你调 API——而是教你怎么想。
🩺 为什么选医学影像？
Python 的科学计算生态
用 C++ 做图像处理，你可能会用 OpenCV：
// C++ OpenCV
Mat img = imread("scan.png", IMREAD_GRAYSCALE);
Mat binary;
threshold(img, binary, 128, 255, THRESH_BINARY);
Python 做同样的事：
# Python scikit-image
import skimage.io as io
from skimage.filters import threshold_otsu
img = io.imread("scan.png", as_gray=True)
thresh = threshold_otsu(img)
binary = img > thresh
Python 的优势不只是代码更短——而是整个生态更适合科学计算： - NumPy：所有数据都是数组，运算直观高效 - matplotlib：可视化调试极其方便 - scikit-image：学术级图像处理算法，文档详细 - pydicom：专门处理医学影像格式 DICOM
医学影像基础概念
CT（计算机断层扫描）：用 X 射线从多角度扫描，重建出身体的三维切片图像。每张切片是一张灰度图，像素值表示 HU（Hounsfield Unit，亨氏单位）：
| 组织类型 | HU 值范围 |
|---|---|
| 空气 | -1000 |
| 肺部（充气） | -600 ~ -400 |
| 脂肪 | -100 ~ -50 |
| 软组织 | 40 ~ 80 |
| 骨骼 | 400 ~ 1000 |
MRI（磁共振成像）：利用磁场和射频脉冲成像，对软组织分辨率更好，但不用于这章的例子。
💡 物理直觉：HU 值本质是 X 射线衰减系数的归一化。水是 0，空气是 -1000，骨骼因为钙含量高所以衰减强，值很大。
🔧 环境准备
验证安装：
import numpy as np
import scipy
import matplotlib
import skimage
import pydicom
print("All good! 🎉")
📂 读取医学影像数据
DICOM 格式简介
医学影像通常用 DICOM（Digital Imaging and Communications in Medicine）格式存储。DICOM 文件不只是图像——它还包含大量元数据：患者信息、扫描参数、像素间距等。
CT_scan.dcm
├── 元数据
│ ├── PatientName: "Anonymous"
│ ├── SliceThickness: 2.5 mm
│ ├── RescaleIntercept: -1024
│ └── RescaleSlope: 1.0
└── 像素数据（PixelData）
└── 512×512 的灰度图
用 pydicom 读取 CT 图像
import pydicom
import numpy as np
import matplotlib.pyplot as plt
def load_dicom(filepath: str) -> tuple[np.ndarray, pydicom.Dataset]:
"""
读取 DICOM 文件，返回 HU 值图像和原始数据集。
为什么要转换成 HU？
原始像素值是整数存储的，需要用 RescaleSlope 和 RescaleIntercept
转换成真实的物理单位（HU），才有实际意义。
"""
ds = pydicom.dcmread(filepath)
# 获取原始像素数组
pixel_array = ds.pixel_array.astype(np.float32)
# 转换为 HU 值
slope = float(getattr(ds, 'RescaleSlope', 1.0))
intercept = float(getattr(ds, 'RescaleIntercept', 0.0))
hu_image = pixel_array * slope + intercept
return hu_image, ds
# 使用示例
hu_image, ds = load_dicom("CT_slice.dcm")
print(f"图像尺寸: {hu_image.shape}")
print(f"HU 值范围: [{hu_image.min():.0f}, {hu_image.max():.0f}]")
窗宽窗位（Window Width/Level）
CT 图像的 HU 范围很宽（-1000 到 3000），但人眼只能分辨约 40 个灰度级。窗宽窗位就是把感兴趣的 HU 范围映射到显示的灰度范围：
常用窗口预设： - 肺窗：WL=-600, WW=1500（看肺部结构） - 纵隔窗：WL=40, WW=400（看软组织） - 骨窗：WL=400, WW=1800（看骨骼）
def apply_window(hu_image: np.ndarray, wl: float, ww: float) -> np.ndarray:
"""
应用窗宽窗位，将 HU 值映射到 [0, 255]。
np.clip 把范围外的值截断，然后线性映射到 [0, 255]。
这和 C++ 里手写一个 clamp + linear_map 是完全等价的。
"""
lower = wl - ww / 2
upper = wl + ww / 2
windowed = np.clip(hu_image, lower, upper)
# 线性归一化到 [0, 255]
windowed = (windowed - lower) / (upper - lower) * 255.0
return windowed.astype(np.uint8)
# 对比不同窗口
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(apply_window(hu_image, wl=-600, ww=1500), cmap='gray')
axes[0].set_title("肺窗 (WL=-600, WW=1500)")
axes[1].imshow(apply_window(hu_image, wl=40, ww=400), cmap='gray')
axes[1].set_title("纵隔窗 (WL=40, WW=400)")
axes[2].imshow(apply_window(hu_image, wl=400, ww=1800), cmap='gray')
axes[2].set_title("骨窗 (WL=400, WW=1800)")
plt.tight_layout()
plt.savefig("windows_comparison.png", dpi=150)
plt.show()
🫁 实战案例一：CT 图像人体轮廓分割
问题定义
目标：给定一张 CT 图像，分割出人体所在的区域（排除床板和背景空气）。
这看起来简单，但实际有几个难点： 1. 背景（空气）和人体之间 HU 差异很大，阈值容易定 2. 但 CT 床板也有较高的 HU 值，会干扰结果 3. 人体内部可能有"空腔"（肺、胃），不应该被排除在外
思路分析
原始 CT 图
↓ 阈值分割（HU > -300 认为是物质）
二值图（前景/背景）
↓ Flood Fill（从四个角填充背景）
只保留从外部能到达的背景
↓ 取反 + 形态学填洞
人体轮廓 mask
↓ 去除小连通域（床板等噪声）
最终结果
💡 关键洞察：直接做阈值分割，空气是背景，但人体内的肺（HU ≈ -600）也会被误判为背景。用 Flood Fill 从图像四角开始填充，只有真正"连通到边缘"的区域才算外部背景，这样肺内部的空气就被正确保留了。
Flood Fill 原理
类似于图像软件中的"油漆桶"工具： 1. 从种子点（四个角）开始 2. 向四个方向扩展，直到遇到前景像素 3. 所有被访问到的像素标记为"外部背景"
from scipy import ndimage
from skimage import morphology, measure
import numpy as np
def flood_fill_background(binary_mask: np.ndarray) -> np.ndarray:
"""
从四个角 Flood Fill，标记外部背景。
scipy.ndimage.label 会找连通域，我们利用它来做 flood fill：
先在图像四周加一圈 0（确保四角连通），然后找到包含角点的连通域。
为什么不用 skimage.segmentation.flood_fill？
那个是单点出发，我们需要同时从四个角出发，这种方法更简洁。
"""
# 在四周添加一圈 0（背景），确保四角一定是背景
padded = np.pad(binary_mask, pad_width=1, mode='constant', constant_values=0)
# 对取反后的 mask 做连通域标记（0=背景，找连通的背景）
inverted = ~padded.astype(bool)
labeled, num_features = ndimage.label(inverted)
# 角点 (0,0) 所在的连通域就是"外部背景"
background_label = labeled[0, 0]
# 外部背景 mask（去掉 padding）
external_bg = (labeled == background_label)[1:-1, 1:-1]
return external_bg
def segment_body(hu_image: np.ndarray) -> np.ndarray:
"""
从 CT 图像中分割人体轮廓。
返回 bool 类型的 mask，True 表示人体区域。
Parameters
----------
hu_image : np.ndarray
HU 值图像，shape (H, W)
Returns
-------
body_mask : np.ndarray
人体轮廓 mask，shape (H, W)，dtype bool
"""
# Step 1: 阈值分割
# HU > -300：空气约 -1000，人体组织最低约 -600（肺），
# 取 -300 是个保守阈值，能保留大部分软组织
binary = hu_image > -300
# Step 2: 初步形态学操作——闭运算
# 闭运算 = 先膨胀后腐蚀，作用：填充小孔洞、连接断裂区域
# disk(5) 表示半径为 5 像素的圆形结构元素
selem = morphology.disk(5)
binary = morphology.binary_closing(binary, selem)
# Step 3: Flood Fill 找外部背景
external_bg = flood_fill_background(binary)
# Step 4: 人体 = 不是外部背景的区域
body_rough = ~external_bg
# Step 5: 填洞——人体内部的空腔（肺、胃等）应该包含在人体内
body_filled = ndimage.binary_fill_holes(body_rough)
# Step 6: 去除小连通域（床板碎片、噪声）
# remove_small_objects 会删除像素数少于 min_size 的连通域
body_clean = morphology.remove_small_objects(body_filled, min_size=10000)
return body_clean
完整可视化代码
def visualize_segmentation(hu_image: np.ndarray, body_mask: np.ndarray,
save_path: str = None):
"""可视化分割结果。"""
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 原始图像（肺窗显示）
display_img = apply_window(hu_image, wl=-600, ww=1500)
axes[0].imshow(display_img, cmap='gray')
axes[0].set_title("原始 CT 图像（肺窗）")
axes[0].axis('off')
# 分割 mask
axes[1].imshow(body_mask, cmap='gray')
axes[1].set_title("人体轮廓 Mask")
axes[1].axis('off')
# 叠加显示
axes[2].imshow(display_img, cmap='gray')
# 用红色轮廓线标出人体边界
contours = measure.find_contours(body_mask.astype(float), 0.5)
for contour in contours:
axes[2].plot(contour[:, 1], contour[:, 0], 'r-', linewidth=2)
axes[2].set_title("分割结果叠加")
axes[2].axis('off')
plt.tight_layout()
if save_path:
plt.savefig(save_path, dpi=150, bbox_inches='tight')
plt.show()
# ── 主流程 ──────────────────────────────────────────────
if __name__ == "__main__":
# 加载数据
hu_image, ds = load_dicom("CT_slice.dcm")
# 人体轮廓分割
body_mask = segment_body(hu_image)
# 可视化
visualize_segmentation