如何检测Matplotlib Quiver绘制的向量场中的涡旋及其位置？

问题背景 #

我正在使用Matplotlib的quiver函数在图像上叠加向量场，通过以下代码实现（使用imgT和imgL两张图像获取水平和垂直向量分量）：

gd=1  # sampling 
sy,sx = imgT.shape[:2]
imgTscaled = resize(imgT, (int(sy/gd), int(sx/gd)))
imgLscaled = resize(imgL, (int(sy/gd), int(sx/gd)))
Y, X = np.mgrid[0:sy:gd, 0:sx:gd]
dX=2*(imgTscaled-minT)/(maxT-minT)-1    
dY=2*(imgLscaled-minL)/(maxT-minT)-1  
ax[0].quiver(X, Y, dX, dY, color='#FE7C25', headwidth=5, scale=30)

请问如何从该向量场中检测是否存在涡旋，并确定其位置？

示例图：
-> 涡旋
-> 无涡旋

解决方案 #

1. 计算向量场的旋度（Curl） #

涡旋的核心特征是向量场的旋转趋势，二维向量场的旋度（标量值，对应三维旋度的z分量）可以量化这种旋转强度：

• 对于向量场$(dX, dY)$，旋度计算公式为：
  $$\text{curl} = \frac{\partial dY}{\partial X} - \frac{\partial dX}{\partial Y}$$
• 用NumPy的np.gradient函数快速计算偏导数并得到旋度：

  # 计算dY和dX的梯度（对应偏导数）
  dY_dX, dY_dY = np.gradient(dY)
  dX_dX, dX_dY = np.gradient(dX)
  # 计算旋度
  curl = dY_dX - dX_dY
  

• 旋度的正负对应旋转方向（正为逆时针，负为顺时针），绝对值越大，旋转趋势越强。

2. 阈值筛选与涡旋定位 #

• 先设定旋度阈值（根据数据分布调整，比如取旋度绝对值的前5%，或固定数值），筛选出旋转强度足够的区域：

  # 示例阈值，需根据实际数据调整
  threshold = 0.5
  curl_mask = np.abs(curl) > threshold
  

• 对筛选后的区域做连通域分析，每个独立的连通域对应一个涡旋，连通域的质心就是涡旋的中心位置。可以用OpenCV实现：

  import cv2
  
  # 连通域分析，返回标签、统计信息和质心
  num_labels, labels, stats, centroids = cv2.connectedComponentsWithStats(curl_mask.astype(np.uint8))
  
  # 遍历所有涡旋（跳过背景标签0）
  for i in range(1, num_labels):
      # 还原到原始图像的坐标（如果采样率gd>1）
      centroid_x = centroids[i][0] * gd
      centroid_y = centroids[i][1] * gd
      print(f"涡旋中心位置：({centroid_x:.2f}, {centroid_y:.2f})")
  

3. 可视化验证 #

在原向量场图上标记涡旋中心，直观验证结果：

# 在ax[0]图层绘制涡旋中心，蓝色圆点突出显示
for i in range(1, num_labels):
    centroid_x = centroids[i][0] * gd
    centroid_y = centroids[i][1] * gd
    ax[0].scatter(centroid_x, centroid_y, color='blue', marker='o', s=60, zorder=10)

关键注意事项 #

• 采样率gd会影响旋度计算精度，若gd>1，必须将质心坐标乘以gd还原到原始图像尺寸。
• 阈值需结合旋度的直方图分布调整，避免误检或漏检。
• 若向量场存在噪声，先对dX和dY做高斯模糊（cv2.GaussianBlur），再计算旋度，能有效提升检测准确性。

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