;+

; $Id: hi_index_peakpix.pro,v 1.1 2008/01/21 16:15:54 bewsher Exp $

;

; Project: STEREO-SECCHI

;

; Name: hi_index_peakpix

;

; Purpose: To find pixels within a 2D data array whose absolute value

;          are above a given threshold, or are next to such a pixel

;

; Category: STEREO, SECCHI, HI

;

; Explanation: This routine is intended to be called from

;              hi_remove_starfield, instead of an HI image, it passes

;              the del^2 data of an image.

;

; Syntax: index=hi_index_peakpix(data)

;

; Example: index=hi_index_peakpix(data,thresh=10)

;

; Inputs: data - 2d image 

;

; Opt. Inputs: thresh - threshold for pixel value, defaults is 1

;

; Outputs: index - 2xN output of (i,j) pixel locations

;                  If no pixels are found, array containing

;                  [-999,-999] is returned

;

; Opt. Outputs: None

;

; Keywords: None

;

; Calls: None

;

; Common: None

;

; Restrictions: None 

;

; Side effects: None

;

; Prev. Hist.: None

;

; History: Written 27 February 2007 by Daniel Brown (DOB)

;          23-Mar-07, documentation added (DOB)

;

; Contact: Daniel Brown (dob@aber.ac.uk)

;

; $Log: hi_index_peakpix.pro,v $
; Revision 1.1  2008/01/21 16:15:54  bewsher
; First release of hi_remove_starfield and associated code
;

;-

function hi_index_peakpix,ddat,thresh=thresh



; set default value for threshold

if (n_elements(thresh) eq 0) then thresh=1.0



; find data array size

sz=size(ddat)

nx=sz(1)

ny=sz(2)



; set a counter for number of pixels found, use long integers

npix=0l



; take initial sweep through data to find peak pixels above threshold

; and their immediate neighbours. This first sweep is just to make a

; count so a suitably sized array can be constructed

for j=1,ny-2 do begin

    for i=1,nx-2 do begin



        if (abs(ddat(i,j)) gt thresh) then begin

            ; data above threshold

            npix=npix+1l

        endif else begin

            ; test to find neighbouring pixels

            tst1 = (abs(ddat(i+1,j)) gt thresh)

            tst2 = (abs(ddat(i-1,j)) gt thresh)

            tst3 = (abs(ddat(i,j+1)) gt thresh)

            tst4 = (abs(ddat(i,j-1)) gt thresh)

            if (tst1 or tst2 or tst3 or tst4)  then begin

                npix=npix+1l

            endif

        endelse



    endfor

endfor



; were any pixels found?

if (npix gt 0l) then begin

    ; create index array

    ind=intarr(2,npix)

    n=0l



    ; add pixel locations to index array

    for j=1,ny-2 do begin

        for i=1,nx-2 do begin



            if (abs(ddat(i,j)) gt thresh) then begin

                ind(*,n)=[i,j]

                n=n+1l

            endif else begin

                tst1 = (abs(ddat(i+1,j)) gt thresh)

                tst2 = (abs(ddat(i-1,j)) gt thresh)

                tst3 = (abs(ddat(i,j+1)) gt thresh)

                tst4 = (abs(ddat(i,j-1)) gt thresh)

                if (tst1 or tst2 or tst3 or tst4)  then begin

                    ind(*,n)=[i,j]

                    n=n+1l

                endif

            endelse



        endfor

    endfor

endif else begin

    ; a default value if there are no pixels

    ind=[-999,-999]

endelse



; return index

return,ind

end