;+
; ZUNWRAP
;
; "Unwrap" a solar image, converting X to angle and Y to radius. 
; Use zra2xy to do the transformation.  
;
; Scanning is done widdershins (CCW) over a specified angular range
; in the original image plane; the widdershins scan is converted
; to a horizontal one in the resulting image plane (ie moving
; in the +X direction on the final image plane moves in the CCW
; direction on the original image plane).  Radius is mapped to the
; Y axis.  This preserves angular sign conventions but mirror-reverses
; the image that is being unwrapped.
;
;;
; CALLING SEQUENCE:
; zunwrap,in=in,ihdr=ihdr,out=out,ohdr=ohdr,[optional keywords]
;
; REQUIRED KEYWORDS
;
; in, ihdr - the fits image and header to unwrap about the origin
; out,ohdr - the destination image and header
;
; OPTIONAL KEYWORDS
;
; RMIN - the minimum radius, in solar radii, of the unwrapped image
; RMAX - the maximum radius, in solar radii, of the unwrapped image
; RSTEP - the radial step, in pixels(!) of the input image
; ARANGE - A vector specifying the range (in degrees clockwise from North)
;          of the start and end of the scan.  Default is [0,360].
; ASTEPS - the number of steps around the specified 'arange'.
;
; SMOOTH - binary keyword.  If set, a slower algorithm is used that averages
; 	pixel values where relevant.  Doesn't work on conformal unwrapping
;
; LERP - Binary keyword.  If set, use bilinear interpolation on the original
;     image plane (rather than sampling) to get the destination pixel values.
;
; CONFORMAL - Binary keyword.  If set, make the transformation conformal
;   (by putting 'r' on a log scale).
;
; RADIUS - floating point.  Used with conformal keyword, specifies 
; (in output image plane pixel widths) the size of the smoothing kernel
; being applied to the input image.  (The smoothing kernel size is
; constant in the output image plane, meaning that it is variable in the
; input image plane).  Specify 0.0 for sampling, 1.0 for quick-n-dirty
; anti-moire averaging, and higher to smooth the output image more.
;
; FAST - Binary keyword - Sets "radius" to zero.
;
; MISSING - If set, contains the value which should signal missing datapoints.
;
; GRID - Binary keyword.  If set, the conformal unwrapper engine 
; 	adds a grid indicating 10s of degrees of azimuth and individual
;	solar radii on the output image.  See also GRMASK.  If this is 
;	DEFINED, then a GRMASK is produces; if it is SET, then the grid is
;	superposed on the output image. 
;
; GRMASK - output keyword. A binary mask containing 255 wherever a gridpoint
;	goes.
;
; EXAMPLE:
;	zunwrap,in=a,ihdr=ahdr,out=b,ohdr=bhdr,/conformal,rmin=1.0,rmax=6.0,$
;		asteps=600,arange=[-180,180]
;
; HISTORY:
; 	Written by Craig DeForest, early 1993.
;
; 	29-Jan-97 - Fixed up the bilinear interpolation in the conformal
; 	mapping code (Was fixing the zra2xy output *then* doing bilinear.
; 	Bad.  Definitely bad.)
;
; 	12-Jan-97 - Added conformal mapping stuff; this is slow but very cool.
; 	It really ought to be vectorized in the opposite direction.
;       Added "arange" option:  you can specify a range of angles 
;	to cover (eg "arange=[90,270]" to get the South pole).
;       Also added "rot" option; "rot" lets you specify the 
;	location of the splice, in degrees.  (Default is zero).  This must be 
;	in the arange.
;
;	20-Aug-97 - Upgraded to use structure, rather than string-array, 
;	headers.
; 	
;	3-Jan-97 - Added differential smoothing option, rather than "mere" 
;	sampling. If you set the /smooth flag, then pixels in the source 
;	plane are averaged over regions roughly the size of the relevant 
;	destination-plane pixels.
;	
;	13-Mar-98 - Fixed up to use CRPIX correctly (off-by-1)
;	16-Jun-98 - Hacked on conformal headers; added 'missing' keyword.
;	26-Mar-99 - Removed references to deprecated "reflect" keyword.
;	20-Apr-99 - Improved header docs slightly; modified rmin default for linear
;			case.
;-

