%PHYSICAL EXTENTs OF IMAGES ...  DISPLAY  0, 1 and 2  MOMENTS of ENERGY
% 
%extents2(G); % where "G" is a matrix  RETURNS NOTHING
%
%see    daugman 1985, opt soc. am. (BUT dx and dy USE same NORMALIZ Gabor, 1945)
%
% may use optional forced epoch location indices,
% for example, extents2(glt,[50 50]) takes mse spread about center of 99by99 glt
%
% extents2(A)  % gives prod dx*dy  ANY VECTOR HAS 0 EXTENT2
% extents2(A,[50 50])  % gives mse spread about indices x=50,y=50
%
%      s   ^
%      i   |
%      x   |
%      a   |
%          |
%      y   |_________________________ 
%                   x axis

function extents2(G,forced_epoch_index);

[M N] = size(G);
[X Y] = meshdom(1:N,1:M);  % physical coordinates (meshdom>>> x like N y like M)

m0 = sum(sum(conj(G).*G)) % sqrenergy to norm and make extent2 indep. of scale
mx = sum(sum(conj(G).*X.*G))/m0   % center of gravity of |G|.^2 in x direction;
my = sum(sum(conj(G).*Y.*G))/m0   % mean in y dir
                                  % [mx,my] marks center of grav. of |G|.^2
if nargin > 1;
  mx = forced_epoch_index(1)
  my = forced_epoch_index(2)
end%if

dxdx = sum(sum(conj(G).*(X-mx).^2..*G))/m0  %%% still needs the sqrt(2pi(...
dydy = sum(sum(conj(G).*(Y-my).^2..*G))/m0  %% still needs the 2pi scaling
dxy  = sum(sum(conj(G).*(X-mx).*(Y-my).*G))/m0  % cross terms...
%%% dxdy = (2*pi*(Ett-Et^2)) ^ (1/2);
dxdy = sqrt(dxdx)*sqrt(dydy)  % NOT SURE ABOUT THE NORMALIZATION...
               % this normaliz gives area of ones(m,n) approx = actual area mn

mxx = sum(sum(conj(G).*X.^2..*G))
myy = sum(sum(conj(G).*Y.^2..*G))
mxy = sum(sum(conj(G).*X.*Y.*G))

%%% Et = m1/m0;
%%% Ett= m2/m0;

nastyness = sum(sum(conj(G).*G))/(sum(sum(abs(G)))).^2

kurtosis = sum(sum(conj(G).*G.*conj(G).*G))/(sum(sum(conj(G).*G)))^2
%%%G2 = conj(G).*G;  % power or energy G
%%%kurtosis = sum(sum(G2.*G2))/(sum(sum(G2)))^2

