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

问题描述

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

I have a synthetic image. I want to do eigenvalue decomposition of local structure tensor (LST) of it for some edge detection purposes. I used the eigenvaluesl1 , l2 and eigenvectors e1 ,e2 of LST to generate an adaptive ellipse for each pixel of image. Unfortunately I get unequal eigenvalues l1 , l2 and so unequal semi-axes length of ellipse for homogeneous regions of my figure:

不过，对于一个简单的测试图像，我得到了很好的响应:

However I get good response for a simple test image:

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

I don't know what is wrong in my code:

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

感谢您的帮助.

这是我的椭圆绘图代码:

This is my ellipse drawing code:

% 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

这是我对高斯导数核的新结果:

This my new result with the Gaussian derivative kernel:

这是使用axis equal的新情节:

推荐答案

我创建了一个类似于您的测试图像(可能不太复杂)，如下所示:

I created a test image similar to yours (probably less complicated) as follows:

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之外，其他都没有改变)，并得到以下结果:

Then ran your LST_eig on it (sigma1=1 and rw=3) and your code to draw ellipses (no change to either, except adding axis equal), and got this result:

我怀疑图像的某些蓝色区域有些不均匀，这会导致出现很小的渐变.使用椭圆时定义椭圆的问题是，对于足够定向的图案，即使看不到该图案，您也会得到一条直线.您可以通过如下定义椭圆轴长度来解决此问题:

I suspect some non-uniformity in some of the blue areas of your image, which cause very small gradients to appear. The problem with the definition of the ellipses as you use them is that, for sufficiently oriented patterns, you'll get a line even if that pattern is imperceptible. You can get around this by defining your ellipse axes lengths as follows:

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仅依赖于能量:

The smallest eigenvalue l1 is high in regions with strong gradients but no clear direction. The above does not take this into account. One option could be to make a depend on both energy and anisotropy measures, and b depend only on energy:

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取决于图像强度，请根据需要进行调整.

This way, the whole ellipse shrinks if there are strong gradients but no clear direction, whereas it stays large and circular when there are only weak or no gradients. The threshold T depends on image intensities, adjust as needed.

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

You should probably also consider taking the square root of the eigenvalues, as they correspond to the variance.

一些建议:

1. 你可以写

1. You can write

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通过定义进行了归一化，因此无需再次对其进行归一化.

without a loop. Note that e1 is normalized by definition, there is no need to normalize it again.

使用高斯梯度代替高斯平滑，然后使用Sobel或Schaar过滤器.有关某些MATLAB实现的详细信息，请参见此处.

Use Gaussian gradients instead of Gaussian smoothing followed by Sobel or Schaar filters. See here for some MATLAB implementation details.

当需要所有特征值时，请使用eig，而不是eigs.尤其对于这么小的矩阵，使用eigs没有优势. eig似乎会产生更一致的结果.不需要取特征值的绝对值(D = abs(D))，因为根据定义它们是非负的.

Use eig, not eigs, when you need all eigenvalues. Especially for such a small matrix, there is no advantage to using eigs. eig seems to produce more consistent results. There is no need to take the absolute value of the eigenvalues (D = abs(D)), as they are non-negative by definition.

您的默认值rw = 0.001太小，该大小的sigma对图像没有影响.平滑的目的是平均局部邻域中的梯度.我使用rw=3效果很好.

Your default value of rw = 0.001 is way too small, a sigma of that size has no effect on the image. The goal of this smoothing is to average gradients in a local neighborhood. I used rw=3 with good results.

使用DIPimage.有一个structuretensor函数，高斯梯度和更多有用的东西. 3.0版本(仍在开发中)是一个主要的重写，在处理矢量和矩阵方面有显着改进值的图像.我可以将您的所有LST_eig编写如下:

Use DIPimage. There is a structuretensor function, Gaussian gradients, and a lot more useful stuff. The 3.0 version (still in development) is a major rewrite that improves significantly on dealing with vector- and matrix-valued images. I can write all of your LST_eig as follows:

I = dip_image(I);
g = gradient(I, sigma1);
H = gaussf(g*g.', rw);
[e,l] = eig(H);
% Equivalences with your outputs:
l1 = l{2};
l2 = l{1};
e1 = e{2,:};
e2 = e{1,:};

这篇关于Matlab中结构张量的特征值分解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！