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

    @compile_opt.pro        ; On error, return to caller

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