OpenCV 各向异性图像分割实战：基于梯度结构张量（GST）的方向估计、相干性度量与阈值分割

梯度结构张量（Gradient Structure Tensor，GST）是图像处理与计算机视觉中分析局部纹理方向的核心工具，广泛用于图像分割、运动检测、自适应滤波与局部特征提取。本文基于 OpenCV 官方教程《Anisotropic image segmentation by a gradient structure tensor》（作者 Karpushin Vladislav，适用于 OpenCV ≥ 3.0）展开，结合 OpenCV 仓库中的 C++ 与 Python 示例源码，系统讲解 GST 的数学原理、方向与相干性的估计算法，以及如何在 OpenCV 中用 Sobel、boxFilter、phase、inRange 等算子实现一张"单一局部方向"各向异性图像的分割。读完本文，你将能够独立复现完整算法流程，理解窗口尺寸 W、方向阈值 LowThr/HighThr 与相干性阈值 C_Thr 四个核心参数如何决定分割效果。

补充说明：本页所述的"各向异性图像"指真实世界图像（real world image）。教程原文与示例代码位于 各向异性图像分割教程，C++ 示例见 anisotropic_image_segmentation.cpp，Python 示例见 anisotropic_image_segmentation.py。

一、什么是梯度结构张量（GST）

在数学上，梯度结构张量（又称 second-moment matrix、二阶矩张量、惯性张量等）是由函数的梯度推导出的一个矩阵。它在一个点附近的指定邻域内，汇总梯度的主导方向（predominant directions）以及这些方向相互一致（相干）的程度。各向异性图像的两个重要特征是局部各向异性的方向（orientation）与相干性（coherency），GST 恰能同时刻画这两者：

• 图像的梯度结构张量是一个 2×2 对称矩阵；
• 其特征向量指示局部方向（orientation）；
• 其特征值给出相干性，即各向异性的度量。

1.1 GST 的数学定义

图像 $Z$ 的梯度结构张量 $J$ 写作：

$$J = \begin{bmatrix} J_{11} & J_{12} \\ J_{12} & J_{22} \end{bmatrix}$$

其中各分量定义为：

• $J_{11} = M[Z_x^{2}]$
• $J_{22} = M[Z_y^{2}]$
• $J_{12} = M[Z_x Z_y]$

这里 $M[\cdot]$ 表示数学期望符号（在实际计算中可视为在一个窗口 $w$ 内的平均），$Z_x$ 与 $Z_y$ 分别是图像 $Z$ 对 $x$、$y$ 的偏导数。

1.2 特征值

张量的特征值由下式求得：

$$\lambda_{1,2} = \frac{1}{2} \left[ J_{11} + J_{22} \pm \sqrt{(J_{11} - J_{22})^{2} + 4J_{12}^{2}} \right]$$

其中 $\lambda_1$ 为最大特征值，$\lambda_2$ 为最小特征值。

二、方向与相干性的估计

2.1 方向（Orientation）

各向异性图像的方向角 $\alpha$ 由下式给出：

$$\alpha = 0.5 \cdot \arctan \frac{2J_{12}}{J_{22} - J_{11}}$$

注意系数 $0.5$：由于张量自身具有对称性（旋转 $180°$ 的方向等价），角度周期是 $180°$ 而非 $360°$，因此需要先对 $2\alpha$ 求反正切再减半。

2.2 相干性（Coherency）

相干性 $C$ 定义为：

$$C = \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2}$$

相干性取值范围为 0 到 1：

• 对理想的单一局部方向（$\lambda_2 = 0$，$\lambda_1 > 0$），$C = 1$；
• 对各向同性的灰度结构（$\lambda_1 = \lambda_2 > 0$），$C = 0$。

直观理解：相干性越高，说明局部窗口内的纹理方向越整齐划一；接近 0 则说明局部结构各方向梯度相当、没有明确的主导方向。

三、算法整体流程

从示例源码可以看出，各向异性图像分割算法由以下四步组成（即源码中 main 部分的核心逻辑）：

