OpenCV combat | Hessian matrix and its application in blood vessel enhancement

Click “Xiaobai learns vision” above, and choose to add “Star” or “Stick“

Heavy dry goods, delivered as soon as possible

Terminology Explanation

– Since the code in this article is based on the OpenCV basic library, the words “OpenCV implementation” are added to the title.

– Due to the two-dimensional nature of the image, all “Hessian matrix” in the following refers specifically to “two-dimensional Hessian matrix”.

Hessian matrix and other related theoretical foundations

There are a lot of basic theories here, you can go through it first, and then go back to deepen your understanding when reading the code, which is better.

1. Origin and definition of Hessian matrix

It can be seen from advanced mathematics knowledge that if the unary function f(x)

In

There is any order derivative in a certain neighborhood of the point, then

in The Taylor expansion at the point is:

Where ,

Binary function in The Taylor expansion at the point is:

in

Write the above expansion in matrix form, then:

that is

in:

is at point The Hessian matrix at . It is created by the function In The square matrix of the second order partial derivatives at the points. We generally express it as:

We can also abbreviate it as:

2. Scale space theory of digital image processing

The basic idea of the scale space theory is to introduce a parameter regarded as a scale in the image information processing model, obtain the scale space representation sequence under multiple scales by continuously changing the scale parameter, extract the main contour of the scale space from these sequences, and The main contour is used as a feature vector to realize edge and corner detection and feature extraction at different resolutions.

The characteristic of the scale space theory is that it incorporates the traditional single-scale image information processing technology into the dynamic analysis framework with changing scales, which makes it easier to obtain the essential characteristics of the image. The blurring degree of each scale image in the scale space gradually increases, which can simulate the formation process of the target on the retina when the distance from the target is close to far.

An intuitive understanding of scale space theory: We humans can recognize the same object at different distances.

The Gaussian convolution kernel is the only linear kernel that implements scale transformation (the Gaussian function can be called the greatest and most important function). The scale space of an image can be defined as:

The symbol “*” represents the convolution operation.

σ is a scale space factor, and the smaller the value, the less the image is smoothed, and the corresponding scale is smaller. The large scale corresponds to the general features of the image, and the small scale corresponds to the detailed features of the image.

3. Hessian simplification algorithm based on scale theory

For a two-dimensional image, the IHessian matrix describes the two-dimensional derivative of each pixel in the main direction as:

According to the scale space theory, the second order derivative can be obtained by convolution of the image with the Gaussian function, for example, at the point (x, y) there is

Where is the scale is the Gaussian function.

If we accept this theory, then we can get the corollary:

Do partial derivative operation directly on the whole image = Convolute the original image with a specific Gaussian kernel

Based on this inference, the Hessian operation on the graph will greatly reduce the amount of calculation while improving the operation speed.

4. The solution method of the eigenvalue of the Hessian matrix

First recall the undergraduate knowledge, and find the eigenvalues of the second-order matrix according to the definition:

By definition, for a matrix A, its eigenvalues satisfy

|λE-A|=0

in

E is a second order diagonal matrix

(1 0)
 (0 1)

We denote A as

|a b|
|c d|

Then λE-A|

= (λ-a)(λ-d)-bc

= λ^2-(a + d)λ + ad-bc = 0
This is a quadratic equation in one variable, and it can be calculated that there are two solutions, which are
λ1=(a + d + √((a-d)^2 + 4bc))/2
λ2=(a + d-√((a-d)^2 + 4bc))/2

From the previous information, we have obtained the definition of the Hessian matrix:

According to the definition of the eigenvalue of a two-dimensional matrix: “Assume that A is a square matrix of order n, if there is a number m and a non-zero n-dimensional column vector x, so that Ax=mx holds true, then m is said to be a characteristic value of the matrix A (characteristic value ) or eigenvalue (eigenvalue)”, the equation can be obtained

And according to the characteristics of the image, we can get

=

Introduce the above equation to get the solution of the eigenvalue of Hessian:

Keep this conclusion in mind, we will refer to it again in the code section.

5. Image properties of eigenvalues of Hessian matrix

A Hessian matrix can be decomposed into two eigenvalues and defined eigenvectors. and

