feat: 简化版增加逐格分割+后处理+斑点统计
现在简化版也具备完整处理链： 网格划线 → 逐格Otsu → keep_largest_object → remove_small_objects → 统计 输出三栏图：网格/分割/后处理结果
This commit is contained in:
+98
-126
@@ -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)
|
||||
|
||||
Reference in New Issue
Block a user