alllocmax.pro

function alllocmax, cubein, indcube = indcube, friends = friends, $
                    specfriends = specfriends, round = round

;+
;
; NAME:
;   ALLLOCMAX()
; PURPOSE:
;   To establish a candidate set of local maxima within a data cube by
;   searching over a F by F by S box of pixels for a point
;   that's larger than all others in the box.  F and S default to 1
;
; CALLING SEQUENCE:
;   local_maxima = ALLLOCMAX(cube, [friends = friends, specfriends =
;   specfriends, indcube = indcube])
;
; INPUTS:
;   CUBE -- data cube to find local maxima in.
;             
; KEYWORD PARAMETERS:
;   INDCUBE -- (optional) if supplied, then the indices returned are
; drawn from this cube instead of indexing the supplied cube itself.
;   FRIENDS -- Sets the search region in the position dimensions to be
;              2*FRIENDS+1 in size.  
;   SPECFRIENDS -- Sets the search region in the velocity dimensions to be
;              2*SPECFRIENDS+1 in size.  
;
; OUTPUTS:
;   LOCAL_MAXIMA -- indices in CUBE (or in INDCUBE if set) that are
;                   local maxima.
;
; MODIFICATION HISTORY:
;
;       Tue Apr 16 15:09:40 2013, erosolo <erosolo@>
;
;		Enabled specfriends = 0 for 3D case.
;
;       Wed Feb 13 20:44:16 2013, erosolo <erosolo@> 
;           Added in the /round keyword to enable circular patches.
;
;       Documented -- Fri Sep 2 15:48:17 2005, Erik Rosolowsky
;                     <erosolow@asgard.cfa.harvard.edu>
;-

; GET THE SIZE AND MAKE AN EMPTY CUBE TO FLAG WHERE LOCAL MAXIMA
; EXIST.
  sz = size(cubein)
  if sz[0] eq 2 then begin
     specfriends = 0
     sz[3] = 1
  endif
; INITIALIZE EVERYTHING TO BE A LOCAL MAX, WE WILL REJECT THEM OVER TIME
  lmaxcube = bytarr(sz[1], sz[2], sz[3]) + 1B

; INITIALIZE THE DEFAULT BOX SIZE TO BE +/- ONE PIXEL (3 x 3 BOX)
  if (n_elements(friends) eq 0) then $
     friends = 1

; SET THE NONSENSICAL INDICES TO NOT-A-NUMBERS
  badind = where(cubein ne cubein, badct)
  cube = cubein
  if (badct gt 0) then begin
     cube[badind] = -!values.f_infinity
     lmaxcube[badind] = 0B
  endif

  if n_elements(specfriends) eq 0 then specfriends = 1 
  
; Round up to integers in case we have fractional pixel sizes make it this far.  Thanks a lot, user.
  specfriends = ceil(specfriends)
  friends = ceil(friends)

; A LOCAL MAXIMA IS A POINT WHICH IS GREATER THAN ALL OF THE POINTS
; AROUND IT.


  if keyword_set(round) then begin
     patch = shift(dist(2*friends+1,2*friends+1),friends,friends) le friends
  endif else patch = bytarr(2*friends+1,2*friends+1)+1b


  if specfriends gt -1 then begin
;;     for k = -specfriends, specfriends do $
;;       for j = -friends, friends do $
;;         for i = -friends, friends do $
;;           if NOT ((i eq 0) AND (j eq 0) AND (k eq 0)) then $
;;             lmaxcube = $
;;       lmaxcube*((cube gt shift(cube, i, j, k)) OR $
;;                 (finite(shift(cube, i, j, k)) eq 0)) 
     for k = -specfriends, specfriends do begin
        for j = -friends, friends do begin
           for i = -friends, friends do begin
              if (NOT ((i eq 0) AND (j eq 0) AND (k eq 0))) AND patch[i+friends,j+friends] then $
                 lmaxcube = $
                 lmaxcube*((cube gt shift(cube, i, j, k)) OR $
                           (finite(shift(cube, i, j, k)) eq 0))
           endfor
        endfor
     endfor
  endif else begin
     for j = -friends, friends do $
        for i = -friends, friends do $
           if (NOT ((i eq 0) AND (j eq 0))) AND patch[i+friends,j+friends] then $
              lmaxcube = $
        lmaxcube*((cube gt shift(cube, i, j)) OR $
                  (finite(shift(cube, i, j)) eq 0))
  endelse

  
  lmaxind = where(lmaxcube eq 1B, num)
  if (num eq 0) then begin
     message, 'No true local max found, defaulting to high point in data.', /con
     dummy = max(lmaxcube, lmaxind, /nan)
  endif

; IF THE INDEX CUBE IS SUPPLIED AND THERE ARE LOCAL MAXIMA THEN
; SUBSTITUTE THE INDICES FROM THE CUBE FOR THE ACTUAL INDICES
  if ((n_elements(indcube) gt 0)) then begin
     lmaxind = indcube[lmaxind]
  endif 

  return, lmaxind
end                             ; OF ALLLOCMAX