function hexpand,h0,hx,hy,hc
; perform one level of expansion of h0,hx,hy,hc; returns image im
; this version works on odd-sized images
;
; this probably could be improved using new IDL ability to specify step size
; in subscripts [i:j:2] -- but until that's been around longer it may not
; be safe to assume everyone has a version that can handle it.

	on_error, 2
	temp=size(h0)
	nx=temp[1]
	ny=temp[2]

	; nx2 = nx*2 or nx*2-1 depending on whether original image was even or odd,
	; and similarly for ny2
	;
	; hx or hy may be undefined if nx=nx2=1

	if n_elements(hx) gt 0 then begin
		temp=size(hx)
		nx2=nx + temp[1]
	endif else begin
		nx2 = nx
	endelse
	if n_elements(hy) gt 0 then begin
		temp=size(hy)
		ny2=ny + temp[2]
	endif else begin
		ny2 = ny
	endelse

	; force im to be float

	zero = h0[0]*0.0
	im=reform(replicate(zero,2,nx,2,ny),2,nx,2,ny)

	; pad with zeros if necessary for odd images
	;
	; The reform is necessary because if ny=1 replicate(zero,nx,ny) turns into
	; a 1-D array.
	; (I consider this an IDL bug, but I suppose it could be called a feature.)

	hhx=reform(replicate(zero,nx,ny),nx,ny)
	if n_elements(hx) gt 0 then hhx[0:nx2-nx-1,0:ny-1]=hx
	hhy=reform(replicate(zero,nx,ny),nx,ny)
	if n_elements(hy) gt 0 then hhy[0:nx-1,0:ny2-ny-1]=hy
	hhc=reform(replicate(zero,nx,ny),nx,ny)
	if n_elements(hc) gt 0 then hhc[0:nx2-nx-1,0:ny2-ny-1]=hc

	im[1,*,1,*]=reform(0.0+h0+hhx+hhy+hhc,nx,ny)
	im[1,*,0,*]=reform(0.0+h0+hhx-hhy-hhc,nx,ny)
	im[0,*,1,*]=reform(0.0+h0-hhx+hhy-hhc,nx,ny)
	im[0,*,0,*]=reform(0.0+h0-hhx-hhy+hhc,nx,ny)
	im=reform(im,2*nx,2*ny)

	; return odd image without padding

	return,reform(im[0:nx2-1,0:ny2-1],nx2,ny2)
end