1. 梯度结构张量计算——得到 $J_{11}$、$J_{22}$、$J_{12}$；
2. 方向计算——由张量分量求出逐像素方向角 $\alpha$；
3. 相干性计算——由特征值组合求出逐像素相干性 $C$；
4. 方向与相干性阈值分割——分别对方向图与相干性图做阈值化，再求交集得到分割掩码。

这四个步骤集中封装在自定义函数 calcGST() 中，其输入参数 w 决定局部平均的窗口大小。下面分别从 C++ 与 Python 两个版本源码逐段拆解实现细节。

四、核心实现剖析（C++ 版本）

完整 C++ 示例位于 samples/cpp/tutorial_code/ImgProc/anisotropic_image_segmentation/anisotropic_image_segmentation.cpp，以下按源码中的标注片段（calcGST、thresholding、combining）进行说明。

4.1 参数定义与图像加载

int W = 52;             // window size is WxW
double C_Thr = 0.43;    // threshold for coherency
int LowThr = 35;        // threshold1 for orientation, it ranges from 0 to 180
int HighThr = 57;       // threshold2 for orientation, it ranges from 0 to 180

samples::addSamplesDataSearchSubDirectory("doc/tutorials/imgproc/anisotropic_image_segmentation/images");
Mat imgIn = imread(samples::findFile("gst_input.jpg"), IMREAD_GRAYSCALE);
if (imgIn.empty()) //check whether the image is loaded or not
{
    cout << "ERROR : Image cannot be loaded..!!" << endl;
    return -1;
}

• W = 52：局部平均窗口取 52×52，属于偏大的平滑窗口，用于聚合足够大的邻域统计量；
• C_Thr = 0.43：相干性阈值；
• LowThr = 35、HighThr = 57：方向阈值，注释中说明方向取值范围为 0~180 度；
• 示例通过 OpenCV 示例数据查找机制 samples::addSamplesDataSearchSubDirectory(...) 与 samples::findFile(...) 定位输入图 gst_input.jpg（实际位于 images 目录）。

4.2 calcGST()：GST 分量计算

void calcGST(const Mat& inputImg, Mat& imgCoherencyOut, Mat& imgOrientationOut, int w)
{
    Mat img;
    inputImg.convertTo(img, CV_32F);

    // GST components calculation (start)
    // J =  (J11 J12; J12 J22) - GST
    Mat imgDiffX, imgDiffY, imgDiffXY;
    Sobel(img, imgDiffX, CV_32F, 1, 0, 3);
    Sobel(img, imgDiffY, CV_32F, 0, 1, 3);
    multiply(imgDiffX, imgDiffY, imgDiffXY);

    Mat imgDiffXX, imgDiffYY;
    multiply(imgDiffX, imgDiffX, imgDiffXX);
    multiply(imgDiffY, imgDiffY, imgDiffYY);

    Mat J11, J22, J12;      // J11, J22 and J12 are GST components
    boxFilter(imgDiffXX, J11, CV_32F, Size(w, w));
    boxFilter(imgDiffYY, J22, CV_32F, Size(w, w));
    boxFilter(imgDiffXY, J12, CV_32F, Size(w, w));
    // GST components calculation (stop)
    ...
}

实现要点：

• 输入图像先 convertTo(CV_32F) 转成浮点型，保证后续梯度与统计计算精度；
• 求导：用 3×3 的 Sobel 分别计算 $Z_x$、$Z_y$，这一步在 imgproc 模块中对应导数算子实现；
• 逐项乘积：multiply 逐元素相乘得到 $Z_xZ_y$、$Z_x^2$、$Z_y^2$ 三个中间图；
• 局部平均：对三个中间图分别做 boxFilter（窗口 Size(w, w)）得到 $J_{11}$、$J_{22}$、$J_{12}$。boxFilter 默认 normalize=true，即每个输出像素等于对应 $w\times w$ 邻域内的算术平均，正好对应数学定义中的 $M[\cdot ]$ 局部期望操作。boxFilter 的实现位于 imgproc 模块的 box_filter.dispatch.cpp。

