如何用Python检测显微图像中类弹坑特征的中心与半径？

处理对象为代表表面高度图的numpy数组，目标是检测类弹坑特征的中心与半径，用于后续生成径向剖面。当前手动阈值二值化效果不稳定，高低阈值均会降低圆/椭圆拟合精度，Canny边缘检测未达预期，以下是针对性优化方案：

1. 自适应阈值+形态学处理替代手动阈值 #

手动阈值无法适配不同图像的高度分布，改用自适应阈值二值化结合形态学操作，精准提取弹坑边缘区域：

from skimage.filters import threshold_local
from skimage.morphology import binary_closing, disk

# 自适应阈值二值化，block_size根据图像尺寸调整（建议为奇数）
local_thresh = threshold_local(image_grey, block_size=35, offset=0.1)
rim_binarized = image_grey >= local_thresh

# 形态学闭操作填补边缘缺口，disk半径根据弹坑实际尺寸调整
rim_binarized = binary_closing(rim_binarized, disk(3))
points = np.argwhere(rim_binarized)

2. 基于梯度边缘的霍夫圆变换 #

跳过二值化区域提取，直接用边缘梯度结合霍夫圆变换检测，鲁棒性更强：

from skimage.feature import canny
from skimage.transform import hough_circle, hough_circle_peaks

# Canny边缘检测，sigma值根据图像噪声程度调整
edges = canny(image_grey, sigma=2)

# 设定弹坑半径范围（需根据实际样本尺寸调整）
hough_radii = np.arange(10, 50, 2)
hough_res = hough_circle(edges, hough_radii)

# 提取最显著的圆特征
accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii, total_num_peaks=1)

# 转换为整数坐标并绘制
cy_int, cx_int, r_int = [int(round(x)) for x in (cy[0], cx[0], radii[0])]
circ_y, circ_x = circle_perimeter(cy_int, cx_int, r_int)

3. 峰值点初始化+径向剖面拟合 #

先定位弹坑的高度峰值作为初始中心，再沿径向找高度突变点（弹坑边缘），拟合更精准的圆参数：

# 找到高度图的峰值点作为初始中心
y_peak, x_peak = np.unravel_index(np.argmax(image_grey), image_grey.shape)

# 沿360度方向采样径向剖面
angles = np.linspace(0, 2*np.pi, 360)
edge_distances = []
for angle in angles:
    # 生成沿角度方向的采样坐标
    x_steps = np.arange(x_peak, image_grey.shape[1], np.cos(angle))
    y_steps = np.arange(y_peak, image_grey.shape[0], np.sin(angle))
    # 截断到图像有效范围
    valid_idx = (x_steps >=0) & (x_steps < image_grey.shape[1]) & (y_steps >=0) & (y_steps < image_grey.shape[0])
    x_steps = x_steps[valid_idx].astype(int)
    y_steps = y_steps[valid_idx].astype(int)
    # 计算高度梯度，找到边缘位置
    profile = image_grey[y_steps, x_steps]
    grad = np.abs(np.diff(profile))
    if len(grad) > 0:
        edge_idx = np.argmax(grad) + 1
        distance = np.sqrt((x_steps[edge_idx]-x_peak)**2 + (y_steps[edge_idx]-y_peak)**2)
        edge_distances.append(distance)

# 用所有边缘点拟合圆
edge_points = []
for d, angle in zip(edge_distances, angles):
    x = int(x_peak + d*np.cos(angle))
    y = int(y_peak + d*np.sin(angle))
    edge_points.append([y, x])

circ = CircleModel()
circ.estimate(np.array(edge_points))
yc, xc, r = (int(round(x)) for x in circ.params)

4. 利用Diplib增强边缘检测 #

Diplib的边缘检测工具更适配高度图特征，替代skimage的Canny：

# 用Diplib进行边缘检测
dipimg = dip.Image(image_grey)
edges = dip.Canny(dipimg, sigma=1.5, lowThreshold=0.1, highThreshold=0.3)
edges_np = np.array(edges)

# 提取边缘点拟合圆
points = np.argwhere(edges_np)
circ = CircleModel()
circ.estimate(points)
yc, xc, r = (int(round(x)) for x in circ.params)

整合优化后的完整代码示例 #

from skimage.measure import EllipseModel, CircleModel
from skimage.draw import ellipse_perimeter, circle_perimeter
from skimage.filters import threshold_local
from skimage.morphology import binary_closing, disk
from skimage.feature import canny
from skimage.transform import hough_circle, hough_circle_peaks
import numpy as np
import diplib as dip
import matplotlib.pyplot as plt
from pathlib import Path

base_path = Path("folder including multiple images")

for file in base_path.rglob("*_f.npz"):
    image_grey: np.ndarray = np.load(file)["arr_0"]

    # 方案1：自适应阈值+形态学处理
    local_thresh = threshold_local(image_grey, block_size=35, offset=0.1)
    rim_binarized = image_grey >= local_thresh
    rim_binarized = binary_closing(rim_binarized, disk(3))
    points = np.argwhere(rim_binarized)

    # 拟合椭圆和圆
    ell = EllipseModel()
    succ = ell.estimate(points)
    ye, xe, a, b = (int(round(x)) for x in ell.params[:-1])
    ey, ex = ellipse_perimeter(ye,xe,a,b,ell.params[-1])

    circ = CircleModel()
    circ.estimate(points)
    yc, xc, r = (int(round(x)) for x in circ.params)
    cy, cx = circle_perimeter(yc,xc,r)

    # 方案2：霍夫圆检测
    edges = canny(image_grey, sigma=2)
    hough_radii = np.arange(10, 50, 2)
    hough_res = hough_circle(edges, hough_radii)
    accums, h_cx, h_cy, h_radii = hough_circle_peaks(hough_res, hough_radii, total_num_peaks=1)
    h_cy_int, h_cx_int, h_r_int = [int(round(x)) for x in (h_cy[0], h_cx[0], h_radii[0])]
    h_cy_draw, h_cx_draw = circle_perimeter(h_cy_int, h_cx_int, h_r_int)

    # 绘图展示
    fig2, (ax1, ax2, ax3) = plt.subplots(ncols=3, nrows=1, figsize=(12, 4))
    ax1.set_title('Original + Detected Circles')
    ax1.imshow(image_grey)
    ax1.plot(cx,cy, "r,", label="Fitted Circle")
    ax1.plot(ex,ey, "m,", label="Fitted Ellipse")
    ax1.plot(h_cx_draw, h_cy_draw, "g-", linewidth=2, label="Hough Circle")
    ax1.legend()

    ax2.set_title('Binarized Rim')
    ax2.imshow(rim_binarized, cmap="Greys")

    ax3.set_title("Radial Profile")
    dipimg = dip.Image(image_grey)
    rad = dip.RadialMean(dipimg, binSize=1, center=(ye,xe), maxRadius=r)
    ax3.plot(rad, label="Ellipse Center")
    rad = dip.RadialMean(dipimg, binSize=1, center=(yc,xc), maxRadius=r)
    ax3.plot(rad, label="Fitted Circle Center")
    rad = dip.RadialMean(dipimg, binSize=1, center=(h_cy_int,h_cx_int), maxRadius=h_r_int)
    ax3.plot(rad, label="Hough Circle Center")
    ax3.legend()

    plt.show()

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