EIGENDECOMP2X2SYM - Test function for GSTDECOMP.

Contents

Description

Check that GSTDECOMP returns proper eigenvalues/vectors. Compare that code with the function EIGS implemented in Matlab.

Syntax

   [l1,l2,e1,e2] = EIGENDECOMP2X2SYM(A, method);

Inputs

*|A|* : symmetric |(2,2)| matrix.

*|method|* : one method among |'koe', 'kro', 'pey'| or |'eig'|; the method
  |'eig'| runs the function EIGS implemented in Matlab.

Outputs

*|[l1,l2,e1,e2]|* : output eigenvalues/vectors.

See also

Ressembles: GSTDECOMP. Requires: EIGS.

Function implementation

function [l1,l2,e1,e2] = eigendecomp2x2sym(A, method)

gx2 = A(1,1);
gy2 = A(2,2);
gxy = A(1,2);

if A(2,1)~=gxy,
    error('eigendecomp2x2sym:inputerror','symmetric matrix expected');
elseif ~any(strcmpi(method,{'koe','kro','pey','eig'}));
    error('eigendecomp2x2sym:inputerror',...
        'choose one method among ''koe'', ''kro'', ''pey'' or ''eig''')
end

% run Matlab internal function
if strcmpi(method,'eig')
    [V,D] = eigs(A);
    l1 = D(1); l2 = D(4);
    e1 = V(:,1)'; e2 = V(:,2)';
    return
end

if strcmpi(method,'pey')
    % another expression for eigenvalues
    t = 0.5 * (gx2+gy2);     % trace/2
    a = gx2 - t;  % a = 0.5 * (gx2-gy2);
    ab2 = sqrt(a.^2 + gxy.^2); % we have root = 2 * ab2; (see above)
    l1 = t + ab2;
    l2 = t - ab2;

elseif any(strcmpi(method,{'koe','kro'}))
    % VIGRA expression: see file tensorutilities.hxx, function
    % tensorEigenRepresentation
    % Kroon expression: see file EigenVectors2D.m
    d1 = gx2 + gy2; %src.getComponent(s,0) + src.getComponent(s,2);
    d2 = gx2 - gy2; %src.getComponent(s,0) - src.getComponent(s,2);
    d3 = 2 * gxy;   %2.0 * src.getComponent(s,1);
    d4 = sqrt(d2.^2 + d3.^2); % hypot(d2, d3)
    % root = d4;
    l1 = 0.5 * (d1+d4); %dest.setComponent(0.5 * (d1 + d4), d, 0);
    l2 = 0.5 * (d1-d4); %dest.setComponent(0.5 * (d1 - d4), d, 1);
end

if strcmpi(method,'koe')
    % double angle representation: theta in [-pi/2,pi/2]
    %    %theta = gstfeature(gx2, gy2, gxy, 'orient');
    theta = 0.5 * atan2(d3, d2);
    %dest.setComponent(0.5 * VIGRA_CSTD::atan2(d3, d2), d, 2);
    % note: atan2(0,0)=0
    % principal eigen vector
    e1(1) = cos(theta);
    e1(2) = sin(theta);

elseif strcmpi(method,'pey')
    % atan2 expression for orientation: atan2(l1-gx2,gxy)
    theta = atan2( ab2-a, gxy );
    e1(1) = cos(theta);
    e1(2) = sin(theta);

elseif strcmpi(method,'kro')
    % note: no theta computed here
    vx = d3;
    vy = d4 - d2;
    mag = sqrt(vx.^2 + vy.^2);
    i = (mag ~= 0);
    vx(i) = vx(i)./mag(i);
    vy(i) = vy(i)./mag(i);

    e1(1) = vx;
    e1(2) = vy;
end

%  eigenvectors are orthogonal
e2(1) = -e1(2);
e2(2) = e1(1);

end % end of eigendecomp2x2sym