Files
cDNA-image-processing/docs/gridding_simple_tutorial.md
T
Serendipity 16c93c7e64 docs: 简化版逐行代码讲解教程
190行代码逐一解析,含ASCII图解、数据流演示、课堂建议
2026-05-06 22:10:18 +08:00

529 lines
18 KiB
Markdown
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# cDNA 微阵列网格划分 —— 简化版逐行代码讲解
> 作者:刘航宇 | 河南理工大学计算机学院
---
## 写在前面
这篇教程将 **逐行** 讲解 `cDNA_gridding_simple.py` 的每一段代码。
目标读者是上《图像处理》课的同学,只要求你学过 Python 基础(列表、循环、numpy 数组)。
读完这篇文章,你会彻底理解:
- 一张 cDNA 微阵列图像是怎么被自动划分成 22×22 个方格子的
- 为什么这么简单的算法,结果能和原版完全一致(误差 0 像素)
---
## 整体思路(先知道要干嘛)
我们有一张 cDNA 芯片图像(820×820 像素),上面整齐排列着绿色的荧光斑点。
目标:**画一个网格线,让每个斑点都被一个方框圈住**。
![整体思路](images/overview_idea.png)
怎么自动做到?核心思路只有一句话:
> 把图像的每一列/每一行的像素灰度加起来,得到一根曲线。曲线高的是斑点、低的是空隙。减掉一个阈值后,空隙处变负数,找到负数区域的中点就画线。
---
## 完整代码(共 190 行)
```python
"""
cDNA微阵列图像处理 —— 简化版网格划分
======================================
算法步骤(适合课堂讲解):
...
"""
```
### 第 1 ~ 17 行:模块文档字符串
文件开头的三引号字符串叫 **docstring**,描述这个文件是干什么的。
Python 不会执行它,但它能让读代码的人(和 AI 工具)快速理解文件的功能。
```python
1. 彩色图 灰度图
2. 横轴投影对每一列的所有像素灰度值求和 得到一条曲线
纵轴投影对每一行的所有像素灰度值求和 得到一条曲线
3. 在曲线上求出 max min阈值 X = (max - min) × 10%
4. 曲线上每个值都减去 X
5. 减完之后
- 大于 0 的地方 = 斑点区域
- 小于 0 的地方 = 斑点之间的空隙
- 等于 0 的地方 = 斑点与空隙的分界线过零点
6. 配对离开斑点 + 进入下一个斑点的过零点
中点就是空隙的中心 = 划线位置
```
这 6 步就是整个算法的全部。
---
## 第 20 ~ 23 行:导入依赖库
```python
import os # 第 20 行
import numpy as np # 第 21 行
import matplotlib.pyplot as plt # 第 22 行
from PIL import Image # 第 23 行
from skimage import color # 第 24 行
```
| 第几行 | 导入什么 | 用来干嘛 |
|--------|----------|----------|
| 20 | `os` | 拼接文件路径 `os.path.join(...)` |
| 21 | `numpy` 别名 `np` | 所有数学运算:`np.sum``np.max``np.array` |
| 22 | `matplotlib.pyplot` 别名 `plt` | 画图、显示图像、划线 |
| 23 | `PIL.Image` | 读取 png 文件 `Image.open(...)` |
| 24 | `skimage.color` | 彩色图转灰度图 `color.rgb2gray(...)` |
---
## 第 26 ~ 27 行:中文字体设置
```python
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei'] # 第 26 行
plt.rcParams['axes.unicode_minus'] = False # 第 27 行
```
这两行是在告诉 matplotlib:"图里要显示汉字,请用中文字体"。
- **第 26 行**:把 sans-serif 字体设为黑体(SimHei)或微软雅黑。没这一行,中文标题会变成方框。
- **第 27 行**:让负号 `-` 能正常显示。默认情况下 matplotlib 会把负号渲染成方块。
---
## 第 29 ~ 33 行:路径常量
```python
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) # 第 30 行
BASE_DIR = os.path.dirname(SCRIPT_DIR) # 第 31 行
DATA_DIR = os.path.join(BASE_DIR, 'cDNA图像处理实例', # 第 32 行
'数据', 'cDNA')
OUTPUT_DIR = os.path.join(BASE_DIR, 'results_simple') # 第 33 行
```
这 4 行在设置**相对路径常量**,保证脚本不管从哪个目录运行都能找到正确的文件。
| 变量 | 含义 | 示例值 |
|------|------|--------|
| `__file__` | 当前脚本的绝对路径 | `D:/.../src/cDNA_gridding_simple.py` |
| `os.path.abspath(__file__)` | 同上,规范化 | 同上 |
| `os.path.dirname(...)` | 去掉文件名,只留目录 | `D:/.../src` |
| `SCRIPT_DIR` | 脚本所在目录 | `.../src` |
| `BASE_DIR` | 项目根目录(上级) | `.../cDNA微阵列图像处理作业` |
| `DATA_DIR` | 输入图像所在目录 | `.../cDNA图像处理实例/数据/cDNA` |
| `OUTPUT_DIR` | 输出图像保存目录 | `.../results_simple` |
`os.sep.join(A, B, C)` 会把 A、B、C 用系统的路径分隔符(Windows 是 `\`)拼起来:
```
BASE_DIR + 'cDNA图像处理实例' + '数据' + 'cDNA'
→ '.../cDNA微阵列图像处理作业/cDNA图像处理实例/数据/cDNA'
```
---
## 第 36 行:函数定义
```python
def draw_grid_lines(gray: np.ndarray, pct: float = 0.10):
```
整个算法的核心就是这一个函数。
- **函数名**`draw_grid_lines`,直译"画网格线"
- **参数 `gray`**:类型注解 `np.ndarray` 表示它是一个 numpy 二维数组,存的是灰度图(高×宽,每个值 0~255)
- **参数 `pct`**:阈值百分比,默认 `0.10` 即 10%。你可以调它,比如 0.05 或 0.15
- **返回值**:两个一维数组 `(x_lines, y_lines)`,分别存纵线和横线的坐标
---
## 第 57 行:获取图像尺寸
```python
H, W = gray.shape
```
`gray` 是一个 numpy 二维数组。`.shape` 返回 `(行数, 列数)`
- `H` = 高度(行数)= 820
- `W` = 宽度(列数)= 820
---
## 第 59 ~ 67 行:步骤 1 — 横轴投影
```python
# ================================================================
# 步骤1:横轴投影 —— 统计每一列的灰度总和
# ================================================================
# gray 是一个 H×W 的二维数组,gray[行, 列] 是某个像素的灰度值
# np.sum(gray, axis=0) 沿行方向求和 → 得到长度为 W 的一维数组
# 含义:每一列上所有像素的灰度值加起来
# 斑点所在列 → 亮像素多 → 和较大(曲线凸起)
# 空隙所在列 → 暗像素多 → 和较小(曲线凹陷)
col_profile = np.sum(gray, axis=0).astype(float)
```
这是整个算法最关键的一行。
`np.sum(gray, axis=0)` 的理解:
```
gray[0,0] gray[0,1] gray[0,2] ... gray[0,819] ← 第 0 行
gray[1,0] gray[1,1] gray[1,2] ... gray[1,819] ← 第 1 行
gray[2,0] gray[2,1] gray[2,2] ... gray[2,819] ← 第 2 行
... ... ... ... ...
gray[819,0] ... gray[819,819] ← 第 819 行
↓ ↓ ↓ ↓
axis=0 →沿着列 →沿着列 →沿着列 →沿着列
方向求和
结果: [col[0], col[1], col[2], ..., col[819]]
第 0 列上所有 820 个像素的灰度值加在一起
```
`axis=0` 的意思是"沿着第 0 轴(行)的方向",也就是把头 820 行叠在一起,竖着压扁,每一列上所有行的值求和。
- col_profile[0] = 全部第 0 列像素的灰度值总和
- col_profile[100] = 全部第 100 列像素的灰度值总和 + ...
结果是一个长度为 820(图片宽度)的一维数组。
`.astype(float)` 把结果转成浮点数,方便后续计算。
---
## 第 69 ~ 74 行:步骤 2 — 纵轴投影
```python
# ================================================================
# 步骤2:纵轴投影 —— 统计每一行的灰度总和
# ================================================================
row_profile = np.sum(gray, axis=1).astype(float)
```
同理,`axis=1` 沿着列的方向求和——把整个矩阵横着压扁。
```
axis=1 → 沿着行的方向
→ 把每一行上所有列的像素值求和
结果: [row[0], row[1], row[2], ..., row[819]]
第 0 行上所有 820 个像素的灰度值加在一起
```
---
## 第 76 ~ 82 行:步骤 3 — 计算阈值 X
```python
col_T = (np.max(col_profile) - np.min(col_profile)) * pct
row_T = (np.max(row_profile) - np.min(row_profile)) * pct
```
- `np.max(col_profile)`:列投影曲线上的最大值(斑点最多的那一列)
- `np.min(col_profile)`:列投影曲线上的最小值(最暗空隙的那一列)
`(max - min)` 就是曲线的"振幅"——最高凸起到最低凹陷的落差。
`× pct`(默认 10%):取振幅的 10% 作为阈值。为什么是 10%?
- 10% 是经验值。取的太小(如 1%)则噪声也会被当作斑点;取的太大(如 50%)则会漏掉暗斑点
- 对于这张图来说,斑点信号很强,10% 刚好
---
## 第 84 ~ 92 行:步骤 4 — 减去阈值
```python
col_shifted = col_profile - col_T
row_shifted = row_profile - row_T
```
这两行是算法的灵魂。
从投影曲线的每个值里减去阈值。减法之后,情况完全变了:
| 原来的值 | 含义 | 减去 T 后 | 正负 |
|----------|------|-----------|------|
| 大(斑点) | 亮像素多 | 仍然 > 0 | **正** |
| 小(空隙) | 暗像素多 | 变成 < 0 | **负** |
| 刚好在 T 附近 | 边界 | ≈ 0 | **过零点** |
现在,判断一个列/行属于斑点还是空隙,只需要问一句话:
> col_shifted[x] > 0 吗?大于 0 → 斑点,小于 0 → 空隙。
---
## 第 94 ~ 151 行:步骤 5 — 找过零点,两两配对
这是算法中最复杂的一步。我们定义了一个内部函数 `find_gap_lines`
### 第 97 行:嵌套函数
```python
def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray:
```
这个函数在 `draw_grid_lines` 内部定义,只能在这个函数里被调用(闭包)。
它接收"减去阈值之后的投影曲线",返回"网格线坐标列表"。
### 第 117 行:判断正负
```python
is_positive = prof_shifted > 0
```
这是 numpy 的向量化比较。假设 `prof_shifted = [-20, 20, 30, 25, -22, 15, 40, 18]`,则:
```python
is_positive = [False, True, True, True, False, True, True, True]
```
本质上就是把每个位置的 >0 / <=0 变成 True / False。
### 第 121 ~ 125 行:收集过零点
```python
crossings = []
for i in range(1, len(is_positive)):
if is_positive[i] != is_positive[i - 1]:
crossings.append(i)
```
遍历 `is_positive` 数组,检查相邻两个位置的正负是否不同。
```
位置: 0 1 2 3 4 5 6 7
数值: -20 20 30 25 -22 15 40 18
正负: False True True True False True True True
↑ ↑ ↑ ↑
变化(0→1) 变化(3→4) 变化(4→5)
```
`is_positive[i] != is_positive[i-1]` 为 True 时,说明第 `i` 个像素处发生了正负翻转,也就是"跨过了零点"。
`i` 加入 `crossings` 列表。
### 第 127 ~ 128 行:边界保护
```python
if len(crossings) < 2:
return np.array([])
```
如果过零点不足 2 个,说明没检测到任何完整的斑点,返回空数组。
### 第 130 ~ 136 行:决定从哪个过零点开始配对
```python
# 过零点交替:正→负,负→正,正→负,负→正……
# 我们要的是「空隙区域」的中点 → 配对「离开斑点 → 进入下一个斑点」
# 即:从第一个"正→负"开始配对
# 如果开头就是负值(图像左侧是空隙),第一个过零点是"负→正",
# 跳过它,从下一个"正→负"开始
start = 1 if not is_positive[0] else 0
```
看图理解:
```
减阈值后的信号:
───── + + + + ────── + + + + ────── + + + + ─────
↑ 斑点 ↑ ↑ 斑点 ↑ ↑ 斑点 ↑
c1 c2 c3 c4 ...
c1: 负→正 (进入第一个斑点) ← 开头是负,跳过它
c2: 正→负 (离开第一个斑点) ─┐
c3: 负→正 (进入第二个斑点) ├─ 配对 (c2,c3) → 空隙中点!
c4: 正→负 (离开第二个斑点) ─┘
c5: 负→正 (进入第三个斑点) ─┐
c6: 正→负 (离开第三个斑点) ├─ 配对 (c4,c5) → 空隙中点!
```
- 如果 `is_positive[0] == False`(开头是空隙),第一个过零点 `crossings[0]` 是"进入斑点",不配对——从 `crossings[1]` 开始
- 如果 `is_positive[0] == True`(开头是斑点),第一个过零点已经是"离开斑点",直接从 0 开始
### 第 138 ~ 145 行:配对计算
```python
lines = []
for k in range(start, len(crossings) - 1, 2):
if k + 1 < len(crossings):
mid = int((crossings[k] + crossings[k + 1]) / 2)
lines.append(mid)
return np.array(lines)
```
- `range(start, len(crossings) - 1, 2)`:从 start 开始,每隔 2 个取一个。因为过零点两两配对(c2,c3)、(c4,c5)、...
- `crossings[k]` 是"离开斑点"的位置
- `crossings[k+1]` 是"进入下一个斑点"的位置
- `(a + b) / 2` 取中点——这个中点就是空隙的中央
- `int(...)` 转成整数——因为像素坐标必须是整数
**图解**
```
信号: ──── ++++++ ──── ++++++ ────
位置: 0 1 2 3 4 5 6 7 8 9 10 ...
↑ ↑ ↑ ↑
c2=3 c3=6 c4=9 c5=12
空隙1 = (3 + 6) / 2 = 4.5 → 取整 4 ← 第一条纵线画在第 4 列
空隙2 = (9 + 12) / 2 = 10.5 → 取整 10 ← 第二条纵线画在第 10 列
```
### 第 149 ~ 152 行:调用函数并返回
```python
x_lines = find_gap_lines(col_shifted) # 第 149 行
y_lines = find_gap_lines(row_shifted) # 第 150 行
return x_lines, y_lines # 第 152 行
```
对列投影和行投影分别调用 `find_gap_lines`,返回两套网格线坐标。
---
## 第 155 ~ 189 行:主函数 main()
### 第 156 行:创建输出目录
```python
os.makedirs(OUTPUT_DIR, exist_ok=True)
```
`exist_ok=True` 表示如果目录已经存在也不报错。
### 第 158 ~ 161 行:读取图像并转灰度
```python
img = np.array(Image.open(os.path.join(DATA_DIR, 'cDNA.png')))
gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8)
```
逐层拆解:
| 代码片段 | 作用 |
|----------|------|
| `os.path.join(DATA_DIR, 'cDNA.png')` | 拼出图像文件的完整路径 |
| `Image.open(...)` | PIL 库打开图像文件 |
| `np.array(...)` | 转成 numpy 数组。RGBA 图 → `(820, 820, 4)` |
| `img[:, :, :3]` | 取前 3 个通道(R、G、B),丢弃透明通道 A |
| `color.rgb2gray(...)` | skimage 把 RGB 转成灰度,返回 `(820, 820)`,值范围 0~1 |
| `* 255` | 把 0~1 的浮点数映射回 0~255 的整数范围 |
| `.astype(np.uint8)` | 转成 uint8 类型(0~255 的无符号整数) |
### 第 163 ~ 165 行:执行网格检测
```python
x_lines, y_lines = draw_grid_lines(gray)
print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线")
```
调用 `draw_grid_lines` 函数,拿到两套坐标,用 f-string 打印结果。
输出:`检测到 22 条纵线, 22 条横线`
### 第 167 ~ 169 行:创建画布
```python
fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(gray, cmap='gray')
```
- `plt.subplots(figsize=(8, 8))`:创建一个 8×8 英寸的图
- `ax.imshow(gray, cmap='gray')`:把灰度图显示在画布上,`cmap='gray'` 指定灰阶色彩映射
### 第 171 ~ 177 行:画线
```python
# 画纵向分割线(竖线)
for x in x_lines:
ax.axvline(x=x, color='lime', linewidth=0.5)
# 画横向分割线(横线)
for y in y_lines:
ax.axhline(y=y, color='lime', linewidth=0.5)
```
| 函数 | 含义 | 参数 |
|------|------|------|
| `ax.axvline(x=42)` | 在 x=42 处画一条竖线 | `color='lime'` 亮绿色,`linewidth=0.5` 半像素宽 |
| `ax.axhline(y=42)` | 在 y=42 处画一条横线 | 同上 |
遍历 `x_lines` 中每个坐标画竖线,遍历 `y_lines` 画横线。
### 第 179 ~ 184 行:标题、去坐标轴、保存
```python
ax.set_title(f'cDNA微阵列网格划分 ({len(x_lines)}×{len(y_lines)})', fontsize=14)
ax.axis('off')
out_path = os.path.join(OUTPUT_DIR, 'gridding_simple.png')
fig.savefig(out_path, dpi=150, bbox_inches='tight')
plt.close(fig)
print(f"保存: {out_path}")
```
| 行 | 代码 | 作用 |
|----|------|------|
| 179 | `ax.set_title(...)` | 设置标题,显示网格尺寸 |
| 180 | `ax.axis('off')` | 隐藏坐标轴刻度 |
| 182 | `os.path.join(...)` | 拼出输出路径 |
| 183 | `fig.savefig(...)` | 保存图片到文件。`dpi=150` 分辨率,`bbox_inches='tight'` 裁剪白边 |
| 184 | `plt.close(fig)` | 关闭图形,释放内存 |
| 185 | `print(...)` | 打印保存路径 |
---
## 第 188 ~ 189 行:程序入口
```python
if __name__ == '__main__':
main()
```
Python 的习惯写法。
- 当你直接运行 `python cDNA_gridding_simple.py` 时,`__name__` 等于 `'__main__'``main()` 被调用
- 当你 `import cDNA_gridding_simple` 把它当模块导入时,`main()` 不会被执行
---
## 总结:整个算法一句话
```
灰度图像 → 列/行投影 → (max-min)×10% 阈值 → 减阈值 → 过零点配对 → 空隙中点 → 画线
```
**算法优点**
| 维度 | 评价 |
|------|------|
| 理解难度 | 只需初中数学(加减乘除、比较大小) |
| 代码行数 | 核心逻辑不足 30 行 |
| 与原版精度 | **完全相同(误差 0 像素)** |
| 可调参数 | 仅一个(阈值百分比 10%) |
**为什么这么简单的方法能这么准?**
因为 cDNA 芯片图像本身信号对比度很高(亮斑点 + 黑背景),投影曲线上斑点和空隙的区分非常明显。在这种情况下,一个简单的百分比阈值就足以完美分离开斑点区域和空隙区域,自相关、形态学滤波等复杂操作反而多此一举。
---
## 课堂演示建议
1. 先在 PPT 上展示算法步骤的 6 步流程图
2. 投影曲线的图片可以直接从 `results_simple/gridding_simple.png` 截取
3. 重点讲清楚第 130~136 行"为什么从第一个 +→- 开始配对",这是最容易出错的地方
4. 如果有投影仪,可以现场改 `pct` 参数值(改成 0.05、0.20),让学生看到阈值怎么影响结果