【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取

这篇具有很好参考价值的文章主要介绍了【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

Hess矩阵是一个多元函数的二阶偏导数构成的方阵，描述了函数的局部曲率。Hess矩阵经常用在牛顿法中求多元函数的极值问题，将目标函数在某点领域内进行二阶泰勒展开，其中的二阶导数就是Hess矩阵。

海森矩阵的意义

应用在图像中，将图像中在某点领域内进行泰勒展开：
$\ F(x_1+\Delta x) = F(x_1) + J(x_1)^\mathrm{T} {\Delta x} + \frac{1}{2} \Delta x^\mathrm{T}H(x_1) \Delta x \,$
其中 $J(x_1)\,$ 是F(x)在 $x_1\,$ 处的一阶导数（梯度）， $\ H(x1) \,$ 是二阶导数，图像 $x_1\,$ 领域内增量是 $\ \Delta x \,$ ；
求图像 $x_1\,$ 点领域的极值，对上述等式右侧等式关于 $\ \Delta x \,$ 求导，并令求导后等于0，得到关系式：
$\ J(x_1) + H(x_1) \Delta x = 0； H(x_1) \Delta x = -J(x_1)；\,$
那么对于图像处理中，图像 $x_1\,$ 领域内增量 $\ \Delta x = [1,1]^\mathrm{T} \,$ 是像素，那么Hess矩阵的特征向量就是反向一阶梯度；特征值正负表示的是局部领域内函数的凹凸性，大小表示凹凸的程度；特征值绝对值最大的特征向量表示一阶梯度的变化快慢，乘以增量 $\ \Delta x = [1,1]^\mathrm{T}\,$ 像素，表示的一阶梯度；

如下图所示，线条处两个特征值表示了图像在两个特征向量所指的方向上图像变化的各向异性。相似的，圆具有各项同性；线性结构越强的越是各项异性。
【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）

海森矩阵的应用（Steger算法）

【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）
其中 $r_{xx} \,$ 表示图像沿x的二阶偏导，其他类似； $\ \lambda \,$ 则表示特征向量的长度。

Hess矩阵按照变量求二阶偏导后按顺序排列，因各二阶导连续故原函数的混合偏导数全相等，导致其是对称的，这时就可以与正负定矩阵知识联系起来。若是正负定矩阵，二阶导为系数的多项式符号确定（前提在该领域内一阶导单调），证明一阶导有零点，故配合二阶导的符号即可确定是极小还是极大。

如何判断一个矩阵是否是正定的，负定的，还是不定的呢？
一个最常用的方法就是顺序主子式。实对称矩阵为正定矩阵的充要条件是的各顺序主子式都大于零。
另一个判定方法：
矩阵正定的充分必要条件是它的特征值都大于零。矩阵负定的充分必要条件是它的特征值都小于零。否则是不定的。

因此在算法中求解处Hess矩阵特征值后，可以依照图像线某一方向极大值点添加约束条件：
两个约束条件（1、起码有一个特征值小于0；2、其中一个特征值的绝对值 >（另一个特征值的绝对值 + 1））；
另一个约束条件，中心点处图像的灰度值要大于某一阈值(结合图像的特性)。

中心点位置亚像素求解，沿着一阶梯度方向，利用下述泰勒展开获得亚像素的值，求极值让下式一阶导为0，得到偏差值距离值；
$\ F(x_0 + tn_x, y_0 + tn_y) = F(x_0, y_0) + J(x_1)^\mathrm{T} [tn_x, tn_y]^\mathrm{T} + \frac{1}{2}[tnx, tny]*H(x1)*[tn_x, tn_y]^\mathrm{T} \,$ ;
【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）

注：在求上述hessian矩阵之前，根据Steger文章中，设计高斯方差 $\ \sigma \leqslant \frac{w}{\sqrt3} \,$ ；核过大，高斯分布顶部满足条件的点变多，选择中心点的精度不够，提取到的中心线会小范围的波动；核过小，中心点在二阶法线长度局部最大值处，非最小点处。

