OpenCV如何通过梯度结构张量进行各向异性图像分割(69)

愚梦者

1973人浏览 · 2024-05-04 07:04:03

愚梦者 · 2024-05-04 07:04:03 发布

目标

在本教程中，您将学习：

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

理论

注意

该解释基于书籍[134]，[27]和[281]。[306]中给出了梯度结构张量的良好物理解释。另外，您可以参考维基百科页面结构张量。

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

什么是梯度结构张量？

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

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

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

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

其中张量的分量m[] 是数学期望的符号（我们可以将此操作视为窗口 w 中的平均值）Zx和 Zy 是图像 Z 相对于x 和y 的偏导数。

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

其中 $\lambda$ 1 最大特征值， $\lambda$ 2 - 最小特征值。

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

各向异性图像的方向：

一致性：

相干性范围从 0 到 1。对于理想的局部方向 $\lambda$ 2= 0， $\lambda$ 1 > 0）它是 1，对于各向同性灰度值结构（ $\lambda$ 1= $\lambda$ 2> 0）它是零。

C++源代码

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

#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; // 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;
 }
 
 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 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)
}

解释

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

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 定义窗口大小：

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 定义方向范围：

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

最后，我们结合阈值结果：

 Mat imgBin;
 imgBin = imgCoherencyBin & imgOrientationBin;

Result

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

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

下面您可以看到细分结果：

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

引用

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

参考文献：

1、《Anisotropic image segmentation by a gradient structure tensor》---Karpushin Vladislav

DAMO开发者矩阵

DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐

卡尔曼滤波及其应用

卡尔曼滤波把不准的运动模型和不准的传感器合成当前最可信的状态，并同时给出不确定度。每拍先预测再更新：预测用模型前推并把协方差P撑大；更新用测量纠偏。增益K由P与测量噪声R算出，不是超参；要调的是过程噪声Q和R。P是此刻有多不确定，Q是模型每走一步还要再添多少不确定。机器人中常见于GPS与IMU融合、动态障碍跟踪。

DAMO开发者矩阵

基于蚁群优化算法的直流电机模糊PID控制（Matlab实现）

直流电机因其调速性能优异、响应速度快等特点，在工业控制、电动汽车、机器人等领域有着广泛应用。然而，传统PID控制存在参数固定、难以应对非线性及环境扰动等问题。模糊PID控制虽通过引入模糊逻辑实现了参数自适应调整，但在复杂场景下仍面临规则库构建与参数优化的难题。本文将蚁群优化算法（Ant Colony Optimization, ACO）引入直流电机模糊PID控制系统中，利用其全局寻优能力对模糊PI