diff --git a/results_simple/gridding_simple.png b/results_simple/gridding_simple.png index f7eb087..2cab982 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 00eb2d1..a34a33b 100644 --- a/src/cDNA_gridding_simple.py +++ b/src/cDNA_gridding_simple.py @@ -5,11 +5,12 @@ cDNA微阵列图像处理 —— 简化版 D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py 算法步骤(划线): + 先用 Otsu 算出像素级最佳阈值 T,百分比 = T / 255 来自数据而非写死 1. 彩色图 → 灰度图 2. 横轴投影:对每一列的所有像素灰度值求和 → 得到一条曲线 纵轴投影:对每一行的所有像素灰度值求和 → 得到一条曲线 -3. 在曲线上,求出 max 和 min,阈值 X = (max - min) × 10% +3. 在曲线上,求出 max 和 min,阈值 X = (max - min) × (T / 255) 4. 曲线上每个值都减去 X 5. 减完之后: - 大于 0 的地方 = 斑点区域 @@ -44,9 +45,33 @@ DATA_DIR = os.path.join(BASE_DIR, 'cDNA图像处理实例', '数据', 'cDNA') OUTPUT_DIR = os.path.join(BASE_DIR, 'results_simple') -def draw_grid_lines(gray: np.ndarray, pct: float = 0.10): +def otsu_threshold_pixels(gray: np.ndarray) -> float: """ - 核心算法:检测网格分割线。 + Otsu 自动阈值——在像素级找到最佳分割灰度值 T。 + + 遍历 0~255,对每个候选 T 计算类内方差,选最小的。 + 返回 T / 255,作为投影曲线阈值的自适应百分比。 + """ + best_T = 0 + best_cost = float('inf') + total = gray.size + for T in range(1, 255): + 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 + 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 + + +def draw_grid_lines(gray: np.ndarray): + """ + 检测网格分割线。先对像素做 Otsu,用 T/255 作为自适百分比。 原理: 灰度图的每一列/行,属于斑点的像素灰度值高,属于背景的灰度值低。 @@ -58,13 +83,11 @@ def draw_grid_lines(gray: np.ndarray, pct: float = 0.10): 过零点的位置就是斑点和空隙的分界线, 两个分界线中点就是划线位置。 - 参数: - gray: 灰度图 (高×宽) - pct: 阈值百分比,默认10% - 返回: - (纵线x坐标列表, 横线y坐标列表) + (纵线x坐标列表, 横线y坐标列表, 自适应百分比) """ + # 0. 先用 Otsu 算出自适应百分比 + pct = otsu_threshold_pixels(gray) H, W = gray.shape # ================================================================ @@ -85,10 +108,10 @@ def draw_grid_lines(gray: np.ndarray, pct: float = 0.10): row_profile = np.sum(gray, axis=1).astype(float) # ================================================================ - # 步骤3:计算阈值 X = (max - min) × 10% + # 步骤3:计算阈值 X = (max - min) × 自适应百分比 # ================================================================ - # max-min 是曲线的"振幅",取10%作为阈值 - # 大于这个阈值的才是真正的斑点信号,小于的是噪声 + # 百分比来自 Otsu(像素级最佳分割线),不是写死的 10% + # max-min 是曲线的"振幅" col_T = (np.max(col_profile) - np.min(col_profile)) * pct row_T = (np.max(row_profile) - np.min(row_profile)) * pct @@ -160,46 +183,7 @@ def draw_grid_lines(gray: np.ndarray, pct: float = 0.10): x_lines = find_gap_lines(col_shifted) y_lines = find_gap_lines(row_shifted) - return x_lines, y_lines - - -def otsu_threshold(gray: np.ndarray) -> tuple[np.ndarray, float]: - """ - Otsu 自动阈值分割。 - - 原理: - 遍历灰度值 0~255,对每个候选 T: - - 将像素分为两组:前景(>T) 和 背景(≤T) - - 计算两组各自的权重 w 和方差 σ² - - 类内方差 = w_bg*σ²_bg + w_fg*σ²_fg - 选使类内方差最小的 T = 最佳分割线。 - """ - best_T = 0 - best_cost = float('inf') - total = gray.size - - for T in range(1, 255): - # 背景(灰度 ≤ T) - bg = gray[gray <= T] - w_bg = len(bg) / total - # 前景(灰度 > T) - fg = gray[gray > T] - w_fg = len(fg) / total - - if w_bg == 0 or w_fg == 0: - continue - - var_bg = np.var(bg) - var_fg = np.var(fg) - # 类内方差 = 加权平均 - cost = w_bg * var_bg + w_fg * var_fg - - if cost < best_cost: - best_cost = cost - best_T = T - - binary = (gray > best_T).astype(np.uint8) - return binary, best_T + return x_lines, y_lines, pct def main(): @@ -207,16 +191,17 @@ def main(): # ---- 读取图像并转为灰度图 ---- img = np.array(Image.open(os.path.join(DATA_DIR, 'cDNA.png'))) - # 原图是 RGBA(红绿蓝+透明通道),取前三个通道转为灰度 gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8) - # ---- 1. 网格划线 ---- - x_lines, y_lines = draw_grid_lines(gray) + # ---- 1. 网格划线(内部自动用 Otsu 算自适应百分比) ---- + x_lines, y_lines, pct = draw_grid_lines(gray) print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线") + print(f"自适应百分比: {pct*100:.1f}%") # ---- 2. Otsu 阈值分割 ---- - binary, T_otsu = otsu_threshold(gray) - print(f"Otsu 最佳阈值: {T_otsu}") + T_otsu = int(pct * 255) + binary = (gray > T_otsu).astype(np.uint8) + print(f"Otsu 阈值: T={T_otsu}") # ---- 3. 输出:左右对比(划线 vs 分割)---- fig, axes = plt.subplots(1, 2, figsize=(14, 7)) @@ -232,7 +217,7 @@ def main(): # 右:Otsu 分割结果 axes[1].imshow(binary, cmap='gray') - axes[1].set_title(f'Otsu 阈值分割 (T={T_otsu})', fontsize=13) + axes[1].set_title(f'Otsu 阈值分割 (T={T_otsu}, pct={pct*100:.1f}%)', fontsize=13) axes[1].axis('off') out_path = os.path.join(OUTPUT_DIR, 'gridding_simple.png')