基于Steger改善算法示例代码

clc; clear; close all;
tic;
%% Extraction of normal binarization center based on Steger algorithm
image = imread(['.\Precision\stair\8\a.bmp']); 
Sz = size(image);

% 自动提取ROI区域阈值
Img_index = zeros(1, 245);
for i = 1:Sz(1)
    for j = 1:Sz(2)
        if image(i, j) > 10
           Img_index(image(i, j) - 10) = Img_index(image(i, j) - 10) + 1;
        end
    end
end
Img_index(245) = 0;
ind_x = Img_index(1:end-1) - Img_index(2:end);
ind_x = imfilter(ind_x, [1 3 6 9 10 9 6 3 1]./48, 'replicate');% ;
for i = 2:245
    if ind_x(i) < 10
       Intensity_threshold = i + 10;%20
       break;
    end
end

% choose the region of interest
BC = roicolor(image, Intensity_threshold, 256);
[r, c] = find(BC == 1);
Range_vertical = lineboundary(r, c);
global num
nn = ones(1, Sz(2));
nn_sz = size(nn(Range_vertical(1, :) ~= 0));
num = sum(Range_vertical(2, :) - Range_vertical(1, :)) + nn_sz(2);

%对所有ROI区域执行平滑操作
[img_filter, Index_c] = filter_gauss_1(image, Range_vertical);
% Guass smoothing： Notices: sigma >= w/sqrt(3);first w = 8; 

DIR = [ [-1; -1], [-1;  0], [-1;  1] ...
        [ 0; -1], [ 0;  0], [ 0;  1]... 
        [ 1; -1], [ 1;  0], [ 1;  1] ];
% 找到对应的中心点与线宽大小
CenterCord = zeros(2, Sz(2)); 
img_2filter = zeros(1, num);
LineW = zeros(1, Sz(2));
for i = 2:Sz(2)-1   
    if i == 1110
        i = 1110;
    end
    % 判断初始检索位置
    y0 = (Range_vertical(2, i-1) + Range_vertical(1, i-1))/2;
    y1 = (Range_vertical(2, i) + Range_vertical(1, i))/2;
    y2 = (Range_vertical(2, i+1) + Range_vertical(1, i+1))/2;
    if (y0==0 || y1==0 || y2==0)|| abs(y1 - y0) > 5 || abs(y1 - y2) > 5
        continue;
    else
        tempInd = find(img_filter(Index_c(i):Index_c(i+1) - 1) > 220);
        if length(tempInd) > 1
            temp = tempInd(1) + floor(length(tempInd)/2);            
        else
            [~, temp] = max(img_filter(Index_c(i):Index_c(i+1) - 1));
        end            
        tempy = Range_vertical(1, i) + temp - 1;        
    end  
    CenterCord(:, i) = [i; tempy];
    
    if image(tempy, i) > 150
       % 计算有限个纵向滤波      
       img_2f_array = Block_filter(img_filter, img_2filter, [i; tempy], Index_c, Range_vertical, gausFilter);    
       for k = [1 2 4 5 6 8]  % 开始检索([1 2 *; 4 5 6; * 8 *]')  
           m = Index_c(i + DIR(1, k)) + tempy + DIR(2, k) - Range_vertical(1, i + DIR(1, k));
           img_2filter(m) = img_2f_array(k);           
       end               
      [norm, lamda, t] = StegerHess(img_2f_array);        
      if lamda(2) < 0 && abs(lamda(2))>(4 + abs(lamda(1)))         
          W = Linewidth([i, tempy], norm, Range_vertical); %    
          if W <= 40
             LineW(i) = W;
          end
      end
   end
end  
se = strel('line', 5, 0);
LineWtemp = imopen(LineW, se);
LineW = round(imfilter(LineWtemp, [1 3 6 9 10 9 6 3 1]./48, 'replicate'));
LineW(LineW <= 5) = 0; LineW(1) = 0;


% 构建高斯二阶偏导数模板
W = 5;
sigma = 2.5;
xxGauKernel = zeros(2*W + 1, W + 1);
yyGauKernel = zeros(2*W + 1, W + 1); 
xyGauKernel = zeros(2*W + 1, W + 1);  
xGauKernel = zeros(2*W + 1, W + 1);
yGauKernel = zeros(2*W + 1, W + 1);
for i = -W:W	
	for j = -W:W		
		yyGauKernel(i + W + 1, j + W + 1) = (1 - (i*i) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * pi*sigma^4));
		xxGauKernel(i + W + 1, j + W + 1) = (1 - (j*j) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * pi*sigma^4));
		xyGauKernel(i + W + 1, j + W + 1) = ((i*j))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * pi*sigma^6));
        yGauKernel(i + W + 1, j + W + 1) = ( - i / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * pi*sigma^2));
        xGauKernel(i + W + 1, j + W + 1) = ( - j / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * pi*sigma^2));
    end
