FUNCTION HWT,image,INVERSE=inverse
;+
;B=HWT(A) returns the Haar wavelet transform of array A.
;A may be a row vector or a 2D array. The number of
;columns must be even (or the result will be truncated
;to an even number). The result B is a
;floating point array with the same dimensions as A
;provided A has an even number of columns.
;
;If A has dimension NxM then B[0:N/2-1,*] contains the
;low-frequency part of the result and B[N/2:2*(N/2)-1,*] contains
;the high-frequency part.
;
;The inverse transform is computed by A=HWT(B,/INVERSE).
;The result is a floating point array. If it is important
;to maintain the array type then that must be done
;as an external step.
;
;EXAMPLE 1:
;A=[1,-3,0,4]
;B=HWT(A)
;PRINT,B
;     -1.41421      2.82843      2.82843     -2.82843
;C=HWT(B,/INVERSE)
;PRINT,C
;      1.00000     -3.00000     0.000000      4.00000
;
;EXAMPLE 2:
;A=DIST(4) & PRINT,A
;
;      0.000000      1.00000      2.00000      1.00000
;      1.00000      1.41421      2.23607      1.41421
;      2.00000      2.23607      2.82843      2.23607
;      1.00000      1.41421      2.23607      1.41421
;
;Compute the HWT along the rows.
;B1=HWT(A) & PRINT,B1
;
;      0.707107      2.12132    -0.707107     0.707107
;      1.70711      2.58114    -0.292893     0.581139
;      2.99535      3.58114    -0.166925     0.418861
;      1.70711      2.58114    -0.292893     0.581139
;
;Transpose B1 and compute the HWT, which does the columns.
;B2=HWT(TRANSPOSE(B1)) & PRINT,B2
;
;      1.70711      3.32514    -0.707107     0.910927
;      3.32514      4.35739    -0.325141     0.707107
;    -0.707107    -0.325141    -0.292893    0.0890728
;     0.910927     0.707107    0.0890728    -0.114748
;
;Compute the inverse HWT of B2
;C1=HWT(B2,/INVERSE) & PRINT,C1
;
;      0.707107      1.70711      2.99535      1.70711
;      2.12132      2.58114      3.58114      2.58114
;    -0.707107    -0.292893    -0.166925    -0.292893
;     0.707107     0.581139     0.418861     0.581139
;
;C1 is seen to be the same as TRANSPOSE(B1). We can
;recover the original array by
;C=HWT(TRANSPOSE(C1),/INVERSE)
;
;      0.000000      1.00000      2.00000      1.00000
;      1.00000      1.41421      2.23607      1.41421
;      2.00000      2.23607      2.82843      2.23607
;      1.00000      1.41421      2.23607      1.41421
;
;To do the HWT on both rows and columns it is necessary that
;the number of rows also be an even number.
;
;;EXAMPLE 3: Odd dimensions.
;A=DIST(5) & PRINT,A
;      0.000000      1.00000      2.00000      2.00000      1.00000
;      1.00000      1.41421      2.23607      2.23607      1.41421
;      2.00000      2.23607      2.82843      2.82843      2.23607
;      2.00000      2.23607      2.82843      2.82843      2.23607
;      1.00000      1.41421      2.23607      2.23607      1.41421
;
;B1=HWT(A) & PRINT,B1
;
;      0.707107      2.82843    -0.707107     0.000000
;      1.70711      3.16228    -0.292893     0.000000
;      2.99535      4.00000    -0.166925     0.000000
;      2.99535      4.00000    -0.166925     0.000000
;      1.70711      3.16228    -0.292893     0.000000
;
;B2=HWT(TRANSPOSE(B1) & PRINT,B2
;
;      1.70711      4.23607    -0.707107     0.000000
;      4.23607      5.65685    -0.236068     0.000000
;    -0.707107    -0.236068    -0.292893     0.000000
;     0.000000     0.000000     0.000000     0.000000
;
;We now have an even number of rows and columns. Now
;see what we get when we transform back.
;
;B3=HWT(B2,/INVERSE) & PRINT,B3
;
;     0.707107      1.70711      2.99535      2.99535
;      2.82843      3.16228      4.00000      4.00000
;    -0.707107    -0.292893    -0.166925    -0.166925
;     0.000000     0.000000     0.000000     0.000000
;
;This is the transpose of B1 with one less column. Now
;invert this.
;
;B4=HWT(TRANSPOSE(B3),/INVERSE) & PRINT,B4
;
;     0.000000      1.00000      2.00000      2.00000
;      1.00000      1.41421      2.23607      2.23607
;      2.00000      2.23607      2.82843      2.82843
;      2.00000      2.23607      2.82843      2.82843

;The result is the same as the original array except for the
;truncation of the number of rows and columns to an even
;number.
;
;
;H. Rhody
;November, 2002
;-

;Set up the filter coefficient vectors
phi=[1,1]/sqrt(2)
psi=[-1,1]/sqrt(2)
;Determine the size of the input array. Somewhat
;different arrangements are needed for a row vector
;and an array.
n=size(image)
CASE n[0] OF
	1: BEGIN
	    nf=n[1]
	    nr=1
	   END
	2: BEGIN
	    nf=n[1]
	    nr=n[2]
	    END
	ELSE: Message,'Incompatible input array.'
END
f=float(image)
;This section does the inverse HWT. The forward
;HWT always produces an even number for nf, so
;this section does not need to worry about odd
;nf behavior.
IF Keyword_Set(inverse) THEN BEGIN
	g0=f[0:nf/2-1,*]; Low Frequency Part
	g1=f[nf/2:nf-1,*]; High Frequency Part
	f0=fltarr(nf,nr) & f0[2*indgen(nf/2),*]=g0
	f1=fltarr(nf,nr) & f1[2*indgen(nf/2),*]=g1
;	Because the scaling function phi=[1,1]/sqrt(2) and
;   the wavelet function psi=[-1,1]/sqrt(2) have the
;   [1,1] and [-1,1] quality, the convolution can be
;   carried out by shifting and adding or subtracting
;   the arrays. First set up the shifted arrays.
	f0a=[f0,fltarr(1,nr)]
	f0b=[fltarr(1,nr),f0]
	f1a=[f1,fltarr(1,nr)]
	f1b=[fltarr(1,nr),f1]
;   Now combine them to get the effect of convolution
	h0=f0a+f0b
	h1=f1a-f1b
	h=(h0+h1)/sqrt(2)
	return,h[0:nf-1,*]
	;
	;The following would be used if you needed to do a
	;true convolution, say with another kind of wavelet.
	;h0=convolve(f,phi)
	;h1=convolve(f,reverse(psi))
	;return,h[-:nf-1,*]

;This section does the forward HWT
ENDIF ELSE BEGIN
    ;construct the sampling index ds=[1,3,5,...]
	ds=2*indgen(nf/2)+1
;	Because the scaling function phi=[1,1]/sqrt(2) and
;   the wavelet function psi=[-1,1]/sqrt(2) have the
;   [1,1] and [-1,1] quality, the convolution can be
;   carried out by shifting and adding or subtracting
;   the arrays. First set up the shifted arrays.

	fa=[f,fltarr(1,nr)]
	fb=[fltarr(1,nr),f]
;   Now combine the shifted arrays to compute the convolutions
	g0=fa+fb
	g1=fb-fa
	;construct the return array with
	;the low-frequency part to the left.
	return,[g0[ds,*],g1[ds,*]]/sqrt(2)

	;
	;The following would be used if you needed to do a
	;true convolution, say with another kind of wavelet.
	;g0=convolve(f,phi)
	;g1=convolve(f,psi)
    ;return,[g0[ds,*],g1[ds,*]]

ENDELSE
end