function IRIS_PREP_GEOWAVE_CORRECT, img, params, $
	inhead = inhead, outhead = outhead, $
	cubic = cubic, boxwarp = boxwarp, update = update, $
	loud = loud, pstop = pstop, check = check, missing = missing, $
	keyword_only = keyword_only, nospread = nospread

;+
;
; Spun off from S. Jaeggli's iris_geowave_correct
;
; INPUT:
;	img	-	the 2D or 3D IRIS image array to be corrected
;	params	-	the structure describing the correction to be applied (see
;				IRIS_PREP_GET_GEOWAVE)
;
; KEYWORD:
;	cubic=	-	an integer specifying how the 2D polynomial interpolation
;				should be done
;	/loud	-	if set, then print some diagnostic information
;	inhead=	-	set to a FITS header (or array for 3D img data). If set, then 
; 				any CRS/ROI fields in the header are updated
;	outhead	-	a FITS header (or structure array) with updated header information
;	/update	-	if set, then the output header has some information added
;				to the HISTORY keyword
;	boxwarp -	Integer code specifying how the ROIs should be squared up following
;				the image warping. 
;					0: (default) the output ROI is the biggest rectangle that 
;						fits inside the transformed ROI;
;					1:	the output ROI is the smallest rectangle that holds
;						all valid data from the input ROI
;					2:	the output ROI is the same as the input ROI (not
;						generally a good idea, since the output range
;						may have moved as well as changed its shape and size)
;	/check	-	If set then the input image is displayed with the output ROI
;				boxes drawn over it for comparison
;	/pstop	-	This program brakes (breaks?) for debugging
;	missing=	Can be set to a value to be used to replace missing data. By 
;				default, NaNs are used.
;	nospread-	Passed through to reduce the spreading of infinite values; defaults
;				to 1, so it must be explicitly turned off if you want spreading
;
; RETURNS:
;	newimg	-	a 2D or 3D iris image array with the correction applied
;
; WRITTEN:
;	P.Boerner 2013/06 boerner ~ at ~ lmsal ~ dot ~ com
;
;-

keyword_only = KEYWORD_SET(keyword_only)
if not KEYWORD_SET(nospread) then nospread = 1

;+++++++++++++++++++++++++++++++++++++++++
; Currently handles array inputs by 
; calling itself recursively; not very 
; efficient, but we'll sort that out later
;-----------------------------------------
imsize = SIZE(img)
if imsize[0] gt 2 then begin
	t0 = SYSTIME(/sec)
	newimg = img
	for i = 0, imsize[3] - 1 do begin
		newimg[*,*,i] = IRIS_PREP_GEOWAVE_CORRECT(img[*,*,i], params, $
			inhead = inhead[i], outhead = thishead, cubic = cubic, boxwarp = boxwarp, $
			keyword_only = keyword_only)
		if i eq 0 then outhead = thishead else outhead = [outhead, thishead]
		PRINT, SYSTIME(/sec) - t0, form = '(f10.2)'
	endfor
	RETURN, newimg
endif

;+++++++++++++++++++++++++++++++++++++++++
; Set up defaults
;-----------------------------------------
if N_ELEMENTS(missing) eq 0 then nan = !values.f_nan else nan = missing
if N_ELEMENTS(cubic) eq 0 then cubic = -0.5

; Decide how to apply the correction based on the tags in the params structure
param_tags = TAG_NAMES(params)
dummy = WHERE(param_tags eq 'FUVS', isfuv)
if isfuv then param_tags = TAG_NAMES(params.FUVS)
poly2d = FIX(GETENV('IRIS_PREP_POLY2D'))

if (keyword_only eq 1) then begin
	newimg = 0.
	goto, warp_roi
endif

if N_ELEMENTS(inhead) eq 0 then begin
	imgsize = SIZE(img)
	if isfuv then maxsptrl = 4144 else maxsptrl = 2072
	sumsptrl = maxsptrl/imgsize[1]
	sumspat = 1096/imgsize[2]
endif else begin
	sumspat = inhead.sumspat
	sumsptrl = inhead.sumsptrl
endelse
numpix = N_ELEMENTS(img)