Among them, the largest absolute eigenvalue indicates the largest local grayscale change, and its The eigenvector represents its direction, which can be considered as the tangent direction; and the smaller one represents the vertical direction, which is the normal direction.

This diagram is a good illustration of the concept of tangents and normals.

These will be utilized in the algorithm below.

6. Gaussian equation and second derivative

The Gaussian function was mentioned earlier, and some knowledge is added here, which is useful below.

The Gaussian distribution is the normal distribution, and its graph is as follows:

Two-dimensional Gaussian distribution, its graph is as follows:

According to the definition, we can obtain the content.

The first partial derivative of a 2D Gaussian function:

The second partial derivative of a two-dimensional Gaussian function

The inference from the original function to the second-order partial derivative here is all undergraduate knowledge. It is recommended that you push it yourself. It is very simple. It is very helpful for an in-depth understanding of this problem. We will mention it again in the code section later.

Principle of “vascular enhancement” algorithm (Frangi algorithm)

The ujkHessian matrix and its eigenvalues can well describe the information of common geometric shapes, and we will use it for blood vessel enhancement; the simplified algorithm of the Hessian matrix will provide us with a possible method for coding. We are mainly based on the most famous “Frangi Filtering” paper.

Frangi A F, Niessen W J, Vincken K L, et al. “Multiscale vessel enhancement filtering” [C]//International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer Berlin Heidelberg, 1998: 130-137.

1. Understanding blood vessels and their enhancement

In the collected pictures, blood vessels should present a “tubular network” structure, which provides a basic basis for our image processing. The “Frangi” algorithm enhanced above is itself used to identify pipelines.

2. Fundamentals of Frangi’s thesis

Based on the “accelerated algorithm” we explained earlier, the blood vessels are first subjected to Gaussian filtering at multiple scales, then the second derivative of each pixel is calculated to construct the Hessian matrix, and two eigenvalues are calculated (this place is implemented in the code time has skill).

Although we have obtained the Hessian matrix and its eigenvalues, we can already see the enhanced effect from the image, but this is not enough. next

The obtained eigenvalues are brought into the pre-established blood vessel similarity function to obtain the filter responses at different scales.

When the scale matches the local structure, the maximum filter response is calculated to determine whether the current pixel belongs to the vascular structure.

In order to obtain the enhanced effect as much as possible, the “multi-scale” superposition method is adopted in the paper. Specifically, different convolution kernels are used to process at the same time to obtain multiple processing effects, and then “colorize” the results. The parts with better effect are superimposed.

3. Advantages and disadvantages of Frangi paper

This method obtains an effective blood vessel enhancement method, but it can be seen that the algorithm introduces many factors that need to be considered and defined; at the same time, the larger and more floating-point operations are difficult to run in real time on the embedded system; about ” The definition of “vessel similarity function” lacks theoretical basis, and is more like the result of the definition. Finally, the algorithm cannot directly segment blood vessels, so this step is often used in the preprocessing stage of blood vessel segmentation.

It’s hard to figure it out just by looking at the theory. Here we continue to understand while explaining the code.

Write code for concrete implementation

Let’s start to explain the specific code, and I will pay for the entire runnable project to the end of the article. In the implementation process, we refer to libfrangi

https://ntnu-bioopt.github.io/software/libfrangi.html

The high-quality code provided is explained, and I made the necessary simplification and annotation in the process.

It must be said that this code largely refers to frangi’s matlab version code in terms of content and style

https://www.mathworks.com/matlabcentral/fileexchange/24409-hessian-based-frangi-vesselness-filter,

If you are also good at matlab, it is recommended to learn by comparison.

1. Project structure

First, let’s start with the structure, so that if you have a certain foundation, you can do your own research first, and then compare my explanation, which is helpful for learning.

This project is very simple, except for main.cpp, the frangi algorithm is encapsulated into frangi.h + frangi.cpp, and provided directly in the form of functions.

