FUNCTION find_2DBlobs, loc, xdim, SIZE=size, N=n $
                     , Fail2DBlobs=fail2db

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Identifies groups of connected data in two dimensions. For example, 
;;	in an image, can identify blobs of adjacent pixels that are all
;;	under a given threshold. 
;;	Blobs are found by starting at known flagged locations (ie. in loc)
;;	and then expanding spacially using a graph search technique until
;;	all the points in the blob are found. Note that points are not
;;	considered adjacent if they border diagonally.
;;
;; CALLING SEQUENCE:
;;	RESULT = find_2dBlobs, loc, dim
;;
;; REQUIRED INPUTS:
;;	loc - vector of positive indices identifing data to block. This is 
;;		usually the result of a call to WHERE.
;;	xdim - the size of the x dimension in the 2D data
;;
;; OUTPUTS:
;;	RETURN - A MxN long array where each [*,i] matrix is a list of
;;		indices that exist in the blob. M is the maximum number
;;		of points in a blob. N is the number of blobs in the data.
;;		Any element in this array that does not have useful data will
;;		be filled with -1. There is no guarantee as to the order of
;;		the elements in the list
;;		
;;
;; OPTIONAL INPUT KEYWORDS:
;;	SIZE - A named variable to receive a 2xN array where each [*,i]
;;		contains [xSize,ySize], the maximum size of each blob
;;	N - named variable to receive the number of blobs found
;;
;; EXAMPLE:
;;      IDL> blobs = find_2DBlobs(loc, xdim, SIZE=sz, N=n)
;;
;; PROCEDURES USED (i.e. called directly!):
;;
;; MODIFICATION HISTORY:
;;   2004-06-03  M. Desnoyer    Created
;;   2006-04-10  M. Desnoyer    Speedup (this line added by B. Carcich)
;;   2006-08-30  B. Carcich     Bumped size of stack by 1
;;
;;-----------------------------------------------------------------------------

fail2db=1b    ;; assume failure
p = N_ELEMENTS(loc)

if p gt 10000L then return,-1b   ;; fail if loc array is too large
                                 ;; - 10k => out will be 100M elements

lc = loc ;A copy of loc that will be emptied
n=0
m=0
out = lonarr(p,p)-1
size = lonarr(2,p)-1

;; Find each blob
while size(lc,/N_DIMENSIONS) GT 0 DO BEGIN
	;; Reset for the next blob
	stack = lonarr(1L+(3*p>4),/NOZERO)   ;;; bumped size by +1 (B. Carcich)
	hd = 0	;Next free slot in the stack 
	curBlob = [-1]
	
	;; Insert a flagged location into the stack
	stack[hd] = lc[0]
	hd = hd+1

	;; Do DFS search
	while hd GT 0 DO BEGIN
		;; Pop the top element
		cur = stack[hd-1]
		hd = hd-1

		;; If this location has already been seen, next
		IF size(where(curBlob EQ cur),/N_DIMENSIONS) GT 0 THEN CONTINUE

		;; If this location is not flagged, next
		IF size(where(lc EQ cur),/N_DIMENSIONS) EQ 0 THEN CONTINUE

		;; Add this location to the blob
		curBlob = [curBlob,cur]

		;; Add all the adjacent locations to the stack but pay 
		;; attention to edges
		IF (cur MOD xdim) NE 0 THEN BEGIN
			stack[hd] = cur-1
			hd = hd+1
		ENDIF
		IF (cur+1 MOD xdim) NE 0 THEN BEGIN
			stack[hd] = cur+1
			hd = hd+1
		ENDIF
		IF cur GE xdim THEN BEGIN
			stack[hd] = cur-xdim
			hd = hd+1
		ENDIF
		stack[hd] = cur+xdim
		hd = hd+1
	ENDWHILE
	curBlob = curBlob[1:*]
	IF(n_elements(curBlob) EQ 1 AND n_elements(lc) GT 1) THEN BEGIN
		lc = lc[where(lc NE curBlob[0])]	;;Shortcut to speedup code
	ENDIF ELSE BEGIN
		lc = setdifference(lc, curBlob)
	ENDELSE
	blbCnt = N_ELEMENTS(curBlob)
	out[0:blbCnt-1,n] = curBlob
	
	;; Determine the size of the blob
	a = max(curBlob MOD xdim,min=b)
	size[0,n] = a-b+1
	a = max(curBlob / xdim,min=b)
	size[1,n] = a-b+1
	
	;; Update counters
	n = n+1
	m = m>blbCnt
ENDWHILE

;; Clean up the output arrays
size = size[*,0:(n-1)>0]
fail2db=0b
RETURN, out[0:(m-1)>0,0:(n-1)>0]
END