勾配構造テンソルによる異方性画像セグメンテーション

目次

• 目的
• Theory
  • 勾配構造テンソルとは何か?
  • 勾配構造テンソルによって異方性画像の方向とコヒーレンシをどのように推定するか?
• ソースコード
• 解説
• 結果
• 参考文献

前のチュートリアル: モーションブラー除去フィルタ
次のチュートリアル: 周期性ノイズ除去フィルタ

原著者	Karpushin Vladislav
互換性	OpenCV >= 3.0

目的

このチュートリアルでは以下を学ぶ:

• 勾配構造テンソルとは何か
• 勾配構造テンソルによって異方性画像の方向とコヒーレンシをどのように推定するか
• 勾配構造テンソルによって単一の局所方向を持つ異方性画像をどのようにセグメンテーション(領域分割)するか

理論

覚え書き

この解説は書籍 [140]、[29]、[289] に基づいている。勾配構造テンソルの優れた物理的な説明は [313] に示されている。また、Wikipediaのページ Structure tensor も参照できる。

このページの異方性画像は実世界の画像である。

勾配構造テンソルとは何か?

数学において、勾配構造テンソル（2次モーメント行列、2次モーメントテンソル、慣性テンソルなどとも呼ばれる）は、関数の勾配から導かれる行列である。これはある点の指定された近傍における勾配の支配的な方向と、それらの方向がどの程度コヒーレント（コヒーレンシ）であるかを要約する。勾配構造テンソルは、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;             // 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)
}

Python

 
import cv2 as cv
import numpy as np
import argparse
 
W = 52          # window size is WxW
C_Thr = 0.43    # threshold for coherency
LowThr = 35     # threshold1 for orientation, it ranges from 0 to 180
HighThr = 57    # threshold2 for orientation, it ranges from 0 to 180
 
 
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
 
 
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);
 
 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 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)
}

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 で計算された。アルゴリズムが単一方向を持つ領域のみを選択していることがわかる。

参考文献

• Structure tensor - Wikipediaの構造テンソルの説明