上一个教程 ： 去动态模糊滤波器

下一个教程 ： 周期性噪声去除滤波器

原作者	Karpushin Vladislav
兼容性	OpenCV >= 3.0

目标

在本教程中，你将学习

• 什么是梯度结构张量
• 如何用梯度结构张量估计各向异性图像的方向和一致性
• 如何通过梯度结构张量分割具有单一局部方向的各向异性图像

理论

注释
本解释基于书籍 [123]、[24] 和 [267]。梯度结构张量的物理解释见 [292]。此外，您还可以参考维基百科页面结构张量。
该页面上的各向异性图像就是现实世界中的图像。

什么是梯度结构张量？

在数学中，梯度结构张量（也称为二阶矩阵、二阶力矩张量、惯性张量等）是由函数梯度导出的矩阵。它概括了某点特定邻域内梯度的主要方向，以及这些方向的一致性（连贯性）程度。梯度结构张量被广泛应用于图像处理和计算机视觉领域，如二维/三维图像分割、运动检测、自适应滤波、局部图像特征检测等。

各向异性图像的重要特征包括局部各向异性的方向和一致性。本文将介绍如何估计方位和一致性，以及如何通过梯度结构张量分割具有单一局部方位的各向异性图像。

图像的梯度结构张量是一个 2x2 对称矩阵。梯度结构张量的特征向量表示局部方向，而特征值则表示一致性（各向异性的度量）。

图像 Z 的梯度结构张量 J 可以写成

$$J=\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{12}& J_{22}\end{array}\right]$$
其中，$$J_{11}=M[Z_x^2],J_{22}=M[Z_y^2],J_{22}=M[Z_y^2]$$ 张量的分量，***M[]***是数学期望的符号（我们可以将此操作视为在窗口 w 中进行平均），Zx 和 Zy 是图像 Z 相对于 $$x和y$$ 的偏导数。

张量的特征值可按下式求得：

$${\lambda }_{1,2}=\frac{1}{2}\left[J_{11}+J_{22}\pm \sqrt{(J_{11}-J_{22}{)}^2+4J_{12}^2}\right]$$
其中 λ1 - 最大特征值，λ2 - 最小特征值。

如何通过梯度结构张量估计各向异性图像的方向和一致性？

各向异性图像的方向：

$$\alpha =0.5arctg\frac{2J_{12}}{J_{22}-J_{11}}$$
相干性：
$$C=\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1+{\lambda }_2}$$
相干性范围为 0 至 1。对于理想的局部方向（λ2 = 0, λ1 > 0），一致性为 1；对于各向同性的灰度值结构（λ1 = λ2 > 0），一致性为 0。

源代码

您可以在 OpenCV 源代码库的 samples/cpp/tutorial_code/ImgProc/anisotropic_image_segmentation/anisotropic_image_segmentation.cpp 中找到源代码。
C++

