function particleIdentify, inpIm, thresh, clusSize=clusSize, britecut=britecut, $ nearneigh=nearneigh, info=info, relative=relative ;+ ;NAME: PARTICLEIDENTIFY ;SYNTAX: RESULT=PARTICLEIDENTIFY(INPIM, THRESH [, CLUSSIZE=CLUSSIZE,BRITECUT = ; BRITECUT,/NEARNEIGH, INFO=INFO, /RELATIVE] ) ;PURPOSE: Given an input image and a threshold, returns a corresponding featured ; image with elements above a threshold & in a "cluster" set to the corresponding ; cluster value, and all others set to 0. ;MANDATORY INPUTS: ; INPIM: The input image to be thresholded. D-dimensional numerical array. ; THRESH: The value to threshold at. Float/double scalar. The image is ; thresholded based on the pix >= thresh. If /RELATIVE is set, then ; THRESH must be on [0,1], and is the image is thresholded based on ; the pix >= thresh*(val-minVal)/(maxVal-minVal). ;OPTIONAL INPUTS: ; CLUSSIZE: 2-element int array, containing the [min, max] number of elements ; in a cluster for it to be considered. Defaults to [600,3000] ; BRITECUT: 2-element float array, containing the [min,max] average brightness ; (i.e. integrated brightness divided by # of voxels) of a cluster for it ; to be considered. Defaults to [0,!values.f_infinity] -- i.e. defaults to ; not mattering. ; /NEARNEIGH: If set to 1, only considers the 4 nearest neighbors in a connected ; cluster analysis. ; /RELAVITE: If set, the threhold is assumed to be on [0,1] and is used as a relative ; value rather than an absolute. ;OUTPUTS: ; RESULT: Integer D-dimensional array of the same size as inpIm, with #s for ; "pixels" in a cluster above threshold and 0s for pixels below threshold. ;OPTIONAL OUTPUTS: ; INFO: 2+2*N-element long array containing [# of points above threshold in image, ; # of original clusters from CCA, *], with info[2+2*a:2+2*a+1] = [mass, # of ; pixels in cluster] of the (a+1)^th returned cluster/particle. ;PROCEDURE: ; 1) Thresholds image, 2) applies CCA (using LABEL_REGION), 3) eliminates clusters ; outside the desired size & mass range, renames remaining clusters. ;MODIFICATION HISTORY: ; -Modified from thresholdimage.pro (written by B. Leahy), 12/01/2011, B. Leahy. ; -Changed default operation to an absolute threshold rather than relative; added ; /RELATIVE as an option, B. Leahy 12/12/2011. ; -Fixed bug with reporting the # of particles above threshold (info), 01/20/2012. ; -Changed massCut to briteCut (from total mass to average mass), 05/03/12, B. Leahy. ; ;- on_error,1 if not (keyword_set(inpIm) and keyword_set(thresh)) then message, $ 'inpIm and thresh must be set!' if not keyword_set(clusSize) then clusSize=[600,3000] if not keyword_set(britecut) then britecut=[0,!values.f_infinity] alln=1-keyword_set(nearneigh) ;for nearest neighbors or no if keyword_set(relative) then begin if (thresh gt 1) or (thresh lt 0) then message, $ 'THRESH must be on [0,1] if /RELATIVE is set!' maxVal=max(inpIm) minVal=min(inpIm) thresh=thresh*(maxval-minval)+minval ;rescaling if /relative is set endif thresIm=reform(inpIm) ge thresh dum=where(thresIm ge 1, numAbove) ;numAbove is the total # of points above threshold ;embedding thresIm in a larger region, since label_region ignores pixels at the edges imInf=size(thresIm) thresIm0=bytarr(imInf[1:imInf[0]]+2) case imInf[0] of 2: thresIm0[1:-2,1:-2]=thresIm 3: thresIm0[1:-2,1:-2,1:-2]=thresIm endcase thresIm=label_region(temporary(thresIm0), all_neighbors=alln) case imInf[0] of ;re-selecting the image of the appropriate size 2: thresIm=thresIm[1:-2,1:-2] 3: thresIm=thresIm[1:-2,1:-2,1:-2] endcase numClusters=max(thresIm);the # of clusters numTossed = 0 ;the number of clusters tossed out b/c they're not in cutoffs for a=1,numClusters do begin w=where(thresIm eq a, count) wrongSize = ( (count lt clusSize[0]) or (count gt clusSize[1])) ;it's outside the # pixel range brite=total(inpIm[w])/float(count) ;the AVERAGE brightness of the cluster wrongBrite = ( (brite lt briteCut[0]) or (brite gt briteCut[1])) ;it's outside the mass range if ( wrongSize or wrongBrite ) then begin ;cluster's got the wrong properties thresIm[w]=0 numTossed++ endif else begin ;cluster looks fine. thresIm[w] = a - numTossed ;renaming the cluster endelse endfor numPoints=total(thresIm gt 0) info=fltarr(2+2*max(thresIm)) ;initializing info info[0]=numAbove ;the # of pixels above threshold originally = # of pixels in all clusters info[1]=numClusters ;the # of original clusters of pixels for a=1,max(thresIm) do begin ;over cluster labels w=where(thresIm eq a,count) ;picking out the cluster labeled 'a' brite=total(inpIm[w])/float(count) ;average brightness of the cluster info[2*a] = brite info[2*a+1]=count;# of pixels in cluster endfor;a return, thresIm end