透過梯度結構張量進行各向異性影像分割

目錄

• 目標
• 理論
  • 什麼是梯度結構張量？
  • 如何透過梯度結構張量估計各向異性影像的方向和相干性？
• 原始碼
• 解釋
• 結果
• 參考文獻

上一個教程: 運動去模糊濾波器
下一個教程: 週期性噪聲去除濾波器

原始作者	Karpushin Vladislav
相容性	OpenCV >= 3.0

目標

在本教程中，您將學習

• 瞭解什麼是梯度結構張量
• 瞭解如何透過梯度結構張量估計各向異性影像的方向和相干性
• 瞭解如何透過梯度結構張量分割具有單一區域性方向的各向異性影像

理論

注意

本解釋基於書籍 [137], [28] 和 [286]。關於梯度結構張量良好的物理解釋可在 [310] 中找到。此外，您還可以參考維基百科頁面 結構張量。

此頁面上的各向異性影像是真實世界的影像。

什麼是梯度結構張量？

在數學中，梯度結構張量（也稱為二階矩矩陣、二階矩張量、慣性張量等）是一個從函式梯度匯出的矩陣。它總結了在點指定鄰域內梯度的主要方向，以及這些方向的一致性程度（相干性）。梯度結構張量廣泛用於影像處理和計算機視覺，用於2D/3D影像分割、運動檢測、自適應濾波、區域性影像特徵檢測等。

各向異性影像的重要特徵包括區域性各向異性的方向和相干性。本文將展示如何估計方向和相干性，以及如何透過梯度結構張量分割具有單一區域性方向的各向異性影像。

影像的梯度結構張量是一個2x2對稱矩陣。梯度結構張量的特徵向量表示區域性方向，而特徵值則表示相干性（各向異性度量）。

影像 \(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[]\) 是數學期望的符號（我們可以將此操作視為在視窗w中的平均），\(Z_{x}\) 和 \(Z_{y}\) 是影像 \(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 ] \]

其中 \(\lambda_1\) 是最大特徵值，\(\lambda_2\) 是最小特徵值。

如何透過梯度結構張量估計各向異性影像的方向和相干性？

各向異性影像的方向

\[\alpha = 0.5arctg\frac{2J_{12}}{J_{22} - J_{11}}\]

相干性

\[C = \frac{\lambda_1 - \lambda_2}{\lambda_1 + \lambda_2}\]

相干性範圍從0到1。對於理想的區域性方向（\(\lambda_2\) = 0, \(\lambda_1\) > 0），它為1；對於各向同性灰度結構（\(\lambda_1\) = \(\lambda_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)
    # 相干性計算（停止）
 
    # 方向角計算（開始）
    # 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='各向異性影像分割教程程式碼。')
parser.add_argument('-i', '--input', help='輸入影像路徑。', required=True)
args = parser.parse_args()
 
imgIn = cv.imread(args.input, cv.IMREAD_GRAYSCALE)
if imgIn is None
print('無法開啟或找到影像: {}'.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);
 
    Mat imgBin;
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分量計算（開始）
    # 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)
    # 相干性計算（停止）
 
    # 方向角計算（開始）
    # 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

以下程式碼將閾值 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 計算的。我們可以看到演算法只選擇了具有單一方向的區域。

參考文獻

• 結構張量 - 維基百科上的結構張量描述