matlab中结构张量的特征值分解

Question

我有一个合成图像。为了达到边缘检测的目的，我想对局部结构张量(LST)进行特征值分解。利用图像的特征值l1、l2和特征向量e1、e2为图像的每个像素生成一个自适应椭圆。不幸的是，对于我的图中的同质区域，我得到了不相等的特征值l1、l2和不相等的椭圆半轴长度：

​

然而，对于一个简单的测试映像，我得到了很好的响应：

​

我不知道我的代码有什么问题：

function [H,e1,e2,l1,l2] = LST_eig(I,sigma1,rw)

%  LST_eig - compute the structure tensor and its eigen
% value decomposition
%
%   H = LST_eig(I,sigma1,rw);
%
%   sigma1 is pre smoothing width (in pixels).
%   rw is filter bandwidth radius for tensor smoothing (in pixels).
%


n = size(I,1);
m = size(I,2);
if nargin<2
    sigma1 = 0.5;
end
if nargin<3
    rw = 0.001;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pre smoothing
J = imgaussfilt(I,sigma1);
% compute gradient using Sobel operator
Sch = [-3 0 3;-10 0 10;-3 0 3];
%h = fspecial('sobel');
gx = imfilter(J,Sch,'replicate');
gy = imfilter(J,Sch','replicate');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute tensors

gx2 = gx.^2;
gy2 = gy.^2;
gxy = gx.*gy;

% smooth
gx2_sm = imgaussfilt(gx2,rw); %rw/sqrt(2*log(2))
gy2_sm = imgaussfilt(gy2,rw);
gxy_sm = imgaussfilt(gxy,rw);
H = zeros(n,m,2,2);
H(:,:,1,1) = gx2_sm; 
H(:,:,2,2) = gy2_sm; 
H(:,:,1,2) = gxy_sm; 
H(:,:,2,1) = gxy_sm; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% eigen decomposition
l1 = zeros(n,m);
l2 = zeros(n,m);
e1 = zeros(n,m,2);
e2 = zeros(n,m,2);
for i = 1:n
    for j = 1:m
        Hmat = zeros(2);
        Hmat(:,:) = H(i,j,:,:);
        [V,D] = eigs(Hmat);
        D = abs(D);
        l1(i,j) = D(1,1); % eigen values
        l2(i,j) = D(2,2); 
        e1(i,j,:) = V(:,1); % eigen vectors
        e2(i,j,:) = V(:,2); 
    end
end

任何帮助都是非常感谢的。

这是我的椭圆绘图代码：

% determining ellipse parameteres from eigen value decomposition of LST

M = input('Enter the maximum allowed semi-major axes length: ');
I = input('Enter the input data: ');

row = size(I,1);
col = size(I,2);
a = zeros(row,col);
b = zeros(row,col);
cos_phi = zeros(row,col);
sin_phi = zeros(row,col);


for m = 1:row
  for n = 1:col

    a(m,n) = (l2(m,n)+eps)/(l1(m,n)+l2(m,n)+2*eps)*M;
    b(m,n) = (l1(m,n)+eps)/(l1(m,n)+l2(m,n)+2*eps)*M;

    cos_phi1 = e1(m,n,1);
    sin_phi1 = e1(m,n,2);
    len = hypot(cos_phi1,sin_phi1);           
    cos_phi(m,n) = cos_phi1/len;
    sin_phi(m,n) = sin_phi1/len;

  end
end

%% plot elliptic structuring elements using parametric equation and superimpose on the image 


figure; imagesc(I); colorbar; hold on

t = linspace(0,2*pi,50);

for i = 10:10:row-10
  for j = 10:10:col-10
    x0 = j;
    y0 = i;

    x = a(i,j)/2*cos(t)*cos_phi(i,j)-b(i,j)/2*sin(t)*sin_phi(i,j)+x0;
    y = a(i,j)/2*cos(t)*sin_phi(i,j)+b(i,j)/2*sin(t)*cos_phi(i,j)+y0;

    plot(x,y,'r','linewidth',1);
    hold on
  end
end 

这是我用高斯导数核得出的新结果：

​

这是axis equal的新情节

​

​

Answer 1

我创建了一个类似于您的测试映像(可能不那么复杂)，如下所示：

pos = yy([400,500]) + 100 * sin(xx(400)/400*2*pi);
img = gaussianlineclip(pos+50,7) + gaussianlineclip(pos-50,7);
I = double(stretch(img));

(这需要运行DIPimage )

然后在其上运行LST_eig (sigma1=1和rw=3)和绘制省略号的代码(除添加axis equal外，两者都不更改)，并得到以下结果：

​

​

我怀疑你的图像中的一些蓝色区域有一些不均匀性，这会导致很小的梯度出现。当您使用椭圆时，其定义的问题是，对于充分定向的模式，即使该模式是不可察觉的，您也会得到一条线。您可以通过定义椭圆轴的长度来解决这个问题，如下所示：

a = repmat(M,size(l2)); % longest axis is always the same
b = M ./ (l2+1); % shortest axis is shorter the more important the largest eigenvalue is 

在梯度强但方向不清晰的区域，最小特征值l1较高。上述情况并没有考虑到这一点。一种选择是使a同时依赖于能量和各向异性度量，而b只依赖于能量：

T = 1000; % some threshold
r = M ./ max(l1+l2-T,1); % circle radius, smaller for higher energy
d = (l2-l1) ./ (l1+l2+eps); % anisotropy measure in range [0,1]
a = M*d + r.*(1-d); % use `M` length for high anisotropy, use `r` length for high isotropy (circle)
b = r; % use `r` width always

这样，当存在强梯度而没有明确的方向时，整个椭圆就会收缩，而当只有弱的或没有梯度的时候，它就会保持大的圆形。阈值T取决于图像强度，根据需要进行调整。

您可能还应该考虑取特征值的平方根，因为它们对应于方差。

几点建议：

1. 你可以写 A=(l2+eps)/(l1+l2+2*eps)* M；b= (l1+eps)./(l1+l2+2*eps) * M；cos_phi = e1(:，:,1)；sin_phi = e1(:，:,2)； 没有循环。请注意，e1是按定义规范化的，因此没有必要再次对其进行规范化。
2. 使用高斯梯度代替高斯平滑，然后是Sobel或Schaar滤波器。有关一些MATLAB实现细节，请参见这里。
3. 当您需要所有的特征值时，请使用eig而不是eigs。特别是对于如此小的矩阵，使用eigs没有任何好处。eig似乎产生了更一致的结果。不需要取特征值(D = abs(D))的绝对值，因为它们在定义上是非负的。
4. 您的默认值rw = 0.001太小了，这个大小的西格玛对图像没有影响。这种平滑的目的是在局部邻域中平均梯度。我使用rw=3取得了很好的效果。
5. 使用DIPimage。有一个structuretensor函数，高斯梯度，还有更多有用的东西。3.0版(仍在开发中)是一种主要的重写，在处理向量和矩阵值图像方面有很大的改进.我可以按以下方式编写您的所有LST_eig： I= dip_image(I)；g=梯度(i，sigma1)；H= gaussf(g*g.'，rw)；e，l= eig(H)；与您的输出相等的百分比: l1 = l{2}；l2 = l{1}；e1 = e{2，:}；e2 = e{1，:}；