1int main(int argc, char *argv[]) {
 2 //Set Frangi with default parameters
 3 frangi2d_opts_t opts;
 4 frangi2d_createopts( &opts);
 5 //Read the picture and process it
 6 Mat input_img = imread("E:/template/51.bmp", 0);
 7 Mat input_img_fl;
 8 //Convert to single channel, floating point operation
 9 input_img.convertTo(input_img_fl, CV_32FC1);
10 //Processing
11 Mat vesselness, scale, angles;
12 frangi2d(input_img_fl, vesselness, scale, angles, opts);
13 //Display the result
14 vesselness.convertTo(vesselness, CV_8UC1, 255);
15 scale.convertTo(scale, CV_8UC1, 255);
16 angles.convertTo(angles, CV_8UC1, 255);
17 imshow("result", vesselness);
18}

And main.cpp is as simple as possible, except for the necessary image format conversion, the main process is encapsulated in

frangi2d(input_img_fl, vesselness, scale, angles, opts);

Open frangi.h, we can see the following:

1/
 2//Frangi filter//
 3/
 4//frangi filtering main process
 5 void frangi2d(const cv::Mat & amp;src, cv::Mat & amp;J, cv::Mat & amp;scale, cv::Mat & amp;directions, frangi2d_opts_t opts);
 6
 7//Helper functions//
 8
 9//Calculate the Hessian matrix with parameter sigma on src, save to Dxx, Dxy and Dyy.
10void frangi2d_hessian(const cv::Mat & amp;src, cv::Mat & amp;Dxx, cv::Mat & amp;Dxy, cv::Mat & amp;Dyy, float sigma);
11//Parameter setting (sigma_start = 3, sigma_end = 7, sigma_step = 1, BetaOne = 1.6, BetaTwo 0.08)
12void frangi2d_createopts(frangi2d_opts_t *opts);
13//Calculate eigenvalues from Dxx, Dxy, Dyy. Save results to lambda1, lambda2, Ix, Iy.
14 void frangi2_eig2image(const cv::Mat & amp;Dxx, const cv::Mat & amp;Dxy, const cv::Mat & amp;Dyy, cv::Mat & amp;lambda1, cv::Mat & amp;lambda2 , cv::Mat & amp;Ix, cv::Mat & amp;Iy);

What remains is the process of one main function + three Helper functions.

2. Calculate the Hessian matrix

Let’s look at the frangi2d_hessian function. As the comment states, it is the specific implementation of the Hessian operation:

1//Calculate Hessian matrix with parameter sigma on src, save to Dxx, Dxy and Dyy.
2void frangi2d_hessian(const cv::Mat & amp;src, cv::Mat & amp;Dxx, cv::Mat & amp;Dxy, cv::Mat & amp;Dyy, float sigma);

One of the more difficult to understand is this paragraph, which seems to have done more numerical calculations

1for (int x = -round(3*sigma); x <= round(3*sigma); x + + ){
 2j=0;
 3 for (int y = -round(3*sigma); y <= round(3*sigma); y + + ){
 4 kern_xx_f[i*n_kern_y + j] = 1.0f/(2.0f*M_PI*sigma*sigma*sigma*sigma) * (x*x/(sigma*sigma) - 1) * exp(-(x*x + y*y)/(2.0f*sigma*sigma));
 5 kern_xy_f[i*n_kern_y + j] = 1.0f/(2.0f*M_PI*sigma*sigma*sigma*sigma*sigma*sigma)*(x*y)*exp(-(x*x + y*y) /(2.0f*sigma*sigma));
 6j++;
 7}
 8i++;
 9    }
10 for (int j=0; j < n_kern_y; j ++ ){
11 for (int i=0; i < n_kern_x; i ++ ){
12 kern_yy_f[j*n_kern_x + i] = kern_xx_f[i*n_kern_x + j];
13}
14 }

It looks rather strange. Here, because the code is relatively long and difficult to understand, it is actually calculating the convolution kernel of Dxx and so on. First of all, let’s recall an important conclusion obtained in the theoretical stage:

Directly perform partial derivative operation on the whole image = perform convolution on the original image with a specific Gaussian kernel

Here we are first to calculate the “specific Gaussian kernel”. Then we recall the second-order partial derivative of the two-dimensional Gaussian function introduced at that time

So what does it look like translated into code? matlab version

