【金山文档 | WPS云文档】 传统影像分割 https://www.kdocs.cn/l/caKCDb2cyYfL
一 阈值分割（如Otsu）
概述
Otsu方法的精髓在于它能够自动找到一个最佳全局阈值，而无需人为指定。这个阈值的标准是：将图像分为前景和背景两类后，使得两类像素的类间方差（Between-class variance）最大。因为方差越大，说明两类差别越大，错分的可能性就越小。
具体流程
- 预处理： 将原始的RGB彩色影像转换为灰度影像。有时会先进行平滑滤波（如高斯模糊）以减少噪声干扰。
- 应用Otsu算法： 算法会遍历所有可能的灰度阈值，计算每个阈值下前景和背景像素的类间方差，并选择使方差最大的那个阈值作为最佳分割阈值。
- 二值化： 使用Otsu计算出的阈值，将灰度图像转换为黑白二值图像。白色部分（通常代表树冠）为潜在的目标区域。
- 后处理： 对二值图像进行操作：
形态学操作： 使用“开运算”（先腐蚀后膨胀）来去除小的噪声点（如阳光斑、石头），并使用“闭运算”（先膨胀后腐蚀）来填充树冠内部可能存在的空洞或细小缝隙。
连通区域分析： 标记二值图像中所有连通的白色区域，每一个连通区域就可以被初步识别为一棵单独的树。
受光照/阴影影响的原因： 树冠内部可能因结构产生阴影（暗），而阳光直射的树叶很亮（亮）。一个全局阈值很难同时正确分割出明亮和阴暗的部分，可能导致同一棵树冠被分割成多个部分，或阴暗的树冠被误判为背景。背景中如果有大片明亮土地，也可能被错误地分割为前景。
实践操作
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 1. 读取图像并转换为灰度图
img = cv2.imread('tree_canopy_image.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 2. 高斯模糊预处理 - 减少噪声和细节干扰
gray_blur = cv2.GaussianBlur(gray, (5, 5), 0) # 内核大小(5,5)可根据图像调整
# 3. 应用Otsu阈值分割
_, thresh = cv2.threshold(gray_blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
print(f"Otsu算法计算出的最佳阈值: {_}")
# 4. 形态学操作 - 优化分割结果
kernel = np.ones((3, 3), np.uint8)
# 先开运算(去噪)后闭运算(填充)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2) # 去除小噪声点
closing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernel, iterations=3) # 填充树冠内部空洞
# 5. 执行连通区域分析
num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(closing, connectivity=8)
# 6. 过滤小区域（可能是噪声而非树冠）并标记结果
min_area = 50 # 最小树冠面积阈值，需根据图像实际尺度调整
tree_count = 0
result_img = img.copy() # 用于绘制结果的副本
# 创建一个彩色图像用于可视化标签
label_hue = np.uint8(179 * labels / np.max(labels))
blank_ch = 255 * np.ones_like(label_hue)
labeled_img = cv2.merge([label_hue, blank_ch, blank_ch])
labeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)
labeled_img[label_hue == 0] = 0 # 将背景设为黑色
for i in range(1, num_labels): # 跳过背景(标签0)
area = stats[i, cv2.CC_STAT_AREA]
if area > min_area:
tree_count += 1
# 获取当前区域的外接矩形坐标
x = stats[i, cv2.CC_STAT_LEFT]
y = stats[i, cv2.CC_STAT_TOP]
w = stats[i, cv2.CC_STAT_WIDTH]
h = stats[i, cv2.CC_STAT_HEIGHT]
# 在原图上绘制矩形框标记树冠
cv2.rectangle(result_img, (x, y), (x + w, y + h), (0, 255, 0), 2)
# 标记树冠编号
cv2.putText(result_img, str(tree_count), (x, y-10),
cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0, 255, 0), 2)
# 7. 显示所有处理阶段的结果
plt.figure(figsize=(18, 12))
plt.subplot(231)
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.title('原始图像')
plt.axis('off')
plt.subplot(232)
plt.imshow(gray, cmap='gray')
plt.title('灰度图像')
plt.axis('off')
plt.subplot(233)
plt.imshow(gray_blur, cmap='gray')
plt.title('高斯模糊后')
plt.axis('off')
plt.subplot(234)
plt.imshow(thresh, cmap='gray')
plt.title('Otsu阈值分割结果')
plt.axis('off')
plt.subplot(235)
plt.imshow(closing, cmap='gray')
plt.title('形态学处理后的二值图像')
plt.axis('off')
plt.subplot(236)
plt.imshow(cv2.cvtColor(result_img, cv2.COLOR_BGR2RGB))
plt.title(f'最终识别结果: 检测到 {tree_count} 棵树')
plt.axis('off')
plt.tight_layout()
plt.show()
# 额外显示连通区域标签图
plt.figure(figsize=(10, 8))
plt.imshow(cv2.cvtColor(labeled_img, cv2.COLOR_BGR2RGB))
plt.title('连通区域标签可视化（不同颜色代表不同树冠）')
plt.axis('off')
plt.show()
print(f"识别出的树木数量: {tree_count}")
在进行Otsu阈值分割之前，先使用高斯模糊对数据进行预处理，目的是平滑图像，减少噪声和细小纹理的干扰，使树冠区域更加均匀；之后通过Otsu阈值分割核心分割步骤，自动确定最佳阈值将图像二值化；直接得出的结果往往比较粗糙会出现噪声和不完整区域，因此我们通过后处理采用开运算和闭运算的方法进行形态学操作；对连通区域进行分析，每个连通区域代表一个潜在树冠；最后进行面积过滤，将面积过小（可能是噪声而不是树冠）的区域过滤掉。
开运算：先腐蚀后膨胀，用于去除小的白噪声（如阳光斑点）和分离粘连不严重的树冠。
闭运算：先膨胀后腐蚀，用于填充树冠内部由于阴影产生的小黑洞，使其成为一个完整的连通区域。
高斯模糊内核大小、形态学操作迭代次数、面积阈值各个步骤参数，要根据具体图像特点进行调整
二 边缘检测（如Canny）
概述
Canny算法是计算机视觉中最经典、最常用的边缘检测算法之一，由John F. Canny于1986年提出。它的目标是一个最优的边缘检测器，不关注区域内部，而是通过检测图像中灰度值急剧变化的位置（即边缘）来勾勒出物体的轮廓。Canny是公认的最优边缘检测算法，它通过高斯滤波、计算梯度、非极大值抑制和双阈值检测等步骤来提取精细、连续的边缘。
最优边缘准则
Canny 的目标是找到一个最优的边缘检测算法，最优边缘检测的含义是：
(1)最优检测：算法能够尽可能多地标识出图像中的实际边缘，漏检真实边缘的概率和误检非边缘的概率都尽可能小；
(2)最优定位准则：检测到的边缘点的位置距离实际边缘点的位置最近，或者是由于噪声影响引起检测出的边缘偏离物体的真实边缘的程度最小；
(3)检测点与边缘点一一对应：算子检测的边缘点与实际边缘点应该是一一对应。
为了满足这些要求 Canny 使用了变分法（calculus of variations），这是一种寻找优化特定功能的函数的方法。最优检测使用四个指数函数项表示，但是它非常近似于高斯函数的一阶导数
具体流程
1 对RGB图像进行灰度化
2 高斯模糊（降噪）
边缘检测对噪声非常敏感，因为噪声会导致图像灰度值的剧烈变化，从而被误认为是边缘。此步骤旨在平滑图像，抑制噪声。
3计算梯度强度和方向
使用一个边缘检测算子（如Sobel、Prewitt）卷积平滑后的图像，分别计算水平方向（Gx）和垂直方向（Gy）的一阶导数。
Sobel算子：
Sobel_X = [[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]]
Sobel_Y = [[-1, -2, -1],
[ 0, 0, 0],
[ 1, 2, 1]]
Gx = Sobel_X * SmoothedSmoothed为高斯模糊处理后的输入图像
Gy = Sobel_Y * Smoothed
计算梯度强度：梯度强度表示了边缘的“强弱”。
G = √(Gx² + Gy²) （为了计算效率，也常用近似值：|Gx| + |Gy|）
计算梯度方向：梯度方向垂直于边缘方向，用于下一步的非极大值抑制。
θ = arctan2(Gy, Gx) （结果通常被四舍五入到0°、45°、90°、135°四个方向之一，以便于后续处理）
4 非极大值抑制（NMS）
细化边缘，在梯度强度图像中，边缘通常会有一定的宽度，NMS通过保留梯度方向上局部最大值点像素细化边缘。
- 根据该点的梯度方向（θ），找到其正方向和反方向上的相邻像素。
例如，如果θ=0°（东/西方向），则比较东和西方向的像素。
如果θ=90°（北/南方向），则比较北和南方向的像素。
如果θ=135°（东北/西南方向），则比较东北和西南方向的像素。
如果θ=45°（西北/东南方向），则比较西北和东南方向的像素。
- 将当前像素的梯度强度（
G）与这两个相邻像素的梯度强度进行比较。 - 如果当前像素的梯度强度是这三个点中最大的，则保留其值；否则，将其抑制（设置为0）。
例：
图中的数字代表了像素点的梯度强度，箭头方向代表了梯度方向。以第二排第三个像素点为例，由于梯度方向向上，则将这一点的强度（7）与其上下两个像素点的强度（5和4）比较，由于这一点强度最大，则保留，最后会输出一个更薄更清晰的边缘图像。
5 双阈值检测
区分“强边缘”、“弱边缘”和“噪声”。
- 设定两个阈值：高阈值（
Thigh） 和 低阈值（Tlow）。通常Thigh是Tlow的2到3倍（例如，Thigh = 200,Tlow = 100）。 - 遍历经过NMS处理的图像中的每一个像素：
如果像素梯度值 ≥ Thigh，则被标记为强边缘像素（确定为真实边缘）。
如果像素梯度值 介于 Tlow 和 Thigh 之间，则被标记为弱边缘像素（可能是边缘，需要进一步验证）。
如果像素梯度值 < Tlow，则被抑制（排除，认为是非边缘）。
6 滞后边界跟踪（通过弱边缘）
通过检查弱边缘像素是否与强边缘像素相连，来决定是否最终保留它们。这确保了边缘的连续性，并减少了孤立的噪声点被保留的可能性。
- 从所有被标记为强边缘的像素开始，沿着8邻域（上、下、左、右、左上、右上、左下、右下）进行追踪。
- 在追踪过程中，如果发现弱边缘像素，则将其提升为强边缘像素（因为它与一个确定的真实边缘相连）。
- 继续以这个新提升的强边缘像素为起点，递归地进行步骤2的追踪，直到没有连通的弱边缘为止。
- 追踪完成后，所有剩余的、未被提升的弱边缘像素则被抑制（丢弃），认为它们是噪声。
最终，输出图像中只包含所有强边缘像素（包括由追踪提升上来的），这就是Canny算法检测到的边缘。
实现代码
纯python代码实现
import matplotlib.pyplot as plt
import numpy as np
import math
import cv2
# 读取图像
img = cv2.imread('./Images/lena.jpg')
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
# 使用OpenCV的Canny实现作为对比
blur = cv2.GaussianBlur(img, (5, 5), 0)
canny = cv2.Canny(blur, 50, 150)
# 手动实现Canny
sigma1 = sigma2 = 1
sum_val = 0
# 生成高斯核
gaussian = np.zeros([5, 5])
for i in range(5):
for j in range(5):
gaussian[i, j] = math.exp(-1/2 * (np.square(i-2)/np.square(sigma1) +
np.square(j-2)/np.square(sigma2))) / (2*math.pi*sigma1*sigma2)
sum_val += gaussian[i, j]
gaussian = gaussian / sum_val
# RGB转灰度函数
def rgb2gray(rgb):
return np.dot(rgb[..., :3], [0.299, 0.587, 0.114])
# Step 1: 高斯滤波
gray = rgb2gray(img)
W, H = gray.shape
new_gray = np.zeros([W-4, H-4])
for i in range(W-4):
for j in range(H-4):
new_gray[i, j] = np.sum(gray[i:i+5, j:j+5] * gaussian)
# Step 2: 计算梯度幅值和方向
W1, H1 = new_gray.shape
dx = np.zeros([W1-1, H1-1])
dy = np.zeros([W1-1, H1-1])
d = np.zeros([W1-1, H1-1])
for i in range(W1-1):
for j in range(H1-1):
dx[i, j] = new_gray[i, j+1] - new_gray[i, j]
dy[i, j] = new_gray[i+1, j] - new_gray[i, j]
d[i, j] = np.sqrt(np.square(dx[i, j]) + np.square(dy[i, j]))
# Step 3: 非极大值抑制
W2, H2 = d.shape
NMS = np.copy(d)
NMS[0, :] = NMS[W2-1, :] = NMS[:, 0] = NMS[:, H2-1] = 0
for i in range(1, W2-1):
for j in range(1, H2-1):
if d[i, j] == 0:
NMS[i, j] = 0
else:
gradX = dx[i, j]
gradY = dy[i, j]
gradTemp = d[i, j]
# 确定梯度方向并插值
if np.abs(gradY) > np.abs(gradX):
weight = np.abs(gradX) / np.abs(gradY)
grad2 = d[i-1, j]
grad4 = d[i+1, j]
if gradX * gradY > 0:
grad1 = d[i-1, j-1]
grad3 = d[i+1, j+1]
else:
grad1 = d[i-1, j+1]
grad3 = d[i+1, j-1]
else:
weight = np.abs(gradY) / np.abs(gradX)
grad2 = d[i, j-1]
grad4 = d[i, j+1]
if gradX * gradY > 0:
grad1 = d[i+1, j-1]
grad3 = d[i-1, j+1]
else:
grad1 = d[i-1, j-1]
grad3 = d[i+1, j+1]
# 双线性插值并比较
gradTemp1 = weight * grad1 + (1 - weight) * grad2
gradTemp2 = weight * grad3 + (1 - weight) * grad4
if gradTemp >= gradTemp1 and gradTemp >= gradTemp2:
NMS[i, j] = gradTemp
else:
NMS[i, j] = 0
# Step 4: 双阈值检测和滞后边界跟踪
def hysteresis_edge_tracking(nms, low_threshold_ratio=0.05, high_threshold_ratio=0.09):
# 计算高低阈值
high_threshold = np.max(nms) * high_threshold_ratio
low_threshold = high_threshold * low_threshold_ratio
# 初始化结果矩阵
result = np.zeros_like(nms, dtype=np.uint8)
# 标记强边缘和弱边缘
strong_edges = (nms > high_threshold)
weak_edges = (nms >= low_threshold) & (nms <= high_threshold)
# 强边缘直接保留
result[strong_edges] = 255
# 对弱边缘进行滞后跟踪
# 找到所有强边缘点的坐标
strong_y, strong_x = np.where(strong_edges)
# 使用队列进行边缘跟踪
from collections import deque
queue = deque(zip(strong_y, strong_x))
# 8邻域方向
directions = [(-1, -1), (-1, 0), (-1, 1),
(0, -1), (0, 1),
(1, -1), (1, 0), (1, 1)]
# 已访问标记
visited = np.zeros_like(nms, dtype=bool)
visited[strong_edges] = True
# BFS遍历
while queue:
y, x = queue.popleft()
for dy, dx in directions:
ny, nx = y + dy, x + dx
# 检查边界
if ny < 0 or ny >= nms.shape[0] or nx < 0 or nx >= nms.shape[1]:
continue
# 如果是弱边缘且未访问过
if weak_edges[ny, nx] and not visited[ny, nx]:
visited[ny, nx] = True
result[ny, nx] = 255 # 标记为边缘
queue.append((ny, nx)) # 加入队列继续搜索
return result
# 应用改进的滞后边界跟踪
my_canny = hysteresis_edge_tracking(NMS)
# 可视化结果
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.imshow(img)
plt.title("Original Image")
plt.axis('off')
plt.subplot(1, 3, 2)
plt.imshow(canny, cmap='gray')
plt.title("OpenCV Canny")
plt.axis('off')
plt.subplot(1, 3, 3)
plt.imshow(my_canny, cmap='gray')
plt.title("My Improved Canny")
plt.axis('off')
plt.tight_layout()
plt.show()
python使用OpenCV
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 1. 读取图像并转换为灰度图
img = cv2.imread('your_image.jpg', cv2.IMREAD_GRAYSCALE)# 直接以灰度模式读取
# 2. 应用Canny边缘检测
# cv2.Canny() 函数内部已经包含了高斯模糊、计算梯度、NMS、双阈值和滞后跟踪所有步骤。
# 参数：
# image: 输入图像
# threshold1: 低阈值 (Tlow)
# threshold2: 高阈值 (Thigh)
# apertureSize: Sobel算子的孔径大小（默认3）
# L2gradient: 是否使用更精确的L2范数计算梯度强度（√(Gx²+Gy²))。默认为False，使用L1范数(|Gx|+|Gy|)
edges = cv2.Canny(image=img, threshold1=100, threshold2=200)
# 3. 显示结果
plt.figure(figsize=(12,6))
plt.subplot(121)
plt.imshow(img, cmap='gray')
plt.title('Original Grayscale Image')
plt.axis('off')
plt.subplot(122)
plt.imshow(edges, cmap='gray')
plt.title('Canny Edge Detection')
plt.axis('off')
plt.tight_layout()
plt.show()
三 区域生长
概述
区域生长是一种经典的图像分割算法，它基于"相似性"原则，从预先定义的种子点开始，逐步将具有相似特性的相邻像素合并到同一区域中。这种方法特别适用于具有均匀区域的图像分割。
操作步骤
1 选择种子点
确定区域生长的起始点，可以手动选择也可以根据某一特定规则自动选择，常见的自动选择方法有
（1）选择图像中灰度值最大/最小的点（2）使用边缘检测结果中的特定点（3）基于图像直方图的极值点
2 定义生长准则
确定哪些像素可以被合并到当前区域，常用准则有
灰度级准则：像素灰度值与区域平均灰度的差值小于阈值T
颜色准则：对于彩色图像，使用颜色距离（如欧氏距离）
纹理准则：基于局部纹理特征的相似性
梯度准则：考虑边缘信息，避免跨边缘生长
数学表达：
对于像素p和区域R，如果满足：|I(p) - μ(R)| ≤ T
其中I(p)是像素p的灰度值，μ(R)是区域R的平均灰度值，T是预设阈值
3 确定生长停止条件
没有更多的像素满足生长准则、区域面积达到预设上限、区域增长速率低于某个阈值
4 实现生长过程
实际执行区域合并过程
- 初始化一个空队列或栈，用于存储待检查的像素
- 将种子点加入队列，并标记为已访问
- 当队列不为空时：
a. 从队列中取出一个像素
b. 检查该像素的所有未访问邻域像素
c. 如果邻域像素满足生长准则，则：
将其标记为当前区域
将其加入队列
更新区域统计信息（如平均灰度）
- 重复步骤3直到队列为空
区域生长的基本思想是将具有相似性质的像素集合起来构成区域。具体先对每个需要分割的区域找一个种子像素作为生长起点，然后将种子像素和周围邻域中与种子像素有相同或相似性质的像素（根据某种事先确定的生长或相似准则来判定）合并到种子像素所在的区域中。将这些新像素当作新的种子继续上面的过程，直到没有满足条件的像素可被包括进来。这样一个区域就生长成了。
例：
图(a)为原始图像，数字表示像素的灰度。以灰度值为8的像素为初始的生长点，记为f(i,j)。在8邻域内，生长准则是待测点灰度值与生长点灰度值相差为1或0。那么图(b)是第一次区域生长后，f(i-1,j)、f(i,j-1)、f(i,j+1)和生长点灰度值相差都是1，因而被合并。图(c)是第二次生长后，f(i+1,j)被合并。图(d)为第三次生长后，f(i+1,j-1)、f(i+2,j)被合并，至此，已经不存在满足生长准则的像素点，生长停止。
代码：
基于灰度图像的简单区域生长算法
import numpy as np
import matplotlib.pyplot as plt
import cv2
from collections import deque
def region_growing(img, seed_point, threshold=10):
"""
区域生长算法实现
参数:
img: 输入灰度图像
seed_point: 种子点坐标 (row, col)
threshold: 生长阈值
返回:
region: 分割出的区域(二值图像)
"""
# 获取图像尺寸
height, width = img.shape
# 初始化输出区域(初始全为0)
region = np.zeros_like(img, dtype=np.uint8)
# 初始化访问标记矩阵
visited = np.zeros_like(img, dtype=bool)
# 定义8邻域方向
directions = [(-1, -1), (-1, 0), (-1, 1),
(0, -1), (0, 1),
(1, -1), (1, 0), (1, 1)]
# 初始化队列
queue = deque()
queue.append(seed_point)
# 标记种子点已访问
visited[seed_point] = True
region[seed_point] = 255
# 初始区域统计
region_mean = float(img[seed_point])
region_pixels = 1
# 开始生长
while queue:
# 取出当前像素
current_point = queue.popleft()
y, x = current_point
# 检查所有邻域像素
for dy, dx in directions:
ny, nx = y + dy, x + dx
# 确保邻域像素在图像范围内
if 0 <= ny < height and 0 <= nx < width:
# 如果邻域像素未被访问
if not visited[ny, nx]:
# 计算当前像素与区域平均值的差异
diff = abs(int(img[ny, nx]) - region_mean)
# 如果差异小于阈值，则添加到区域
if diff < threshold:
# 标记为已访问和区域部分
visited[ny, nx] = True
region[ny, nx] = 255
# 更新区域统计
region_mean = (region_mean * region_pixels + img[ny, nx]) / (region_pixels + 1)
region_pixels += 1
# 将新像素加入队列
queue.append((ny, nx))
return region
# 示例使用
def main():
# 读取图像并转换为灰度
img = cv2.imread('./Images/lena.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 手动选择种子点 (需要根据图像内容调整)
seed_point = (200, 200) # (行, 列)
# 应用区域生长
segmented_region = region_growing(gray, seed_point, threshold=15)
# 可视化结果
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.imshow(gray, cmap='gray')
plt.title('Original Image')
plt.plot(seed_point[1], seed_point[0], 'ro') # 标记种子点
plt.axis('off')
plt.subplot(1, 3, 2)
plt.imshow(segmented_region, cmap='gray')
plt.title('Region Growing Result')
plt.axis('off')
# 将分割结果叠加到原图
colored_region = img.copy()
colored_region[segmented_region == 255] = [255, 0, 0] # 将分割区域标记为红色
plt.subplot(1, 3, 3)
plt.imshow(cv2.cvtColor(colored_region, cv2.COLOR_BGR2RGB))
plt.title('Segmentation Overlay')
plt.axis('off')
plt.tight_layout()
plt.show()
if __name__ == "__main__":
main()
改进代码
- 多种子点区域生长
import numpy as np
from collections import deque
def multi_seed_region_growing(img, seed_points, threshold=10, max_area=10000, min_area=10):
"""
多种子点区域生长 - 改进版本
参数:
img: 输入灰度图像 (numpy数组)
seed_points: 种子点列表 [(row1, col1), (row2, col2), ...]
threshold: 生长阈值 (像素值与区域平均值的最大允许差异)
max_area: 单个区域的最大允许面积 (防止过度生长)
min_area: 单个区域的最小允许面积 (过滤太小的区域)
返回:
segmented_img: 分割后的图像，不同区域用不同标签值表示
region_stats: 每个区域的统计信息字典
"""
# 参数验证
if not isinstance(img, np.ndarray) or len(img.shape) != 2:
raise ValueError("输入图像必须是二维numpy数组")
if not seed_points:
raise ValueError("必须提供至少一个种子点")
height, width = img.shape
segmented_img = np.zeros_like(img, dtype=np.uint16) # 使用16位以支持更多区域
visited = np.zeros_like(img, dtype=bool)
# 存储每个区域的统计信息
region_stats = {}
# 为每个区域分配唯一的标签值
region_label = 1 # 从1开始，0保留给背景
for seed in seed_points:
# 检查种子点是否有效
if (seed[0] < 0 or seed[0] >= height or
seed[1] < 0 or seed[1] >= width or
visited[seed]):
continue
# 初始化当前区域
queue = deque([seed])
visited[seed] = True
segmented_img[seed] = region_label
# 初始化区域统计
region_pixels = [img[seed]]
region_sum = float(img[seed])
region_pixel_count = 1
# 区域边界像素（用于后续分析）
boundary_pixels = set()
# 使用BFS进行区域生长
while queue and region_pixel_count < max_area:
y, x = queue.popleft()
# 检查8邻域
for dy in [-1, 0, 1]:
for dx in [-1, 0, 1]:
if dy == 0 and dx == 0:
continue
ny, nx = y + dy, x + dx
# 检查边界和访问状态
if (0 <= ny < height and 0 <= nx < width and not visited[ny, nx]):
# 计算与区域平均值的差异
current_mean = region_sum / region_pixel_count
diff = abs(int(img[ny, nx]) - current_mean)
if diff < threshold:
# 添加到当前区域
visited[ny, nx] = True
segmented_img[ny, nx] = region_label
queue.append((ny, nx))
# 更新区域统计
region_pixels.append(img[ny, nx])
region_sum += img[ny, nx]
region_pixel_count += 1
else:
# 记录边界像素
boundary_pixels.add((ny, nx))
# 过滤太小的区域
if region_pixel_count < min_area:
# 移除太小的区域
segmented_img[segmented_img == region_label] = 0
visited[segmented_img == region_label] = False
continue
# 记录区域统计信息
region_stats[region_label] = {
'area': region_pixel_count,
'mean_intensity': region_sum / region_pixel_count,
'min_intensity': np.min(region_pixels),
'max_intensity': np.max(region_pixels),
'std_intensity': np.std(region_pixels),
'boundary_pixels': list(boundary_pixels),
'seed_point': seed
}
# 增加标签值
region_label += 1
return segmented_img, region_stats
# 辅助函数：可视化分割结果
def visualize_segmentation(original_img, segmented_img, region_stats):
"""
可视化分割结果
参数:
original_img: 原始图像
segmented_img: 分割结果图像
region_stats: 区域统计信息
返回:
显示分割结果的matplotlib图像
"""
import matplotlib.pyplot as plt
from matplotlib import colors
# 创建彩色映射
cmap = plt.cm.get_cmap('tab20', np.max(segmented_img) + 1)
plt.figure(figsize=(15, 5))
# 原始图像
plt.subplot(1, 3, 1)
plt.imshow(original_img, cmap='gray')
plt.title('Original Image')
plt.axis('off')
# 分割结果
plt.subplot(1, 3, 2)
plt.imshow(segmented_img, cmap=cmap)
plt.title('Segmentation Result')
plt.axis('off')
# 叠加显示
plt.subplot(1, 3, 3)
plt.imshow(original_img, cmap='gray')
# 为每个区域绘制边界
for label, stats in region_stats.items():
if label == 0: # 跳过背景
continue
# 绘制边界点
boundary_y, boundary_x = zip(*stats['boundary_pixels'])
plt.scatter(boundary_x, boundary_y, s=1, c='red', alpha=0.5)
# 标记种子点
seed_y, seed_x = stats['seed_point']
plt.scatter(seed_x, seed_y, s=30, c='blue', marker='x')
plt.title('Segmentation Overlay')
plt.axis('off')
plt.tight_layout()
plt.show()
return plt.gcf()
# 示例使用
def example_usage():
# 生成示例图像（在实际应用中应从文件读取）
height, width = 200, 200
img = np.random.rand(height, width) * 255
img = img.astype(np.uint8)
# 手动创建一些种子点（在实际应用中应自动检测或手动选择）
seed_points = [(50, 50), (100, 100), (150, 150)]
# 应用多种子点区域生长
segmented_img, region_stats = multi_seed_region_growing(
img, seed_points, threshold=20, max_area=1000, min_area=50
)
# 可视化结果
visualize_segmentation(img, segmented_img, region_stats)
# 打印区域统计信息
print("Region Statistics:")
for label, stats in region_stats.items():
print(f"Region {label}:")
print(f" Area: {stats['area']} pixels")
print(f" Mean Intensity: {stats['mean_intensity']:.2f}")
print(f" Seed Point: {stats['seed_point']}")
if __name__ == "__main__":
example_usage()
- 自适应阈值区域生长
import numpy as np
import cv2
from collections import deque
import matplotlib.pyplot as plt
# 您提供的自适应区域生长函数
def adaptive_region_growing(img, seed_point, initial_threshold=10, adaptive_factor=0.1):
"""
自适应阈值区域生长
参数:
img: 输入灰度图像
seed_point: 种子点
initial_threshold: 初始阈值
adaptive_factor: 自适应因子，控制阈值如何随区域标准差变化
返回:
region: 分割出的区域
"""
height, width = img.shape
region = np.zeros_like(img, dtype=np.uint8)
visited = np.zeros_like(img, dtype=bool)
queue = deque([seed_point])
visited[seed_point] = True
region[seed_point] = 255
# 初始化区域统计
region_pixels = [img[seed_point]]
region_mean = np.mean(region_pixels)
region_std = np.std(region_pixels)
directions = [(-1, -1), (-1, 0), (-1, 1),
(0, -1), (0, 1),
(1, -1), (1, 0), (1, 1)]
while queue:
y, x = queue.popleft()
for dy, dx in directions:
ny, nx = y + dy, x + dx
if 0 <= ny < height and 0 <= nx < width and not visited[ny, nx]:
# 计算自适应阈值
adaptive_threshold = initial_threshold + adaptive_factor * region_std
diff = abs(int(img[ny, nx]) - region_mean)
if diff < adaptive_threshold:
visited[ny, nx] = True
region[ny, nx] = 255
queue.append((ny, nx))
# 更新区域统计
region_pixels.append(img[ny, nx])
region_mean = np.mean(region_pixels)
region_std = np.std(region_pixels)
return region
# 主程序
def main():
# 1. 读取图像
img = cv2.imread('./Images/lena.jpg')
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
# 2. 选择种子点（这里使用手动选择方法）
def select_seed_point_interactive(image):
clone = image.copy()
seed_point = [None]
def mouse_callback(event, x, y, flags, param):
if event == cv2.EVENT_LBUTTONDOWN:
seed_point[0] = (y, x) # 转换为(row, col)格式
cv2.circle(clone, (x, y), 5, (0, 255, 0), -1)
cv2.imshow("Select Seed Point", clone)
print(f"Selected seed point: {seed_point[0]}")
cv2.namedWindow("Select Seed Point")
cv2.setMouseCallback("Select Seed Point", mouse_callback)
while True:
cv2.imshow("Select Seed Point", clone)
key = cv2.waitKey(1) & 0xFF
if key == 13 or seed_point[0] is not None: # 按Enter键或选择点后退出
break
cv2.destroyAllWindows()
return seed_point[0]
seed_point = select_seed_point_interactive(img)
if seed_point is None:
print("No seed point selected. Using default.")
seed_point = (100, 100) # 默认种子点
# 3. 应用自适应区域生长
segmented_region = adaptive_region_growing(gray, seed_point,
initial_threshold=15,
adaptive_factor=0.2)
# 4. 可视化结果
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
plt.title('Original Image')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.imshow(gray, cmap='gray')
plt.title('Grayscale Image')
plt.plot(seed_point[1], seed_point[0], 'ro') # 标记种子点
plt.axis('off')
plt.subplot(1, 3, 3)
plt.imshow(segmented_region, cmap='gray')
plt.title('Segmented Region')
plt.axis('off')
plt.tight_layout()
plt.show()
# 5. 可选：将分割结果叠加到原图
result_img = img.copy()
result_img[segmented_region == 255] = [255, 0, 0] # 将分割区域标记为红色
cv2.imshow("Segmentation Result", result_img)
cv2.waitKey(0)
cv2.destroyAllWindows()
if __name__ == "__main__":
main()
四 分水岭算法
概述
传统分水岭算法（watershed）的基本思想是将图像看做为地理学的拓扑地貌，将图像想象成一个立体的地形表面。图像中的每一个像素的灰度值作为该点的海拔高度，每一个局部极小值及它影响的区域被称为集水盆，则集水盆的边界就形成了分水岭。分水岭的概念和形成能够通过模拟浸入的过程来说明，假设在每个区域极小值的位置刺穿一个小孔，从孔中向该区域不断注水，则地势较低的区域会首先被淹没，在浸入加深的过程中，每一个局部极小值的影响域都处在慢慢扩展过程中，当水淹没到一定高度时两个集水盆会相连合并为同一个区域，此时应在两个集水盆汇合处构筑“大坝”来阻止区域合并，最终形成分水岭，这些“大坝”的边界对应于分水岭的分割线，同时也是由分水岭算法提取出的连续边界轮廓。要实现分水岭方法，可分为以下几步进行：读入图像，求取图像的边界，应用分割算法。
但是由于图像中的噪声或任何其他不规则性，这种方法会造成过度分割的结果。所以OpenCV实现了一种基于标记的分水岭算法，可以指定哪些是所有要合并的山谷点，哪些不是。这是一种交互式图像分割。我们所做的是为我们知道的对象给出不同的标签。用一种颜色（或强度）标记我们确定是前景或对象的区域，用另一种颜色标记我们确定是背景或非对象的区域，最后用0标记我们不确定的区域。这就是我们的标记。然后应用分水岭算法。然后我们的标记将使用我们给出的标签进行更新，并且对象的边界将具有-1的值。
流程
- 图像预处理（灰度化 + 去噪）。
- 二值化（确定前景区域）。
- 形态学操作（去噪、增强前景）。
- 距离变换（计算前景区域的中心）。
- 标记前景、背景、未知区域。
- 应用分水岭算法，获取分割边界。
- 可视化结果。
代码
import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读取图像
image = cv2.imread('D:\\resource\\filter\\coins.jpeg')
# 转换为灰度图像
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# 进行高斯模糊以减少噪声
blurred = cv2.GaussianBlur(gray, (5, 5), 0)
# 使用OTSU方法进行二值化
_, thresh = cv2.threshold(blurred, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
# 进行形态学开运算以去除小的噪声点
kernel = np.ones((3, 3), np.uint8)
opening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel, iterations=2)
# 膨胀操作得到背景区域
sure_bg = cv2.dilate(opening, kernel, iterations=3)
# 距离变换得到前景区域
dist_transform = cv2.distanceTransform(opening, cv2.DIST_L2, 5)
_, sure_fg = cv2.threshold(dist_transform, 0.7 * dist_transform.max(), 255, 0)
# 找到未知区域
sure_fg = np.uint8(sure_fg)
unknown = cv2.subtract(sure_bg, sure_fg)
# 标记连通区域
_, markers = cv2.connectedComponents(sure_fg)
# 所有标记加1，确保背景标记为1而不是0
markers = markers + 1
# 未知区域标记为0
markers[unknown == 255] = 0
# 使用分水岭算法进行分割
markers = cv2.watershed(image, markers)
image[markers == -1] = [0, 255, 0]
# 显示结果
plt.subplot(121), plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
plt.title('Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(markers, cmap='jet')
plt.title('Markers'), plt.xticks([]), plt.yticks([])
plt.show()
代码说明
读取图像并转换为灰度图像：使用cv2.imread读取图像，然后使用cv2.cvtColor将其转换为灰度图像。
高斯模糊和二值化：使用cv2.GaussianBlur进行高斯模糊以减少噪声，然后使用 Otsu's 方法进行二值化。
形态学开运算和膨胀操作：使用cv2.morphologyEx进行形态学开运算以去除小的噪声点，然后使用cv2.dilate进行膨胀操作得到背景区域。
距离变换和前景区域提取：使用cv2.distanceTransform进行距离变换，然后使用cv2.threshold提取前景区域。
找到未知区域：通过背景区域减去前景区域得到未知区域。
标记连通区域：使用cv2.connectedComponents标记连通区域，并将所有标记加 1，确保背景标记为 1 而不是 0。
使用分水岭算法进行分割：使用cv2.watershed进行分割，并将分割边界标记为 - 1。
显示结果：使用matplotlib显示分割后的图像和标记图
注意：分水岭算法容易受到噪声的影响，因此在进行分割之前需要进行适当的预处理，如高斯模糊或者形态学操作。