FUNCTION GroupPixels, Locations, Positions, sizeframe=sizeframe, dimension=dimension,	$
	ellipsoid=ellipsoid, range=Range, ngroup=nGroup, pGroup=pGroup, lgroup=lGroup

;+
; NAME:
;	GroupPixels
; PURPOSE:
;	Divide a group of pixels in separate groups based on nearness of the pixels in a group
; CATEGORY:
;	Tools
; CALLING SEQUENCE:
;	Loc = GroupPixels(Locations, Positions, range=range, ngroup=nGroup, pgroup=pGroup, lgroup=lGroup)
; INPUTS:
;	Locations		1D array, type integer
;						array of linear indices of pixel locations
;	Positions		3D array, type: float
;						x,y,z positions of pixels in 'Locations' array
;						(keywords sizeframe or dimension provide alternative
;						ways to specify x,y,z positions)
; OPTIONAL INPUT ARRAY:
;	sizeframe=sizeframe	size vector of image cube
;	dimension=dimension	dimensions of image cube
;						if either of these is specified, then 'Positions' is ignored
;						Instead the dimensions of the image cube are used to
;						to get x,y,z pixel coordinates from 'Locations' using
;						href=ArrayLocation=.		
; OPTIONAL INPUT PARAMETERS:
;	range=Range		scalar, or vector with as many elements as Frame has dimensions
;						default: Range = 1
;						defines the distance in each dimension of Frame defining
;						how close the members of a group are.
;	/ellipsoid		if set, Range is used to define an ellipsoidal region
;					if not set, Range defines a square region.
; OUTPUTS:
;	Loc				1D array, type integer
;						same as input array Locations but with elements reordered
;					if Frame is supplied as an array (i.e. if /sizeframe is not set),
;					then groups will be ordered with the values of the highest
;					pixel in each group in decreasing order
; OPTIONAL OUTPUT PARAMETERS:
;	ngroup=nGroup	scalar, type long integer
;						number of groups of pixels in Loc
;	pgroup=pGroup	1D array[nGroup], type long integer
;						indicates the position in Loc where each group begins
;	lgroup=lGroup	1D array[nGroup], type long integer
;						indicates number of pixels in each group
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsType, ArrayLocation
; PROCEDURE:
;	This procedure was written to group pixels in a 3D image cube, but I think it works
;	with cubes of any number of dimensions.
; MODIFICATION HISTORY:
;	NOV-1998, Paul Hick (UCSD/CASS)
;	FEB-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added possibility to specify x,y,z coordinates for pixels explicitly
;		(in addition to specification of size vector of image cube only)
;-

InitVar, ellipsoid, /key

; Get # dimensions in image cube and 3D pixel locations of pixels in image cube

Loc0 = Locations

IF IsType(dimension, /defined) THEN BEGIN					; Dimensions of image cube specified
	ndim = n_elements(dimension)							; # dimensions
	Pos0 = ArrayLocation(Loc0,dimension=dimension)			; Multi-dim pixel locations
ENDIF ELSE IF IsType(sizeframe, /defined) NE 0 THEN BEGIN	; Size vector of image cube specified
	ndim = Positions[0]										; # dimensions
	Pos0 = ArrayLocation(Loc0,size=sizeframe)				; Multi-dim pixel locations
ENDIF ELSE BEGIN
	ndim = (size(Positions,/dim))[0]						; # dimensions
	Pos0 = Positions										; Multi-dim pixel locations
ENDELSE

InitVar, Range, 1
InitVar, Range, count=1, replicate(Range,ndim)
Range = float(Range)

;Loc = 0*Loc0-1
;Pos = 0*Pos0-1
Loc = make_array(size=size(Loc0), value=-1)
Pos = make_array(size=size(Pos0), value=-1)

nLoc = 0
Loc[  nLoc] = Loc0[  0]
Pos[*,nLoc] = Pos0[*,0]
Loc0[0] = -1

iLoc = -1
pGroup = 0

i0 = where(Loc0 NE -1,n0)
WHILE i0[0] NE -1 DO BEGIN
	iLoc = iLoc+1
	qLoc = Pos[*,iLoc]
	CASE ellipsoid OF
	0: BEGIN
		One = replicate(1L,n0)
		i = where( total( abs(Pos0[*,i0]-qLoc#One) LE Range#One, 1) EQ ndim, n )
	END
	1: BEGIN
		i1 = where(Range NE 0,n1)
		i2 = where(Range EQ 0,n2)
		One = replicate(1L,n0)
		IF n1 EQ 0 THEN			$
			i = where( total( Pos0[*,i0]-qLoc#One EQ 0, 1) EQ ndim, n )	$
		ELSE IF n2 EQ 0 THEN	$
			i = where( total(((Pos0[*,i0]-qLoc#One)/Range#One)^2,1) LE 1, n )		$
		ELSE BEGIN
			p1 = Pos0[*,i0]
			i = where(	total(((p1[i1,*]-qLoc[i1]#One)/Range[i1]#One)^2,1) LE 1 AND	$
						total(  p1[i2,*]-qLoc[i2]#One EQ 0, 1) EQ ndim, n )
		ENDELSE
	END
	ENDCASE

	IF i[0] NE -1 THEN		$
		i = i0[i]			$
	ELSE IF iLoc EQ nLoc THEN BEGIN
		pGroup = [pGroup,nLoc+1]
		i = i0[0]
		n = 1
	ENDIF

	IF i[0] NE -1 THEN BEGIN
		Loc[  nLoc+1:nLoc+n] = Loc0[  i]
		Pos[*,nLoc+1:nLoc+n] = Pos0[*,i]
		nLoc = nLoc+n
		Loc0[i] = -1
	ENDIF
	i0 = where(Loc0 NE -1,n0)
ENDWHILE

nGroup = n_elements(pGroup)				; # groups
lGroup = lonarr(nGroup)					; # pixels in group
n = [pGroup,nLoc+1]
FOR i=0,nGroup-1 DO lGroup[i] = n[i+1]-n[i]

RETURN, Loc  &  END