【opencv 450 Image Processing】Anisotropic image segmentation by a gradient structure tensor

Anisotropic image segmentation by a gradient structure tensor

梯度结构张量的各向异性图像分割

Goal

在本教程中，您将学习：

梯度结构张量是什么

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

如何通过梯度结构张量分割具有单个局部方向的各向异性图像

Theory

笔记

解释基于书籍 [120]、[22] 和 [260]。 [284] 中给出了梯度结构张量的良好物理解释。 此外，您可以参考维基百科页面结构张量https://en.wikipedia.org/wiki/Structure_tensor 。

此页面上的各向异性图像是真实世界的图像。

What is the gradient structure tensor?

在数学上，梯度结构张量（也称为二阶矩矩阵、二阶矩张量、惯性张量等）是从函数的梯度导出的矩阵。 它总结了点的指定邻域中梯度的主要方向，以及这些方向的相干程度（coherency相干性）。 梯度结构张量广泛应用于图像处理和计算机视觉中，用于 2D/3D 图像分割、运动检测、自适应滤波、局部图像特征检测等。

各向异性图像的重要特征包括局部各向异性的方向和相干性。 在本文中，我们将展示如何估计方向和相干性，以及如何通过梯度结构张量分割具有单个局部方向的各向异性图像。

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

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

 -张量的分量

M[] 是数学期望的符号（我们可以将此操作视为在窗口 w 中求平均），Zx 和 Zy 是图像 Z 关于 x 和 y 的偏导数。

张量的特征值可以在下面的公式中找到：

 其中 λ1 - 最大特征值，λ2 - 最小特征值。

How to estimate orientation and coherency of an anisotropic image by gradient structure tensor?

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

各向异性图像的方向：

 相干性：Coherency:

 相干性范围从 0 到 1。对于理想的局部方向 (λ2 = 0, λ1 > 0)，它是 1，对于各向同性灰度值结构 (λ1 = λ2 > 0)，它是零。

Source code

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

/**
* @brief You will learn how to segment an anisotropic image with a single local orientation by a gradient structure tensor (GST)
* 您将学习如何通过梯度结构张量 (GST) 分割具有单个局部方向的各向异性图像
* @author Karpushin Vladislav, karpushin@ngs.ru, https://github.com/VladKarpushin
*/

#include <iostream>
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"

using namespace cv;
using namespace std;

//! [calcGST_proto]
void calcGST(const Mat& inputImg, Mat& imgCoherencyOut, Mat& imgOrientationOut, int w);
//! [calcGST_proto]