1//Generate convolution kernel
2 //DGaussxx = 1/(2*pi*Sigma^4) * (X.^2/Sigma^2 - 1) .* exp(-(X.^2 + Y.^2)/(2*Sigma ^2));
3 //DGaussxy = 1/(2*pi*Sigma^6) * (X .* Y) .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
4 //DGaussyy = DGaussxx';

c++ (opencv) version

1cv::Mat tmp1 = 1/(2*PI*Sigma*Sigma*Sigma*Sigma)*(EWM(X,X)/(Sigma*Sigma)-1);
2cv::Mat tmp2 = -1*((EWM(X,X) + EWM(Y,Y))/(2*Sigma*Sigma));
3exp(tmp2,tmp2);

So you can understand the meaning of this code. Next we need to use these 3 specific convolution kernels for convolution

flip(Mat(n_kern_y, n_kern_x, CV_32FC1, kern_xx_f), kern_xx, -1);

In this place, about how OpenCV’s filter2D function implements convolution, and what is the difference between it and Matlab, you can refer to “Actually Comparing the Relationship Between Filter2D and Imfilter”

3. Calculation of Hessian eigenvalues

Let’s recall the conclusion we got earlier:

Then you can understand the code here:

1//calculate eigenvectors of J, v1 and v2
 2Mat tmp, tmp2;
 3tmp2 = Dxx - Dyy;
 4sqrt(tmp2.mul(tmp2) + 4*Dxy.mul(Dxy), tmp);
 5Mat v2x = 2*Dxy;
 6Mat v2y = Dyy - Dxx + tmp;
 7//normalize
 8Mat mag;
 9sqrt((v2x.mul(v2x) + v2y.mul(v2y)), mag);
10Mat v2xtmp = v2x.mul(1.0f/mag);
11v2xtmp.copyTo(v2x, mag != 0);
12Mat v2ytmp = v2y.mul(1.0f/mag);
13v2ytmp.copyTo(v2y, mag != 0);
14//eigenvectors are orthogonal
15Mat v1x, v1y;
16v2y.copyTo(v1x);
17v1x = -1*v1x;
18v2x.copyTo(v1y);
19//compute eigenvalues
20Mat mu1 = 0.5*(Dxx + Dyy + tmp);
21Mat mu2 = 0.5*(Dxx + Dyy - tmp);

Basically, it translates mathematical formulas into codes. Note that .mul is the “dot multiplication” method in OpenCV.

Notice:

1tmp2 = Dxx - Dyy;
2sqrt(tmp2.mul(tmp2) + 4*Dxy.mul(Dxy), tmp);

temp =

4.frangi2d main process

Below we focus on explaining

void frangi2d(const Mat & amp;src, Mat & amp;maxVals, Mat & amp;whatScale, Mat & amp;outAngles, frangi2d_opts_t opts)

function, which is the main process of the whole Frangi algorithm

1for (float sigma = opts.sigma_start; sigma <= opts.sigma_end; sigma + = opts.sigma_step){……}
2
3//Multi-scale superposition
4for (int i=1; i < ALLfiltered. size(); i ++ ){
5 maxVals = max(maxVals, ALLfiltered[i]);
6 whatScale.setTo(sigma, ALLfiltered[i] == maxVals);
7 ALLangles[i].copyTo(outAngles, ALLfiltered[i] == maxVals);
8 sigma += opts.sigma_step;
9}

The outermost frame of the program tells us that the whole program is a superposition of multiple operations (scales).

1//According to the definition of the paper, correct the results
2Dxx = Dxx*sigma*sigma;
3Dyy = Dyy*sigma*sigma;
4Dxy = Dxy*sigma*sigma;
5
6//calculate (abs sorted) eigenvalues and vectors
7Mat lambda1, lambda2, Ix, Iy;
8frangi2_eig2image(Dxx, Dxy, Dyy, lambda1, lambda2, Ix, Iy);

In each calculation, first calculate the value of , the code for The variables are called lambda1, lambda2.

