function htrans,im

; return the 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]

	; force h to be float

	h=im+0.0

	; now do log2(n) reductions
	;
	; do first order outside loop (no divide by 2 in first step)

	hreduce,h,h0,hx,hy,hc

	; odd values get rounded up

	sx2=nx-1
	sy2=ny-1
	sx=(nx+1)/2
	sy=(ny+1)/2
	h[sx:sx2, 0:sy-1] = hx
	h[0:sx-1, sy:sy2] = hy
	h[sx:sx2, sy:sy2] = hc
	while (sx gt 1) or (sy gt 1) do begin
		hreduce,h0,nh0,hx,hy,hc
		sx2=sx-1
		sy2=sy-1
		sx=(sx+1)/2
		sy=(sy+1)/2
		if sx2 ge sx then h[sx:sx2, 0:sy-1] = hx/2
		if sy2 ge sy then h[0:sx-1, sy:sy2] = hy/2
		if sx2 ge sx and (sy2 ge sy) then h[sx:sx2, sy:sy2] = hc/2
		h0=reform(nh0/2,sx,sy)
	endwhile
	h[0,0] = h0
	return,h
end