refactor: 恢复简化版所有详细行内注释
每函数、每关键行都有中文注释说明原理
This commit is contained in:
+147
-43
@@ -20,7 +20,7 @@ D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py
|
|||||||
|
|
||||||
2. Otsu 自动阈值
|
2. Otsu 自动阈值
|
||||||
遍历灰度 0~255,每个候选 T 将像素分为前景(>T)和背景(≤T),
|
遍历灰度 0~255,每个候选 T 将像素分为前景(>T)和背景(≤T),
|
||||||
计算类内方差 w_bg×σ²_bg + w_fg×σ²_fg,选使类内方差最小的 T。
|
计算类内方差 w_bg×σ²_bg + w_fg×σ²_fg,选使方差最小的 T。
|
||||||
|
|
||||||
3. 投影
|
3. 投影
|
||||||
横轴:np.sum(每列) → 曲线,高点=斑点列,低点=空隙列
|
横轴:np.sum(每列) → 曲线,高点=斑点列,低点=空隙列
|
||||||
@@ -32,6 +32,7 @@ D:\ProgramData\anaconda3\envs\my_env\python.exe src/cDNA_gridding_simple.py
|
|||||||
过零点 = 斑点和空隙的分界线
|
过零点 = 斑点和空隙的分界线
|
||||||
|
|
||||||
6. 过零点配对
|
6. 过零点配对
|
||||||
|
过零点交替:正→负(离开斑点)、负→正(进入下一斑点)
|
||||||
配对「离开斑点 + 进入下一斑点」,中点 = 空隙中央 = 划线位置
|
配对「离开斑点 + 进入下一斑点」,中点 = 空隙中央 = 划线位置
|
||||||
|
|
||||||
7. 逐格分割 + 后处理
|
7. 逐格分割 + 后处理
|
||||||
@@ -48,82 +49,175 @@ from PIL import Image
|
|||||||
from skimage import color
|
from skimage import color
|
||||||
from scipy import ndimage
|
from scipy import ndimage
|
||||||
|
|
||||||
|
# matplotlib 中文字体设置
|
||||||
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) -> int:
|
def otsu_threshold_pixels(gray: np.ndarray) -> int:
|
||||||
"""像素级 Otsu 自动阈值——遍历0~255找最小类内方差,返回 T"""
|
"""
|
||||||
best_T = 0
|
对图像像素做 Otsu 自动阈值检测。
|
||||||
best_cost = float('inf')
|
|
||||||
total = gray.size
|
遍历灰度值 0~255,对每个候选 T:
|
||||||
|
- 将像素分为两组:前景(>T) 和 背景(≤T)
|
||||||
|
- 计算类内方差 = w_bg × σ²_bg + w_fg × σ²_fg
|
||||||
|
- 选使类内方差最小的 T
|
||||||
|
|
||||||
|
返回 T(0~255 整数)。
|
||||||
|
"""
|
||||||
|
best_T = 0 # 当前最佳阈值
|
||||||
|
best_cost = float('inf') # 当前最小类内方差
|
||||||
|
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 # 某组为空(T 太极端),跳过
|
||||||
|
|
||||||
|
# 类内方差 = 加权平均方差
|
||||||
|
# 方差小 = 组内像素灰度接近 = 分组效果好
|
||||||
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
|
return best_T
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 函数2:网格划线
|
||||||
|
# ================================================================
|
||||||
def draw_grid_lines(gray: np.ndarray):
|
def draw_grid_lines(gray: np.ndarray):
|
||||||
"""网格划线:投影 → 减阈值 → 过零点配对 → 空隙中点"""
|
"""
|
||||||
T = otsu_threshold_pixels(gray)
|
检测网格分割线。
|
||||||
pct = T / 255.0
|
|
||||||
|
流程:
|
||||||
|
Otsu 求自适应百分比 → 列/行投影 → 减阈值 → 过零点配对 → 空隙中点
|
||||||
|
|
||||||
|
返回 (纵线x列表, 横线y列表, T, 自适应百分比)。
|
||||||
|
"""
|
||||||
|
T = otsu_threshold_pixels(gray) # 像素级最佳阈值
|
||||||
|
pct = T / 255.0 # 归一化为百分比
|
||||||
H, W = gray.shape
|
H, W = gray.shape
|
||||||
# 列/行投影
|
|
||||||
|
# ---- 1. 横轴投影 ----
|
||||||
|
# 对每一列上所有行的灰度求和 → 长度=宽度的数组
|
||||||
|
# 斑点所在列亮像素多 → 和较大(曲线凸起)
|
||||||
|
# 空隙所在列暗像素多 → 和较小(曲线凹陷)
|
||||||
col_profile = np.sum(gray, axis=0).astype(float)
|
col_profile = np.sum(gray, axis=0).astype(float)
|
||||||
|
|
||||||
|
# ---- 2. 纵轴投影 ----
|
||||||
|
# 对每一行上所有列的灰度求和 → 长度=高度的数组
|
||||||
row_profile = np.sum(gray, axis=1).astype(float)
|
row_profile = np.sum(gray, axis=1).astype(float)
|
||||||
|
|
||||||
|
# ---- 3. 投影阈值 ----
|
||||||
|
# 振幅 × 自适应百分比 = 动态阈值(不写死,每张图自调整)
|
||||||
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. 曲线减去阈值 ----
|
||||||
|
# 减完后:正 = 斑点区域,负 = 空隙,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
|
||||||
|
|
||||||
def find_gap_lines(prof):
|
# ---- 5. 过零点配对 → 空隙中线 ----
|
||||||
is_pos = prof > 0
|
def find_gap_lines(prof_shifted: np.ndarray) -> np.ndarray:
|
||||||
|
"""
|
||||||
|
在减去阈值后的曲线上,配对过零点,取空隙中央。
|
||||||
|
|
||||||
|
原理图解:
|
||||||
|
信号: ----++++----++++----++++
|
||||||
|
↑ ↑ ↑ ↑
|
||||||
|
过零点配对:离开斑点 + 进入下一个斑点
|
||||||
|
→ 中点 = 空隙中央 = 划线位置
|
||||||
|
"""
|
||||||
|
# 每个位置是正(斑点)还是负(空隙)
|
||||||
|
is_positive = prof_shifted > 0
|
||||||
|
|
||||||
|
# 收集符号变化位置(过零点)
|
||||||
crossings = []
|
crossings = []
|
||||||
for i in range(1, len(is_pos)):
|
for i in range(1, len(is_positive)):
|
||||||
if is_pos[i] != is_pos[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([])
|
||||||
start = 1 if not is_pos[0] else 0
|
|
||||||
|
# 过零点交替:正→负(离开斑点), 负→正(进入下一斑点)
|
||||||
|
# 要配对的是"离开斑点 → 进入下一斑点",即空隙的两端
|
||||||
|
# 如果信号开头是负,跳过第一个 crossing
|
||||||
|
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):
|
||||||
lines.append(int((crossings[k] + crossings[k + 1]) / 2))
|
# crossings[k]: 正→负(离开斑点)
|
||||||
|
# crossings[k+1]: 负→正(进入下一斑点)
|
||||||
|
# 中点 = 空隙中央 = 划线位置
|
||||||
|
mid = int((crossings[k] + crossings[k + 1]) / 2)
|
||||||
|
lines.append(mid)
|
||||||
|
|
||||||
return np.array(lines)
|
return np.array(lines)
|
||||||
|
|
||||||
return find_gap_lines(col_shifted), find_gap_lines(row_shifted), T, pct
|
x_lines = find_gap_lines(col_shifted)
|
||||||
|
y_lines = find_gap_lines(row_shifted)
|
||||||
|
|
||||||
|
return x_lines, y_lines, T, pct
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 函数3:后处理(完全自动,无需人工设定阈值)
|
||||||
|
# ================================================================
|
||||||
def keep_largest_object(binary: np.ndarray) -> np.ndarray:
|
def keep_largest_object(binary: np.ndarray) -> np.ndarray:
|
||||||
"""每格只保留最大连通域"""
|
"""
|
||||||
|
每个格子里只保留面积最大的连通域。
|
||||||
|
|
||||||
|
ndimage.label 给每个白色连通域编号 → 算面积 → 只留最大那块。
|
||||||
|
不需要设定任何阈值。
|
||||||
|
"""
|
||||||
labeled, num = ndimage.label(binary)
|
labeled, num = ndimage.label(binary)
|
||||||
if num == 0:
|
if num == 0:
|
||||||
return np.zeros_like(binary)
|
return np.zeros_like(binary) # 全黑,直接返回
|
||||||
|
# 统计每个连通域的像素数
|
||||||
areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
||||||
return (labeled == (int(np.argmax(areas)) + 1)).astype(np.uint8)
|
# 找面积最大的编号
|
||||||
|
max_idx = int(np.argmax(areas)) + 1
|
||||||
|
return (labeled == max_idx).astype(np.uint8)
|
||||||
|
|
||||||
|
|
||||||
def remove_small_objects(binary: np.ndarray) -> np.ndarray:
|
def remove_small_objects(binary: np.ndarray) -> np.ndarray:
|
||||||
"""自动去噪——连通域面积 < 中位数25% 的剔除"""
|
"""
|
||||||
|
自动去除小连通域(噪声)。
|
||||||
|
|
||||||
|
统计所有连通域面积的中位数,
|
||||||
|
小于中位数 25% 的视为噪声,直接剔除。
|
||||||
|
换图换分辨率都自动适应,不需要手动调参。
|
||||||
|
"""
|
||||||
labeled, num = ndimage.label(binary)
|
labeled, num = ndimage.label(binary)
|
||||||
if num == 0:
|
if num == 0:
|
||||||
return binary
|
return binary # 全黑,直接返回
|
||||||
|
|
||||||
|
# 收集所有连通域的面积
|
||||||
areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
areas = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
||||||
median = np.median(areas)
|
median = np.median(areas) # 面积中位数
|
||||||
min_size = max(1, int(median * 0.25))
|
min_size = max(1, int(median * 0.25)) # 中位数的25%,最少1像素
|
||||||
|
|
||||||
|
# 面积不达标的连通域整块置0
|
||||||
result = binary.copy()
|
result = binary.copy()
|
||||||
for i in range(1, num + 1):
|
for i in range(1, num + 1):
|
||||||
if areas[i - 1] < min_size:
|
if areas[i - 1] < min_size:
|
||||||
@@ -131,43 +225,51 @@ def remove_small_objects(binary: np.ndarray) -> np.ndarray:
|
|||||||
return result
|
return result
|
||||||
|
|
||||||
|
|
||||||
|
# ================================================================
|
||||||
|
# 主流程
|
||||||
|
# ================================================================
|
||||||
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')))
|
||||||
|
# 原图 RGBA,取前三个通道转为 0~255 灰度图
|
||||||
gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8)
|
gray = (color.rgb2gray(img[:, :, :3]) * 255).astype(np.uint8)
|
||||||
|
|
||||||
# 1. 网格划线
|
# ---- 1. 网格划线 ----
|
||||||
x_lines, y_lines, T_otsu, pct = draw_grid_lines(gray)
|
x_lines, y_lines, T_otsu, pct = draw_grid_lines(gray)
|
||||||
print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线")
|
print(f"检测到 {len(x_lines)} 条纵线, {len(y_lines)} 条横线")
|
||||||
print(f"Otsu 阈值: T={T_otsu}, 自适应百分比: {pct*100:.1f}%")
|
print(f"Otsu 阈值: T={T_otsu}, 自适应百分比: {pct*100:.1f}%")
|
||||||
|
|
||||||
# 2. 逐格分割 + 后处理
|
# ---- 2. 逐格分割 + 后处理 ----
|
||||||
bw_full = np.zeros_like(gray)
|
bw_full = np.zeros_like(gray) # 初始化全黑底图
|
||||||
for i in range(len(y_lines) - 1):
|
for i in range(len(y_lines) - 1): # 遍历每一行格子
|
||||||
for j in range(len(x_lines) - 1):
|
for j in range(len(x_lines) - 1): # 遍历每一列格子
|
||||||
r1, r2 = y_lines[i], y_lines[i + 1]
|
r1, r2 = y_lines[i], y_lines[i + 1] # 格子行范围
|
||||||
c1, c2 = x_lines[j], x_lines[j + 1]
|
c1, c2 = x_lines[j], x_lines[j + 1] # 格子列范围
|
||||||
blk = gray[r1:r2, c1:c2]
|
blk = gray[r1:r2, c1:c2] # 提取该格子
|
||||||
if blk.size == 0:
|
if blk.size == 0:
|
||||||
continue
|
continue
|
||||||
|
# 对单个格子做 Otsu 分割
|
||||||
T = otsu_threshold_pixels(blk)
|
T = otsu_threshold_pixels(blk)
|
||||||
bw_blk = (blk > T).astype(np.uint8)
|
bw_blk = (blk > T).astype(np.uint8) # 大于 T 为前景
|
||||||
|
# 每格只保留最大连通域(去噪点)
|
||||||
bw_blk = keep_largest_object(bw_blk)
|
bw_blk = keep_largest_object(bw_blk)
|
||||||
bw_full[r1:r2, c1:c2] = bw_blk
|
bw_full[r1:r2, c1:c2] = bw_blk # 拼回全局图
|
||||||
|
|
||||||
|
# 全局去除小连通域(自动判别噪声阈值)
|
||||||
bw_clean = remove_small_objects(bw_full)
|
bw_clean = remove_small_objects(bw_full)
|
||||||
|
|
||||||
# 3. 统计斑点
|
# ---- 3. 统计斑点 ----
|
||||||
labeled, num = ndimage.label(bw_clean)
|
labeled, num = ndimage.label(bw_clean)
|
||||||
spot_sizes = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
spot_sizes = [int(np.sum(labeled == i)) for i in range(1, num + 1)]
|
||||||
valid = [s for s in spot_sizes if s >= 10]
|
valid = [s for s in spot_sizes if s >= 10] # 过滤极小噪声
|
||||||
print(f"检测到 {len(valid)} 个斑点")
|
print(f"检测到 {len(valid)} 个斑点")
|
||||||
|
|
||||||
# 4. 输出三栏图
|
# ---- 4. 输出三栏图 ----
|
||||||
fig, axes = plt.subplots(1, 3, figsize=(20, 7))
|
fig, axes = plt.subplots(1, 3, figsize=(20, 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)
|
||||||
@@ -176,10 +278,12 @@ def main():
|
|||||||
axes[0].set_title(f'网格划分 ({len(x_lines)}x{len(y_lines)})', fontsize=13)
|
axes[0].set_title(f'网格划分 ({len(x_lines)}x{len(y_lines)})', fontsize=13)
|
||||||
axes[0].axis('off')
|
axes[0].axis('off')
|
||||||
|
|
||||||
|
# 中:逐格 Otsu 分割结果(后处理前)
|
||||||
axes[1].imshow(bw_full, cmap='gray')
|
axes[1].imshow(bw_full, cmap='gray')
|
||||||
axes[1].set_title('逐格 Otsu 分割', fontsize=13)
|
axes[1].set_title('逐格 Otsu 分割', fontsize=13)
|
||||||
axes[1].axis('off')
|
axes[1].axis('off')
|
||||||
|
|
||||||
|
# 右:后处理之后的最终二值图
|
||||||
axes[2].imshow(bw_clean, cmap='gray')
|
axes[2].imshow(bw_clean, cmap='gray')
|
||||||
axes[2].set_title(f'后处理结果 ({len(valid)}个斑点)', fontsize=13)
|
axes[2].set_title(f'后处理结果 ({len(valid)}个斑点)', fontsize=13)
|
||||||
axes[2].axis('off')
|
axes[2].axis('off')
|
||||||
|
|||||||
Reference in New Issue
Block a user