pro hreduce,im,h0,hx,hy,hc

; perform one level of reduction of image im
; works on odd-sized images

	on_error, 2
	temp = size(im)
	nx = temp[1]
	if temp[0] eq 1 then ny = 1L else ny = temp[2]

	; round nx/2,ny/2 up if nx,ny odd

	nx2 = (nx+1)/2
	ny2 = (ny+1)/2

	; add row or column pad to make images even in size
	; force h to be float

	h = replicate(im[0]*0.0,2*nx2,2*ny2)
	h[0:nx-1,0:ny-1] = im
	if 2*nx2 ne nx then h[nx,*] = h[nx-1,*]
	if 2*ny2 ne ny then h[*,ny] = h[*,ny-1]

	h = reform(h,2,nx2,2,ny2)
	h0 = reform(h[1,*,1,*]+h[1,*,0,*]+h[0,*,1,*]+h[0,*,0,*],nx2,ny2)
	hx = reform(h[1,*,1,*]+h[1,*,0,*]-h[0,*,1,*]-h[0,*,0,*],nx2,ny2)
	hy = reform(h[1,*,1,*]-h[1,*,0,*]+h[0,*,1,*]-h[0,*,0,*],nx2,ny2)
	hc = reform(h[1,*,1,*]-h[1,*,0,*]-h[0,*,1,*]+h[0,*,0,*],nx2,ny2)

	; convert back to odd sizes if necessary
	; don't bother if size=1 (no difference coefficients required then)

	nnx=nx2
	nny=ny2
	if (2*nx2 ne nx) and (nx2 gt 1) then begin
		nnx = nx2-1
		hx = reform( hx[0:nnx-1,*], nnx, ny2)
		hc = reform( hc[0:nnx-1,*], nnx, nny)
	endif
	if (2*ny2 ne ny) and (ny2 gt 1) then begin
		nny = ny2-1
		hy = reform( hy[*,0:nny-1], nx2, nny)
		hc = reform( hc[*,0:nny-1], nnx, nny)
	endif
	return
end