#include <iostream>
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
using namespace cv;
using namespace std;
void calcGST(const Mat& inputImg, Mat& imgCoherencyOut, Mat& imgOrientationOut, int w);
int main()
{
 int W = 52; // 窗口大小为 WxW
 double C_Thr = 0.43; // 一致性阈值
 int LowThr = 35; // 方向阈值1，范围为 0 至 180
 int HighThr = 57; // 方向阈值 2，范围为 0 至 180
 samples::addSamplesDataSearchSubDirectory("doc/tutorials/imgproc/anisotropic_image_segmentation/images");
 Mat imgIn = imread(samples::findFile("gst_input.jpg"), IMREAD_GRAYSCALE);
 if (imgIn.empty()) //检查图像是否已加载
 {
 cout << "ERROR : Image cannot be loaded..!!" << endl;
 return -1;
 }
 Mat imgCoherency, imgOrientation;
 calcGST(imgIn, imgCoherency, imgOrientation, W);
 Mat imgCoherencyBin;
 imgCoherencyBin = imgCoherency > C_Thr;
 Mat imgOrientationBin;
 inRange(imgOrientation, Scalar(LowThr), Scalar(HighThr), imgOrientationBin);
 Mat imgBin;
 imgBin = imgCoherencyBin & imgOrientationBin;
 normalize(imgCoherency, imgCoherency, 0, 255, NORM_MINMAX, CV_8U);
 normalize(imgOrientation, imgOrientation, 0, 255, NORM_MINMAX, CV_8U);
 imshow("Original", imgIn);
 imshow("Result", 0.5 * (imgIn + imgBin));
 imshow("Coherency", imgCoherency);
 imshow("Orientation", imgOrientation);
 imwrite("result.jpg", 0.5*(imgIn + imgBin));
 imwrite("Coherency.jpg", imgCoherency);
 imwrite("Orientation.jpg", imgOrientation);
 waitKey(0);
 return 0;
}
void calcGST(const Mat& inputImg, Mat& imgCoherencyOut, Mat& imgOrientationOut, int w)
{
 Mat img;
 inputImg.convertTo(img, CV_32F);
 // GST 分量计算（开始）
 // 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 和 J12 是 GST 的组成部分
 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 成分计算（停止）
 // 特征值计算（开始）
 // 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; // 最大特征值
 lambda2 = tmp1 - tmp4；
 lambda2 = 0.5*lambda2; // 最小特征值
 // 特征值计算（停止）
 // 一致性计算（开始）
 // 一致性 = (lambda1 - lambda2)/(lambda1 + lambda2)) - 各向异性的度量
 // 一致性是各向异性程度（局部方向的一致性）
 divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut)；
 // 相干性计算（停止）
 // 方向角计算（开始）
 // 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；
 // 方向角计算（停止）
}

Python

import cv2 as cv
import numpy as np
import argparse
W = 52 # 窗口大小为 WxW
C_Thr = 0.43 # 一致性阈值
LowThr = 35 # 方向阈值 1，范围从 0 到 180
HighThr = 57 # 方向阈值 2，范围在 0 到 180 之间
def calcGST(inputIMG, w)：
 img = inputIMG.astype(np.float32)
 # GST 分量计算（开始）
 # j = (j11 j12; j12 j22) - gst
 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))
 # GST 成分计算（停止）
 # 特征值计算（开始）
 # lambda1 = 0.5*(J11 + J22 + sqrt((J11-J22)^2 + 4*J12^2))
 # lambda2 = 0.5*(J11 + J22 - sqrt((J11-J22)^2 + 4*J12^2))
 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) # 最大特征值
 lambda2 = 0.5*(tmp1 - tmp4) # 最小特征值
 # 特征值计算（停止）
 # 一致性计算（开始）
 # 一致性 = (lambda1 - lambda2)/(lambda1 + lambda2)) - 各向异性的度量
 # 一致性是各向异性程度（局部方向的一致性）
 imgCoherencyOut = cv.divide(lambda1 - lambda2, lambda1 + lambda2)
 # 一致性计算（停止）
 # 方向角计算 (start)
 # tan(2*Alpha) = 2*J12/(J22 - J11)
 # Alpha = 0.5 atan2(2*J12/(J22 - J11))
 imgOrientationOut = cv.phase(J22 - J11, 2.0 * J12, angleInDegrees = True)
 imgOrientationOut = 0.5 * imgOrientationOut
 # 方向角计算（停止）
 return imgCoherencyOut, imgOrientationOut
parser = argparse.ArgumentParser(description='Code for Anisotropic image segmentation tutorial.')
parser.add_argument('-i', '--input', help='Path to input image.', required=True)
args = parser.parse_args()
imgIn = cv.imread(args.input, cv.IMREAD_GRAYSCALE)
if imgIn is None:
 print('Could not open or find the image: {}'.format(args.input))
 exit(0)
imgCoherency, imgOrientation = calcGST(imgIn, W)
_, imgCoherencyBin = cv.threshold(imgCoherency, C_Thr, 255, cv.THRESH_BINARY)
_, imgOrientationBin = cv.threshold(imgOrientation, LowThr, HighThr, cv.THRESH_BINARY)
imgBin = cv.bitwise_and(imgCoherencyBin, imgOrientationBin)
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)

