Previous tutorial: Demotion blur filter
Next tutorial: Periodic noise removal filter
| Original author | Karpushin Vladislav |
|---|---|
| Compatibility | OpenCV >= 3.0 |
Goals
In this tutorial you will learn
- What is a gradient structure tensor
- How to estimate the direction and consistency of anisotropic images with gradient structure tensors
- How to segment anisotropic image with single local direction by gradient structure tensor
Theory
Comments
This explanation is based on the books [123], [24] and [267]. See [292] for the physical interpretation of the gradient structure tensor. Also, you can refer to the Wikipedia page Structure Tensors.
The anisotropic images on this page are real-world images.
What is a gradient structure tensor?
In mathematics, a gradient structure tensor (also called a second-order matrix, second-order moment tensor, inertia tensor, etc.) is a matrix derived from the gradient of a function. It summarizes the main directions of gradients within a specific neighborhood of a point, and how consistent (coherent) these directions are. Gradient structure tensors are widely used in image processing and computer vision fields, such as two-dimensional/three-dimensional image segmentation, motion detection, adaptive filtering, local image feature detection, etc.
Important features of anisotropic images include the direction and consistency of local anisotropy. This article will show how to estimate orientation and coherence, and how to segment anisotropic images with a single local orientation via gradient structure tensors.
The image’s gradient structure tensor is a 2×2 symmetric matrix. The eigenvectors of a gradient structure tensor represent local directions, while the eigenvalues represent consistency (a measure of anisotropy).
The gradient structure tensor J of image Z can be written as
J
=
[
J
11
J
12
J
12
J
twenty two
]
J = \begin{bmatrix} J_{11} & amp; J_{12} \ J_{12} & amp; J_{22} \end{bmatrix}
J=[J11?J12J12?J22]
in,
J
11
=
M
[
Z
x
2
]
,
J
twenty two
=
M
[
Z
y
2
]
,
J
twenty two
=
M
[
Z
y
2
]
J_{11} = M[Z_{x}^{2}],J_{22} = M[Z_{y}^{2}],J_{22} = M[Z_{y}^{2}]
J11?=M[Zx2?],J22?=M[Zy2?],J22?=M[Zy2?] The components of the tensor, ***M[]*** are the symbols of mathematical expectations (we can use this The operation is considered as averaging in the window w), Zx and Zy are images Zrelative to
x
and
y
x and y
Partial derivatives of x and y.
The eigenvalues of the tensor can be obtained as follows:
λ
1
,
2
=
1
2
[
J
11
+
J
twenty two
±
(
J
11
?
J
twenty two
)
2
+
4
J
12
2
]
\lambda_{1,2} = \frac{1}{2} \left [ J_{11} + J_{22} \pm \sqrt{(J_{11} – J_{22})^ {2} + 4J_{12}^{2}} \right ]
λ1,2?=21?[J11? + J22?±(J11J22?)2 + 4J122?
How to estimate the direction and consistency of an anisotropic image from a gradient structure tensor?
Orientation of anisotropic images:
α
=
0.5
a
r
c
t
g
2
J
12
J
twenty two
?
J
11
\alpha = 0.5arctg\frac{2J_{12}}{J_{22} – J_{11}}
α=0.5arctgJ22J11?2J12
Relevance:
C
=
λ
1
?
λ
2
λ
1
+
λ
2
C = \frac{\lambda_1 – \lambda_2}{\lambda_1 + \lambda_2}
C=λ1? + λ2?λ1λ2
Coherence ranges from 0 to 1. The consistency is 1 for ideal local directions (λ2 = 0, λ1 > 0) and 0 for isotropic gray value structures (λ1 = λ2 > 0).
Source code
You can find the source code in the OpenCV source code repository at 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 & amp; inputImg, Mat & amp; imgCoherencyOut, Mat & amp; imgOrientationOut, int w);
int main()
{<!-- -->
int W = 52; // Window size is WxW
double C_Thr = 0.43; // Consistency threshold
int LowThr = 35; // Direction threshold 1, ranging from 0 to 180
int HighThr = 57; // Direction threshold 2, ranging 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 has been loaded
{<!-- -->
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 & amp; 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 & amp; inputImg, Mat & amp; imgCoherencyOut, Mat & amp; imgOrientationOut, int w)
{<!-- -->
Mat img;
inputImg.convertTo(img, CV_32F);
//GST component 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 components of 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 component calculation (stopped)
// 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; // maximum eigenvalue
lambda2 = tmp1 - tmp4;
lambda2 = 0.5*lambda2; // minimum eigenvalue
// Eigenvalue calculation (stop)
// Consistency calculation (start)
// Consistency = (lambda1 - lambda2)/(lambda1 + lambda2)) - measure of anisotropy
// Consistency is the degree of anisotropy (consistency in local directions)
divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut);
//Coherence calculation (stop)
// direction 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;
// direction angle calculation (stop)
}
Python
import cv2 as cv
import numpy as np
import argparse
W = 52 # The window size is WxW
C_Thr = 0.43 # Consistency threshold
LowThr = 35 # Direction threshold 1, ranging from 0 to 180
HighThr = 57 # Direction threshold 2, ranging from 0 to 180
def calcGST(inputIMG, w)：
img = inputIMG.astype(np.float32)
#GST component 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 composition calculation (stopped)
# 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) # Maximum eigenvalue
lambda2 = 0.5*(tmp1 - tmp4) # Minimum eigenvalue
# Eigenvalue calculation (stop)
# Consistency calculation (start)
# Consistency = (lambda1 - lambda2)/(lambda1 + lambda2)) - measure of anisotropy
# Consistency is the degree of anisotropy (consistency in local directions)
imgCoherencyOut = cv.divide(lambda1 - lambda2, lambda1 + lambda2)
# Consistency calculation (stop)
# Direction 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
# Direction 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)
Description
The anisotropic image segmentation algorithm includes gradient structure tensor calculation, direction calculation, consistency calculation, and direction and consistency threshold calculation:
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 & amp; 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)
The function calcGST() calculates direction and consistency using gradient structure tensors. The input parameter w defines the window size:
C++
void calcGST(const Mat & amp; inputImg, Mat & amp; imgCoherencyOut, Mat & amp; imgOrientationOut, int w)
{<!-- -->
Mat img;
inputImg.convertTo(img, CV_32F);
//GST component 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 components of 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 component calculation (stopped)
// 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; // maximum eigenvalue
lambda2 = tmp1 - tmp4;
lambda2 = 0.5*lambda2; // minimum eigenvalue
// Eigenvalue calculation (stop)
// Consistency calculation (start)
// Consistency = (lambda1 - lambda2)/(lambda1 + lambda2)) - measure of anisotropy
// Consistency is the degree of anisotropy (consistency in local directions)
divide(lambda1 - lambda2, lambda1 + lambda2, imgCoherencyOut);
//Coherence calculation (stop)
// direction 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;
// direction 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
The following code applies thresholds LowThr and HighThr to the image direction and a threshold C_Thr to the image coherence, which are calculated by the previous function. LowThr and HighThr define the direction range:
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)
Finally, we merge the thresholding results:
C++
Mat imgBin; imgBin = imgCoherencyBin & amp; imgOrientationBin;
Python
imgBin = cv.bitwise_and(imgCoherencyBin, imgOrientationBin)
Results
Below is a real unidirectional image of an anisotropic image:
Anisotropic image in one direction
Direction
Coherence
The following is the segmentation result:
Segmentation results
The calculation results are w = 52, C_Thr = 0.43, LowThr = 35, HighThr = 57. We can see that the algorithm only selects areas in a single direction.
References
- Structure Tensor – Structure tensor description on Wikipedia