function theta=steger(EinRDR)
%光条中心
 im=double(EinRDR);
im = im/max(im(:)); 
%imshow(im)
bw=im2bw(im);
[x0,y0]=find(bw==1);
sigma=1;
% 求hessian矩阵
[Dx,Dy,Dxx,Dxy,Dyy] = Hessian2D(im,sigma);
 
[eigenvalue1, eigenvalue2, eigenvectorx, eigenvectory]=eig2image(Dxx, Dxy, Dyy);  
  
% 
t = -(Dx.*eigenvectorx + Dy .* eigenvectory) ./...  
    (Dxx .* eigenvectorx.^2 + 2*Dxy.*eigenvectorx.*eigenvectory + Dyy.*eigenvectory.^2 );  
  
px = t.*eigenvectorx;  
py = t.*eigenvectory;  
  
[candidateX1, candidateY1] = find(px >= -0.5 & px <= 0.5 & py >= -0.5 & py <= 0.5 & bw==1);  
  
linePixel = [candidateX1, candidateY1];


for i=1:size(candidateX1,1)
    m1=candidateX1(i,1);
    n1=candidateY1(i,1);
    px1(i,1)=px(m1,n1);
    py1(i,1)=py(m1,n1);
    x1(i,1)=m1+px(m1,n1);
    x2(i,1)=py(m1,n1)+n1;
end

p=fitPolynomialRANSAC([x2,x1],1,0.5);

if isempty(p)
    
   xfc=cov(x2,x1);
   [v,d]=eig(xfc);
   [u,q]=find(d==max(d(:)));
   theta=atan(v(2,u)/v(1,u))*180/pi;
   
else
    theta=atan(p(1))*180/pi;
end
if theta<0
    theta=180+theta;
end
%w1=theta1-InRandomdegs8
% p2=polyfit(y0,x0,1);
% theta2=atan(p2(1))*180/pi;
% w2=theta2-InRandomdegs8
  
% figure(1), imshow(im);  
% hold on  
% plot(linePixel(:,2), linePixel(:,1), 'g.');  
% hold off;  

function [Dx,Dy,Dxx,Dxy,Dyy] = Hessian2D(I,Sigma)
%  This function Hessian2 Filters the image with 2nd derivatives of a 
%  Gaussian with parameter Sigma.
% 
% [Dxx,Dxy,Dyy] = Hessian2(I,Sigma);
% 
% inputs,
%   I : The image, class preferable double or single
%   Sigma : The sigma of the gaussian kernel used
%
% outputs,
%   Dxx, Dxy, Dyy: The 2nd derivatives
%
% example,
%   I = im2double(imread('moon.tif'));
%   [Dxx,Dxy,Dyy] = Hessian2(I,2);
%   figure, imshow(Dxx,[]);
%
% Function is written by D.Kroon University of Twente (June 2009)

if nargin < 2, Sigma = 1; end

% Make kernel coordinates
[X,Y]   = ndgrid(-round(3*Sigma):round(3*Sigma));

% Build the gaussian 2nd derivatives filters
DGaussx  = 1/(2*pi*Sigma^4)*(-X).* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussy  = 1/(2*pi*Sigma^4)*(-Y).* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussxx = 1/(2*pi*Sigma^4) * (X.^2/Sigma^2 - 1) .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussxy = 1/(2*pi*Sigma^6) * (X .* Y)           .* exp(-(X.^2 + Y.^2)/(2*Sigma^2));
DGaussyy = DGaussxx';

Dx  = imfilter(I,DGaussx,'conv');
Dy  = imfilter(I,DGaussy,'conv');
Dxx = imfilter(I,DGaussxx,'conv');
Dxy = imfilter(I,DGaussxy,'conv');
Dyy = imfilter(I,DGaussyy,'conv');

function [Lambda1,Lambda2,Ix,Iy]=eig2image(Dxx,Dxy,Dyy)
% This function eig2image calculates the eigen values from the
% hessian matrix, sorted by abs value. And gives the direction
% of the ridge (eigenvector smallest eigenvalue) .
% 
% [Lambda1,Lambda2,Ix,Iy]=eig2image(Dxx,Dxy,Dyy)
%

%
% | Dxx  Dxy |
% |          |
% | Dxy  Dyy |


% Compute the eigenvectors of J, v1 and v2
tmp = sqrt((Dxx - Dyy).^2 + 4*Dxy.^2);
v2x = 2*Dxy; v2y = Dyy - Dxx + tmp;

% Normalize
mag = sqrt(v2x.^2 + v2y.^2); i = (mag ~= 0);
v2x(i) = v2x(i)./mag(i);
v2y(i) = v2y(i)./mag(i);

% The eigenvectors are orthogonal
v1x = -v2y; 
v1y = v2x;

% Compute the eigenvalues
mu1 = 0.5*(Dxx + Dyy + tmp);
mu2 = 0.5*(Dxx + Dyy - tmp);

% Sort eigen values by absolute value abs(Lambda1)<abs(Lambda2)
check=abs(mu1)>abs(mu2);

Lambda1=mu1; Lambda1(check)=mu2(check);
Lambda2=mu2; Lambda2(check)=mu1(check);

Ix=v1x; Ix(check)=v2x(check);
Iy=v1y; Iy(check)=v2y(check);