eigenvalue decomposition of structure tensor in ma

Question

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

Any help is appreciated.

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:

This is the new plot with axis equal:

Answer 1

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));

(This requires DIPimage to run)

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 

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

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.

Some suggestions:

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);
   

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

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

3. 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.

4. 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.

5. 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,:};