int main()
{
    int W = 52;             // 窗口大小为 WxW
    double C_Thr = 0.43;    // 相干性阈值 threshold for coherency
    int LowThr = 35;        // 方向阈值1threshold1 for orientation, it ranges from 0 to 180
    int HighThr = 57;       // 方向阈值2threshold2 for orientation, it ranges from 0 to 180

    Mat imgIn = imread("input.jpg", IMREAD_GRAYSCALE);
    if (imgIn.empty()) //检查图像是否加载 check whether the image is loaded or not
    {
        cout << "ERROR : Image cannot be loaded..!!" << endl;
        return -1;
    }

    //! [main_extra]
    //! [main]
    Mat imgCoherency, imgOrientation;//相干性图，角度图
    calcGST(imgIn, imgCoherency, imgOrientation, W);//gradient structure tensor 

    //! [thresholding]
    Mat imgCoherencyBin;
    imgCoherencyBin = imgCoherency > C_Thr;//相干性图的阈值处理
    Mat imgOrientationBin;
    inRange(imgOrientation, Scalar(LowThr), Scalar(HighThr), imgOrientationBin);//方向图阈值处理
    //! [thresholding]

    //! [combining]
    Mat imgBin;
    imgBin = imgCoherencyBin & imgOrientationBin;//组合相干性图和方向图
    //! [combining]
    //! [main]

    normalize(imgCoherency, imgCoherency, 0, 255, NORM_MINMAX);//归一化相干性图  用于显示
    normalize(imgOrientation, imgOrientation, 0, 255, NORM_MINMAX);//归一化方向图  用于显示

    imwrite("result.jpg", 0.5*(imgIn + imgBin));//显示原始图和分割结果的加权组合
    imwrite("Coherency.jpg", imgCoherency);
    imwrite("Orientation.jpg", imgOrientation);
    //! [main_extra]
    return 0;
}
//! [calcGST]
//! [calcJ_header]
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);//https://zhuanlan.zhihu.com/p/40491339 
    Sobel(img, imgDiffY, CV_32F, 0, 1, 3);//图像在X方法与Y方向梯度图像
    multiply(imgDiffX, imgDiffY, imgDiffXY);// ZxZy    最终图像梯度
    //! [calcJ_header]

    Mat imgDiffXX, imgDiffYY;
    multiply(imgDiffX, imgDiffX, imgDiffXX);// ZxZx
    multiply(imgDiffY, imgDiffY, imgDiffYY);// ZyZy

    Mat J11, J22, J12;      // J11, J22 and J12 are GST components
    boxFilter(imgDiffXX, J11, CV_32F, Size(w, w));//在窗口 w 中求和
    boxFilter(imgDiffYY, J22, CV_32F, Size(w, w));
    boxFilter(imgDiffXY, J12, CV_32F, Size(w, w));
    // GST components calculation (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))
    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;      // 梯度结构张量 J的最大特征值biggest eigenvalue
    lambda2 = tmp1 - tmp4;	
    lambda2 = 0.5*lambda2;      // 梯度结构张量 J的最小特征值smallest eigenvalue
    // eigenvalue calculation (stop)

    //相干性计算 Coherency calculation (start)
    // Coherency = (lambda1 - lambda2)/(lambda1 + lambda2)) - 各向异性测量measure of anisotropism
    //Coherency 是各向异性程度（局部方向的一致性） Coherency is anisotropy degree (consistency of local orientation)
    divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut);
    // Coherency calculation (stop)

    // 方位角计算 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)
}
//! [calcGST]

Explanation

各向异性图像分割算法由梯度结构张量计算、方向计算、相干性计算以及方向和相干性阈值处理组成：

An anisotropic image segmentation algorithm consists of a gradient structure tensor calculation, an orientation calculation, a coherency calculation and an orientation and coherency thresholding:

    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;

函数 calcGST() 通过使用梯度结构张量计算方向和相干性。 输入参数 w 定义窗口大小：

A function calcGST() calculates orientation and coherency by using a gradient structure tensor. An input parameter w defines a window size:

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)
    // 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
    // Coherency is anisotropy degree (consistency of local orientation)
    divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut);
    // Coherency calculation (stop)
    // 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)
}

下面的代码将阈值 LowThr 和 HighThr 应用于方向图像，并将阈值 C_Thr 应用于由前一个函数计算的图像相干性。 LowThr 和 HighThr 定义方向范围：

The below code applies a thresholds LowThr and HighThr to image orientation and a threshold C_Thr to image coherency calculated by the previous function. LowThr and HighThr define orientation range:

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

最后我们结合阈值结果：

And finally we combine thresholding results:

Mat imgBin;
imgBin = imgCoherencyBin & imgOrientationBin;

Result

下面你可以看到真实的单向各向异性图像：

Below you can see the real anisotropic image with single direction:

具有单一方向的各向异性图像

Anisotropic image with the single direction

 您可以在下面看到各向异性图像的方向和相干性：

Orientation

 

Coherency

下面你可以看到分割结果：

Segmentation result

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

References

• Structure tensor - structure tensor description on the wikipedia

参考：

图像处理（各项异性）_caojinpei123的博客-CSDN博客_图像各向异性

对图像来说各向异性就是在每个像素点周围四个方向上梯度变化都不一样，滤波的时候我们要考虑图像的各向异性对图像的影响，而各向同性显然是说各个方向的值都一致，常见的图像均值或者高斯均值滤波可以看成是各向同性滤波。

OpenCV图像处理|1.17 Sobel算子 - 知乎 (zhihu.com)

Sobel算子：是离散微分算子（discrete differentiation operator），用来计算图像灰度的近似梯度，梯度越大越有可能是边缘。