4.3 特征值与相干性计算

    // eigenvalue calculation (start)
    // lambda1 = 0.5*(J11 + J22 + sqrt((J11-J22)^2 + 4*J12^2))
    // lambda2 = 0.5*(J11 + J22 - sqrt((J11-J22)^2 + 4*J12^2))
    Mat tmp1, tmp2, tmp3, tmp4;
    tmp1 = J11 + J22;
    tmp2 = J11 - J22;
    multiply(tmp2, tmp2, tmp2);
    multiply(J12, J12, tmp3);
    sqrt(tmp2 + 4.0 * tmp3, tmp4);

    Mat lambda1, lambda2;
    lambda1 = tmp1 + tmp4;
    lambda1 = 0.5*lambda1;      // biggest eigenvalue
    lambda2 = tmp1 - tmp4;
    lambda2 = 0.5*lambda2;      // smallest eigenvalue
    // eigenvalue calculation (stop)

    // Coherency calculation (start)
    // Coherency = (lambda1 - lambda2)/(lambda1 + lambda2)) - measure of anisotropism
    divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut);
    // Coherency calculation (stop)

这段代码直接用逐像素矩阵运算复现了前文的解析公式：先算 tmp1 = J11+J22、判别式根号项 tmp4 = sqrt((J11-J22)^2 + 4*J12^2)，进而得到 $\lambda_1$ 与 $\lambda_2$，最后按 $(\lambda_1-\lambda_2)/(\lambda_1+\lambda_2)$ 求出相干性图。需要留意的是，在完全平坦（梯度为零）的区域，特征值均为零，相干性会出现 0/0 的不定情况——从实现结构看，这些像素在后续阈值比较中天然被排除，不会误选为分割目标。

4.4 方向角计算

    // orientation angle calculation (start)
    // tan(2*Alpha) = 2*J12/(J22 - J11)
    // Alpha = 0.5 atan2(2*J12/(J22 - J11))
    phase(J22 - J11, 2.0*J12, imgOrientationOut, true);
    imgOrientationOut = 0.5*imgOrientationOut;
    // orientation angle calculation (stop)

方向计算巧妙地复用了 OpenCV 核心模块的 phase() 函数：它按 x = J22 - J11、y = 2*J12 计算 $\arctan 2(y,x)$，恰好对应公式 $2\alpha =\arctan \frac{2J_{12}}{J_{22}-J_{11}}$。phase 的第四个参数 angleInDegrees = true 表示输出角度单位为度；随后乘以 0.5 还原出方向角 $\alpha$。phase 的函数签名声明位于 core 模块的 core.hpp，默认 angleInDegrees=false 时返回弧度。

4.5 阈值化（thresholding）

Mat imgCoherencyBin;
imgCoherencyBin = imgCoherency > C_Thr;
Mat imgOrientationBin;
inRange(imgOrientation, Scalar(LowThr), Scalar(HighThr), imgOrientationBin);

• 相干性阈值：imgCoherency > C_Thr 做逐像素比较，选出相干性高于 0.43 的"方向高度一致"区域，输出为 0/255 二值掩码；
• 方向阈值：inRange(imgOrientation, LowThr, HighThr) 取出方向角落在 $[35°, 57°]$ 区间内的像素，两个阈值共同框定用户关心的目标方向范围。

4.6 结果合并（combining）

Mat imgBin;
imgBin = imgCoherencyBin & imgOrientationBin;

最终分割掩码是两幅二值图的逐像素逻辑与（AND）：既要方向落在目标区间，又要局部相干性足够高，二者缺一不可。这也是本算法能够"只选出具有单一方向的区域"的关键——仅靠方向阈值或仅靠相干性阈值都无法同时抑制方向不符区域与方向杂乱区域。

合并后，示例将相干性图与方向图 normalize 到 0~255 便于显示，并用 0.5 * (imgIn + imgBin) 将分割掩码半透明叠加回原图，输出 result.jpg，同时通过 imshow 展示原图、方向图、相干性图与叠加结果。

五、Python 版本实现要点

