function htrans_var,var

; return the variance of the H-transform of a 2-d image
; works on images of any size
; uses floating point arithmetic

	on_error, 2
	temp = size(var)
	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 = var+0.0

	; now do log2(n) reductions
	;
	; do first order outside loop so we can do rounding differently
	; (no divide by 2 in first step)

	hreduce,h,h0,hx,hy,hxy

	; odd values get rounded up

	sx2 = nx-1
	sy2 = ny-1
	sx = (nx+1)/2
	sy = (ny+1)/2

	; increase variances at edge if section size is odd because edge pixels
	; are used in sum twice

	if sx*2 ne sx2+1 then h0[sx-1, *] = h0[sx-1, *]*2
	if sy*2 ne sy2+1 then h0[*, sy-1] = h0[*, sy-1]*2

	h[sx:sx2, 0:sy-1] = h0[0:sx2-sx, *]
	h[0:sx-1, sy:sy2] = h0[*, 0:sy2-sy]
	h[sx:sx2, sy:sy2] = h0[0:sx2-sx, 0:sy2-sy]
	while (sx gt 1) or (sy gt 1) do begin
		hreduce,h0,nh0,hx,hy,hxy
		sx2 = sx-1
		sy2 = sy-1
		sx = (sx+1)/2
		sy = (sy+1)/2

		; increase variances at edge if section size is odd

		if sx*2 ne sx2+1 then nh0[sx-1, *] = nh0[sx-1, *]*2
		if sy*2 ne sy2+1 then nh0[*, sy-1] = nh0[*, sy-1]*2

		if sx2 ge sx then h[sx:sx2, 0:sy-1] = nh0[0:sx2-sx, *]/4
		if sy2 ge sy then h[0:sx-1, sy:sy2] = nh0[*, 0:sy2-sy]/4
		if sx2 ge sx and sy2 ge sy then $
			h[sx:sx2, sy:sy2] = nh0[0:sx2-sx, 0:sy2-sy]/4
		h0 = nh0/4
	endwhile
	h[0,0] = h0
	return,h
end