;+
; NAME: ircal_dewarp 
; PURPOSE:
; 	Remove distortion from an IRCAL image
;
; NOTES:
; 	Right now this is -incredibly- basic and just handles the anisomorphic
; 	magnification term, and nothing higher order than that. 
;
; INPUTS:
; 		image0			an image or stack of images.
; KEYWORDS:
; 		/NOFLUXSCALE	Don't rescale image counts to preserve total flux.
; 						(default is to do scaling)
; 		PLATESCALE=		new platescale to use. Default is 0.040 arcsec
; 		CUBIC=			use cubic resampling. See docs for INTERPOLATE.
; 		/BYTE			return a byte image, not float. Useful for dewarping
; 						masks
; 		/CROPTOP		Drop the top 5 rows of the IRCAL detector. These are the
; 						rows full of many bad pixels and sometimes it's just
; 						easiest to get rid of them.
; OUTPUTS:
;
; HISTORY:
; 	Began 2005-11-07 15:24:51 by Marshall Perrin 
; 	2006-05-09	Actually got it working finally. Major code updates.
; 	2006-06-09	Added /croptop option
;-

function ircal_dewarp,image0,nofluxscale=nofluxscale,platescale=platescale,$
	cubic=cubic,byte=byte,croptop=croptop


	@instruments.include
	scalex = ircal.platescale
	scaley = ircal.platescale_y

	; what should be the new (regridded) platescale?
	;if ~(keyword_set(platescale)) then platescale=scalex
	if ~(keyword_set(platescale)) then platescale=0.040 ; Nyquist at J for Lick.
							; TODO would this be better at twice Nyquist?
	

	image = image0
	sz = size(image)

	; because pixel coordinates are the CENTERS of pixels, an image with
	; N pixels is actually (N-1)*platescale units wide. 
	;  (well, not counting the half-pixel border on either side)
	
	; new image pixel counts
	nsx = floor((sz[1]-1)*scalex/platescale)+1
	nsy = floor((sz[2]-1)*scaley/platescale)+1
	; should we just drop the top 5 rows (with tons of bad pixels?)
	ncrop=5
	if keyword_set(croptop) then  nsy = floor((sz[2]-ncrop-1)*scaley/platescale)+1
		
	; handle 3D arrays if necessary via recursion
	if sz[0] eq 3 then begin
		out = dblarr(nsx,nsy,sz[3])
		for i=0L,sz[3]-1 do begin
			statusline, "Dewarping "+strc(i+1)+"/"+strc(sz[3])
			out[0,0,i] = ircal_dewarp(image[*,*,i],nofluxscale=nofluxscale,$
				platescale=platescale,cubic=cubic,croptop=croptop)
		endfor 
		return,out
	endif

	; now actually handle 2D arrays.
	;
	; indices of new pixels relative to old
	ix = findgen(nsx)*platescale/scalex 
	iy = findgen(nsy)*platescale/scaley

	; can't use interpolate with NANs, so fix if needed. 
	wnans = where(~finite(image),nanct)
	if nanct gt 0 then begin
		nanmask = finite(image)
		;TODO add an option for more clever nan repair here. 
		image[wnans] = median(image)
		ix2 = ix ; needed because bilinear stomps on its input IX and IY arrays.
		iy2 = iy
		newnanmask = bilinear(nanmask,ix2,iy2)
	endif

	; On tests with linear ramp data, the following code works very well on all
	; central pixels, but has slight errors on pixels within 2 rows of the edge.
	; This is presumably due to how the cubic interpolater deals with edge
	; effects. For now I am going to IGNORE this.

	;TODO how well will this deal with IRCAL's undersampled data? Is there something
	; we can do which will improve this?
	; ANS: it definitely rings off of bad pixels. But it does OK on other stuff,
	; as far as I can tell. So do the best job of removing hot/cold/bad pixels
	; first and this should be OK.

	;im_out = interpolate(image,ix,iy,cubic=-0.5,/grid)
	im_out = interpolate(image,ix,iy,cubic=cubic,/grid)

if ~(keyword_set(nofluxscale)) then begin
	; Rescale the image to preserve the total flux
	; 
	; Pixel area conversion factor, to preserve flux
	pixelarearatio = platescale^2 / (scalex*scaley)

	im_out *= pixelarearatio

endif 

; TODO James points out that for just rescaling the axes, with no spatially
; varying distortions, then one can use the Jacobian to compensate for the
; change of variables between the different coordinate systems.
; See eqn 6 on http://mathworld.wolfram.com/Jacobian.html
; But in practice here - we distort all pixels equally, so it's just some
; scalar coefficient, and then we measure our photometry on standards
; reduced using the code, so the fixed scalar coefficients should cancel out
; in the end. - MDP 2006-06-13

	if nanct gt 0 then begin
		newwnan = where(newnanmask eq 0)
		im_out[newwnan] = !values.f_nan
	endif

	if keyword_set(byte) then im_out = byte(round(im_out))

	return,im_out

end