function [L1,L2,phi,aniso,pol,con,window_sizes]=wismm3(V);
% [L1,L2,phi,aniso,pol,con,window_sizes]=wismm3(V);
% Calculate windowed image second moment matrices for image V and return
% the following values:
%
% L1 is the larger eigenvalue (lambda_1)
% L2 is the smaller eigenvalue (lambda_2)
% phi is the angle of the 1st eigenvector (phi)
%

% note "natural" smoothing is done at a scale proportional to the "artificial" scale, as Lindeberg suggests

[Gx,Gy]=gradient(V);

GGTxx=Gx.^2;
GGTxy=Gx.*Gy;
GGTyy=Gy.^2;

[r,c]=size(V);

window_sizes=[1 2 5 10 17 26 37 50]; % true values
s=sqrt(window_sizes-1)/2;
max_count=length(window_sizes);

L1=zeros(r,c,max_count);
L2=zeros(r,c,max_count);
phi=zeros(r,c,max_count);
pol=zeros(r,c,max_count);
con=zeros(r,c,max_count);

%fprintf(1,'Integration window size: ');
counter=1;
for n=window_sizes
   %fprintf(1,'%d ',n);
   Txx=gaussN(GGTxx,n);
   Txy=gaussN(GGTxy,n);
   Tyy=gaussN(GGTyy,n);
   
   tr=Txx+Tyy;
   V1=0.5*sqrt(tr.^2-4*(Txx.*Tyy-Txy.^2));
   V1=real(V1);
   
   L1(:,:,counter)=0.5*tr+V1;
   L2(:,:,counter)=0.5*tr-V1;
   phi(:,:,counter)=0.5*atan2(2*Txy,Txx-Tyy);

   % do polarity stuff here
   [P,angle_vector]=polarity(Gx,Gy,n);
   quant_bound=angle_vector(2)/2;
   % quantize angle and pull corresponding polarity out of P
   % (perhaps use set-theoretic functions for masking out P???)
   % (can this be made faster by using an incremental method?)
   for m=1:length(angle_vector);
      a=angle_vector(end-m+1);
      old_pol=pol(:,:,counter);
      Pmask=abs(cos(phi(:,:,counter)-a))>=cos(quant_bound);
      Pmask=Pmask&(old_pol==0); % prevent pileup on quant. boundaries
      pol(:,:,counter)=old_pol+(Pmask.*P(:,:,m));
   end
   
   % contrast calculation
   con(:,:,counter)=tr;
   counter=counter+1;
end
%fprintf(1,'\n')

aniso=1-L2./(L1+eps);
con=con/0.25; % normalize contrast to [0,1]
con=min(con,1);
con=sqrt(con);

return


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


function [P,angle_vector]=polarity(Gx,Gy,N);
% [P,angle_vector]=polarity(Gx,Gy,N);
% 
% prepare pre-convolved polarity measures as several discrete angles from
% 0 to pi
%
% uses integration window of size N

angle_spacing=pi/8;
num_angles=pi/angle_spacing;
angle_vector=(0:num_angles-1)*angle_spacing;
P=zeros(size(Gx,1),size(Gx,2),num_angles);

G=Gx+i*Gy;

for n=1:num_angles
   Eall=real(G.*exp(i*angle_spacing*(n-1)));
   Eplus=Eall.*(Eall>=0);
   Eminus=-Eall.*(Eall<0);
   EplusW=gaussN(Eplus,N)+eps;
   EminusW=gaussN(Eminus,N)+eps;
   P(:,:,n)=abs(EplusW-EminusW)./(EplusW+EminusW);
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

