针对破损颗粒的Watershed算法优化及参数自动化设置问询

问题背景 #

正在分析颗粒图像，部分颗粒存在破损，现有基于Watershed的分析程序无法正确识别这类颗粒。示例图像如下：

颗粒原始图像 #

带轮廓的图像（显示识别错误） #

带颜色标记的图像（显示识别错误） #

当前分析代码 #

import cv2
import numpy as np
from IPython.display import Image, display
from matplotlib import pyplot as plt
from skimage import io, img_as_float, color, measure, img_as_ubyte
from skimage.segmentation import clear_border
from skimage.restoration import denoise_nl_means, estimate_sigma
path = 'xxx'

# import image and sigma factor for denoising
img_sigma = img_as_float(io.imread(path))
img = img_sigma[0:1024, :]
sigma_est = np.mean(estimate_sigma(img_sigma))
img_ori = cv2.imread(path)
img_ori = img_ori[:1024, :]
sigma_est

# denoise the image
sigma_fak = 2
denoise_img = denoise_nl_means(img, h = sigma_fak*sigma_est, fast_mode = True, patch_size = 5, patch_distance = 5)
denoise_img_ubyte = img_as_ubyte(denoise_img)
denoise_img = cv2.medianBlur(denoise_img_ubyte, 5)
#denoise_img = denoise_nl_means(denoise_img, fast_mode = True, patch_size = 3, patch_distance = 3)
#denoise_img_ubyte = img_as_ubyte(denoise_img)

# turn into gray scale
gray = denoise_img_ubyte

# tresholding
ret, bin_img = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)

# smoothing the image
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
#kernel = np.ones((5,5), np.uint8)
#bin_img = cv2.morphologyEx(bin_img, cv2.MORPH_CLOSE, kernel, iterations=3)
bin_img = cv2.morphologyEx(bin_img, cv2.MORPH_OPEN, kernel, iterations=2)

# removal of the objects that touch the edge of the image
bin_img = clear_border(bin_img)

# sure background area
sure_bg = cv2.dilate(bin_img, kernel, iterations=15)

# Distance transform
dist = cv2.distanceTransform(bin_img, cv2.DIST_L2, 5)

# Make the distance transform normal. 
dist = cv2.normalize(dist, None, 0, 1.0, cv2.NORM_MINMAX) 

# foreground area
ret, sure_fg = cv2.threshold(dist, 0.15 * dist.max(), 255, cv2.THRESH_BINARY)
sure_fg = sure_fg.astype(np.uint8)  

# unknown area
unknown = cv2.subtract(sure_bg, sure_fg)

# Marker labelling sure foreground 
ret, markers = cv2.connectedComponents(sure_fg)
markers += 255
markers[unknown == 255] = 0

# watershed Algorithm
gray = cv2.cvtColor(gray, cv2.COLOR_GRAY2RGB)
markers = cv2.watershed(gray, markers)

# drawing contours from color marks
labels = np.unique(markers)

coins = []
for label in labels[3:]:    
    target = np.where(markers == label, 255, 0).astype(np.uint8)
    contours, hierarchy = cv2.findContours(target, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    
    coins.append(contours[0])

#img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
img = cv2.drawContours(img_ori, coins, -1, color=(0, 255, 0), thickness=2)

img2 = color.label2rgb(markers, bg_label=255)
img2[img2 == 0] = 0

当前手动调整的参数 #

目前需手动调整多类参数优化结果，包括：

• 添加Median模糊或修改denoise_nl_means的patch_size、patch_distance参数：

  denoise_img = denoise_nl_means(img, h = sigma_fak*sigma_est, fast_mode = True, patch_size = 5, patch_distance = 5)
  denoise_img = cv2.medianBlur(denoise_img_ubyte, 5)
  

• 调整确定背景区域的迭代次数：

  sure_bg = cv2.dilate(bin_img, kernel, iterations=15)
  

• 设置距离变换中确定前景的比例：

  ret, sure_fg = cv2.threshold(dist, 0.15 * dist.max(), 255, cv2.THRESH_BINARY)
  sure_fg = sure_fg.astype(np.uint8)
  

• 修改开运算的迭代次数：

  bin_img = cv2.morphologyEx(bin_img, cv2.MORPH_OPEN, kernel, iterations=2)
  

咨询问题 #

1. 如何确保算法模糊/忽略颗粒裂缝，但不影响颗粒间边界？
2. 是否有自动化方案让程序依据图像自动选择参数？

一、忽略颗粒裂缝但保留颗粒边界的方法 #

1. 针对性形态学操作替换开运算 #

当前开运算会同时磨蚀颗粒边界和裂缝，建议优先用闭运算填补裂缝，再用小尺寸开运算清理噪点：

# 用小尺寸椭圆核闭运算填充颗粒内部裂缝（迭代次数根据裂缝宽度调整）
kernel_close = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3))
bin_img = cv2.morphologyEx(bin_img, cv2.MORPH_CLOSE, kernel_close, iterations=2)
# 用极小尺寸开运算清理微小噪点，避免破坏颗粒间边界
kernel_open = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2,2))
bin_img = cv2.morphologyEx(bin_img, cv2.MORPH_OPEN, kernel_open, iterations=1)

