function hinv,im

; Return the inverse H-transform of a 2-d image.
; Works on images of any size.
; Uses floating point arithmetic rather than integer.

	on_error, 2
	temp = size(im)
	if temp[0] ne 2 then message,'Only works on 2-D images'
	nx = temp[1]
	ny = temp[2]

	; have to calculate sequence of sx,sy first because it must be computed
	; from the biggest values down
	;
	; odd values get rounded up

	sx2 = [nx-1]
	sy2 = [ny-1]
	sx = [(nx+1)/2]
	sy = [(ny+1)/2]
	ns = 1
	while (sx[0] gt 1) or (sy[0] gt 1) do begin
		sx2 = [sx[0]-1 , sx2]
		sy2 = [sy[0]-1 , sy2]
		sx = [(sx[0]+1)/2 , sx]
		sy = [(sy[0]+1)/2 , sy]
		ns = ns+1
	endwhile

	; now do log2(n) expansions
	; force h to be float

	h = im+0.0
	for i = 0, ns-2 do begin
		h0 = reform(h[0:sx[i]-1,0:sy[i]-1], sx[i], sy[i])
		if sx2[i] ge sx[i] then $
			hx = reform(h[sx[i]:sx2[i],0:sy[i]-1], sx2[i]-sx[i]+1, sy[i])
		if sy2[i] ge sy[i] then $
			hy = reform(h[0:sx[i]-1,sy[i]:sy2[i]], sx[i], sy2[i]-sy[i]+1)
		if sx2[i] ge sx[i] and sy2[i] ge sy[i] then $
			hc= reform(h[sx[i]:sx2[i],sy[i]:sy2[i]], $
				sx2[i]-sx[i]+1, sy2[i]-sy[i]+1)

		h[0:sx2[i],0:sy2[i]] = hexpand(h0,hx,hy,hc)/2
	endfor

	; do last expansion separately to change scaling

	i = ns-1
	h0 = reform(h[0:sx[i]-1,0:sy[i]-1], sx[i], sy[i])
	if sx2[i] ge sx[i] then $
		hx = reform(h[sx[i]:sx2[i],0:sy[i]-1], sx2[i]-sx[i]+1, sy[i])
	if sy2[i] ge sy[i] then $
		hy = reform(h[0:sx[i]-1,sy[i]:sy2[i]], sx[i], sy2[i]-sy[i]+1)
	if sx2[i] ge sx[i] and sy2[i] ge sy[i] then $
		hc = reform(h[sx[i]:sx2[i],sy[i]:sy2[i]], $
			sx2[i]-sx[i]+1, sy2[i]-sy[i]+1)
	h[0:sx2[i],0:sy2[i]] = hexpand(h0,hx,hy,hc)/4
	return,h
end