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

三维成像魔法

已于 2023-03-17 11:18:49 修改

阅读量3.9k

点赞数 4

分类专栏：机器视觉文章标签： 3d成像视觉检测

于 2022-07-06 19:34:44 首次发布

本文链接：https://blog.csdn.net/weixin_40204808/article/details/125210303

版权

机器视觉专栏收录该内容

6 篇文章 4 订阅

订阅专栏

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}\,$ 像素，表示的一阶梯度；

如下图所示，线条处两个特征值表示了图像在两个特征向量所指的方向上图像变化的各向异性。相似的，圆具有各项同性；线性结构越强的越是各项异性。
在这里插入图片描述

海森矩阵的应用（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矩阵之前，根据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矩阵）。

三维成像魔法

关注

4
点赞
踩
30

收藏

觉得还不错? 一键收藏
4
评论
【图像处理】海森矩阵（Hessian Matrix）及用例（基于Steger的中心提取_含代码）

线结构光中心提取算法是线扫描三维重建精度高低以及光条轮廓准确性的重要因素。本方法可以改善提取的光条中心有断线；当光条较宽时，线条值局部会出现上下波动；当有弧光反射时光条中心线局部上下跳变等情况。...
复制链接

扫一扫