pro zunwrap,in=in,ihdr=iihdr,out=out,ohdr=ohdr,rmin=rmin,rmax=rmax,rstep=rstep $
	,asteps=asteps,smooth=smooth,monitor=monitor,monout=monout,lerp=lerp $
	,foo_out=foo_out,conformal=conformal,radius=smoothing,arange=arange $
	,rot=rot,grid=grid,grmask=grmask,fast=fast,missing=missing

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Cheezy fake-vectorization recursion goes here:  If we get more than
;; one image, we call ourselves multiple times.  This could be cleaned
;; up, but what the heck -- it works for now, and the overhead isn't
;; too humongous...
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if(nlm(iihdr) gt 1) then begin
	ifoohdr = iihdr(0)
	zunwrap,in=in(*,*,0),ihdr=ifoohdr,out=foo,ohdr=foohdr,rmin=rmin,rmax=rmax,rstep=rstep,asteps=asteps,smooth=smooth,monitor=monitor,monout=monout,lerp=lerp,foo_out=foo_out,conformal=conformal,radius=smoothing,arange=arange,rot=rot,grid=grid,grmask=grmask,fast=fast
	out = replicate(foo(0),nlm(foo(*,0)),nlm(foo(0,*)),nlm(iihdr))
	out(*,*,0) = foo
	ohdr = replicate(foohdr,nlm(iihdr))

	for i=1,nlm(iihdr)-1 do begin
		ifoohdr = iihdr(i)
		zunwrap,in=in(*,*,i),ihdr=ifoohdr,out=foo,ohdr=foohdr,rmin=rmin,rmax=rmax,rstep=rstep,asteps=asteps,smooth=smooth,monitor=monitor,monout=monout,lerp=lerp,foo_out=foo_out,conformal=conformal,radius=smoothing,arange=arange,rot=rot,grid=grid,grmask=grmask,fast=fast
		out(*,*,i) = foo
		ohdr(i) = foohdr
	end
	return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

ihdr = iihdr				; Copy ihdr so as not to mess it up
v = keyword_set(verbose) 
if(v) then openu,44,'/dev/tty',error=err
	
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Initial setup ... check and generate defaults for the various keywords
;;

if ((not isvalid(in)) or (not isvalid(ihdr))) then begin
	print,"You must specify a valid input image and header!"
	return
end

w = ihdr.NAXIS1
h = ihdr.NAXIS2

xy0 = [(ihdr.CRPIX1-1), (ihdr.CRPIX2-1)]

zhdrunit,ihdr,'solar-radii'
pix2radii = ihdr.CDELT1
radii2pix = 1.0 / pix2radii

if (not isvalid(arange)) then arange = [0.0,360.0] else arange = float(arange)
if (not isvalid(asteps)) then asteps=720
if (not isvalid(rmin)) then begin
	if keyword_set(conformal) then rmin = 0.5 else rmin = 0.0
end

;;;;;;  Default rmax is the largest radius in the image plane
if (not isvalid(rmax)) then begin
	xys =  transpose([ [0,0,w,w] , $
	 	 [0,h,0,h] ])
	corner_rs = fltarr(4)
	for i=0,3 do corner_rs(i) = cartdist(xy0,xys(*,i))

	rmax = max(corner_rs) * pix2radii
end

;;;;;  rstep defaults to one pixel in the input plane...
;;;;;  rpix is rstep in pixels on the input plane.
if (not isvalid(rstep)) then rstep=pix2radii * 1.0
rpix = fix(0.5 + rstep*radii2pix)

if (not isvalid(smoothing)) then smoothing = (1.0+sqrt(2))/2.0
if keyword_set(fast) then smoothing = 0

if(keyword_set(monitor)) then begin
	monout = fltarr(ihdr.naxis1,ihdr.naxis2)
	wmatch,15,monout

	foosiz = 500
	foo_out = fltarr(foosiz,foosiz)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; conformal case:  we do some fancy footwork...
; To make the transformation conformal, we have to put the radial scale
; on a log scale.  
;
; In this case, we almost certainly want to do some averaging over larger
; and larger areas, to take advantage of the image degradation we're doing
; toward the edges of the image.  Since the mapping is conformal, we 
; know that small shapes on the final image will be the same shapes (but
; rotated) on the initial image.  We sample over small circles, which have
; the advantage of rotational symmetry -- so that we don't have to rotate
; (as well as scale) the sampling zone.
;

