# 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），让学生看到阈值怎么影响结果