1//According to the definition, calculate the "vessel enhancement response function"
 2lambda2.setTo(nextafterf(0, 1), lambda2 == 0);
 3Mat Rb = lambda1.mul(1.0/lambda2);
 4Rb = Rb.mul(Rb);
 5Mat S2 = lambda1.mul(lambda1) + lambda2.mul(lambda2);
 6Mat tmp1, tmp2;
 7exp(-Rb/beta, tmp1);
 8exp(-S2/c, tmp2);
 9//Save the single-scale result
10Mat Ifiltered = tmp1.mul(Mat::ones(src.rows, src.cols, src.type()) - tmp2);

Here is the most valuable part of the Frangi calculation: “vascular response function”, its mathematical representation is:

Among them, and C are hyperparameters, they are all in front of the code has been based on The input parameters are transformed.

float beta = 2*opts.BetaOne*opts.BetaOne;
float c = 2*opts.BetaTwo*opts.BetaTwo;

The reference values are given in the paper, and the variable names in the code are all corresponding. You can also view it with modifications.

In the later literature, other “vascular enhancement” functions also appeared, such as LiQiang

and GVF

These are easy to implement based on the previously calculated eigenvalues.

4. References:

11.Hessian matrix and its application in image
 2https://blog.csdn.net/lwzkiller/article/details/55050275
 32. Literature review of blood vessel segmentation technology
 4https://blog.csdn.net/u013102349/article/details/53819314
 53. "Crack Segmentation Method of CT Sequence Images Based on Hessian Matrix and Entropy", Wang Huiqian et al., Journal of Instrumentation, August 2016
 64. Baidu Encyclopedia
 7https://baike.baidu.com/item/Hessian Matrix/2248782?fr=aladdin
 85. Scale Space Theory of Digital Image Processing https://blog.csdn.net/TJC3014583925/article/details/88836485>
 96. Implementation and Research of CLAHE
10 https://www.cnblogs.com/jsxyhelu/p/6435601.html
117. "Actually compare the relationship between filter2D and imfilter"
12 https://www.cnblogs.com/jsxyhelu/p/6597544.html
138. "Visual Image: Gaussian Equation and Its Derivatives"
14https://jingyan.baidu.com/album/5bbb5a1bedf94413eba179ba.html?picindex=2

All code:

1 link: https://pan.baidu.com/s/1Q0bbqkJJSHqhN7Htg3dhzQ
2 Extraction code: gzmg

Finally, thank you for studying so far. If you have gained one-tenth of the insight in understanding the content I expressed when I wrote it, then you must be a lucky reader.

Good news!

Xiaobai learns visual knowledge planet

Open to the outside world

Download 1: OpenCV-Contrib extension module Chinese version tutorial

Reply in the backstage of the "Xiaobai Xue Vision" public account: Chinese tutorial of extension module, you can download the first Chinese version of OpenCV extension module tutorial on the whole network, covering extension module installation, SFM algorithm, stereo vision, target tracking, biological vision, super Resolution processing and more than 20 chapters.


Download 2: Python visual combat project 52 lectures
Reply in the background of the "Xiaobai Xue Vision" public account: Python visual combat project, you can download including image segmentation, mask detection, lane line detection, vehicle counting, adding eyeliner, license plate recognition, character recognition, emotion detection, text content extraction, Facial recognition and other 31 visual combat projects help fast school computer vision.


Download 3: OpenCV Practical Project 20 Lectures
Reply in the background of the "Xiaobai Xue Vision" public account: 20 lectures on OpenCV practical projects, you can download 20 practical projects based on OpenCV to achieve advanced OpenCV learning.


Exchange group

Welcome to join the reader group of the official account to communicate with peers. Currently, there are WeChat groups such as SLAM, 3D vision, sensors, automatic driving, computational photography, detection, segmentation, recognition, medical imaging, GAN, algorithm competitions (will be gradually subdivided in the future), Please scan the WeChat ID below to join the group, and remark: "nickname + school/company + research direction", for example: "Zhang San + Shanghai Jiaotong University + visual SLAM". Please make notes according to the format, otherwise it will not be approved. After the addition is successful, it will be invited to enter the relevant WeChat group according to the research direction. Please do not send advertisements in the group, otherwise you will be asked to leave the group, thank you for understanding~