if(keyword_set(conformal)) then begin
	; For conformal mapping, 'rstep' is the step in log-r space.  We
	; adjust it to match the "asteps" argument.  At r=1, a=0, then
	; dr/dy = 1 and da/dx = 1/2pi.  So rstep should be 2pi * astep
	; there.   (Note that 'astep' is the (sort-of) reciprocal of 'asteps',
	; corrected for the arc length of the "arange".)
	rstep =  (arange(1) - arange(0)) / !radeg / asteps

	if(rmin le 0) then rmin = rstep   ; Make sure we got no singularity...

	ow = asteps
	oh = fix(alog(rmax/rmin)/rstep) + 1

	print,"Output image:  width = ",fug(asteps),"; height = ",fug(oh)

	ohdr = ihdr
	ohdr.CTYPE1=	'degrees'
	ohdr.CDELT1=	(arange(1)-arange(0)) / asteps
	ohdr.CTYPE2=	'alog solar-radii'
	ohdr.CDELT2=	rstep
	ohdr.CRPIX1=	(-arange(0) * asteps / (arange(1) - arange(0))) + 1
	ohdr.CRPIX2=	1
	ohdr.NAXIS1=	ow
	ohdr.NAXIS2=	oh
	ohdr.CRVAL1 = 0
	ohdr.CRVAL2 = alog(rmin)

	out = fltarr(ohdr.naxis1,ohdr.naxis2)
	if isvalid(missing) then out(*) = missing

	
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;; Coordinate arrays.  "as" will be filled with copies of
	;; the current angle for each loop iteration; "rs" and 
	;; "sizes" hold the 'r' location of each pixel in the 
	;; source image plane (in solar radii) and size of each 
	;; region to be snarfed (in pixels in the source image
	;; plane).
	;;
	;; "Weights" holds the weighting information for boundary
	;; pixels, in squares with the lower left corner at (0,0) in
	;; each layer of the cube.
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	sizes = fltarr(oh)
	rs = exp(findgen(oh)*rstep)*rmin
	sizes(1:oh-2) =	(((rs(2:oh-1) - rs(0:oh-3))/4.0) * radii2pix * smoothing) > 0 
	sizes(0) = 0
        sizes(oh-1) = sizes(oh-2)

	if(keyword_set(monitor)) then begin
		window,14,xsiz=400,ysiz=300
		plot,sizes
		wset,15
	end

	msiz = fix(max(sizes)+1)
	if(max(fix(sizes)) eq 0) then msiz = 0

;	sz = fix(sizes)+1

	; 
	; On a higher-resolution grid, draw a circle in single bits.
	; Then rebin the high-resolution grid back down to weighted
	; pixel averages.
	;
	wn = 5

	foo = 2*msiz+1
	weights = fltarr(foo,foo,oh)
	weights2 = fltarr(foo,foo,oh)

	sketchpad = bytarr(wn*foo,wn*foo)
	skxy = [ transpose(indgen(n_elements(sketchpad)) /   (wn*foo)) ,$ 
	         transpose(indgen(n_elements(sketchpad)) mod (wn*foo)) ]
	skc = skxy
	
	if(v) then writeu,44,"Calculating kernels..."
	
	sizes = sizes > 1.0
	weights(msiz,msiz,*) = 1
	weights2(msiz,msiz,*) = 1

	for i=max(where(sizes eq 1.0)),oh-1 do begin
		skc(*) = msiz*wn+fix(wn/2)
		sketchpad(skxy(0,*),skxy(1,*)) = 100 * (cartdist(skxy,skc) le sizes(i)*wn)
		weights2(*,*,i) = rebin(sketchpad,foo,foo)/100.0
		weights (*,*,i) = weights2(*,*,i) / total(weights2(*,*,i))

		if keyword_set(monitor) then begin
			tvscl,rebin(weights(*,*,i),5*wn*foo,5*wn*foo,/sample)
		end
	end
	if(v) then writeu,44,"done -- kernel size is ",string(foo,format="(I3.3)"),byte(13),byte(10)

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;; Begin of actual transformation...  we make it bilinear
	;; by pixel.  The innermost loop is over the "weights" array.
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	xs = intarr(oh)
	ys = intarr(oh)
	msk = bytarr(oh)

;	in00 = in(0,0)
;	in(0,0) = 0

	for i=0,asteps - 1 do begin

		;; Calculate the angle we're doing now, and then convert to
		;; (x,y) coordinates.
		a = ( arange(0) + i * (arange(1) - arange(0)) / asteps )
		xys = zra2xy(rs,a,ihdr) - msiz

		for xx = 0,foo-1 do for yy = 0,foo-1 do begin
				 
		    n = where ( (weights(xx,yy,*) gt 0) and  $
				((xys(*,0) ge 0) and         $
			         (xys(*,0) lt w) and         $
			         (xys(*,1) ge 0) and         $
		 	         (xys(*,1) lt h)             $
				) )
		    ni = n_elements(n)-1
		    if(ni ge 1) then $
			out(i,n(*)) = out(i,n(*)) + weights(xx,yy,n(*))*interpolate(in,xys(n,0)+xx,xys(n,1)+yy)
		
		end

		if(i mod 10 eq 0 or i eq asteps-1) then begin
			if(v) then writeu,44,byte(13),fug(fix(100*(i/(asteps - 1.0)))),"%, ",fug(i)," rasters completed ..."
			if(keyword_set(monitor)) then begin
				wset,15
				tvscl,monout<10
				wset,14
				tvscl,out
			end
		end
	end		

	if(isvalid(grid)) then begin
		if(where((size(grmask))(0:2) ne (size(out))(0:2)))(0) ne -1 then $
			grmask = bytarr(nlm(out(*,0)),nlm(out(0,*)))

		for i=fix(rmin)+1,fix(rmax) do begin
			y0 = alog(float(i)/rmin)/rstep
			for j=0,ow/10-1 do grmask(j*10,y0) = 1
		end

		for i=0,fix(arange(1) - arange(0))/10-1 do begin
			x0 = (i*10)*float(asteps)/(arange(1) - arange(0))
			for j=0,fix(oh/10)-1 do grmask(x0,j*10) = 1
		end

		if keyword_set(grid) then begin
			m = max(out)+(256.0/255.0)*(max(out) - min(out))
			out(where(grmask)) = m
		end
	end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;  End of conformal map section
