19 KiB
19 KiB
cDNA 微阵列网格划分 —— 简化版逐行代码讲解
作者：刘航宇 | 河南理工大学计算机学院
写在前面
这篇教程将 逐行 讲解 cDNA_gridding_simple.py 的每一段代码。
目标读者是上《图像处理》课的同学，只要求你学过 Python 基础（列表、循环、numpy 数组）。
读完这篇文章，你会彻底理解：
- 一张 cDNA 微阵列图像是怎么被自动划分成 22×22 个方格子的
- 为什么这么简单的算法，结果能和原版完全一致（误差 0 像素）
整体思路（先知道要干嘛）
我们有一张 cDNA 芯片图像（820×820 像素），上面整齐排列着绿色的荧光斑点。
目标：画一个网格线，让每个斑点都被一个方框圈住。
怎么自动做到？核心思路只有一句话：
把图像的每一列/每一行的像素灰度加起来，得到一根曲线。曲线高的是斑点、低的是空隙。减掉一个阈值后，空隙处变负数，找到负数区域的中点就画线。
完整代码（共 190 行）
"""
cDNA微阵列图像处理 —— 简化版网格划分
======================================
算法步骤（适合课堂讲解）：
...
"""
第 1 ~ 17 行：模块文档字符串
文件开头的三引号字符串叫 docstring，描述这个文件是干什么的。 Python 不会执行它，但它能让读代码的人（和 AI 工具）快速理解文件的功能。
1. 彩色图 → 灰度图
2. 横轴投影：对每一列的所有像素灰度值求和 → 得到一条曲线
纵轴投影：对每一行的所有像素灰度值求和 → 得到一条曲线
3. 在曲线上，求出 max 和 min，阈值 X = (max - min) × 10%
4. 曲线上每个值都减去 X
5. 减完之后：
- 大于 0 的地方 = 斑点区域
- 小于 0 的地方 = 斑点之间的空隙
- 等于 0 的地方 = 斑点与空隙的分界线（过零点）
6. 配对「离开斑点 + 进入下一个斑点」的过零点，
中点就是空隙的中心 = 划线位置
这 6 步就是整个算法的全部。
第 20 ~ 23 行：导入依赖库
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 行：中文字体设置
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 行：路径常量
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 行：函数定义
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 行：获取图像尺寸
H, W = gray.shape
gray 是一个 numpy 二维数组。.shape 返回 (行数, 列数)。
H= 高度（行数）= 820W= 宽度（列数）= 820
第 59 ~ 67 行：步骤 1 — 横轴投影
# ================================================================
# 步骤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 — 纵轴投影
# ================================================================
# 步骤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
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 — 减去阈值
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 行：嵌套函数
def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray:
这个函数在 draw_grid_lines 内部定义，只能在这个函数里被调用（闭包）。
它接收"减去阈值之后的投影曲线"，返回"网格线坐标列表"。
第 117 行：判断正负
is_positive = prof_shifted > 0
这是 numpy 的向量化比较。假设 prof_shifted = [-20, 20, 30, 25, -22, 15, 40, 18]，则：
is_positive = [False, True, True, True, False, True, True, True]
本质上就是把每个位置的 >0 / <=0 变成 True / False。
第 121 ~ 125 行：收集过零点
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 行：边界保护
if len(crossings) < 2:
return np.array([])
如果过零点不足 2 个，说明没检测到任何完整的斑点，返回空数组。
第 130 ~ 136 行：决定从哪个过零点开始配对
# 过零点交替：正→负，负→正，正→负，负→正……
# 我们要的是「空隙区域」的中点 → 配对「离开斑点 → 进入下一个斑点」
# 即：从第一个"正→负"开始配对
# 如果开头就是负值（图像左侧是空隙），第一个过零点是"负→正"，
# 跳过它，从下一个"正→负"开始
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 行：配对计算
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 行：调用函数并返回
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 行：创建输出目录
os.makedirs(OUTPUT_DIR, exist_ok=True)
exist_ok=True 表示如果目录已经存在也不报错。
第 158 ~ 161 行：读取图像并转灰度
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 |
.astype(np.uint8) |
转成 uint8 类型（0~255 的无符号整数） |
第 163 ~ 165 行：执行网格检测
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 行：创建画布
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 行：画线
# 画纵向分割线（竖线）
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 行：标题、去坐标轴、保存
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 行：程序入口
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 芯片图像本身信号对比度很高（亮斑点 + 黑背景），投影曲线上斑点和空隙的区分非常明显。在这种情况下，一个简单的百分比阈值就足以完美分离开斑点区域和空隙区域，自相关、形态学滤波等复杂操作反而多此一举。
课堂演示建议
- 先在 PPT 上展示算法步骤的 6 步流程图
- 投影曲线的图片可以直接从
results_simple/gridding_simple.png截取 - 重点讲清楚第 130~136 行"为什么从第一个 +→- 开始配对"，这是最容易出错的地方
- 如果有投影仪，可以现场改
pct参数值（改成 0.05、0.20），让学生看到阈值怎么影响结果