%PHYSICAL EXTENT OF IMAGES about mean epoch in terms of the 0, 1 and 2  MOMENTS
% 
%dt = extent2(G); % where "G" is a matrix
%
%extent is given by: sqrt(2pi*MSE) and is independent of relative amplitude
%
%extent is padding invariant:  extent([1 1]) = extent([0 0 0 0 1 1])
%
%see    daugman 1985, opt soc. am.  BUT USES same NORMALIZ as Gabor, 1945)
%
% may use optional forced epoch location indices,
% for example, extent2(glt,[50 50]) takes mse spread about center of 99by99 glt
%
% [dx, dy] = extent2(A)
% [dxdy] = extent2(A)  % gives prod dx*dy  ANY VECTOR HAS 0 EXTENT2
% [dx, dy] = extent(A,[50 50])  % gives mse spread about indices x=50,y=50
%
% [dx, dy, dxy2] = ...  gives measure of skewness: 0 means not skew (manhatten)
%                      >0 for ellipse along X=Y (up and to right)
%                      <0 for ellipse going from top left to lower right (X=-Y)
%          not sqrt since its often <0    (dxy would be cmpx)

function [dx, dy, dxy2] = extent2(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;   % mean in x direction
my = sum(sum(conj(G).*Y.*G))/m0;   % actual mean in y dir

if nargin > 1;
  mx = forced_epoch_index(1);
  my = forced_epoch_index(2);
end%if

dx = sqrt(2*pi*sum(sum(conj(G).*(X-mx).^2..*G))/m0);
dy = sqrt(2*pi*sum(sum(conj(G).*(Y-my).^2..*G))/m0);
dxy2 = (2*pi*sum(sum(conj(G).*(X-mx).*(Y-my).*G))/m0); % not sqrt since
%        dxy2 is often negative  dxy would be cmpx


%%% 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;

%%% dxdy = (2*pi*(Ett-Et^2)) ^ (1/2);
dxdy = dx*dy;  % NOT SURE ABOUT THE NORMALIZATION...

if nargout < 2  % nargout 0 or 1   (0 if command typed to displ as "ans")
  dx = dxdy;
end%if