if (poly2d eq 0) then begin

	; Use INTERPOLATE
	if KEYWORD_SET(loud) then BOX_MESSAGE, ['IRIS_PREP_GEOWAVE_CORRECT:', $
		'Correcting image geometry using a distortion map']
	if not isfuv then begin
			dspec = 4144 / sumsptrl / 2
			dspat = 1096 / sumspat
			dmapx = REBIN(params.dmapx/sumsptrl, dspec, dspat) - (sumsptrl-1)/2./sumsptrl
			dmapy = REBIN(params.dmapy/sumspat, dspec, dspat) - (sumspat-1)/2./sumspat
			newimg = INFTERPOLATE(img, dmapx, dmapy, missing = nan, $
				cubic = cubic, inf=2e4, nospread = nospread)
	endif else begin
		; FUV images have to be split up
		dspec = 4144 / sumsptrl / 2
		dspat = 1096 / sumspat
		dmapx_s = REBIN(params.fuvs.dmapx/sumsptrl, dspec, dspat) - (sumsptrl-1)/2./sumsptrl
		dmapy_s = REBIN(params.fuvs.dmapy/sumspat, dspec, dspat) - (sumspat-1)/2./sumspat
		newimg_s = INFTERPOLATE(img[0:dspec-1,*], dmapx_s, dmapy_s, $
			cubic = cubic, missing = nan, inf=2e4, nospread = nospread)
		dmapx_l = REBIN(params.fuvl.dmapx/sumsptrl, dspec, dspat) - (sumsptrl-1)/2./sumsptrl
		dmapy_l = REBIN(params.fuvl.dmapy/sumspat, dspec, dspat) - (sumspat-1)/2./sumspat
		newimg_l = INFTERPOLATE(img[dspec:*,*], dmapx_l, dmapy_l, $
			cubic = cubic, missing = nan, inf=2e4, nospread = nospread)
		newimg = [newimg_s, newimg_l]
	endelse

endif else begin

	; Use POLY_2D
	interp = 2		;	Seems not to matter (cubic overrides it?)
	pixcen = 0.0
	if KEYWORD_SET(loud) then BOX_MESSAGE, ['IRIS_PREP_GEOWAVE_CORRECT:', $
		'Correcting image geometry using a 2D polynomial']
	if not isfuv then begin
		kx = params.kx
		ky = params.ky
		polysize = SQRT(N_ELEMENTS(kx))
		for i = 0., polysize - 1 do for j = 0., polysize - 1 do begin
			kx[i,j] = kx[i,j] * sumspat^i * sumsptrl^(j-1.)
			ky[i,j] = ky[i,j] * sumspat^(i-1.) * sumsptrl^j 
		endfor
		newimg = INF_POLY_2D(img, kx, ky, interp, inf=2e4, $
			cubic = cubic, missing = nan, pixel = pixcen, nospread = nospread)
	endif else begin
		; FUV images have to be split up
		kxs = params.fuvs.kx
		kys = params.fuvs.ky
		kxl = params.fuvl.kx
		kyl = params.fuvl.ky
		polysize = SQRT(N_ELEMENTS(kxs))
		for i = 0., polysize - 1 do for j = 0., polysize - 1 do begin
			kxs[i,j] = kxs[i,j] * sumspat^i * sumsptrl^(j-1.)
			kys[i,j] = kys[i,j] * sumspat^(i-1.) * sumsptrl^j 
			kxl[i,j] = kxl[i,j] * sumspat^i * sumsptrl^(j-1.)
			kyl[i,j] = kyl[i,j] * sumspat^(i-1.) * sumsptrl^j 
		endfor
		dspec = 4144 / sumsptrl / 2
		dspat = 1096 / sumspat
		newimg_s = INF_POLY_2D(img[0:dspec-1,*], kxs, kys, inf=2e4, $
			interp, cubic = cubic, missing = nan, pixel = pixcen, nospread = nospread)
		newimg_l = INF_POLY_2D(img[dspec:*,*], kxl, kyl, inf=2e4, $
			interp, cubic = cubic, missing = nan, pixel = pixcen, nospread = nospread)
		newimg = [newimg_s, newimg_l]
	endelse

endelse

if KEYWORD_SET(pstop) then STOP
warp_roi:

; Warp the header, and apply the resulting ROI mask
if N_ELEMENTS(inhead) gt 0 then begin
	if KEYWORD_SET(check) then begin
		if isnuv then begin
			WDEF, check + 10, inhead.naxis1, inhead.naxis2, title = 'IRIS_PREP_GEOWAVE_CORRECT CHECK'
			TVSCL, newimg>0
		endif else begin
			WDEF, check + 20, inhead.naxis1/2, inhead.naxis2/2, title = 'IRIS_PREP_GEOWAVE_CORRECT CHECK'
			TVSCL, REBIN(newimg>0, inhead.naxis1/2, inhead.naxis2/2)
		endelse
	endif
	outhead = IRIS_PREP_GEOWAVE_ROI(inhead, params, outmask = outmask, $
		boxwarp = boxwarp, cubic = cubic, check = check, update = update)
	maskedpoints = WHERE(outmask eq 0, nummasked)
	if nummasked gt 0 then newimg[maskedpoints] = nan
endif

RETURN, newimg

end