PRO dical_interp, img, loc, VERT=vert, MAXSIZE=maxsize, FLAGS=flags

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Performs thin plate spline interpolation across segments of the given
;;	image. At the edges, pixels are repeated, instead of extrapolated.
;;
;; CALLING SEQUENCE:
;;	DICal_interp, img, loc, /VERT
;;
;; REQUIRED INPUTS:
;;	img - The image to perform the interpolation on
;;	loc - result of a WHERE call specifiying where to interpolate over
;;	
;; OUTPUTS:
;;	The image with the areas interpolated over
;;
;; OPTIONAL INPUT KEYWORDS:
;;	VERT - only interpolate in the vertical direction. Used for IR images.
;;	MAXSIZE - the maximum size of a region to interpolate over in 2D interpolation
;;	FLAGS - flags identifying where the invalid pixels are
;;
;; EXAMPLE:
;;      IDL> dical_interp, in, loc
;;
;; PROCEDURES USED (i.e. called directly!):
;;    find_2Dblobs
;;
;; MODIFICATION HISTORY:
;;   2004-08-19  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

@dical_flags

img = double(img)
dim = size(img, /dimensions)

IF NOT keyword_set(maxSize) THEN maxSize = 9999

IF size(loc, /n_dimension) EQ 0 THEN RETURN ;; Nothing to interpolate over

;; Determine where the invalid data resides
badData = where(flags NE 0)

IF keyword_set(VERT) THEN BEGIN ;; 1D interpolation along the columns
	img = dical_irinterp(img,loc,badData)

	;; Modify the flags
	flags[loc] = flags[loc] OR FLAG_INTERP	
ENDIF ELSE BEGIN ;; 2D Interpolation

	;; Group the gaps by size
	gaps = find_2DBlobs(loc, dim[0], SIZE=gapSize)

	;; Remove the gaps that are too big
	justRight = where((gapSize[0,*]<gapSize[1,*]) LE maxSize, cnt)
	IF cnt EQ 0 THEN RETURN ;; Nothing to interpolate over

	toInterp = gaps[*,justRight]
	loc = toInterp[where(toInterp GE 0)]
	loc = reform(loc, n_elements(loc))
	
	;; Figure out where the edges are
	X = loc mod dim[0]
	Y = loc / dim[0]
	edge = where((X EQ 0) OR (Y EQ 0) OR (X EQ dim[0]-1) OR (Y EQ dim[1]-1), cnt,$
		complement=w, ncomplement=compCnt)
	IF cnt GT 0 THEN BEGIN
		;; Fill in the edges, for now, just use the ir routine because we
		;; there will only be holes that can be filled using the vertical and
		;; development time is limited. THIS MAY BE REVISITED
		img = dical_irinterp(img,loc[edge],badData)

		;; Update the flag map
		flags[loc[edge]] = flags[loc[edge]] OR FLAG_INTERP
		
		;; Update the pixels that need interpolation
		IF compCnt EQ 0 THEN RETURN
		loc = loc[w]
	ENDIF
							

	;; Interpolate over the gaps
	IF size(toInterp, /N_DIMENSIONS) GT 1 THEN $
		gapCnt = (size(toInterp, /DIMENSIONS))[1] $
	ELSE $
		gapCnt = 1

	FOR i=0, gapCnt-1 DO BEGIN
		gloTodo = toInterp[*,i]	;; Pixels to interpolate in global coordinates
		gloTodo = gloTodo[where(gloTodo GE 0)]
		gloTodo = SetIntersection(gloTodo, loc)

		;; Make sure the gap still exists
		IF (size(gloTodo,/n_dimensions))[0] EQ 0 THEN CONTINUE
	
		;; Move the gap so that it's top corner is at (1,1)
		todo = [[gloTodo MOD dim[0]], [gloTodo/dim[0]]]
		corner = [min(todo[*,0]), min(todo[*,1])]-2
		dataSize = gapSize[*,i]+4
	
		;; Adjust the corner in case we're on an edge
		IF corner[0] LT 0 THEN BEGIN
			corner[0] = 0
			dataSize[0] = dataSize[0]-1
		ENDIF
		IF corner[1] LT 0 THEN BEGIN
			corner[1] = 0
			dataSize[1] = dataSize[1]-1
		ENDIF
		IF corner[0]+dataSize[0] GT dim[0] THEN $
			dataSize[0] = dim[0] - corner[0]
		IF corner[1]+dataSize[1] GT dim[1] THEN $
			dataSize[1] = dim[1] - corner[1]
	
		todo[*,0] = todo[*,0]-corner[0]
		todo[*,1] = todo[*,1]-corner[1]
		todo = (dataSize[0])*todo[*,1]+todo[*,0]
	
		;; Get the points around the gap that have valid values
		known = setdifference(indgen(dataSize[0]*dataSize[1],/LONG),todo)
		window = img[corner[0]:corner[0]+dataSize[0]-1, $
			corner[1]:corner[1]+dataSize[1]-1]
		flagWnd = flags[corner[0]:corner[0]+dataSize[0]-1, $
			corner[1]:corner[1]+dataSize[1]-1]
		flagWnd = flagWnd[known]
		known = known[where((flagWnd AND (FLAG_GAP OR FLAG_BAD)) EQ 0)]
		X = (known MOD (dataSize[0]))
		Y = (known / (dataSize[0]))
	
		;; Do the actual interpolation
		IF max(gapSize[*,i]) EQ 1 THEN BEGIN
		;; Do a weighted average of the pixels. This is done to save 
		;; time because 1 pixel holes are probably just bad pixels and
		;; there could be a lot of them
			sq2 = sqrt(2)
			diag = ((X+Y) MOD 2) EQ 0 ; 1 if the element is on a diagonal
			tot = total(diag)/sq2 - total(diag-1)
			weight = diag/sq2 - (diag-1)
			tmp = total(img[X+corner[0],Y+corner[1]]*weight, /nan)$
				/tot
		ENDIF ELSE $
			tmp= grid_tps(X,Y,img[X+corner[0],Y+corner[1]],$
				NGRID = dataSize, START=0, DELTA=1)
		img[gloTodo] = tmp[todo]
		flags[gloTodo] = flags[gloTodo] OR FLAG_INTERP
	ENDFOR
ENDELSE

RETURN

END