PRO MESHGRID,N,M,X,Y
;+
; PRO MESHGRID,N,M,X,Y calculates two N Column x M Row arrays, X
; and Y, that can be used to calculate and plot a 2D function.
;
; USAGE:
; MESHGRID,31,31,X,Y
; X=(X-15)/2.
; Y=(Y-11)/3.
; Z=EXP(-(X^2+Y^2))
; SURFACE,Z
;
; HISTORY:
; Originated by Harvey Rhody September 17, 1997
;-

X=INDGEN(N)#(1+INTARR(M))
Y=INDGEN(M)##(1+INTARR(N))

END

FUNCTION Moments2D,A,p,q
;+
;m=Moments2D(A,p,q) computes the Mpq moment of A. p and q may be
;vectors with a list of moment exponents. In that case, the result
;is a vector of values for the items in the p and q lists. The meshgrid
;procedure is used to construct X and Y arrays. The pq moment is
;computed by summing (X^p)*(Y^q)*A with the TOTAL function.
;
;USES meshgrid to compute the X and Y fields.
;
;HISTORY
;H. Rhody, Sept. 1999
;-

;Construct the X and Y grids of the same
;size as the image.

aSiz=SIZE(A)
meshGrid,aSiz[1],aSiz[2],X,Y
X=FLOAT(X) & Y=FLOAT(Y)

n=N_ELEMENTS(p)

;If p and q have one element then one calculation is needed.

IF n EQ 1 THEN RETURN,TOTAL(X^p*Y^q*A)

;If p and q are vectors, then compute a result for
;each value pair in the list.

IF n GT 1 THEN BEGIN
	m=FLTARR(n)
	FOR k=0,n-1 DO m[k]=TOTAL(X^p[k]*Y^q[k]*A)
ENDIF
RETURN,m
END

;==============================================

FUNCTION CentralMoments2D,A,p,q
;+
;m=CentralMoments2D(A,p,q) computes the Mpq central moment of A.
;p and q may be vectors with a list of moment exponents. In
;that case, the result is a vector of values for the items in
;the p and q lists.
;
;USES meshgrid to compute the X and Y fields.
;
;
;HISTORY
;H. Rhody, Sept. 1999
;-

;Construct the X and Y arrays the same size as A.
aSiz=SIZE(A)
meshGrid,aSiz[1],aSiz[2],X,Y
X=FLOAT(X) & Y=FLOAT(Y)

;Calculate the normalization factor and the centroid values.
m00=TOTAL(A)
xbar=TOTAL(X*A)/m00
ybar=TOTAL(Y*A)/m00

n=N_ELEMENTS(p)

;If p and q are scalar values, then return the related result
IF n EQ 1 THEN RETURN,TOTAL((X-xbar)^p*(Y-ybar)^q*A)

;If p and q are vectors, then return a vector with values for
;each (p,q) in the list.
IF n GT 1 THEN BEGIN
	m=FLTARR(n)
	FOR k=0,n-1 DO m[k]=TOTAL((X-xbar)^p[k]*(Y-ybar)^q[k]*A)
ENDIF
RETURN,m
END

;==============================================

FUNCTION NormalMoments2D,A,p,q
;+
;m=NormalMoments2D(A,p,q) computes the Mpq normalized central moment of A.
;p and q may be vectors with a list of moment exponents. In
;that case, the result is a vector of values for the items in
;the p and q lists.
;
;USES meshgrid to compute the X and Y fields.
;
;
;HISTORY
;H. Rhody, Sept. 1999
;-

;Construct the X and Y arrays the same size as A.

aSiz=SIZE(A)
meshGrid,aSiz[1],aSiz[2],X,Y
X=FLOAT(X) & Y=FLOAT(Y)

;Calculate the normalization factor and the centroid values.
m00=TOTAL(A)
xbar=TOTAL(X*A)/m00
ybar=TOTAL(Y*A)/m00

;Calculate the exponent for the normalization term.
n=N_ELEMENTS(p)
gamma=FLOAT(p+q)/2+1

;If p and q are scalar values, then return the related result
IF n EQ 1 THEN RETURN,TOTAL((X-xbar)^p*(Y-ybar)^q*A)/m00^gamma

;If p and q are vectors, then return a vector with values for
;each (p,q) in the list.
IF n GT 1 THEN BEGIN
	m=FLTARR(n)
	FOR k=0,n-1 DO m[k]=TOTAL((X-xbar)^p[k]*(Y-ybar)^q[k]*A)
ENDIF
RETURN,m/m00^gamma
END

;====================================================

FUNCTION InvariantMoments2D,A
;+
;v=InvariantMoments2D(A) computes the seven invariant moment
;features described in Gonzalez and Woods, page 516.
;
;v[k]=phi[k+1], k=0,1,...,6 where phi[n] is element n in
;Gonzalez and Woods notation.
;
;USES NormalMoments2D.
;
;HISTORY
;H. Rhody, Sept. 1999
;-

;Set up the p and q vectors to create the vector of normal moment 
;results to support the calculation of the invariants.
p=[0,0,1,1,2,2,3]
q=[2,3,1,2,0,1,0]
m=NormalMoments2D(A,p,q)

;Compute the invariants. Note that m[0]=m02, m[1]=m03, etc. from the 
;order of coefficients in the p and q vectors. The invariants phi_1,
;phi_2,...,phi_7 of Gonzalez and Woods are the elements v[],v[1], etc.
;of the v vector.
v=FLTARR(7)
v[0]=m[0]+m[4]
v[1]=(m[4]-m[0])^2+4*m[2]^2
v[2]=(m[6]-3*m[3])^2+(3*m[5]-m[1])^2
v[3]=(m[6]+m[3])^2+(m[5]+m[1])^2
v[4]=(m[6]-3*m[3])*(m[6]+m[3])*((m[6]+m[3])^2- $
	3*(m[5]+m[1])^2) +(3*m[5]-m[1])*(m[5]+m[1])* $
	(3*(m[6]+m[3])^2-(m[5]+m[1])^2)
v[5]=(m[4]-m[0])*((m[6]+m[3])^2-(m[5]+m[1])^2)+ $
	4*m[2]*(m[6]+m[3])*(m[5]+m[1])
v[6]=(3*m[5]-m[6])*(m[6]+m[3])*((m[6]+m[3])^2 $
	-3*(m[5]+m[1])^2)+(3*m[3]-m[6])*(m[5]+m[1])* $
	(3*(m[6]+m[3])^2-(m[3]+m[1])^2)

RETURN,v
END