function elliptrack_sh(Name,level,a_lb, b_hb, area_lb, area_hb)
%This code identifies ellipsoidal features in a dark background
clc

Files=dir([Name '*.tif']);

Result = [];

for k=1:1:length(Files) 

img=imread(Files(k).name); 
clims = [min(min(img)) 3*min(min(img))];
imagesc(img,clims), colormap('gray'), daspect([1 1 1]), axis tight; drawnow; hold on

bw=imbinarize(img,level);
imshow(bw); colormap('gray'); daspect([1 1 1]);
st = regionprops(bw, 'BoundingBox' ); stc = struct2cell(st); bw=transpose(bw);

m00=0; m01=0; m10=0; m11=0; m02=0; m20=0; 
for i=1:1:size(stc,2)
    for X=ceil(stc{i}(1)):1:ceil(stc{i}(1))+stc{i}(3)-1
        for Y=ceil(stc{i}(2)):1:ceil(stc{i}(2))+stc{i}(4)-1
        
        m00=m00+bw(X,Y);    m11=m11+X*Y*bw(X,Y);
        m10=m10+X*bw(X,Y);  m01=m01+Y*bw(X,Y);   
        m20=m20+X^2*bw(X,Y); m02=m02+Y^2*bw(X,Y);
        
        end
    end
    
    u00=m00; u11=m11-m10*m01/m00; u20=m20-m10^2/m00; u02=m02-m01^2/m00; %Calculating the central moments.
  
    if u00 > area_lb
      
    
        if u11 > 0 && u20 > u02
             delta=sqrt(4*u11^2+(u20-u02)^2); phi=0.5*atan(2*u11/(u20-u02));
             a=sqrt(2*(u20+u02+delta)/u11); b=sqrt(2*(u20+u02-delta)/u11); % principal axes, a is the major, b is the minor.
       
        elseif u11 > 0 && u20 < u02
             delta=sqrt(4*u11^2+(u20-u02)^2); phi=pi/2+0.5*atan(2*u11/(u20-u02));
             a=sqrt(2*(u20+u02+delta)/u11); b=sqrt(2*(u20+u02-delta)/u11); % principal axes, a is the major, b is the minor.
         
        elseif u11 < 0 && u20 > u02
                 delta=sqrt(4*u11^2+(u20-u02)^2); phi=0.5*atan(2*u11/(u20-u02));
                 a=sqrt(2*(u20+u02+delta)/u11); b=sqrt(2*(u20+u02-delta)/u11); % principal axes, a is the major, b is the minor.
                
        else
                    delta=sqrt(4*u11^2+(u20-u02)^2); phi=pi/2+0.5*atan(2*u11/(u20-u02));
                    a=abs(sqrt(2*(u20+u02+delta)/u11)); b=abs(sqrt(2*(u20+u02-delta)/u11));
        end
     
         norm=sqrt(u00/(pi*a*b)); a=abs(norm*a); b=abs(norm*b);
         temp = [m10/m00 m01/m00 atan(tan(-phi)) k u00  a b a/b];
         
         if a > a_lb && b < b_hb && u00 < area_hb
             rectangle('position',[ceil(stc{i}(1)) ceil(stc{i}(2)) stc{i}(3)-1 stc{i}(4)-1],'EdgeColor','b')
             fimplicit(@(x,y) ((x-m10/m00)*cos(phi) + (y-m01/m00)*sin(phi)).^2/a^2+((y-m01/m00)*cos(phi) - (x-m10/m00)*sin(phi)).^2/b^2-1, [1 size(bw,1) 1 size(bw,2)],':','Linewidth',2,'Color','y'); drawnow
         else
             rectangle('position',[ceil(stc{i}(1)) ceil(stc{i}(2)) stc{i}(3)-1 stc{i}(4)-1],'EdgeColor','r')
             temp = [];
         end
         % x, y, angle (rad), time (in frame #), a (major), b (minor), a/b.
         Result = [Result ; temp]; %#ok<AGROW>
    else
        
    end
    
    
    m00=0; m01=0; m10=0; m11=0; m02=0; m20=0;
end
 if isempty(Result) == 1
     n=0;
     
 elseif isempty(Result((Result(:,4) == k))) == 1
        n=0;
 else
        n=length(Result((Result(:,4) == k)));
 end

    disp([num2str(n) ' features found in ' num2str(k) ' / ' num2str(length(Files)) '  frame.' ])
  
    hold off
end

csvwrite(['elliptrack_' Name '_a_' num2str(a_lb) '_b_' num2str(b_hb) '_area_' num2str(area_lb) '_' num2str(area_hb) '_full.csv'],Result);

Result(:,5:end)=[];
Result(:,3)=[];
csvwrite(['elliptrack_' Name '_a_' num2str(a_lb) '_b_' num2str(b_hb) '_area_' num2str(area_lb) '_' num2str(area_hb) '_tr.csv'],Result);



end