function [oriented,nonoriented,oriented_size_inv,nonoriented_size_inv,...
       AB,B,zero]=shape(region,K)

Nrows=size(region,1);

% get outline
[rowS,colS,dl,dp]=trace_shape(region,0); %% S=0 ==> don't smooth

% get original coefficients
[a,b,c,d,zero]=fourier_ellipse(colS,Nrows-rowS+1,dl,K);

% convert original coefficients to elliptical coefficients
% (remove starting point variation)
desc=get_ellipse(a,b,c,d);

% convert back to Fourier coefficients (after removing starting point)
[ap,bp,cp,dp]=invert_ellipse(desc);
oriented=[ap ; bp ; cp ; dp];

% convert to Fourier coefficients (size-invariant)
desc2=desc;
desc2(1)=1;
[ap,bp,cp,dp]=invert_ellipse(desc2);
oriented_size_inv=[ap ; bp ; cp ; dp];

% convert to Fourier coefficients (orientation-invariant)
desc(3)=0;
[ap,bp,cp,dp]=invert_ellipse(desc);
nonoriented=[ap ; bp ; cp ; dp];

% convert to Fourier coefficients (orientation-invariant, size-invariant)
desc2=desc;
desc2(1)=1;
desc2(3)=0; % (should already be 0)
[ap,bp,cp,dp]=invert_ellipse(desc2);
nonoriented_size_inv=[ap ; bp ; cp ; dp];

AB=abs(desc(1)*(desc(2)*desc(1)));  % since desc(2) is normalized by desc(1)
B=abs(desc(2));

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [a,b,c,d,zero]=fourier_ellipse(x,y,dl,N)
% follows Kuhl & Giardina 1982

M=length(x);

lcum=dl;
for j=2:M
   lcum(j)=lcum(j-1)+dl(j);
end
L=lcum(M); %% = sum(dl)

dx=[x(2:end) x(1)]-x;
dy=[y(2:end) y(1)]-y;
arg1=2*pi*lcum/L;
arg2=[0 2*pi*lcum(1:end-1)/L];

for n=1:N
   a(n) = L/(2*n^2*pi^2) * sum(dx./dl.*(cos(arg1*n)-cos(arg2*n)));
   b(n) = L/(2*n^2*pi^2) * sum(dx./dl.*(sin(arg1*n)-sin(arg2*n)));
   c(n) = L/(2*n^2*pi^2) * sum(dy./dl.*(cos(arg1*n)-cos(arg2*n)));
   d(n) = L/(2*n^2*pi^2) * sum(dy./dl.*(sin(arg1*n)-sin(arg2*n)));
end

for p=1:M
   c1(p)=sum(dx(1:p-1)) - dx(p)/dl(p)*sum(dl(1:p-1));
   c2(p)=sum(dy(1:p-1)) - dy(p)/dl(p)*sum(dl(1:p-1));
end

zero(1)=1/L * sum(...
    dx./(2*dl).*(lcum.^2-[0 lcum(1:end-1)].^2) + ...
    c1.*(lcum-[0 lcum(1:end-1)]));
zero(2)=1/L * sum(...
    dy./(2*dl).*(lcum.^2-[0 lcum(1:end-1)].^2) + ...
    c2.*(lcum-[0 lcum(1:end-1)]));
zero=zero+[x(1) y(1)];

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function desc=get_ellipse(an,bn,cn,dn);
% follows Staib & Duncan 1992

N=length(an);
lasttheta=0;
lastphi=0;

for n=1:N
   a=an(n);b=bn(n);c=cn(n);d=dn(n);
   alpha=a^2+b^2+c^2+d^2;
   beta=a*d-b*c;
   A=sqrt((alpha+sqrt(alpha^2-4*beta^2))/2);
   B=sqrt(2*beta^2/(alpha+sqrt(alpha^2-4*beta^2)))*sign(beta);
   theta=mod(atan((A*c+B*b)/(A*a-B*d)),pi);
   phi=mod(atan2((B*a-A*d),(A*c+B*b)),2*pi);
%%   theta=atan2((A*c+B*b),(A*a-B*d));
%%   phi=atan((B*a-A*d)/(A*c+B*b));
   thetaprime=theta-lasttheta;
   lasttheta=theta;
   if n==1
      phi1=phi;
      phiprime=0;
   else
      phi=phi-n*phi1;
      phiprime=phi-lastphi;
      lastphi=phi;
   end
   thetaprime=mod(thetaprime,2*pi);
   phiprime=mod(phiprime,2*pi);
   if n==1
      A1=A;
      desc=[A B/A1 thetaprime phiprime]; % note phiprime=0
   else
      desc=[desc A/A1 B/A1 thetaprime phiprime];
   end
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [a,b,c,d]=invert_ellipse(desc)
% follows Staib & Duncan 1992

N=length(desc)/4;

cumtheta=0;
cumphi=0;
A1=desc(1);

for n=1:N
   curr=desc((n-1)*4+1:n*4);
   if n==1
      A=curr(1);
   else
      A=curr(1)*A1;
   end
   B=curr(2)*A1;
   th=curr(3)+cumtheta;
   cumtheta=th;
   ph=curr(4)+cumphi;
   cumphi=ph;
   
   a(n)=  A*cos(th)*cos(ph) - B*sin(th)*sin(ph);
   b(n)= -A*cos(th)*sin(ph) - B*sin(th)*cos(ph);
   c(n)=  A*sin(th)*cos(ph) + B*cos(th)*sin(ph);
   d(n)= -A*sin(th)*sin(ph) + B*cos(th)*cos(ph);
end

return




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1997-1999. The Regents of the University
% of California (Regents). All Rights Reserved.
% 
% Permission to use, copy, modify, and distribute this software and its
% documentation for educational, research, and not-for-profit purposes,
% without fee and without a signed licensing agreement, is hereby
% granted, provided that the above copyright notice, this paragraph and
% the following two paragraphs appear in all copies, modifications, and
% distributions. Contact The Office of Technology Licensing, UC Berkeley,
% 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, (510)
% 643-7201, for commercial licensing opportunities. Created by
% Chad Carson, Serge Belongie, and Jitendra Malik, Department of Electrical
% Engineering and Computer Sciences, University of California, Berkeley.
% For technical information, contact carson@eecs.berkeley.edu.
% 
% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF
% ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION
% TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

