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