refactor: 重写注释,统一风格
- 顶部docstring改为算法流程总览+各步骤详解 - 两个函数各配职责明确的注释 - 主流程三个步骤注释简洁
This commit is contained in:
+104
-112
@@ -4,29 +4,44 @@ cDNA微阵列图像处理 —— 简化版
|
|||||||
|
|
||||||
D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py
|
D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py
|
||||||
|
|
||||||
算法步骤(划线):
|
一、算法流程总览
|
||||||
先用 Otsu 算出像素级最佳阈值 T,百分比 = T / 255 来自数据而非写死
|
|
||||||
|
|
||||||
1. 彩色图 → 灰度图
|
灰度图 ──→ Otsu求像素最佳阈值 T ──→ 百分比 = T/255(自适应,不写死)
|
||||||
2. 横轴投影:对每一列的所有像素灰度值求和 → 得到一条曲线
|
│
|
||||||
纵轴投影:对每一行的所有像素灰度值求和 → 得到一条曲线
|
├─→ 横轴投影/纵轴投影 ──→ X = (max-min) × 百分比 ──→ 减阈值 ──→
|
||||||
3. 在曲线上,求出 max 和 min,阈值 X = (max - min) × (T / 255)
|
│ 过零点配对 ──→ 空隙中点 ──→ 网格线
|
||||||
4. 曲线上每个值都减去 X
|
│
|
||||||
5. 减完之后:
|
└─→ gray > T ──→ 二值图(分割结果)
|
||||||
- 大于 0 的地方 = 斑点区域
|
|
||||||
- 小于 0 的地方 = 斑点之间的空隙
|
|
||||||
- 等于 0 的地方 = 斑点与空隙的分界线(过零点)
|
|
||||||
6. 配对相邻的过零点(离开斑点 + 进入下一个斑点),
|
|
||||||
中点就是空隙的中心 = 划线位置
|
|
||||||
|
|
||||||
算法步骤(Otsu阈值分割):
|
二、各步骤详解
|
||||||
|
|
||||||
1. 遍历所有可能的 T (0~255)
|
1. 彩色图 → 灰度图
|
||||||
2. 计算前景(>T)的方差 + 背景(≤T)的方差 → 类内方差
|
|
||||||
3. 计算前景均值 vs 背景均值的差距 → 类间方差
|
|
||||||
4. 选 T 使 类内方差最小 / 类间方差最大
|
|
||||||
5. 相当于在灰度直方图上找一个"谷底",两座山之间的最低处就是最佳分割线。
|
|
||||||
|
|
||||||
|
2. Otsu 自动阈值
|
||||||
|
遍历灰度 0~255,每个候选 T 将像素分为前景(>T)和背景(≤T),
|
||||||
|
计算类内方差 w_bg×σ²_bg + w_fg×σ²_fg,
|
||||||
|
选使类内方差最小的 T 作为最佳分割线。
|
||||||
|
百分比 = T / 255,取代原来的固定 10%。
|
||||||
|
|
||||||
|
3. 投影
|
||||||
|
横轴:np.sum(每列) → 长度=宽度的曲线,高点=斑点列,低点=空隙列
|
||||||
|
纵轴:np.sum(每行) → 长度=高度的曲线,高点=斑点行,低点=空隙行
|
||||||
|
|
||||||
|
4. 阈值 X = (max-min) × 百分比
|
||||||
|
|
||||||
|
5. 曲线减 X
|
||||||
|
- 大于 0 = 斑点区域
|
||||||
|
- 小于 0 = 斑点之间的空隙
|
||||||
|
- 等于 0 处 = 过零点 = 斑点和空隙的分界线
|
||||||
|
|
||||||
|
6. 过零点配对
|
||||||
|
过零点交替出现:正→负(离开斑点)、负→正(进入下一斑点)
|
||||||
|
配对「离开斑点 + 进入下一斑点」,中点 = 空隙中央 = 划线位置
|
||||||
|
|
||||||
|
7. 分割
|
||||||
|
gray > T 为前景(斑点),≤T 为背景
|
||||||
|
|
||||||
|
8. 输出左右对比图:左=网格划线,右=Otsu分割
|
||||||
"""
|
"""
|
||||||
|
|
||||||
import os
|
import os
|
||||||
@@ -38,143 +53,116 @@ from skimage import color
|
|||||||
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
|
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']
|
||||||
plt.rcParams['axes.unicode_minus'] = False
|
plt.rcParams['axes.unicode_minus'] = False
|
||||||
|
|
||||||
# 路径设置
|
# 路径设置(从脚本位置动态推导,禁止硬编码绝对路径)
|
||||||
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
|
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
|
||||||
BASE_DIR = os.path.dirname(SCRIPT_DIR)
|
BASE_DIR = os.path.dirname(SCRIPT_DIR)
|
||||||
DATA_DIR = os.path.join(BASE_DIR, 'cDNA图像处理实例', '数据', 'cDNA')
|
DATA_DIR = os.path.join(BASE_DIR, 'cDNA图像处理实例', '数据', 'cDNA')
|
||||||
OUTPUT_DIR = os.path.join(BASE_DIR, 'results_simple')
|
OUTPUT_DIR = os.path.join(BASE_DIR, 'results_simple')
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 函数1:Otsu 像素级阈值
|
||||||
|
# ================================================================
|
||||||
def otsu_threshold_pixels(gray: np.ndarray) -> float:
|
def otsu_threshold_pixels(gray: np.ndarray) -> float:
|
||||||
"""
|
"""
|
||||||
Otsu 自动阈值——在像素级找到最佳分割灰度值 T。
|
对图像像素做 Otsu 自动阈值检测。
|
||||||
|
|
||||||
遍历 0~255,对每个候选 T 计算类内方差,选最小的。
|
遍历灰度值 0~255,找到使"类内方差"最小的阈值 T。
|
||||||
返回 T / 255,作为投影曲线阈值的自适应百分比。
|
类内方差 = w_bg × σ²_bg + w_fg × σ²_fg
|
||||||
|
(背景权重 × 背景方差 + 前景权重 × 前景方差)
|
||||||
|
|
||||||
|
返回 T/255,即自适应百分比,供投影曲线使用。
|
||||||
"""
|
"""
|
||||||
best_T = 0
|
best_T = 0 # 当前最佳阈值
|
||||||
best_cost = float('inf')
|
best_cost = float('inf') # 当前最小类内方差
|
||||||
total = gray.size
|
total = gray.size # 总像素数
|
||||||
|
|
||||||
for T in range(1, 255):
|
for T in range(1, 255):
|
||||||
bg = gray[gray <= T]
|
# 将像素按 T 分为两组
|
||||||
fg = gray[gray > T]
|
bg = gray[gray <= T] # 背景
|
||||||
w_bg = len(bg) / total
|
fg = gray[gray > T] # 前景
|
||||||
w_fg = len(fg) / total
|
w_bg = len(bg) / total # 背景权重
|
||||||
|
w_fg = len(fg) / total # 前景权重
|
||||||
|
|
||||||
if w_bg == 0 or w_fg == 0:
|
if w_bg == 0 or w_fg == 0:
|
||||||
continue
|
continue # 某组为空,跳过
|
||||||
|
|
||||||
|
# 类内方差 = 加权平均方差
|
||||||
cost = w_bg * np.var(bg) + w_fg * np.var(fg)
|
cost = w_bg * np.var(bg) + w_fg * np.var(fg)
|
||||||
|
|
||||||
if cost < best_cost:
|
if cost < best_cost:
|
||||||
best_cost = cost
|
best_cost = cost
|
||||||
best_T = T
|
best_T = T
|
||||||
return best_T / 255.0
|
|
||||||
|
return best_T / 255.0 # 归一化为百分比
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 函数2:网格划线
|
||||||
|
# ================================================================
|
||||||
def draw_grid_lines(gray: np.ndarray):
|
def draw_grid_lines(gray: np.ndarray):
|
||||||
"""
|
"""
|
||||||
检测网格分割线。先对像素做 Otsu,用 T/255 作为自适百分比。
|
检测网格分割线。
|
||||||
|
|
||||||
原理:
|
先用 Otsu 算出自适应百分比,
|
||||||
灰度图的每一列/行,属于斑点的像素灰度值高,属于背景的灰度值低。
|
再对列投影和行投影分别处理:
|
||||||
把每列/行的灰度值加起来,就能得到一条曲线:
|
投影 → 减阈值 → 过零点配对 → 空隙中点 = 网格线
|
||||||
——曲线凸起的地方 = 斑点所在位置
|
|
||||||
——曲线凹陷的地方 = 斑点之间的空隙
|
|
||||||
|
|
||||||
去掉一个阈值后,曲线在空隙处会变成负数,
|
返回 (纵线x列表, 横线y列表, 自适应百分比)。
|
||||||
过零点的位置就是斑点和空隙的分界线,
|
|
||||||
两个分界线中点就是划线位置。
|
|
||||||
|
|
||||||
返回:
|
|
||||||
(纵线x坐标列表, 横线y坐标列表, 自适应百分比)
|
|
||||||
"""
|
"""
|
||||||
# 0. 先用 Otsu 算出自适应百分比
|
# ---- 0. 自适应百分比 ----
|
||||||
pct = otsu_threshold_pixels(gray)
|
pct = otsu_threshold_pixels(gray)
|
||||||
|
|
||||||
H, W = gray.shape
|
H, W = gray.shape
|
||||||
|
|
||||||
# ================================================================
|
# ---- 1. 横轴投影(列方向)----
|
||||||
# 步骤1:横轴投影 —— 统计每一列的灰度总和
|
# 对每一列上所有行的像素灰度求和 → 长度为 W 的数组
|
||||||
# ================================================================
|
|
||||||
# gray 是一个 H×W 的二维数组,gray[行, 列] 是某个像素的灰度值
|
|
||||||
# np.sum(gray, axis=0) 沿行方向求和 → 得到长度为 W 的一维数组
|
|
||||||
# 含义:每一列上所有像素的灰度值加起来
|
|
||||||
# 斑点所在列 → 亮像素多 → 和较大(曲线凸起)
|
|
||||||
# 空隙所在列 → 暗像素多 → 和较小(曲线凹陷)
|
|
||||||
col_profile = np.sum(gray, axis=0).astype(float)
|
col_profile = np.sum(gray, axis=0).astype(float)
|
||||||
|
|
||||||
# ================================================================
|
# ---- 2. 纵轴投影(行方向)----
|
||||||
# 步骤2:纵轴投影 —— 统计每一行的灰度总和
|
# 对每一行上所有列的像素灰度求和 → 长度为 H 的数组
|
||||||
# ================================================================
|
|
||||||
# np.sum(gray, axis=1) 沿列方向求和 → 得到长度为 H 的一维数组
|
|
||||||
# 含义:每一行上所有像素的灰度值加起来
|
|
||||||
row_profile = np.sum(gray, axis=1).astype(float)
|
row_profile = np.sum(gray, axis=1).astype(float)
|
||||||
|
|
||||||
# ================================================================
|
# ---- 3. 投影阈值 ----
|
||||||
# 步骤3:计算阈值 X = (max - min) × 自适应百分比
|
|
||||||
# ================================================================
|
|
||||||
# 百分比来自 Otsu(像素级最佳分割线),不是写死的 10%
|
|
||||||
# max-min 是曲线的"振幅"
|
|
||||||
col_T = (np.max(col_profile) - np.min(col_profile)) * pct
|
col_T = (np.max(col_profile) - np.min(col_profile)) * pct
|
||||||
row_T = (np.max(row_profile) - np.min(row_profile)) * pct
|
row_T = (np.max(row_profile) - np.min(row_profile)) * pct
|
||||||
|
|
||||||
# ================================================================
|
# ---- 4. 曲线减去阈值 ----
|
||||||
# 步骤4:曲线上所有值减去阈值
|
# 减完:正 = 斑点区域,负 = 空隙
|
||||||
# ================================================================
|
|
||||||
# 减去阈值后:
|
|
||||||
# 原本在空隙处的值(本来就小)→ 变成负数
|
|
||||||
# 原本在斑点处的值(本来就大)→ 仍然为正数
|
|
||||||
# 等于0的位置 = 斑点与空隙的分界线 = 过零点
|
|
||||||
col_shifted = col_profile - col_T
|
col_shifted = col_profile - col_T
|
||||||
row_shifted = row_profile - row_T
|
row_shifted = row_profile - row_T
|
||||||
|
|
||||||
# ================================================================
|
# ---- 5. 过零点配对 → 空隙中线 ----
|
||||||
# 步骤5:找过零点,两两配对,中间点即划线位置
|
|
||||||
# ================================================================
|
|
||||||
def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray:
|
def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray:
|
||||||
"""
|
"""
|
||||||
在减去阈值后的投影曲线上,找到每两个斑点之间的空隙中线。
|
在减去阈值后的曲线上,找到空隙中线 = 网格线位置。
|
||||||
|
|
||||||
图解说明(设阈值为30):
|
原理图解:
|
||||||
原始曲线: ___/‾‾‾\___/‾‾‾‾\___/‾‾‾\___
|
信号: ----++++----++++----++++
|
||||||
数值: 10 50 60 55 8 45 70 48 5 55 50 12
|
|
||||||
减阈值后: -20 20 30 25 -22 15 40 18 -25 25 20 -18
|
|
||||||
正负: ─ + + + ─ + + + ─ + + ─
|
|
||||||
↑ ↑ ↑ ↑
|
↑ ↑ ↑ ↑
|
||||||
过零点 过零点 过零点 过零点
|
过零点配对:第1个+→- 与 第1个-→+
|
||||||
|
→ 中点 = 空隙中央 = 划线位置
|
||||||
过零点之间的区域:
|
|
||||||
第一个 + 区域 = 一个斑点 (过零点1 → 过零点2)
|
|
||||||
第一个 - 区域 = 空隙 (过零点2 → 过零点3) ← 我们要的!
|
|
||||||
第二个 + 区域 = 下一个斑点 (过零点3 → 过零点4)
|
|
||||||
|
|
||||||
划线位置 = 空隙(负数区域)的中点
|
|
||||||
= (离开斑点的过零点 + 进入下一个斑点的过零点) / 2
|
|
||||||
"""
|
"""
|
||||||
# 判断每个位置是正(斑点)还是负(空隙)
|
# 每个位置是正(斑点)还是负(空隙)
|
||||||
is_positive = prof_shifted > 0
|
is_positive = prof_shifted > 0
|
||||||
|
|
||||||
# 收集所有符号变化的位置(过零点)
|
# 收集符号变化处(过零点)
|
||||||
crossings = []
|
crossings = [] # 存过零点的位置
|
||||||
for i in range(1, len(is_positive)):
|
for i in range(1, len(is_positive)):
|
||||||
# 如果当前和前一个正负不同 → 发生了跨越零点
|
if is_positive[i] != is_positive[i - 1]: # 符号变化
|
||||||
if is_positive[i] != is_positive[i - 1]:
|
|
||||||
crossings.append(i)
|
crossings.append(i)
|
||||||
|
|
||||||
if len(crossings) < 2:
|
if len(crossings) < 2: # 过零点不足 → 放弃
|
||||||
return np.array([])
|
return np.array([])
|
||||||
|
|
||||||
# 过零点交替:正→负,负→正,正→负,负→正……
|
# 过零点交替:正→负,负→正,正→负,负→正 ...
|
||||||
# 我们要的是「空隙区域」的中点 → 配对「离开斑点 → 进入下一个斑点」
|
# 要配对的是「离开斑点(正→负)」+「进入下一斑点(负→正)」
|
||||||
# 即:从第一个"正→负"开始配对
|
# 如果信号开头是负,跳过第一个 crossing
|
||||||
|
|
||||||
# 如果开头就是负值(图像左侧是空隙),第一个过零点是"负→正",
|
|
||||||
# 跳过它,从下一个"正→负"开始
|
|
||||||
start = 1 if not is_positive[0] else 0
|
start = 1 if not is_positive[0] else 0
|
||||||
|
|
||||||
lines = []
|
lines = []
|
||||||
for k in range(start, len(crossings) - 1, 2):
|
for k in range(start, len(crossings) - 1, 2):
|
||||||
if k + 1 < len(crossings):
|
if k + 1 < len(crossings):
|
||||||
# crossings[k]: 正→负(离开斑点)
|
|
||||||
# crossings[k+1]: 负→正(进入下一个斑点)
|
|
||||||
# 中点 = 空隙中央 = 划线位置
|
|
||||||
mid = int((crossings[k] + crossings[k + 1]) / 2)
|
mid = int((crossings[k] + crossings[k + 1]) / 2)
|
||||||
lines.append(mid)
|
lines.append(mid)
|
||||||
|
|
||||||
@@ -186,27 +174,30 @@ def draw_grid_lines(gray: np.ndarray):
|
|||||||
return x_lines, y_lines, pct
|
return x_lines, y_lines, pct
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 主流程
|
||||||
|
# ================================================================
|
||||||
def main():
|
def main():
|
||||||
os.makedirs(OUTPUT_DIR, exist_ok=True)
|
os.makedirs(OUTPUT_DIR, exist_ok=True)
|
||||||
|
|
||||||
# ---- 读取图像并转为灰度图 ----
|
# ---- 读取图像,转为灰度 ----
|
||||||
img = np.array(Image.open(os.path.join(DATA_DIR, 'cDNA.png')))
|
img = np.array(Image.open(os.path.join(DATA_DIR, 'cDNA.png')))
|
||||||
gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8)
|
gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8)
|
||||||
|
|
||||||
# ---- 1. 网格划线(内部自动用 Otsu 算自适应百分比) ----
|
# ---- 1. 网格划线 ----
|
||||||
x_lines, y_lines, pct = draw_grid_lines(gray)
|
x_lines, y_lines, pct = draw_grid_lines(gray)
|
||||||
print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线")
|
print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线")
|
||||||
print(f"自适应百分比: {pct*100:.1f}%")
|
print(f"自适应百分比: {pct*100:.1f}%")
|
||||||
|
|
||||||
# ---- 2. Otsu 阈值分割 ----
|
# ---- 2. Otsu 分割 ----
|
||||||
T_otsu = int(pct * 255)
|
T_otsu = int(pct * 255) # 百分比还原为阈值
|
||||||
binary = (gray > T_otsu).astype(np.uint8)
|
binary = (gray > T_otsu).astype(np.uint8) # 灰度>T 为斑点
|
||||||
print(f"Otsu 阈值: T={T_otsu}")
|
print(f"Otsu 阈值: T={T_otsu}")
|
||||||
|
|
||||||
# ---- 3. 输出:左右对比(划线 vs 分割)----
|
# ---- 3. 输出左右对比图 ----
|
||||||
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
|
fig, axes = plt.subplots(1, 2, figsize=(14, 7))
|
||||||
|
|
||||||
# 左:网格划线
|
# 左图:网格线叠加在灰度图上
|
||||||
axes[0].imshow(gray, cmap='gray')
|
axes[0].imshow(gray, cmap='gray')
|
||||||
for x in x_lines:
|
for x in x_lines:
|
||||||
axes[0].axvline(x=x, color='lime', linewidth=0.5)
|
axes[0].axvline(x=x, color='lime', linewidth=0.5)
|
||||||
@@ -215,9 +206,10 @@ def main():
|
|||||||
axes[0].set_title(f'网格划分 ({len(x_lines)}×{len(y_lines)})', fontsize=13)
|
axes[0].set_title(f'网格划分 ({len(x_lines)}×{len(y_lines)})', fontsize=13)
|
||||||
axes[0].axis('off')
|
axes[0].axis('off')
|
||||||
|
|
||||||
# 右:Otsu 分割结果
|
# 右图:Otsu 二值分割结果
|
||||||
axes[1].imshow(binary, cmap='gray')
|
axes[1].imshow(binary, cmap='gray')
|
||||||
axes[1].set_title(f'Otsu 阈值分割 (T={T_otsu}, pct={pct*100:.1f}%)', fontsize=13)
|
axes[1].set_title(f'Otsu 阈值分割 (T={T_otsu}, pct={pct*100:.1f}%)',
|
||||||
|
fontsize=13)
|
||||||
axes[1].axis('off')
|
axes[1].axis('off')
|
||||||
|
|
||||||
out_path = os.path.join(OUTPUT_DIR, 'gridding_simple.png')
|
out_path = os.path.join(OUTPUT_DIR, 'gridding_simple.png')
|
||||||
|
|||||||
Reference in New Issue
Block a user