Python 示例位于 samples/python/tutorial_code/imgProc/anisotropic_image_segmentation/anisotropic_image_segmentation.py，算法逻辑与 C++ 完全一致，仅 API 形式不同：

def calcGST(inputIMG, w):
    img = inputIMG.astype(np.float32)

    # GST components calculation (start)
    imgDiffX = cv.Sobel(img, cv.CV_32F, 1, 0, 3)
    imgDiffY = cv.Sobel(img, cv.CV_32F, 0, 1, 3)
    imgDiffXY = cv.multiply(imgDiffX, imgDiffY)

    imgDiffXX = cv.multiply(imgDiffX, imgDiffX)
    imgDiffYY = cv.multiply(imgDiffY, imgDiffY)

    J11 = cv.boxFilter(imgDiffXX, cv.CV_32F, (w, w))
    J22 = cv.boxFilter(imgDiffYY, cv.CV_32F, (w, w))
    J12 = cv.boxFilter(imgDiffXY, cv.CV_32F, (w, w))

    # eigenvalue calculation
    tmp1 = J11 + J22
    tmp2 = J11 - J22
    tmp2 = cv.multiply(tmp2, tmp2)
    tmp3 = cv.multiply(J12, J12)
    tmp4 = np.sqrt(tmp2 + 4.0 * tmp3)

    lambda1 = 0.5*(tmp1 + tmp4)    # biggest eigenvalue
    lambda2 = 0.5*(tmp1 - tmp4)    # smallest eigenvalue

    # Coherency calculation
    imgCoherencyOut = cv.divide(lambda1 - lambda2, lambda1 + lambda2)

    # orientation angle calculation
    imgOrientationOut = cv.phase(J22 - J11, 2.0 * J12, angleInDegrees=True)
    imgOrientationOut = 0.5 * imgOrientationOut

    return imgCoherencyOut, imgOrientationOut

阈值化与合并部分使用 OpenCV Python 封装：

# thresholding
_, imgCoherencyBin = cv.threshold(imgCoherency, C_Thr, 255, cv.THRESH_BINARY)
imgOrientationBin = cv.inRange(imgOrientation, LowThr, HighThr)

# combining
imgBin = cv.bitwise_and(imgCoherencyBin, imgOrientationBin)

Python 版与 C++ 版有三处 API 形态差异需要注意：

• C++ 用表达式 imgCoherency > C_Thr 直接生成二值图，Python 用 cv.threshold(..., cv.THRESH_BINARY) 等价实现；
• 二值合并 C++ 用 & 运算符，Python 用 cv.bitwise_and；
• Python 版通过 argparse 接收 -i/--input 命令行参数指定输入图像路径，因此运行前需要先安装 OpenCV Python 绑定与 numpy。

主流程调用链为：calcGST(imgIn, W) → 阈值化 → 合并 → 归一化显示：

imgCoherency, imgOrientation = calcGST(imgIn, W)
# ... thresholding ...
# ... combining ...
imgCoherency = cv.normalize(imgCoherency, None, alpha=0, beta=1,
                            norm_type=cv.NORM_MINMAX, dtype=cv.CV_32F)
imgOrientation = cv.normalize(imgOrientation, None, alpha=0, beta=1,
                              norm_type=cv.NORM_MINMAX, dtype=cv.CV_32F)
cv.imshow('result.jpg', np.uint8(0.5*(imgIn + imgBin)))
cv.imshow('Coherency.jpg', imgCoherency)
cv.imshow('Orientation.jpg', imgOrientation)
cv.waitKey(0)

Python 归一化时把相干性与方向图缩放到 0~1 区间（与 C++ 缩放到 0~255 等价），直接以浮点图显示。

六、核心参数语义速查表

下表汇总了四个核心参数的含义、推荐初始值与取值范围（均来自教程与示例代码中的注释和默认值）：

参数	含义	示例默认值	取值范围说明
W	GST 局部平均的窗口边长（W×W），用于 boxFilter	52	越大越强调大尺度纹理一致性；过小则受噪点干扰，过大则边界被平滑
C_Thr	相干性阈值，高于该值判定为方向一致区域	0.43	理论范围 [0, 1]（1 为理想单一方向，0 为各向同性）
LowThr	方向角区间下限（度）	35	注释标注方向角范围 0~180 度
HighThr	方向角区间上限（度）	57	与 LowThr 一起框定目标方向 $[LowThr, HighThr]$

