function IRIS_PREP_GEOWAVE_ROI, inhead, params, $
	outmask = outmask, cubic = cubic, boxwarp = boxwarp, $
	check = check, update = update
;+
;
; When the IRIS images are warped to correct for geometric/wavelength distortion,
; the square sub-regions read out from the camera map to non-square regions in the
; level 1.5 image. This routine re-rectangularizes the regions of interest (ROIs)
;
; Applys the geometric/wavelength correction to the header CRS regions and returns
; a FITS header with the updated region specifications. Also generates
; a mask that can be applied to the image itself to enforce the output ROIs
;
; INPUT:
;	inhead	-	Input FITS header structure, including the CRS definitions
;				(Does not get modified)
;	params	-	The structure describing the geometric/wavelength correction
;				to be applied to the image (see IRIS_PREP_GEOWAVE_CORRECT)
;
; KEYWORDS:
;	outmask=	-	An integer image with 1 for all pixels included in the updated
;					rectangular ROIs and 0 everywhere else
;	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 output regions are drawn on the input image
;				(assumes that it was called from IRIS_PREP_GEOWAVE_CORRECT(/check)
;	cubic=	-	passed through to IRIS_GEOWAVE_CORRECT	
;	/update	-	If set, then the HISTORY keyword is updated 
;
; RETURNS:
;	outhead	-	a FITS header structure with the CRS keywords updated
;
; WRITTEN:
;	P.Boerner 2013/06 boerner ~ at ~ lmsal ~ dot ~ com
;
;-

; Currently only handles scalar structures for inhead and params

; Create dummy outputs
outhead = inhead
inmask = INTARR(inhead.naxis1, inhead.naxis2)
outmask = inmask
this_inmask = inmask
this_outmask = inmask

; Set up common block used to invert warping
common IRIS_GEOWAVE_ROI_CB, cb_params, target_x, target_y
cb_params = params

; Set up defaults and constants
numreg = inhead.crs_nreg
if numreg eq 0 then begin
	; Don't do anything if this is a full-frame image
	outmask = outmask + 1
	RETURN, outhead
endif
outreg = INTARR(numreg, 4)
corners = ['LL', 'UL', 'LR', 'UR']
if ANYTIM2TAI(inhead.t_obs) lt ANYTIM2TAI('2013-07-01') then begin
	isold = 1
	tagpart = 'CRS_'
	fits_tag_parts = ['SR', 'SC', 'ER', 'EC']
endif else begin
	isold = 0
	tagpart = 'T'
	fits_tag_parts = ['SC', 'SR', 'EC', 'ER']
endelse
minoutbox = [0,0,0,0]
squashbox = 0
maxoutbox = [inhead.naxis1, inhead.naxis2, inhead.naxis1, inhead.naxis2] - 1
tighten_outbox = [0, 1, 0, -1]
xs = DBLARR(4)
ys = DBLARR(4)
tagval = INTARR(4)
if not KEYWORD_SET(cubic) then cubic = -0.5
if N_ELEMENTS(boxwarp) eq 0 then boxwarp = 1

; Check whether the image is FUV (in which case it uses two different transforms)
if inhead.img_path eq 'FUV' then isfuv = 1 else isfuv = 0
if inhead.img_path eq 'NUV' then isnuv = 1 else isnuv = 0

; Loop through the populated CRS regions and add them to the mask
for i = 1, numreg do begin

	for j = 0, 3 do $
		tagval[j] = GT_TAGVAL(inhead, tagpart + fits_tag_parts[j] + STRTRIM(i,2))
	if isnuv and isold then tagval = tagval - 2072 * [1,0,1,0]
	tagval = tagval - 1 	;	CRS entries count from 1, not 0

	; Find the rectangular box that spans the transformed region
	
	; samebox is trivial; just use the input CRS parameters
	if (boxwarp eq 2) then begin
		outbox = tagval
	endif

	; minbox is a little trickier; find it by looking for the corners
	; This may not always work for very small regions and/or very distorted images!
	if (boxwarp eq 0) then begin
		for j = 0, 3 do begin
			target_x = tagval[(j/2)*2]
			target_y = tagval[(j mod 2) * 2 + 1]
			flipcoords = IRIS_PREP_INTERPOLATE_NEWTFUNC([target_x, target_y])
			xs[j] = target_x - flipcoords[0]
			ys[j] = target_y - flipcoords[1]
;			outcoords = NEWTON([target_x, target_y], 'iris_prep_interpolate_newtfunc')
;			xs[j] = outcoords[0]
;			ys[j] = outcoords[1]
		endfor
		outbox = [CEIL(xs[1]>xs[0]), CEIL(ys[0]>ys[2]), FLOOR(xs[2]<xs[3]), FLOOR(ys[1]<ys[3])]
		outbox = outbox + tighten_outbox
		if (outbox[2] le (outbox[0]+1) ) or (outbox[3] le (outbox[1]+1) ) then begin
			BOX_MESSAGE, ['IRIS_PREP: ROI too small/distorted to use /minbox...', $
				'Using /maxbox for this region (consider specifying it for whole image)']
			squashbox = 1
		endif
	endif

	; maxbox is easy; just take all rows and columns with nonzero values
	if (boxwarp eq 1) or (squashbox eq 1) then begin
		; Add the region to the input mask image
		this_inmask[tagval[0]:tagval[2], tagval[1]:tagval[3]] = 1
		inmask = inmask + this_inmask
		
		; Transform the inmask to make an outmask
		this_outmask = IRIS_PREP_GEOWAVE_CORRECT(this_inmask, params, cubic = cubic, missing = 0)
		this_inmask = this_inmask * 0
		outmask_profile1 = TOTAL(this_outmask, 1)
		outmask_profile2 = TOTAL(this_outmask, 2)
		goodcols = WHERE(outmask_profile1 gt 0, numgoodcols)
		goodrows = WHERE(outmask_profile2 gt 0, numgoodrows)
		outbox = [goodrows[0], goodcols[0], goodrows[numgoodrows-1], goodcols[numgoodcols-1]]
		squashbox = 0
	endif

	; Update the outmask and the CRS FITS keywords to match the selected outbox
	outbox = minoutbox > outbox < maxoutbox
	outmask[outbox[0]:outbox[2], outbox[1]:outbox[3]] = 1
	if KEYWORD_SET(check) then begin
		if isfuv then plot_outbox = outbox / 2 else plot_outbox = outbox
		TVLCT, rr, gg, bb, /get
		PB_SET_LINE_COLOR
		PLOTS, /dev, plot_outbox[[0,0,2,2,0]], plot_outbox[[1,3,3,1,1]], color = 2
		TVLCT, rr, gg, bb
	endif
	if isnuv and isold then outbox = outbox + 2072 * [1,0,1,0]
	for j = 0, 3 do $
		CHANGE_TAG_VALUE, outhead, outbox[j] + 1, tagpart + fits_tag_parts[j] + STRTRIM(i,2)
endfor

if KEYWORD_SET(update) then begin
	if isfuv then begin
		recnumstr = STRTRIM(params.fuvs.rec_num, 2) + ' (FUVS), ' + STRTRIM(params.fuvl.rec_num, 2) + ' (FUVL)'
	endif else begin
		recnumstr = STRTRIM(params.rec_num, 2)
	endelse
	UPDATE_HISTORY, outhead, 'ran with rec_num ' + recnumstr
	UPDATE_HISTORY, outhead, 'boxwarp set to ' + STRTRIM(boxwarp, 2)
	IRIS_ISP2WCS, outhead, /level1point5
	UPDATE_HISTORY, outhead, 'updated WCS parameters with iris_isp2wcs'
endif

RETURN, outhead

end