diff --git a/results_simple/gridding_simple.png b/results_simple/gridding_simple.png index 2cab982..b1c0246 100644 Binary files a/results_simple/gridding_simple.png and b/results_simple/gridding_simple.png differ diff --git a/src/cDNA_gridding_simple.py b/src/cDNA_gridding_simple.py index bf18420..067fc70 100644 --- a/src/cDNA_gridding_simple.py +++ b/src/cDNA_gridding_simple.py @@ -6,12 +6,13 @@ D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py 一、算法流程总览 - 灰度图 ──→ Otsu求像素最佳阈值 T ──→ 百分比 = T/255（自适应，不写死） + 灰度图 ──→ Otsu求像素最佳阈值 T ──→ 百分比 = T/255（自适应） │ - ├─→ 横轴投影/纵轴投影 ──→ X = (max-min) × 百分比 ──→ 减阈值 ──→ - │ 过零点配对 ──→ 空隙中点 ──→ 网格线 + ├─→ 投影/减阈值/过零点配对 ──→ 网格线 │ - └─→ gray > T ──→ 二值图（分割结果） + ├─→ 逐格 Otsu 分割 ──→ keep_largest_object（每格留最大块） + │ + └─→ remove_small_objects（中位数25%以下判为噪声）──→ 统计斑点数 二、各步骤详解 @@ -19,29 +20,25 @@ D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py 2. Otsu 自动阈值 遍历灰度 0~255，每个候选 T 将像素分为前景(>T)和背景(≤T)， - 计算类内方差 w_bg×σ²_bg + w_fg×σ²_fg， - 选使类内方差最小的 T 作为最佳分割线。 - 百分比 = T / 255，取代原来的固定 10%。 + 计算类内方差 w_bg×σ²_bg + w_fg×σ²_fg，选使类内方差最小的 T。 3. 投影 - 横轴：np.sum(每列) → 长度=宽度的曲线，高点=斑点列，低点=空隙列 - 纵轴：np.sum(每行) → 长度=高度的曲线，高点=斑点行，低点=空隙行 + 横轴：np.sum(每列) → 曲线，高点=斑点列，低点=空隙列 + 纵轴：np.sum(每行) → 曲线，高点=斑点行，低点=空隙行 - 4. 阈值 X = (max-min) × 百分比 + 4. 阈值 X = (max-min) × (T/255) - 5. 曲线减 X - - 大于 0 = 斑点区域 - - 小于 0 = 斑点之间的空隙 - - 等于 0 处 = 过零点 = 斑点和空隙的分界线 + 5. 曲线减 X → 大于 0 = 斑点区域，小于 0 = 空隙 + 过零点 = 斑点和空隙的分界线 6. 过零点配对 - 过零点交替出现：正→负（离开斑点）、负→正（进入下一斑点） 配对「离开斑点 + 进入下一斑点」，中点 = 空隙中央 = 划线位置 - 7. 分割 - gray > T 为前景（斑点），≤T 为背景 + 7. 逐格分割 + 后处理 + 对每个格子独立做 Otsu → keep_largest_object（留最大块） + → 全局 remove_small_objects（自动去噪）→ 统计斑点数 - 8. 输出左右对比图：左=网格划线，右=Otsu分割 + 8. 输出三栏图：左=网格，中=分割，右=后处理结果 """ import os @@ -49,169 +46,144 @@ import numpy as np import matplotlib.pyplot as plt from PIL import Image from skimage import color +from scipy import ndimage plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei'] plt.rcParams['axes.unicode_minus'] = False -# 路径设置（从脚本位置动态推导，禁止硬编码绝对路径） SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) BASE_DIR = os.path.dirname(SCRIPT_DIR) DATA_DIR = os.path.join(BASE_DIR, 'cDNA图像处理实例', '数据', 'cDNA') OUTPUT_DIR = os.path.join(BASE_DIR, 'results_simple') -# ================================================================ -# 函数1：Otsu 像素级阈值 -# ================================================================ -def otsu_threshold_pixels(gray: np.ndarray) -> float: - """ - 对图像像素做 Otsu 自动阈值检测。 - - 遍历灰度值 0~255，找到使"类内方差"最小的阈值 T。 - 类内方差 = w_bg × σ²_bg + w_fg × σ²_fg - （背景权重 × 背景方差 + 前景权重 × 前景方差） - - 返回 T/255，即自适应百分比，供投影曲线使用。 - """ - best_T = 0 # 当前最佳阈值 - best_cost = float('inf') # 当前最小类内方差 - total = gray.size # 总像素数 - +def otsu_threshold_pixels(gray: np.ndarray) -> int: + """像素级 Otsu 自动阈值——遍历0~255找最小类内方差，返回 T""" + best_T = 0 + best_cost = float('inf') + total = gray.size for T in range(1, 255): - # 将像素按 T 分为两组 - bg = gray[gray <= T] # 背景 - fg = gray[gray > T] # 前景 - w_bg = len(bg) / total # 背景权重 - w_fg = len(fg) / total # 前景权重 - + bg = gray[gray <= T] + fg = gray[gray > T] + w_bg = len(bg) / total + w_fg = len(fg) / total if w_bg == 0 or w_fg == 0: - continue # 某组为空，跳过 - - # 类内方差 = 加权平均方差 + continue cost = w_bg * np.var(bg) + w_fg * np.var(fg) - if cost < best_cost: best_cost = cost best_T = T - - return best_T / 255.0 # 归一化为百分比 + return best_T -# ================================================================ -# 函数2：网格划线 -# ================================================================ def draw_grid_lines(gray: np.ndarray): - """ - 检测网格分割线。 - - 先用 Otsu 算出自适应百分比， - 再对列投影和行投影分别处理： - 投影 → 减阈值 → 过零点配对 → 空隙中点 = 网格线 - - 返回 (纵线x列表, 横线y列表, 自适应百分比)。 - """ - # ---- 0. 自适应百分比 ---- - pct = otsu_threshold_pixels(gray) - + """网格划线：投影 → 减阈值 → 过零点配对 → 空隙中点""" + T = otsu_threshold_pixels(gray) + pct = T / 255.0 H, W = gray.shape - - # ---- 1. 横轴投影（列方向）---- - # 对每一列上所有行的像素灰度求和 → 长度为 W 的数组 + # 列/行投影 col_profile = np.sum(gray, axis=0).astype(float) - - # ---- 2. 纵轴投影（行方向）---- - # 对每一行上所有列的像素灰度求和 → 长度为 H 的数组 row_profile = np.sum(gray, axis=1).astype(float) - - # ---- 3. 投影阈值 ---- col_T = (np.max(col_profile) - np.min(col_profile)) * pct row_T = (np.max(row_profile) - np.min(row_profile)) * pct - - # ---- 4. 曲线减去阈值 ---- - # 减完：正 = 斑点区域，负 = 空隙 col_shifted = col_profile - col_T row_shifted = row_profile - row_T - # ---- 5. 过零点配对 → 空隙中线 ---- - def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray: - """ - 在减去阈值后的曲线上，找到空隙中线 = 网格线位置。 - - 原理图解： - 信号: ----++++----++++----++++ - ↑ ↑ ↑ ↑ - 过零点配对：第1个+→- 与 第1个-→+ - → 中点 = 空隙中央 = 划线位置 - """ - # 每个位置是正(斑点)还是负(空隙) - is_positive = prof_shifted > 0 - - # 收集符号变化处（过零点） - crossings = [] # 存过零点的位置 - for i in range(1, len(is_positive)): - if is_positive[i] != is_positive[i - 1]: # 符号变化 + def find_gap_lines(prof): + is_pos = prof > 0 + crossings = [] + for i in range(1, len(is_pos)): + if is_pos[i] != is_pos[i - 1]: crossings.append(i) - - if len(crossings) < 2: # 过零点不足 → 放弃 + if len(crossings) < 2: return np.array([]) - - # 过零点交替：正→负，负→正，正→负，负→正 ... - # 要配对的是「离开斑点(正→负)」+「进入下一斑点(负→正)」 - # 如果信号开头是负，跳过第一个 crossing - start = 1 if not is_positive[0] else 0 - + start = 1 if not is_pos[0] else 0 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) - + lines.append(int((crossings[k] + crossings[k + 1]) / 2)) return np.array(lines) - x_lines = find_gap_lines(col_shifted) - y_lines = find_gap_lines(row_shifted) - - return x_lines, y_lines, pct + return find_gap_lines(col_shifted), find_gap_lines(row_shifted), T, pct + + +def keep_largest_object(binary: np.ndarray) -> np.ndarray: + """每格只保留最大连通域""" + labeled, num = ndimage.label(binary) + if num == 0: + return np.zeros_like(binary) + areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)] + return (labeled == (int(np.argmax(areas)) + 1)).astype(np.uint8) + + +def remove_small_objects(binary: np.ndarray) -> np.ndarray: + """自动去噪——连通域面积 < 中位数25% 的剔除""" + labeled, num = ndimage.label(binary) + if num == 0: + return binary + areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)] + median = np.median(areas) + min_size = max(1, int(median * 0.25)) + result = binary.copy() + for i in range(1, num + 1): + if areas[i - 1] < min_size: + result[labeled == i] = 0 + return result -# ================================================================ -# 主流程 -# ================================================================ def main(): os.makedirs(OUTPUT_DIR, exist_ok=True) - # ---- 读取图像，转为灰度 ---- + # 读取图像，转灰度 img = np.array(Image.open(os.path.join(DATA_DIR, 'cDNA.png'))) gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8) - # ---- 1. 网格划线 ---- - x_lines, y_lines, pct = draw_grid_lines(gray) + # 1. 网格划线 + x_lines, y_lines, T_otsu, pct = draw_grid_lines(gray) print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线") - print(f"自适应百分比: {pct*100:.1f}%") + print(f"Otsu 阈值: T={T_otsu}, 自适应百分比: {pct*100:.1f}%") - # ---- 2. Otsu 分割 ---- - T_otsu = int(pct * 255) # 百分比还原为阈值 - binary = (gray > T_otsu).astype(np.uint8) # 灰度>T 为斑点 - print(f"Otsu 阈值: T={T_otsu}") + # 2. 逐格分割 + 后处理 + bw_full = np.zeros_like(gray) + for i in range(len(y_lines) - 1): + for j in range(len(x_lines) - 1): + r1, r2 = y_lines[i], y_lines[i + 1] + c1, c2 = x_lines[j], x_lines[j + 1] + blk = gray[r1:r2, c1:c2] + if blk.size == 0: + continue + T = otsu_threshold_pixels(blk) + bw_blk = (blk > T).astype(np.uint8) + bw_blk = keep_largest_object(bw_blk) + bw_full[r1:r2, c1:c2] = bw_blk - # ---- 3. 输出左右对比图 ---- - fig, axes = plt.subplots(1, 2, figsize=(14, 7)) + bw_clean = remove_small_objects(bw_full) + + # 3. 统计斑点 + labeled, num = ndimage.label(bw_clean) + spot_sizes = [int(np.sum(labeled == i)) for i in range(1, num + 1)] + valid = [s for s in spot_sizes if s >= 10] + print(f"检测到 {len(valid)} 个斑点") + + # 4. 输出三栏图 + fig, axes = plt.subplots(1, 3, figsize=(20, 7)) - # 左图：网格线叠加在灰度图上 axes[0].imshow(gray, cmap='gray') for x in x_lines: axes[0].axvline(x=x, color='lime', linewidth=0.5) for y in y_lines: axes[0].axhline(y=y, color='lime', linewidth=0.5) - axes[0].set_title(f'网格划分 ({len(x_lines)}×{len(y_lines)})', fontsize=13) + axes[0].set_title(f'网格划分 ({len(x_lines)}x{len(y_lines)})', fontsize=13) axes[0].axis('off') - # 右图：Otsu 二值分割结果 - axes[1].imshow(binary, cmap='gray') - axes[1].set_title(f'Otsu 阈值分割 (T={T_otsu}, pct={pct*100:.1f}%)', - fontsize=13) + axes[1].imshow(bw_full, cmap='gray') + axes[1].set_title('逐格 Otsu 分割', fontsize=13) axes[1].axis('off') + axes[2].imshow(bw_clean, cmap='gray') + axes[2].set_title(f'后处理结果 ({len(valid)}个斑点)', fontsize=13) + axes[2].axis('off') + out_path = os.path.join(OUTPUT_DIR, 'gridding_simple.png') fig.savefig(out_path, dpi=150, bbox_inches='tight') plt.close(fig)