闭运算会填充颗粒内部的细小裂缝，小尺寸开运算仅清理孤立噪点，不会影响颗粒间的清晰边界。

2. 优化距离变换与标记生成 #

当前全局阈值易把裂缝误判为前景，建议用局部最大值定位颗粒核心作为标记：

from skimage.feature import peak_local_max

# 生成距离变换后，提取局部最大值作为颗粒核心标记
local_max = peak_local_max(dist, min_distance=20, exclude_border=False)
markers = np.zeros_like(dist, dtype=np.int32)
for i, (y, x) in enumerate(local_max):
    markers[y, x] = i + 1
# 执行Watershed分割
markers = cv2.watershed(gray, markers)

这种方式能精准定位每个颗粒的核心区域，避免裂缝干扰标记生成。

3. 调整降噪策略 #

调换降噪顺序，先中值模糊过滤裂缝类高频噪声，再用非局部均值保留颗粒边界细节：

# 先做中值模糊，针对性过滤裂缝这类小尺寸噪声
denoise_img = cv2.medianBlur(img_as_ubyte(img), 3)
# 再用非局部均值降噪，保留颗粒边界的细节
denoise_img = denoise_nl_means(img_as_float(denoise_img), h=sigma_fak*sigma_est, fast_mode=True, patch_size=3, patch_distance=5)

二、参数自动化选择方案 #

1. 基于图像特征的自适应参数 #

• 形态学参数：通过初始轮廓计算颗粒平均尺寸，自动调整核大小和迭代次数：

  # 检测初始轮廓，计算颗粒平均半径
  contours, _ = cv2.findContours(bin_img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
  avg_area = np.mean([cv2.contourArea(c) for c in contours])
  avg_radius = np.sqrt(avg_area / np.pi)
  # 闭运算核尺寸设为预估裂缝宽度（如平均半径的1/10）
  crack_width = int(avg_radius / 10)
  kernel_close = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (crack_width, crack_width))
  

• 距离变换阈值：基于距离变换直方图自动选择阈值，比如取累积概率90%对应的数值：

  hist, bins = np.histogram(dist.flatten(), bins=256)
  cdf = hist.cumsum()
  cdf_normalized = cdf / cdf[-1]
  threshold_idx = np.where(cdf_normalized >= 0.9)[0][0]
  sure_fg_thresh = bins[threshold_idx]
  ret, sure_fg = cv2.threshold(dist, sure_fg_thresh, 255, cv2.THRESH_BINARY)
  

2. 基于网格搜索的参数优化 #

结合分割质量评估指标（如Dice系数），用网格搜索自动寻找最优参数：

from sklearn.model_selection import ParameterGrid

# 定义待优化的参数网格
param_grid = {
    'patch_size': [3,5,7],
    'median_kernel': [3,5],
    'close_iter': [1,2,3],
    'fg_thresh_ratio': [0.1,0.15,0.2]
}

# 评估函数：计算分割结果与标注的Dice系数
def evaluate(params, ground_truth):
    # 代入当前参数执行分割流程（省略重复代码）
    # ...
    # 计算Dice系数
    intersection = np.sum(markers[ground_truth>0] == ground_truth[ground_truth>0])
    dice = 2 * intersection / (np.sum(markers) + np.sum(ground_truth))
    return dice

# 遍历参数网格，记录最优结果
best_score = 0
best_params = {}
for params in ParameterGrid(param_grid):
    score = evaluate(params, ground_truth_img)
    if score > best_score:
        best_score = score
        best_params = params

若无标注数据，可改用无监督指标（如分割区域的圆形度、面积分布合理性）评估。

内容的提问来源于stack exchange，提问作者Tokar