从砾石照片到颗粒参数:Python 自动提取粒径、圆度和形状分布(含代码)

import osimport numpy as npimport pandas as pdimport cv2import matplotlibmatplotlib.use("Agg")import matplotlib.pyplot as pltfrom scipy import ndimage as ndifrom skimage import morphology, measure, segmentation, feature, colorfrom skimage.filters import rank# 当前目录下的图片名称# 改成你的砾石图片文件名即可input_image = "example.jpg"out_dir = "gravel_image_analysis_results"os.makedirs(out_dir, exist_ok=True)np.random.seed(2026)plt.rcParams["font.family"] = "Times New Roman"plt.rcParams["axes.unicode_minus"] = Falseplt.rcParams["figure.dpi"] = 160plt.rcParams["savefig.dpi"] = 500# 如果有比例尺，在这里填写每个像素对应多少 mm# 例如：MM_PER_PIXEL = 0.25# 如果不知道，就保持 None，输出单位为 pxMM_PER_PIXEL = None# 颗粒过滤参数MIN_AREA_PX = 80MAX_AREA_PX = 50000MIN_EQ_DIAM_PX = 6MAX_EQ_DIAM_PX = 260# 如果过滤后颗粒太少，代码会自动放宽筛选MIN_VALID_PARTICLES = 30# 去掉贴边颗粒，避免边缘半颗粒影响统计BORDER_MARGIN = 4# Canny 边缘参数# 如果分割太碎：适当增大 CANNY_LOW / CANNY_HIGH# 如果颗粒粘连严重：适当减小 CANNY_LOW / CANNY_HIGHCANNY_LOW = 45CANNY_HIGH = 115# 距离变换局部峰值参数# 如果颗粒被切太碎：增大 PEAK_MIN_DISTANCE# 如果相邻颗粒分不开：减小 PEAK_MIN_DISTANCEPEAK_MIN_DISTANCE = 8PEAK_THRESHOLD_ABS = 3.0# 梯度分水岭参数GRADIENT_DISK_SIZE = 2EDGE_DILATE_RADIUS = 1WATERSHED_COMPACTNESS = 0.01# 过暗区域剔除阈值# 如果很多深色砾石漏掉，可以调低，例如 20 或 25OBJECT_MIN_GRAY = 25# 粒径累计曲线模式# "area"：按颗粒投影面积加权，更接近二维面积级配# "number"：按颗粒数量加权DISTRIBUTION_MODE = "area"# 读取图片if not os.path.exists(input_image):    raise FileNotFoundError(        f"没有找到图片：{input_image}\n"        f"请确认图片在当前目录下，或者把 input_image 改成完整路径。"    )img_bgr = cv2.imread(input_image)if img_bgr is None:    raise ValueError("图片读取失败，请检查文件格式。")img_rgb = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB)H, W = img_rgb.shape[:2]# 3. 天然灰色砾石图像分割#    CLAHE增强 + Canny边缘 + 距离变换 + 分水岭gray = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2GRAY)# 局部对比度增强，突出颗粒边界clahe = cv2.createCLAHE(    clipLimit=2.5,    tileGridSize=(8, 8))gray_eq = clahe.apply(gray)# 中值滤波，降低纹理噪声gray_blur = cv2.medianBlur(gray_eq, 5)# Canny 边缘检测edges = cv2.Canny(gray_blur, CANNY_LOW, CANNY_HIGH)# 边界膨胀，形成颗粒之间的分隔线edge_bool = morphology.binary_dilation(    edges > 0,    morphology.disk(EDGE_DILATE_RADIUS))# 候选颗粒内部区域：非边缘 + 非极暗区域interior = (~edge_bool) & (gray_blur > OBJECT_MIN_GRAY)interior = morphology.remove_small_objects(interior, min_size=25)interior = morphology.binary_closing(interior, morphology.disk(1))# 距离变换，用局部极大值作为颗粒中心标记distance = ndi.distance_transform_edt(interior)coords = feature.peak_local_max(    distance,    min_distance=PEAK_MIN_DISTANCE,    threshold_abs=PEAK_THRESHOLD_ABS,    labels=interior,    exclude_border=BORDER_MARGIN)# 如果标记点太少，自动放宽一次if len(coords) < 20:    coords = feature.peak_local_max(        distance,        min_distance=max(4, PEAK_MIN_DISTANCE - 3),        threshold_abs=max(1.5, PEAK_THRESHOLD_ABS - 1.5),        labels=interior,        exclude_border=BORDER_MARGIN    )markers = np.zeros_like(distance, dtype=np.int32)for i, (r, c) in enumerate(coords, start=1):    markers[r, c] = imarkers = ndi.label(markers > 0)[0]if markers.max() == 0:    raise RuntimeError(        "没有生成有效分水岭标记。建议降低 CANNY_LOW/CANNY_HIGH，"        "或降低 PEAK_THRESHOLD_ABS。"    )# 梯度图作为分水岭地形gradient = rank.gradient(    gray_blur,    morphology.disk(GRADIENT_DISK_SIZE))# 当前图基本全是砾石，mask 只剔除极暗区域和边界区域mask_clean = gray_blur > OBJECT_MIN_GRAYmask_clean[:BORDER_MARGIN, :] = Falsemask_clean[-BORDER_MARGIN:, :] = Falsemask_clean[:, :BORDER_MARGIN] = Falsemask_clean[:, -BORDER_MARGIN:] = Falselabels_ws = segmentation.watershed(    gradient,    markers,    mask=mask_clean,    compactness=WATERSHED_COMPACTNESS)# 去掉贴边区域labels_ws = segmentation.clear_border(labels_ws)# 提取单颗粒形状与粒径参数def px_to_unit(x):    if MM_PER_PIXEL is None:        return x    return x * MM_PER_PIXELdef px2_to_unit2(x):    if MM_PER_PIXEL is None:        return x    return x * MM_PER_PIXEL * MM_PER_PIXELdef get_unit_name():    return "px" if MM_PER_PIXEL is None else "mm"def get_equivalent_diameter(region):    try:        return region.equivalent_diameter_area    except AttributeError:        return region.equivalent_diameterunit_name = get_unit_name()candidate_records = []regions = measure.regionprops(labels_ws)for region in regions:    source_label = region.label    area = region.area    minr, minc, maxr, maxc = region.bbox    # 去掉明显贴边半颗粒    if (        minr <= BORDER_MARGIN        or minc <= BORDER_MARGIN        or maxr >= H - BORDER_MARGIN        or maxc >= W - BORDER_MARGIN    ):        continue    if area <= 10:        continue    equiv_d = get_equivalent_diameter(region)    perimeter = region.perimeter    major_axis = region.major_axis_length    minor_axis = region.minor_axis_length    solidity = region.solidity    eccentricity = region.eccentricity    extent = region.extent    orientation = np.degrees(region.orientation)    if perimeter <= 0 or minor_axis <= 0 or major_axis <= 0:        continue    circularity = 4 * np.pi * area / (perimeter ** 2 + 1e-12)    aspect_ratio = major_axis / (minor_axis + 1e-12)    roundness = 4 * area / (np.pi * major_axis ** 2 + 1e-12)    coords_region = region.coords    pts = coords_region[:, ::-1].astype(np.float32)    rect = cv2.minAreaRect(pts)    rect_w, rect_h = rect[1]    rect_long = max(rect_w, rect_h)    rect_short = min(rect_w, rect_h)    if rect_short <= 0:        continue    rect_aspect = rect_long / rect_short    cy, cx = region.centroid    candidate_records.append({        "source_label": source_label,        "centroid_x_px": cx,        "centroid_y_px": cy,        "area_px2": area,        f"area_{unit_name}2": px2_to_unit2(area),        "perimeter_px": perimeter,        f"perimeter_{unit_name}": px_to_unit(perimeter),        f"equivalent_diameter_{unit_name}": px_to_unit(equiv_d),        f"major_axis_{unit_name}": px_to_unit(major_axis),        f"minor_axis_{unit_name}": px_to_unit(minor_axis),        f"rect_long_axis_{unit_name}": px_to_unit(rect_long),        f"rect_short_axis_{unit_name}": px_to_unit(rect_short),        "aspect_ratio": aspect_ratio,        "rect_aspect_ratio": rect_aspect,        "circularity": circularity,        "roundness": roundness,        "solidity": solidity,        "extent": extent,        "eccentricity": eccentricity,        "orientation_deg": orientation    })raw_df = pd.DataFrame(candidate_records)if len(raw_df) == 0:    raise RuntimeError(        "分水岭完成后没有得到候选颗粒。建议降低 OBJECT_MIN_GRAY，"        "或减小 CANNY_LOW/CANNY_HIGH。"    )diam_col = f"equivalent_diameter_{unit_name}"# 先按正常标准筛选valid_mask = (    (raw_df["area_px2"] >= MIN_AREA_PX)    & (raw_df["area_px2"] <= MAX_AREA_PX)    & (raw_df[diam_col] >= MIN_EQ_DIAM_PX)    & (raw_df[diam_col] <= MAX_EQ_DIAM_PX)    & (raw_df["aspect_ratio"] <= 12.0)    & (raw_df["circularity"] > 0.03)    & (raw_df["solidity"] > 0.35))metrics_df = raw_df[valid_mask].copy()# 如果筛选后颗粒太少，自动放宽筛选if len(metrics_df) < MIN_VALID_PARTICLES:    print("提示：正常筛选后颗粒数量偏少，已自动放宽筛选条件。")    valid_mask_relaxed = (        (raw_df["area_px2"] >= 30)        & (raw_df["area_px2"] <= max(MAX_AREA_PX, raw_df["area_px2"].quantile(0.995)))        & (raw_df[diam_col] >= 4)        & (raw_df[diam_col] <= max(MAX_EQ_DIAM_PX, raw_df[diam_col].quantile(0.995)))        & (raw_df["aspect_ratio"] <= 20.0)        & (raw_df["solidity"] > 0.15)    )    metrics_df = raw_df[valid_mask_relaxed].copy()# 重新编号metrics_df = metrics_df.reset_index(drop=True)metrics_df.insert(0, "label", np.arange(1, len(metrics_df) + 1))if len(metrics_df) == 0:    raise RuntimeError(        "仍然没有识别到有效砾石颗粒。建议检查图片名称是否正确，"        "或进一步降低 MIN_AREA_PX、OBJECT_MIN_GRAY。"    )# 生成有效标签图valid_label_mask = np.zeros_like(labels_ws, dtype=np.int32)for _, row in metrics_df.iterrows():    source_label = int(row["source_label"])    new_label = int(row["label"])    valid_label_mask[labels_ws == source_label] = new_labelmetrics_path = os.path.join(out_dir, "gravel_particle_metrics.csv")metrics_df.to_csv(metrics_path, index=False, encoding="utf-8-sig")raw_path = os.path.join(out_dir, "gravel_raw_candidate_regions.csv")raw_df.to_csv(raw_path, index=False, encoding="utf-8-sig")# 5. 计算粒径分布与 D10/D30/D50/D60diameters = metrics_df[diam_col].valuesarea_weights = metrics_df["area_px2"].valuesorder = np.argsort(diameters)diam_sorted = diameters[order]if DISTRIBUTION_MODE == "area":    weights_sorted = area_weights[order]    cum_percent = np.cumsum(weights_sorted) / np.sum(weights_sorted) * 100.0    cum_label = "Cumulative area percent (%)"else:    cum_percent = np.arange(1, len(diam_sorted) + 1) / len(diam_sorted) * 100.0    cum_label = "Cumulative number percent (%)"def interpolate_D_value(d_sorted, cum_percent_arr, target_percent):    d_sorted = np.asarray(d_sorted)    cum_percent_arr = np.asarray(cum_percent_arr)    if target_percent < cum_percent_arr.min() or target_percent > cum_percent_arr.max():        return np.nan    log_d = np.log10(d_sorted + 1e-12)    log_D = np.interp(target_percent, cum_percent_arr, log_d)    return float(10 ** log_D)D10 = interpolate_D_value(diam_sorted, cum_percent, 10)D30 = interpolate_D_value(diam_sorted, cum_percent, 30)D50 = interpolate_D_value(diam_sorted, cum_percent, 50)D60 = interpolate_D_value(diam_sorted, cum_percent, 60)Cu = D60 / D10 if D10 > 0 else np.nanCc = D30 ** 2 / (D10 * D60) if D10 > 0 and D60 > 0 else np.nansummary_df = pd.DataFrame([{    "distribution_mode": DISTRIBUTION_MODE,    "particle_count": len(metrics_df),    f"D10_{unit_name}": D10,    f"D30_{unit_name}": D30,    f"D50_{unit_name}": D50,    f"D60_{unit_name}": D60,    "Cu_D60_over_D10": Cu,    "Cc_D30_square_over_D10_D60": Cc,    f"mean_diameter_{unit_name}": np.mean(diameters),    f"median_diameter_{unit_name}": np.median(diameters),    f"std_diameter_{unit_name}": np.std(diameters),    "mean_aspect_ratio": metrics_df["aspect_ratio"].mean(),    "mean_rect_aspect_ratio": metrics_df["rect_aspect_ratio"].mean(),    "mean_circularity": metrics_df["circularity"].mean(),    "mean_solidity": metrics_df["solidity"].mean(),    "mean_roundness": metrics_df["roundness"].mean()}])summary_path = os.path.join(out_dir, "gravel_size_shape_summary.csv")summary_df.to_csv(summary_path, index=False, encoding="utf-8-sig")# 图1：分割结果检查图fig, axes = plt.subplots(1, 3, figsize=(15, 6))axes[0].imshow(img_rgb)axes[0].set_title("Original gravel image", fontsize=13)axes[0].axis("off")label_vis = color.label2rgb(    valid_label_mask,    image=img_rgb,    bg_label=0,    alpha=0.45)axes[1].imshow(label_vis)axes[1].set_title("Separated particle labels", fontsize=13)axes[1].axis("off")overlay = img_rgb.copy()boundaries = segmentation.find_boundaries(valid_label_mask, mode="outer")overlay[boundaries] = [255, 0, 0]axes[2].imshow(overlay)axes[2].set_title(f"Detected particles: n = {len(metrics_df)}", fontsize=13)axes[2].axis("off")plt.tight_layout()fig1_path = os.path.join(out_dir, "fig1_segmentation_check.png")plt.savefig(fig1_path, bbox_inches="tight", facecolor="white")plt.close(fig)# 图2：粒径分布直方图 + 累计曲线fig, ax1 = plt.subplots(figsize=(8.8, 6.4))bins = np.logspace(    np.log10(max(diameters.min() * 0.85, 1e-3)),    np.log10(diameters.max() * 1.15),    28)ax1.hist(    diameters,    bins=bins,    color="#A6CEE3",    edgecolor="white",    alpha=0.90)ax1.set_xscale("log")ax1.set_xlabel(f"Equivalent particle diameter ({unit_name})", fontsize=12)ax1.set_ylabel("Particle count", fontsize=12)ax1.grid(True, which="both", alpha=0.18)ax2 = ax1.twinx()ax2.plot(    diam_sorted,    cum_percent,    color="#D55E00",    linewidth=2.4,    label=cum_label)ax2.set_ylabel(cum_label, fontsize=12)ax2.set_ylim(0, 100)for D_val, P_val, name in [    (D10, 10, "D10"),    (D30, 30, "D30"),    (D60, 60, "D60")]:    ax1.axvline(        D_val,        color="gray",        linestyle="--",        linewidth=1.0,        alpha=0.65    )    ax2.scatter(        D_val,        P_val,        s=58,        color="red",        edgecolors="white",        linewidths=0.8,        zorder=5    )    ax2.text(        D_val * 1.05,        P_val + 3.0,        f"{name}={D_val:.2f}",        fontsize=10,        ha="left",        va="bottom"    )txt = (    f"n = {len(metrics_df)}\n"    f"D10 = {D10:.2f}{unit_name}\n"    f"D30 = {D30:.2f}{unit_name}\n"    f"D60 = {D60:.2f}{unit_name}\n"    f"Cu = {Cu:.2f}\n"    f"Cc = {Cc:.2f}")ax1.text(    0.04,    0.96,    txt,    transform=ax1.transAxes,    ha="left",    va="top",    fontsize=10,    bbox=dict(        boxstyle="round,pad=0.35",        facecolor="white",        edgecolor="#999999",        alpha=0.92    ))ax1.set_title("Particle size distribution from gravel image", fontsize=15)ax1.spines["top"].set_visible(False)ax2.spines["top"].set_visible(False)line_handles, line_labels = ax2.get_legend_handles_labels()ax2.legend(    line_handles,    line_labels,    loc="lower right",    frameon=True,    fontsize=10)plt.tight_layout()fig2_path = os.path.join(out_dir, "fig2_particle_size_distribution.png")plt.savefig(fig2_path, bbox_inches="tight", facecolor="white")plt.close(fig)# 图3：形状参数分布图fig, axes = plt.subplots(2, 2, figsize=(11.2, 8.4))axes = axes.ravel()shape_items = [    ("aspect_ratio", "Aspect ratio", "#4C72B0"),    ("circularity", "Circularity", "#55A868"),    ("roundness", "Roundness", "#C44E52"),    ("solidity", "Solidity", "#8172B2")]for ax, (col, title, color_val) in zip(axes, shape_items):    vals = metrics_df[col].replace([np.inf, -np.inf], np.nan).dropna().values    ax.hist(        vals,        bins=24,        color=color_val,        edgecolor="white",        alpha=0.88    )    mean_val = np.mean(vals)    ax.axvline(        mean_val,        color="black",        linestyle="--",        linewidth=1.3,        label=f"Mean = {mean_val:.2f}"    )    ax.set_title(title, fontsize=13)    ax.set_xlabel(title, fontsize=11)    ax.set_ylabel("Count", fontsize=11)    ax.grid(alpha=0.18)    ax.legend(frameon=True, fontsize=9)    ax.spines["top"].set_visible(False)    ax.spines["right"].set_visible(False)fig.suptitle("Shape parameter distributions of detected gravel particles", fontsize=16, y=0.98)plt.tight_layout(rect=[0, 0, 1, 0.96])fig3_path = os.path.join(out_dir, "fig3_shape_parameter_distribution.png")plt.savefig(fig3_path, bbox_inches="tight", facecolor="white")plt.close(fig)# 图4：粒径-形状参数关系图fig, axes = plt.subplots(1, 3, figsize=(14.8, 4.8))scatter_items = [    ("aspect_ratio", "Aspect ratio"),    ("circularity", "Circularity"),    ("solidity", "Solidity")]sc = Nonefor ax, (col, ylabel) in zip(axes, scatter_items):    sc = ax.scatter(        metrics_df[diam_col],        metrics_df[col],        c=metrics_df["area_px2"],        s=28,        cmap="viridis",        alpha=0.78,        edgecolors="none"    )    ax.set_xscale("log")    ax.set_xlabel(f"Equivalent diameter ({unit_name})", fontsize=11)    ax.set_ylabel(ylabel, fontsize=11)    ax.set_title(f"Size vs. {ylabel}", fontsize=13)    ax.grid(True, which="both", alpha=0.18)    ax.spines["top"].set_visible(False)    ax.spines["right"].set_visible(False)fig.suptitle(    "Relationship between particle size and shape parameters",    fontsize=16,    y=1.02)plt.tight_layout(rect=[0, 0, 0.88, 0.95])cax = fig.add_axes([0.91, 0.18, 0.018, 0.62])cbar = fig.colorbar(sc, cax=cax)cbar.set_label("Projected area (px²)", fontsize=11)fig4_path = os.path.join(out_dir, "fig4_size_shape_scatter.png")plt.savefig(fig4_path, bbox_inches="tight", facecolor="white")plt.close(fig)

如有其他需要，欢迎关注我的咸鱼号:pfc小姐姐