【影象處理】海森矩陣（Hessian Matrix）及一個用例（影象增強）
前言
Hessian Matrix(海森矩陣)在影象處理中有廣泛的應用，比如邊緣檢測、特徵點檢測等。而海森矩陣本身也包含了大量的數學知識，例如泰勒展開、多元函式求導、矩陣、特徵值等。寫這篇部落格的目的，就是想從原理和實現上講一講Hessian Matrix，肯定有不足的地方，希望大家批評指正。
泰勒展開及海森矩陣
將一個一元函式f(x)在x0處進行泰勒展開，可以得到以下公式。
其中餘項為皮亞諾餘項。
其中二階導數的部分對映到二維以及多維空間就是Hessian Matrix。在二維影象中，假設影象畫素值關於座標（x, y）的函式是f(x, y)，那麼將f(x+dx，y+dy)在f(x0, y0)處展開,得到如下式子；
如果將這個式子用矩陣表示，並且捨去餘項，則式子會是下面這個樣子。
上面等式右邊的第三項中的第二個矩陣就是二維空間中的海森矩陣了；從而有了一個結論，海森矩陣就是空間中一點處的二階導數。進而推廣開來，多維空間中的海森矩陣可以表示為。
海森矩陣的意義
眾所周知，二階導數表示的導數的變化規律，如果函式是一條曲線，且曲線存在二階導數，那麼二階導數表示的是曲線的曲率，曲率越大，曲線越是彎曲。以此類推，多維空間中的一個點的二階導數就表示該點梯度下降的快慢。以二維影象為例，一階導數是影象灰度變化即灰度梯度，二階導數就是灰度梯度變化程度，二階導數越大灰度變化越不具有線性性（這裡有一點繞口了，意思就是這裡灰度梯度改變越大，不是線性的梯度）。
但是在二維影象中，海森矩陣是二維正定矩陣，有兩個特徵值和對應的兩個特徵向量。兩個特徵值表示出了影象在兩個特徵向量所指方向上影象變化的各向異性。如果利用特徵向量與特徵值構成一個橢圓，那麼這個橢圓就標註出了影象變化的各向異性。那麼在二維影象中，什麼樣的結構最具各項同性，又是什麼樣的結構的各向異性更強呢？很顯然，圓具有最強的各項同性，線性越強的結構越具有各向異性。如下圖；
注：圖中箭頭的方向不一定正確，我只是隨意標註
且特徵值應該具有如下特性；
|
λ1 |
λ2 |
影象特徵 |
|
-High |
-High |
斑點結構（前景為亮） |
|
+High |
+High |
斑點結構（前景為暗） |
|
Low |
-High |
線性結構（前景為亮） |
|
Low |
+High |
線性結構（前景為暗） |
有一個常見的表格也把不同空間特徵與特徵值之間的對應關係描述的很直觀；
海森矩陣的應用
海森矩陣的應用很多，我在此不一一列舉。上文中提到了矩陣的特徵值與特徵向量構成的橢圓表現出了影象的各向異性，這種各項異性影象處理就得到了應用。以二維影象為例，影象中的點性結構具有各項同性，而線性結構具有各向異性。因此我們可以利用海森矩陣對影象中的線性結構進行增強，濾去點狀的結構和噪聲點。同樣，也可以用於找出影象中的點狀結構，濾除其他資訊。
當然，我們在使用海森矩陣時，不需要把影象進行泰勒展開，我們只需要直接求取矩陣中的元素即可。一般，對數字影象進行二階求導使用的是以下方法；
但是這種方法魯棒性很差，容易受到影象中區域性訊號的干擾，甚至可以說，這種求導方式是存在爭議的。因為這一點的二階導數也可以採用如下方法表示；
除了上面兩種表示方法外，二階導數還可以採用其他方式表示，而且往往不同的方法求得的值不同，因為這種方法只把包含其自身在內的三個點的資訊囊括進去，資訊量不足。因此，提倡大家換一種方法。根據線性尺度空間理論（LOG），對一個函式求導，等於函式與高斯函式導數的卷積。式子如下；
由於高斯模板可以將周圍一矩形範圍內所有的點的資訊都包含進來，這樣就不會有誤差。所以利用影象求取hessian矩陣中的元素時，將影象與高斯函式的二階導數做卷積即可。在此，為大家提供高斯函式的二階偏導。
在編寫程式時，我們只需要事先將影象分別於三個模板進行卷積，生成三種偏導數的“圖”，然後每次根據需要索引對應位置的偏導數即可。
程式碼
需要說明的是，我在程式碼中運用了一些OpenCV的函式；
///---------------------------///
//---海森矩陣二維影象增強
//---潘正宇---2018.03.27
#include <iostream>
#include <vector>
#include<opencv2/opencv.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<opencv2/imgproc/imgproc.hpp>
#include <map>
#define STEP 6
#define ABS(X) ((X)>0? X:(-(X)))
#define PI 3.1415926
using namespace std;
using namespace cv;
int main()
{
Mat srcImage = imread("G:\\部落格\\影象處理\\hessian\\hessian_matrix\\Multiscale vessel enhancement filtering1.bmp");
if (srcImage.empty())
{
cout << "影象未被讀入";
system("pause");
return 0;
}
if (srcImage.channels() != 1)
{
cvtColor(srcImage, srcImage, CV_RGB2GRAY);
}
int width = srcImage.cols;
int height = srcImage.rows;
Mat outImage(height, width, CV_8UC1,Scalar::all(0));
int W = 5;
float sigma = 01;
Mat xxGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
Mat xyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
Mat yyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
//構建高斯二階偏導數模板
for (int i = -W; i <= W;i++)
{
for (int j = -W; j <= W; j++)
{
xxGauKernel.at<float>(i + W, j + W) = (1 - (i*i) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
yyGauKernel.at<float>(i + W, j + W) = (1 - (j*j) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
xyGauKernel.at<float>(i + W, j + W) = ((i*j))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * PI*pow(sigma, 6)));
}
}
for (int i = 0; i < (2 * W + 1); i++)
{
for (int j = 0; j < (2 * W + 1); j++)
{
cout << xxGauKernel.at<float>(i, j) << " ";
}
cout << endl;
}
Mat xxDerivae(height, width, CV_32FC1, Scalar::all(0));
Mat yyDerivae(height, width, CV_32FC1, Scalar::all(0));
Mat xyDerivae(height, width, CV_32FC1, Scalar::all(0));
//影象與高斯二階偏導數模板進行卷積
filter2D(srcImage, xxDerivae, xxDerivae.depth(), xxGauKernel);
filter2D(srcImage, yyDerivae, yyDerivae.depth(), yyGauKernel);
filter2D(srcImage, xyDerivae, xyDerivae.depth(), xyGauKernel);
for (int h = 0; h < height; h++)
{
for (int w = 0; w < width; w++)
{
//map<int, float> best_step;
/* int HLx = h - STEP; if (HLx < 0){ HLx = 0; }
int HUx = h + STEP; if (HUx >= height){ HUx = height - 1; }
int WLy = w - STEP; if (WLy < 0){ WLy = 0; }
int WUy = w + STEP; if (WUy >= width){ WUy = width - 1; }
float fxx = srcImage.at<uchar>(h, WUy) + srcImage.at<uchar>(h, WLy) - 2 * srcImage.at<uchar>(h, w);
float fyy = srcImage.at<uchar>(HLx, w) + srcImage.at<uchar>(HUx, w) - 2 * srcImage.at<uchar>(h, w);
float fxy = 0.25*(srcImage.at<uchar>(HUx, WUy) + srcImage.at<uchar>(HLx, WLy) - srcImage.at<uchar>(HUx, WLy) - srcImage.at<uchar>(HLx, WUy));*/
float fxx = xxDerivae.at<float>(h, w);
float fyy = yyDerivae.at<float>(h, w);
float fxy = xyDerivae.at<float>(h, w);
float myArray[2][2] = { { fxx, fxy }, { fxy, fyy } }; //構建矩陣，求取特徵值
Mat Array(2, 2, CV_32FC1, myArray);
Mat eValue;
Mat eVector;
eigen(Array, eValue, eVector); //矩陣是降序排列的
float a1 = eValue.at<float>(0, 0);
float a2 = eValue.at<float>(1, 0);
if ((a1>0) && (ABS(a1)>(1+ ABS(a2)))) //根據特徵向量判斷線性結構
{
outImage.at<uchar>(h, w) = pow((ABS(a1) - ABS(a2)), 4);
//outImage.at<uchar>(h, w) = pow((ABS(a1) / ABS(a2))*(ABS(a1) - ABS(a2)), 1.5);
}
}
}
//----------做一個閉操作
Mat element = getStructuringElement(MORPH_RECT, Size(3, 2));
morphologyEx(outImage, outImage, MORPH_CLOSE, element);
imwrite("temp.bmp", outImage);
imshow("[原始圖]", outImage);
waitKey(0);
system("pause");
return 0;
}輸出結果
左側是原圖，中間是增強後的結果，右側是將增強後的結果做了二值化的結果圖；
結果分析
從結果來看，影象中的大部分的線性結構都被增強了，但是有一些細微的結構並未被增強太多，而且有些粗的結構中出現了空洞，其實這都與求導視窗的大小有關，求導視窗太小，很多粗的結構會出現中空的現象，因為中心區域被認為是點結構了；求導視窗太大，就容易出現細微結構丟失的情況。此外，高斯模板的方差選取也影響了偏導數的大小。其實，這種方式是使用一個方差為s 的高斯核的二階導數生成一個探測核，用於測量導數方向範圍內（-s，s）內外區域之間的對比度。
但是同一影象中，線性結構的粗細肯定是不同的，同樣的視窗大小是無法全部適用的，針對上面的問題，有人提出了多模板的方法。即對一個點用多種尺度的高斯模板進行卷積，然後選擇各向異性最強的結果作為該點的輸出。值得一試。
回過頭來看，根據海森矩陣，還可以確定一張影象中的角點部分，即前面表格中提到的兩個特徵值的絕對值都較大的情況。其實這就是Harris 角點檢測的主要思想。
最後
這篇部落格準備了一段時間了，主要是公式的原因，我想使用markdown編輯器以及LaTex公式編輯器編輯公式，但是沒搞定，還需要琢磨一下，等我琢磨好了，就把這篇文章中的圖片格式的公式換掉。
參考文獻：
Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A. (1998) Multiscale vessel enhancement filtering. In: Wells W.M., Colchester A., Delp S. (eds) Medical Image Computing and Computer-Assisted Intervention — MICCAI’98. MICCAI 1998. Lecture Notes in Computer Science, vol 1496. Springer, Berlin, Heidelberg