;;  Beginning of linear-radial map section
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

end else begin
	out = fltarr(asteps,((rmax - rmin)/rstep))
	if isvalid(missing) then out(*)=missing

	ohdr = ihdr
	ohdr.CTYPE1 = 'degrees'
	ohdr.CDELT1 = (arange(1)-arange(0)) / asteps
	ohdr.CDELT2 = rstep	; (ohdr.CTYPE2 is already 'solar-radii')
	ohdr.CRPIX1 = (-arange(0) * asteps / (arange(1) - arange(0))) + 1
	ohdr.CRPIX2 = -(rmin/rstep) + 1
	ohdr.NAXIS1 = asteps
	ohdr.NAXIS2 = fix((rmax-rmin)/rstep)
	
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	if(keyword_set(smooth)) then begin
	;
	; Smoothing case:  We unwrap at greater than full resolution, then 
	; average down to scale...
	;
		circumference = rmaxpix * 3.1415926 * 2
		multiplier = fix(circumference / asteps) + 1
		zunwrap,in=in,ihdr=ihdr,out=foo,ohdr=foohdr,rstep = rstep / 2.0, asteps=asteps * multiplier,monitor=monitor,monout=monout,smooth=0,rmin=rmin,rmax=rmax
		hgt = (size(foo))(2)
		hgt2 = hgt - (hgt mod 2)
		out = rebin(foo(*,0:(hgt2 - 1)),asteps,hgt2/2)
		ohdr = foohdr
		ohdr.naxis1 = asteps
		ohdr.cdelt1 = ohdr.cdelt1 * multiplier
		ohdr.naxis2 = (size(foo))(2)/2.0
		ohdr.cdelt2 = ohdr.cdelt2 * 2.0
	
	end else begin
	;
	; Non-smoothing case:  Just do the sampling.  If the interpolation
	; flag is set, we do a bilinear.
	;
	
		; Sampling unwrap algorithm
		rs = float(indgen((rmax-rmin)/rstep) * rstep) + rmin
		as = float(rs)
		
		for i=0,asteps-1 do begin
			as(*) = ( arange(0) + i * (arange(1) - arange(0)) / asteps )
			xys = zra2xy([transpose(rs),transpose(as)],ihdr);
		
			foo = where(((xys(*,0) ge 0) and $
				  	 (xys(*,0) lt w-1) and  $
			  	  	 (xys(*,1) ge 0) and $
				  	 (xys(*,1) lt h-1) ))
			if(foo(0) ge 0) then begin
		;;;;;;;;;;;;;;;;;;;;;;;;;
				if(keyword_set(lerp)) then $
					out(i,foo) = interpolate(in,xys(foo,0),xys(foo,1)) $
				else $
					out(i,foo) = in(xys(foo,0),xys(foo,1))
		;;;;;;;;;;;;;;;;;;;;;;
		
				if(i mod 10 eq 0 or i eq asteps-1) then begin
					if(v) then writeu,44,byte(13),fug(fix(100*(i/(asteps - 1.0)))),"%..."
					if(keyword_set(monitor)) then tvscl,monout<4
				end
			end
	
		end
		if(v) then printf,44,byte(13),"done."
		if(v) then close,44
	
		if(keyword_set(monitor)) then begin
			wmatch,16,foo_out
			tvscl,foo_out
		end
	end
end

;sxaddpar,ohdr,'HISTORY','Unwrapped into radial coords'
;sxaddpar,ohdr,'HISTORY','  - rstep='+fug(sxpar(ohdr,'CDELT2')*zunits(sxpar(ohdr,'CTYPE2'),'solar-radii'))+' solar radii; '+fug(asteps)+' angular steps'

end