end

image = double(image);
xxG = imfilter(image, xxGauKernel, 'replicate');% ;
xyG = imfilter(image, xyGauKernel, 'replicate');% ;
yyG = imfilter(image, yyGauKernel, 'replicate');% ;
xG = imfilter(image, xGauKernel, 'replicate');% ;
yG = imfilter(image, yGauKernel, 'replicate');% 

line = zeros(2, Sz(2));
% 精细中心提取
for i = 2:Sz(2)-1   
    if i == 37
       i = 37; 
    end
    y0 = (Range_vertical(2, i-1) + Range_vertical(1, i-1))/2;
    y1 = (Range_vertical(2, i) + Range_vertical(1, i))/2;
    y2 = (Range_vertical(2, i+1) + Range_vertical(1, i+1))/2;
    if (y0==0 || y1==0 || y2==0)|| abs(y1 - y0) > 5 || abs(y1 - y2) > 5 || LineW(i) == 0
        continue;
    else  
        tempx = round(CenterCord(1, i)); 
        tempy = round(CenterCord(2, i));
    if tempx == 0 && tempy == 0 
        continue;
    end
        
    mx_1order = xG(tempy, tempx);
    my_1order = yG(tempy, tempx);
    % 二阶偏导
    mx2_2order = xxG(tempy, tempx);
    mxy_2order = xyG(tempy, tempx);
    my2_2order = yyG(tempy, tempx); 
    hess = [mx2_2order mxy_2order;
            mxy_2order my2_2order;];
    [V, D] = eig(hess);   
    if abs(D(1, 1)) >= abs(D(2, 2))
       norm = V(:, 1) ;
       lamda1 = D(2, 2);        
       lamda2 = D(1, 1);     
    else
       norm = V(:, 2) ; 
       lamda1 = D(1, 1);     
       lamda2 = D(2, 2);     
    end
    lamda = [lamda1; lamda2];  
    t = -(norm(1)*mx_1order + norm(2)*my_1order)/ ...
         (norm(1)^2*mx2_2order + 2*norm(1)*norm(2)*mxy_2order + norm(2)^2*my2_2order); 
               
    if lamda(2) < 0 && abs(lamda(2))>(1 + abs(lamda(1)))        
          line(:, i) = [tempx - t*norm(1); tempy - t*norm(2)];  
    end              
    end
end 

% 检查
toc;
figure;
imshow(image, []);
hold on
plot(line(1, :), line(2, :));

主要参考有两个，一个是综述，一个是法线法的算法实现代码；都具有一定的参考意义，这里也贴出来仅供大家参考。
1：线结构光中心提取综述。
2：光条中心线提取-Steger算法（基于Hessian矩阵）。文章来源地址https://www.toymoban.com/news/detail-401144.html

到了这里，关于【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！