;+
; NAME:
;		Feature	
;
; PURPOSE:
;		Finds and measures roughly circular 'features' within 
;		an image.
;
; CATEGORY:
;		Image Processing
;
; CALLING SEQUENCE:
;		f = feature( image, diameter, separation, min=min, masscut=masscut,
;			iterate=iterate, maxits=maxits, pickn=pickn, /field, /quiet)
;
; INPUTS:
;		image:	(nx,ny) array which presumably contains some
;			features worth finding
;		diameter: a parameter which should be a little greater than
;			the diameter of the largest features in the image.
;			Diameter MUST BE ODD valued.
;
; OPTIONAL INPUTS:
;		separation: an optional parameter which specifies the 
;			minimum allowable separation between feature 
;			centers. The default value is diameter+1.
;
; KEYWORD PARAMETERS:
;		min: 	Minimum pixel intensity considered to be a
;                       candidate local maximum (default: 64th
;                       percentile).
;                       Setting min appropriately can lead to much
;                       faster operation.
;		masscut: Discard candidate features with integrated
;                       brightness below this value.
;                       Setting masscut appropriately can lead to much
;                       faster operation.
;               pickn:  Return only the PICKN brightest candidate
;                       features.
;		field: 	Set this keyword if image is actually just one field
;			of an interlaced (e.g. video) image. All the masks
;			will then be constructed with a 2:1 aspect ratio.
;		iterate: if the refined centroid position is too far from
;			the initial estimate, iteratively recalc. the centroid
;			using the last cetroid to position the mask.  This 
;			can be useful for really noisy data, or data with
;			flat (e.g. saturated) peaks.  Use with caution- it
;			may 'climb' hills and give you multiple hits.
;               maxits: Maximum number of iterations (default = 10).
;		quiet:	Supress printing of informational messages.
;
; OUTPUTS:
;		f[0,*]:	x centroid positions, in pixels.
;		f[1,*[: y centroid positions, in pixels. 
;		f[2,*]: integrated brightness of the features.
;		f[3,*]: square of the radius of gyration
;
; SIDE EFFECTS:
;		Optionally prints diagnostic messages to the terminal.
;
; RESTRICTIONS:
;    This program finds the centers of brightness of bright, circularly
;    symmetric, non-overlapping features on a dark background.  To
;    identify dark features on a bright background, first invert the
;    image.  
;
;    Performance is degraded by non-zero background values,
;    particularly if the background varies across the image.  
;    High-frequency noise also degrades performance.  Both the
;    background and noise can be suppressed by
;    filtering the image with BPASS before running FEATHRE.
;
; PROCEDURE:
;    First, identify the positions of all the local maxima in the
;    image.  Local maxima are the brightest pixels in a circular
;    neighborhood whose diameter is set by DIAMETER.  Within each
;    circular region, the feature's centroid is calculated as
;    the brightness-weighted center of brightness.  The integrated
;    brightness ("mass") and the brightness-weighted radius
;    ("radius of gyration") also are calculated within each circle.
;    If the centroid is found to be more than 0.5 pixels from the
;    original local maximum, the mask can be moved and the centroid
;    recalculated.  This is controlled by the ITERATE keyword.
;    Finally, centroids closer than SEPARATION are merged. 
;
;    This procedure can yield centroid positions with errors of 0.1
;    pixel or better in each dimension for features larger than about
;    5 pixels across.  Achieving this accuracy requires meeting the
;    conditions described in RESTRICTIONS.  Choosing an inappropriate
;    value for DIAMETER or improper settings in BPASS can degrade
;    performance to single-pixel accuracy.
;
;    Finally, FEATURE can select features based on their peak and
;    integrated brightness.  The former is set with the MIN keyword,
;    and the latter with MASSCUT.  Setting these values appropriately
;    greatly increases FEATURE's speed and yields much better
;    rejection of spurious features.  The PICKN keyword also is useful
;    when the number of objects is known.
;
;    Setting DIAMETER, SEPARATION, MIN and MASSCUT is greatly
;    facilitated with the companion program FEATURETOOL.
;
;    The algorithms used in FEATURE and BPASS and a quantitative
;    analysis of their performance are described in
;
;    J.C. Crocker and D.G. Grier, J. Colloid Interface Sci. 179, 298 (1996).
;
; MODIFICATION HISTORY:
;		Original version: feature_stats2: David G. Grier, U of Chicago 1992.
;		Rewritten by John C. Crocker, U of Chicago, optimizing 
;			runtime and measurement error JCC 10/93.
;		Added field keyword JCC 4/94. 
;		Added eccentricity parameter JCC 5/95.
;		Added quiet keyword JCC 12/95.
;		Added iteration, fixed up the radius/diameter fiasco and
;		debugging to improve non-centroid data JCC 4/96.
;               Memory and run-time optimizations DGG 8/99.
;               Added PICKN DGG and E.R. Dufresne 8/99.
;               Fixed error for spheres very near edges DGG 3/2000.
;   11/12/2005. David G. Grier, New York University
;     Major overhaul ...
;     FEATURE now selects each subarray only once -- major speed-up
;     Removed LMX to eliminate redundant code -- another speed-up.
;     Eliminate for loop in RSQD.
;     Removed eccentricity, (and thus THETARR) -- rarely used.
;     Updated syntax for IDL 6.X
;   5/20/2006: DGG, improved fidelity to SEPARATION parameter.
;     Overhauled documentation.
;   6/13/2006: DGG: Better respect for QUIET.
;     Don't allow SEP < EXTENT.
;
;   Copyright (c) 2006 John C. Crocker and David G. Grier.
; 
; LICENSE
;
;    This program is free software; you can redistribute it and/or
;    modify it under the terms of the GNU General Public License as
;    published by the Free Software Foundation; either version 2 of the
;    License, or (at your option) any later version.
;
;    This program is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;    General Public License for more details.
;
;    You should have received a copy of the GNU General Public License
;    along with this program; if not, write to the Free Software
;    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;    02111-1307 USA
;
;    If the Internet and WWW are still functional when you are using
;    this, you shold be able to access the GPL here: 
;    http://www.gnu.org/copyleft/gpl.html
;-

;;;
;
;	RSQD: produce a parabolic mask
;
function rsqd, w, h

if n_params() eq 1 then h = w
xc = float(w-1) / 2.
yc = float(h-1) / 2.
x2 = (findgen(w) - xc)^2
y2 = (findgen(h)- yc)^2

return,  x2 # replicate(1., h) + replicate(1., w) # y2
end

;;;
;
;	FRACSHIFT: barrel shifts a floating point array by a 
;               fractional pixel amount using bilinear interpolation.
;
function fracshift, im, shiftx, shifty

if n_elements( im ) eq 1 then return, im

ipx = floor( shiftx )           ; integer part of x shift
ipy = floor( shifty )           ; integer part of y shift
fpx = shiftx - ipx              ; fractional part of x shift
fpy = shifty - ipy              ; fractional part of y shift
if fpx lt 0 then begin
   fpx++ & ipx--                ; make sure fractional parts are positive
endif		
if fpy lt 0 then begin
    fpy++ & ipy--
endif

image = im                      ; preserve input data: use local copy

imagex  = shift( image, ipx+1, ipy   )
imagey  = shift( image, ipx  , ipy+1 )
imagexy = shift( image, ipx+1, ipy+1 )
image   = shift( image, ipx  , ipy   )

res   = ( (1. - fpx) * (1. - fpy) * image   ) + $
	( (     fpx) * (1. - fpy) * imagex  ) + $
	( (1. - fpx) * (     fpy) * imagey  ) + $
	( (     fpx) * (     fpy) * imagexy ) 
	
return,res	
end

;;;
;
;	FEATURE:
;
function feature, a, extent, sep,  $
                  min=min, masscut=masscut, pickn=pickn, $
                  field=field, quiet=quiet, $
                  iterate=iterate, maxits=maxits, $
                  lmax=lmax, count=count

; set flags
report = not keyword_set(quiet) ; report normal output as well as errors
field = keyword_set(field)      ; work on a field rather than a frame
iterate = keyword_set(iterate)  ; iterate to improve centroid estimates

; check and process inputs
sz = size( a )
nx = sz[1]                      ; width of image
ny = sz[2]                      ; height of image

extent = floor(extent)          ; diameter of a particle
if (extent mod 2) eq 0 and report then begin
    message,'Requires an odd extent.  Adding 1...',/inf
    extent++
endif

if (n_params() lt 3) then $
  sep = extent + 1              ; separation between particles

if sep le extent then $
  sep++

; derived parameters
radius = float(extent)/2.       ; estimated particle radius
range = floor(sep/2)            ; range over which to search for neighboring features
yrange = field ? range/2. : range
yscale = field ? 2. : 1.

; numerical options
if n_elements(maxits) ne 1 then maxits = 10 ; maximum number of iterations
if n_elements(masscut) ne 1 then masscut = 0 ; minimum acceptable integrated brightness

if not keyword_set( min ) then begin ; minimum acceptable pixel intensity
    h = histogram( a )
    goal = 0.64 * total( h )
    min = 0
    val = h[min]
    while val le goal do begin
        min++
        val += h[min]
    endwhile
    if report then message, "setting threshold to "+strcompress(min), /inf
endif

; find local maxima
mmask = rsqd(sep) lt range^2
if field then $
  mmask = fieldof(mmask, /odd)
b = byte(a)
c = dilate(b, mmask, /gray)
r = where( b eq c and b ge min, ngood )
if ngood lt 1 then begin
    message, "No local maxima were brighter than MIN", /inf
    return, -1
endif

; local maxima provide initial estimates for particle positions
x  = float( r mod nx )
y  = float( floor(r / nx) ) 

; some local maxima will be too close to the edge -- eliminate them
good = where( x ge range and x lt (nx-range) and y ge yrange and y lt (ny-yrange), lmax )
if lmax lt 1 then begin
    message, "All local maxima were too close to edge", /inf
    return, -1
endif
x = x[good]
y = y[good]
if report then message,  strcompress( lmax ) + ' local maxima found.', /inf

; corners of regions around each local maximum
xl = x - floor(extent/2) 
xh = xl + extent - 1
yl = field ? y - floor(extent/4) : y - floor(extent/2) 
yh = field ? yl + floor(extent/2) - 1 : yl + extent - 1

; set up some masks
rsq = rsqd( extent )
mask = rsq le radius^2
xmask = findgen(extent, extent) mod extent
xmask -= (extent-1)/2.                  ; center mask at center!
ymask = transpose( xmask )
xmask *= mask
ymask *= mask
rmask = rsq * mask + 1./6.

; extract fields of the masks for field-based analysis
if field then begin
    mask = fieldof(mask, /odd)
    xmask = fieldof(xmask, /odd)
    ymask = fieldof(ymask, /odd)
    rmask = fieldof(rmask, /odd)
endif

; the array of results
res = fltarr(4, lmax)

; calculate sub-pixel centroid and other characteristics for each feature
for i = 0L, lmax-1L do begin
    xi =  x[i] & yi = y[i]
    xli = xl[i] & xhi = xh[i]
    yli = yl[i] & yhi = yh[i]

    suba = a[xli:xhi,yli:yhi]   ; sub-image around candidate feature
    m = total( suba * mask )    ; integrated brightness

    if m lt masscut then continue ; too small: m = 0 and rg = 0 at loop's end
                                ; doing the test here eliminates
                                ; calculations for spurious features.

                                ; displacement of centroid from center of mask
    xc = total( suba * xmask ) / m
    yc = total( suba * ymask ) / (m * yscale)
        
    if iterate then begin       ; iterate to improve position estimate
        its = 0
        repeat begin
            shifted = 0
            if abs(xc) gt 0.6 then begin
                shifted = 1
                dx = round(xc)
                xli = (xli + dx) > 0
                xhi = (xhi + dx) <  (nx-1) >  (xli+1)
                xi = (xi + dx) > xli <  xhi
            endif
            if abs(yc) gt 0.6 then begin
                shifted = 1
                dy = round(yc)
                yli = (yli + dy) > 0
                yhi = (yhi + dy) <  (ny-1) >  (yli+1)
                yi = (yi + dy) > yli <  yhi
            endif
            if shifted then begin     
                suba = a[xli:xhi,yli:yhi]
                m = total( suba * mask )
                xc = total( suba * xmask ) / m
                yc = total( suba * ymask ) / (m * yscale)
            endif
            its++
        endrep until not shifted or (its eq maxits)
     endif


    rg = total( suba * rmask ) / m ; radius of gyration
    res[*, i] = [xi+xc, (yi+yc)*yscale, m, rg]
endfor

; keep only non-trivial features
w = where(res[2, *] gt 0,  ngood)
if ngood gt 0 then $
   res = res[*, w]

; some features might have converged to the same point
; eliminate duplicates
hash =  floor(res[0, *]/sep) + nx * floor(res[1, *]/sep)
ndx = uniq(hash, sort(hash))
count = n_elements(ndx)
res = res[*, ndx]
if report then message, strcompress(count)+" unique features above threshold", /inf

; select the brightest features
if keyword_set( pickn ) then begin
    order = sort( res[2, *] )   ; sort by integrated brightness
    if count lt pickn then $
      message, "PICKN: Ignored: Fewer than "+strtrim(pickn,2)+ $
      " features to choose from.",/inf $
    else begin
        good = order[count-pickn:*]    
        res = res[*, good]
        if report then message, strcompress(pickn)+" features selected", /inf
    endelse
endif

return, res
end