function IRIS_DITHER_UNDUST, img1, img2, $
	offset = offset, pedestal = pedestal, unslit = unslit, sigma = sigma, $
	pstop = pstop, fuv = fuv, crop = crop, $
	xmin = xmin, ymin = ymin, nx = nx, ny = ny
	
;+
;
; Turns a pair of dithered images into a composite image with the worst of the
; dust spots removed.
;
; INPUTS: 
;	img1/img2	-	The two images. Should be offset by something like 10"
;
; KEYWORDS:
;	offset	-	OPTIONAL, set to a 2-element vector giving estimated offset
;				between the images in pixels X,Y
;
;	pedestal-	OPTIONAL, the amount of pedestal to subtract (try something like 90?)
;
;	/unslit	-	If set, then the routine makes an effort to remove a region around
;				the slit. (Can also set to a width in pixels for the slit tapering
;				region)
;
;	sigma	-	If set, then points where the difference between the images exceeds
;				the specified amount of error, then that pixel is treated as 
;				dust and removed (by default, it merely takes the maximum value
;				for each pixel)
;
;	/fuv	-	set if the image is an FUV slit-jaw image (i.e. you only want to
;				look at the right half of the image).
;
;	/crop	-	set to crop the image to only include the overlapping region (by
;				default, a full-sized image is returned)
;
; WRITTEN: P. Boerner 7/19/2013
;			PB 8/14/2013 Added xmin/nx keywords
;
;-

if N_ELEMENTS(offset) lt 2 then offset = [0,0]
if N_ELEMENTS(pedestal) eq 0 then pedestal = 0
if N_ELEMENTS(xmin) eq 0 then xmin = 100
if N_ELEMENTS(ymin) eq 0 then ymin = 100
if N_ELEMENTS(nx) eq 0 then nx = 400
if N_ELEMENTS(ny) eq 0 then ny = 400

imgs =  ([[[img1]], [[SHIFT(img2, offset)]]]>0 - pedestal)>0
if KEYWORD_SET(fuv) then begin
	imgs[*,*,0] = ROTATE(imgs[*,*,0], 5)
	imgs[*,*,1] = ROTATE(imgs[*,*,1], 5)
endif
imsize = SIZE(imgs)

subimg = imgs[(xmin>0):(xmin+1)>(xmin+nx)<(imsize[1]-1),$
	ymin>0:(ymin+1)>(ymin+ny)<(imsize[2]-1),*] 
corr_disp = TR_GET_DISP(subimg)

if KEYWORD_SET(pstop) then STOP

if KEYWORD_SET(unslit) then begin
	if unslit eq 1 then slitshift = ABS(FIX(corr_disp[0,1]/4)) else slitshift = unslit
	profile = TOTAL(imgs, 2)
	for i = 0, 1 do begin
		slitpoint = WHERE(DERIV(profile[50:-50,i]) eq MAX(DERIV(profile[50:-50,i])))
		for j = 0-slitshift, slitshift do imgs[50+slitpoint+j,*,i] = imgs[50+slitpoint+j,*,i] - 200 + 200*ABS(j/slitshift)^(1/2.)
	endfor
endif

imgs[*,*,1] = SHIFT(imgs[*,*,1], corr_disp[*,1])
maximg = MAX(imgs, dim = 3)
result = maximg

if KEYWORD_SET(sigma) then begin
	noiseimg = SQRT(maximg)
	diffimg = imgs[*,*,1] - imgs[*,*,0]
	badpts = WHERE(ABS(diffimg) gt sigma*noiseimg, numbad)
	maskimg = maximg * 0
	maskimg[badpts] = 1
	meanimg = MEAN(imgs, dim = 3)
	meanimg[badpts] = maximg[badpts]
	result = meanimg	
endif

if KEYWORD_SET(crop) then begin
	corr_disp = ABS(corr_disp + REBIN(offset,2,2))
	tmpresult = result[corr_disp[0,1] : 1036 - corr_disp[0,1],*]
	result = tmpresult[*,corr_disp[1,1] : 1096 - corr_disp[1,1],*]	
endif

if KEYWORD_SET(fuv) then begin
	result = ROTATE(result, 5)
endif

if KEYWORD_SET(pstop) then STOP

RETURN, result

end