调参思路：LowThr/HighThr 描述"期望分割出来的纹理朝哪个方向"；C_Thr 描述"这个方向要多纯粹才值得分割"；W 描述"在多大的局部邻域内统计方向"。教程中给出的 W = 52, C_Thr = 0.43, LowThr = 35, HighThr = 57 是一组针对下方示例图可直接复现结果的组合。

七、运行与复现

7.1 准备示例图像

输入图位于 doc/tutorials/imgproc/anisotropic_image_segmentation/images/gst_input.jpg。C++ 示例启动时会向 OpenCV 示例数据搜索路径注册上述目录，若在自定义环境运行发现图片加载失败，可将该图片复制到当前目录或配置 OPENCV_SAMPLES_DATA_PATH 指向仓库根目录后重试。

7.2 运行 C++ 示例

C++ 示例依赖 opencv_core、opencv_imgproc、opencv_highgui、opencv_imgcodecs 四个模块。一种典型编译方式（假设仓库已按官方流程构建出 OpenCV 库文件并导出 OpenCVConfig.cmake）：

g++ anisotropic_image_segmentation.cpp -o gst_seg \
    $(pkg-config --cflags --libs opencv4)
./gst_seg

程序启动后会依次弹出 "Original"（原图）、"Result"（分割叠加）、"Coherency"（相干性图）、"Orientation"（方向图）四个窗口，并在当前目录写出 result.jpg、Coherency.jpg、Orientation.jpg。

7.3 运行 Python 示例

Python 示例依赖 cv2 与 numpy，通过命令行指定输入图：

python3 anisotropic_image_segmentation.py \
    --input doc/tutorials/imgproc/anisotropic_image_segmentation/images/gst_input.jpg

按下任意键可切换查看叠加结果、相干性图与方向图窗口。

八、运行结果与解读

下面四幅图依次为：输入的各向异性图像（具有单一方向的纹理结构）、算法估计出的方向图、相干性图，以及最终分割结果（以半透明白色叠加在原图上）：

该结果使用 W = 52、C_Thr = 0.43、LowThr = 35、HighThr = 57 计算得到。可以直观看到，算法只保留了图像中具备单一主导方向的区域，而方向偏离目标区间或方向杂乱、不具各向异性的区域被正确排除——这正是 GST 作为"方向+一致性联合描述子"的典型应用价值。

九、扩展讨论

• 应用到 2D/3D 场景：教程所述公式在图像处理与计算机视觉领域有广泛用途，包括 2D/3D 图像分割、运动检测、自适应滤波与局部图像特征检测等。原理基础可进一步参考 Jahne《Computer Vision and Applications》、Bigun《Vision with Direction》、van den Boogaard 关于估计器的专著，以及 Yang 关于结构张量物理含义的论文，教程的引用条目与 @cite 标签一并收录在 OpenCV 文档体系（doc 目录）中。
• 方向周期性的处理：GST 方向角周期为 $180°$（非 $360°$），因此方向阈值需要按 $[0°, 180°)$ 的周期语义设计，示例中 LowThr=35、HighThr=57 正是落在这一区间内。
• 阈值策略推广：若目标图像包含多个局部方向，可将"单一 inRange 方向区间 + 单相干性阈值"的组合推广为多区间、多阈值的逻辑并/或运算；若需软分割，也可将二值化替换为对相干性的连续加权。

参考

• 本教程原文：Anisotropic image segmentation by a gradient structure tensor
• C++ 示例源码：anisotropic_image_segmentation.cpp
• Python 示例源码：anisotropic_image_segmentation.py
• 示例输入与结果图像：images 目录
• 相关 API 声明与实现：cv::phase 见 modules/core/include/opencv2/core.hpp，cv::boxFilter 实现见 modules/imgproc/src/box_filter.dispatch.cpp