Soble算子的功能集合了高斯平滑和微分求导，又被称为一阶微分算子，求导算子，在水平和垂直两个方向上求导，得到的是图像在X方法与Y方向梯度图像。

缺点：比较敏感，容易受影响，要通过高斯模糊（平滑）来降噪。

算子是通过权重不同来扩大差异。

梯度计算：（在两个方向求导，假设被作用图像为 I）

水平变化: 将 I 与一个奇数大小的内核 Gx进行卷积。比如，当内核大小为3时, Gx的计算结果为:

垂直变化: 将:math:I 与一个奇数大小的内核 Gy进行卷积。比如，当内核大小为3时, Gy的计算结果为:

在图像的每一点，结合以上两个结果求出近似梯度：

有时也用下面更简单公式代替，计算速度快：（最终图像梯度）。

盒子(方框)滤波(BoxFilter)原理及C++及Matlab实现_小武~~的博客-CSDN博客_box filter

盒子滤波是一种非常有用的线性滤波，也叫方框滤波，最简单的均值滤波就是盒子滤波归一化的情况。

应用：可以说，一切需要求某个邻域内像素之和的场合，都有盒子滤波的用武之地，比如：均值滤波、引导滤波、计算Haar特征等等。

Opencv之图像滤波：3.方框滤波（cv2.boxFilter）_Justth.的博客-CSDN博客_opencv方框滤波函数

OpenCV还提供了方框滤波方式，与均值滤波的不同在于，方框滤波不会计算像素均值。在均值滤波中，滤波结果的像素值是任意一个点的邻域平均值，等于各邻域像素值之 和除以邻域面积。而在方框滤波中，可以自由选择是否对均值滤波的结果进行归一化，即 可以自由选择滤波结果是邻域像素值之和的平均值，还是邻域像素值之和。

在OpenCV中，实现方框滤波的函数是cv2.boxFilter（），其语法格式为：

        dst=cv2.boxFilter（src,ddepth,ksize,anchor,normalize,borderType）

        式中：

        ● dst是返回值，表示进行方框滤波后得到的处理结果。

        ● src 是需要处理的图像，即原始图像。它能够有任意数量的通道，并能对各个通道独立处理。图像深度应该是CV_8U、CV_16U、CV_16S、CV_32F 或者 CV_64F中的一种。

         ● ddepth是处理结果图像的图像深度，一般使用-1表示与原始图像使用相同的图像深度。

         ● ksize 是滤波核的大小。滤波核大小是指在滤波处理过程中所选择的邻域图像的高 度和宽度。

        ● anchor 是锚点，其默认值是（-1,-1），表示当前计算均值的点位于核的中心点位 置。该值使用默认值即可，在特殊情况下可以指定不同的点作为锚点。

        ● normalize 表示在滤波时是否进行归一化（这里指将计算结果规范化为当前像素值范围内的值）处理，该参数是一个逻辑值，可能为真（值为1）或假（值为0）：

                1.当参数normalize=1时，表示要进行归一化处理，要用邻域像素值的和除以面积。此时方框滤波与均值滤波效果相同。

                2.当参数normalize=0时，表示不需要进行归一化处理，直接使用邻域像素值的和。当 normalize=0时，因为不进行归一化处理，因此滤波得到的值很可能超过当前像素值范围的最大值，从而被截断为最大值。这样，就会得到一幅纯白色的图像。

        ● borderType是边界样式，该值决定了以何种方式处理边界。

        通常情况下，在使用方框滤波函数时，对于参数anchor、normalize和borderType，直接采用其默认值即可。因此，函数cv2.boxFilter（）的常用形式为：

        dst=cv2.boxFilter（src,ddepth,ksize）

phase函数——opencv

phase函数——opencv - 修行君 - 博客园 (cnblogs.com)

cvtColor(left, left, CV_BGR2GRAY);

Sobel(left, grad1, CV_64FC1, 1, 0); //求梯度

Sobel(left, grad2, CV_64FC1, 0, 1);

phase(grad1, grad2, angle, true); //求角度

double min, max;

minMaxLoc(angle, &min, &max);

cout << "min:" << min << " max:" << max << endl;

normalize(angle, angle, 0, 255, NORM_MINMAX); //归一化