说明

各向异性图像分割算法包括梯度结构张量计算、方向计算、一致性计算以及方向和一致性阈值计算：
C++

 Mat imgCoherency, imgOrientation；
 calcGST(imgIn, imgCoherency, imgOrientation, W)；
 Mat imgCoherencyBin；
 imgCoherencyBin = imgCoherency > C_Thr；
 Mat imgOrientationBin；
 inRange(imgOrientation, Scalar(LowThr), Scalar(HighThr), imgOrientationBin)；
 imgBin.Mat
 imgBin = imgCoherencyBin & imgOrientationBin；

Python

imgCoherency, imgOrientation = calcGST(imgIn, W)
_, imgCoherencyBin = cv.threshold(imgCoherency, C_Thr, 255, cv.THRESH_BINARY)
_, imgOrientationBin = cv.threshold(imgOrientation, LowThr, HighThr, cv.THRESH_BINARY)
imgBin = cv.bitwise_and(imgCoherencyBin, imgOrientationBin)

函数 calcGST() 使用梯度结构张量计算方向和一致性。输入参数 w 定义了窗口大小：
C++

void calcGST(const Mat& inputImg, Mat& imgCoherencyOut, Mat& imgOrientationOut, int w)
{
 Mat img；
 inputImg.convertTo(img, CV_32F)；
 // GST 分量计算（开始）
 // 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 和 J12 是 GST 的组成部分
 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 成分计算（停止）
 // 特征值计算（开始）
 // 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; // 最大特征值
 lambda2 = tmp1 - tmp4；
 lambda2 = 0.5*lambda2; // 最小特征值
 // 特征值计算（停止）
 // 一致性计算（开始）
 // 一致性 = (lambda1 - lambda2)/(lambda1 + lambda2)) - 各向异性的度量
 // 一致性是各向异性程度（局部方向的一致性）
 divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut)；
 // 相干性计算（停止）
 // 方向角计算（开始）
 // 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；
 // 方向角计算（停止）
}

Python

def calcGST(inputIMG, w):
 img = inputIMG.astype(np.float32)
 # GST components calculation (start)
 # J = (J11 J12; J12 J22) - GST
 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))
 # GST components calculations (stop)
 # 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))
 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
 # eigenvalue calculation (stop)
 # Coherency calculation (start)
 # Coherency = (lambda1 - lambda2)/(lambda1 + lambda2)) - measure of anisotropism
 # Coherency is anisotropy degree (consistency of local orientation)
 imgCoherencyOut = cv.divide(lambda1 - lambda2, lambda1 + lambda2)
 # Coherency calculation (stop)
 # orientation angle calculation (start)
 # tan(2*Alpha) = 2*J12/(J22 - J11)
 # Alpha = 0.5 atan2(2*J12/(J22 - J11))
 imgOrientationOut = cv.phase(J22 - J11, 2.0 * J12, angleInDegrees = True)
 imgOrientationOut = 0.5 * imgOrientationOut
 # orientation angle calculation (stop)
 return imgCoherencyOut, imgOrientationOut

下面的代码对图像方向应用了阈值 LowThr 和 HighThr，对图像连贯性应用了阈值 C_Thr，这两个阈值是由前面的函数计算得出的。LowThr 和 HighThr 定义了方向范围：
C++

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

Python

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

最后，我们合并阈值处理结果：
C++

 Mat imgBin；
 imgBin = imgCoherencyBin & imgOrientationBin；

Python

imgBin = cv.bitwise_and(imgCoherencyBin, imgOrientationBin)

结果

下面是各向异性图像的真实单向图像：

单一方向的各向异性图像 下面是各向异性图像的方向和一致性：

方向

相干性

下面是分割结果：

分割结果
计算结果为 w = 52，C_Thr = 0.43，LowThr = 35，HighThr = 57。我们可以看到，算法只选择了单一方向的区域。

参考文献

• 结构张量 - 维基百科上的结构张量描述