如何通过OpenCV获取晶粒图像的像素值与多阶交界点坐标？

需求概述 #

• 识别相场模拟生成的晶粒结构图像中的每个晶粒，获取其全部像素坐标
• 提取晶界的像素坐标，同时检测triple junctions（三重交界）及更高阶交界点

现有代码问题分析 #

你提供的代码存在几个核心问题，导致无法准确达成目标：

1. 错误地以灰度模式读取原始图像，丢失了区分不同晶粒的颜色信息
2. 晶界检测逻辑绕远路，未直接利用晶界特征或原始图像的颜色突变
3. 交界点检测的形态学操作精度不足，易出现误检、漏检

改进解决方案 #

核心思路 #

1. 基于颜色聚类的晶粒识别：利用原始彩色图像的颜色特征，通过聚类区分不同晶粒，结合连通分量过滤噪声
2. 直接提取晶界像素：借助高亮晶界的图像快速获取晶界坐标，或从原始图像通过颜色突变检测
3. 邻域分析检测交界点：对晶界骨架做邻域连通数统计，判定三重及更高阶交界

完整实现代码 #

import cv2
import numpy as np
import matplotlib.pyplot as plt

# -------------------------- 1. 图像读取与预处理 --------------------------
# 读取原始彩色晶粒图像（替换为你的图像路径）
original_img = cv2.imread('original_grain.png')
original_img_rgb = cv2.cvtColor(original_img, cv2.COLOR_BGR2RGB)

# 读取高亮晶界的图像（替换为你的图像路径）
grain_boundary_img = cv2.imread('grain_boundary.png', cv2.IMREAD_GRAYSCALE)

# -------------------------- 2. 晶粒识别与像素提取 --------------------------
# 将图像像素转换为二维数组，用于颜色聚类
pixel_values = original_img_rgb.reshape((-1, 3))
pixel_values = np.float32(pixel_values)

# KMeans颜色聚类，k值可根据实际晶粒数量调整
k = 50
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.2)
_, labels, centers = cv2.kmeans(pixel_values, k, None, criteria, 10, cv2.KMEANS_RANDOM_CENTERS)

# 重塑聚类结果为图像格式
centers = np.uint8(centers)
clustered_img = centers[labels.flatten()].reshape(original_img_rgb.shape)

# 过滤小噪声区域，存储每个有效晶粒的像素坐标
grain_data = {}
for label in np.unique(labels):
    mask = (labels == label).reshape(original_img_rgb.shape[:2])
    coords = np.column_stack(np.where(mask))
    # 过滤像素数少于30的小区域（噪声）
    if len(coords) > 30:
        grain_data[label] = {
            'pixel_coords': coords,
            'color': centers[label]
        }

print(f"识别到有效晶粒数量：{len(grain_data)}")

# -------------------------- 3. 晶界像素提取 --------------------------
# 从高亮晶界图像中提取所有非白色像素（晶界）
boundary_coords = np.column_stack(np.where(grain_boundary_img < 255))
print(f"提取到晶界像素数量：{len(boundary_coords)}")

# -------------------------- 4. 交界点检测 --------------------------
# 对晶界图像做骨架提取，细化晶界
skeleton = cv2.ximgproc.thinning(grain_boundary_img)

# 定义3x3邻域卷积核，统计每个像素的邻域连通数
kernel = np.array([[1, 1, 1],
                   [1, 0, 1],
                   [1, 1, 1]], dtype=np.uint8)
neighbor_count = cv2.filter2D(skeleton, -1, kernel)

# 交界点判定：骨架像素且邻域连通数≥3
junction_coords = np.column_stack(np.where((skeleton == 255) & (neighbor_count >= 3)))
print(f"检测到交界点数量：{len(junction_coords)}")

# -------------------------- 5. 结果可视化 --------------------------
plt.figure(figsize=(15, 10))

# 原始晶粒图像
plt.subplot(221)
plt.imshow(original_img_rgb)
plt.title('原始晶粒结构')
plt.axis('off')

# 聚类后的晶粒识别结果
plt.subplot(222)
plt.imshow(clustered_img)
plt.title('聚类后的晶粒识别')
plt.axis('off')

# 标记晶界的图像
boundary_vis = original_img_rgb.copy()
boundary_vis[boundary_coords[:, 0], boundary_coords[:, 1]] = [255, 0, 0]
plt.subplot(223)
plt.imshow(boundary_vis)
plt.title('标记晶界的晶粒结构')
plt.axis('off')

# 标记交界点的图像
junction_vis = original_img_rgb.copy()
for (y, x) in junction_coords:
    cv2.circle(junction_vis, (x, y), 3, (0, 0, 255), -1)
plt.subplot(224)
plt.imshow(junction_vis)
plt.title('标记交界点的晶粒结构')
plt.axis('off')

plt.tight_layout()
plt.show()

# -------------------------- 数据导出（可选） --------------------------
# 导出每个晶粒的像素坐标
for grain_id, data in grain_data.items():
    np.save(f'grain_{grain_id}_coords.npy', data['pixel_coords'])

# 导出晶界和交界点坐标
np.save('grain_boundary_coords.npy', boundary_coords)
np.save('junctions_coords.npy', junction_coords)

关键说明 #

• 颜色聚类保留了原始图像的晶粒区分特征，避免灰度化导致的信息丢失
• 直接利用高亮晶界图像提取晶界，比通过形态学操作推导更精准高效
• 邻域连通数统计能准确识别三重及更高阶交界点，减少误检
• 提供了数据导出功能，方便后续的量化分析

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