; Copyright (c) 1998-2000 A.P. Hitchcock  All rights reserved 
;+ 
;NAME: 
;	AX_PARTICLES 
; 
;LAST CHANGED: ----------------------------------- 19-dec-00 
; 
;PURPOSE: 
;	This FUNCTION derives particle sizes from masked images 
;  A histogram of the particles sizes is generated. 
; 
;CATEGORY: 
;	STAND ALONE: image processing 
; 
;CALLING SEQUENCE: 
;	Result = AX_PARTICLES( [ ]) 
; 
;CALLED FROM AXIS: 
;	->Images->Particle analyze 
; 
;ROUTINES 
;	add_part(i,j) - function that adds (i,j) pixel to list 
;	ax_p_adj(i,j) - function that checks 4-adjacent pixels to see if they are 1 
; 
;INPUTS: 
;	INPUT - POSITIONAL INPUTS only (not keywords) 
;	All input parameters are passed as keywords. 
; 
;KEYWORDS: 
;      SIMAGE =  image in AXIS structure format 
;                (if not included; user is prompted for an AXIS format image) 
;      FILE  =  name of file to read (AXIS binary (*.axb)) 
;      AXIS_ON = flag indicating AX_PARTICLES is being run from AXIS 
; 
;OUTPUT: 
;      array indicating size of all particles 
; 
;OPTIONAL OUTPUT 
;     when run stand alone, ; 
;          the image is plotted with particles indicated 
;          a histogram of the particles sizes is plotted 
; 
;COMMON BLOCKS: 
;	@AXIS_COM	standard set of common blocks 
;  ax_part, image, particles, nx, ny, p, np, p_now, npix 
;               - common for ax_particles only 
; 
;PROCEDURE: 
;	This function is meant to apply to MASKED data (use Generate_mask 
; to convert image to (1/0) map based on user-selected threshold. 
; AX_PARTICLES then loops through all pixels. When it finds a 1 
; it searches all 4 adjacent pixels (N,S,E,W), marking those alos 1, until 
; a local search no longer finds 1's. The number of pixels defines 
; the size of that particle 
; It returns an ARRAY p(np) listing the size (in pixels) of each particle. 
; The array for currently analysed particle is p_now(npix,2) 
; where npix is the number of pixels in particle np 
;       np is the number of particles 
;       2nd index of p_now: 0 = x, 1 = y index co-ordinate in image 
; 
;MODIFICATION HISTORY: 
; (04-oct-99 aph) first developed 
; (31-dec-99 aph) AXIS standard documentation 
; (16-jan-00 aph) add lower, upper bound and display 
; (09-dec-00 cgz) modified names of color indices (blue and green)
; (19-dec-00 cgz) substituted ax_white_color_index with ax_plot_axes_color_index
;- 
 
PRO add_part,i,j 
; add this pixel to current particle - assumes user has tested & not used yet 
 
COMMON ax_part, image, particles, nx, ny, p, np, p_now, npix 
 
 
; marks the (i.j) pixel as part of current particle 
 
; ------- update indices of THIS particle 
	p_now = [p_now,intarr(1,2)]  	; increment pixel track of np-th particle 
	p_now(npix,0) = i 
	p_now(npix,1) = j 
	npix = npix + 1 
; ---------- mark this pixel as identified in particles 
    particles(i,j) = 1 
END 
 
; ******************************************************** 
 
function ax_p_adj,i,j 
; -------------- last changed: 4-Oct-99 
; checks the four adjacent pixels which share edges and identifies 1's 
 
COMMON ax_part, image, particles, nx, ny, p, np, p_now, npix 
 
test = 0 
 
ii= i+1 & jj = j 
if ii NE nx then begin 
 if image(ii,jj) EQ 1 then begin 
   if particles(ii,jj) EQ 0 then begin 
      add_part,ii,jj 
      test = 1 
   endif 
 endif 
endif 
 
ii= i & jj = j+1 
if jj NE ny then begin 
  if image(ii,jj) EQ 1 then begin 
   if particles(ii,jj) EQ 0 then begin 
      add_part,ii,jj 
      test = 1 
   endif 
 endif 
endif 
 
ii= i-1 & jj = j 
if ii NE -1 then begin 
  if image(ii,jj) EQ 1 then begin 
   if particles(ii,jj) EQ 0 then begin 
      add_part,ii,jj 
      test = 1 
   endif 
 endif 
endif 
 
ii= i & jj = j-1 
if jj NE -1 then begin 
  if image(ii,jj) EQ 1 then begin 
   if particles(ii,jj) EQ 0 then begin 
      add_part,ii,jj 
      test = 1 
   endif 
 endif 
endif 
 
return, test 
end 
 
; ******************************************************** 
 
FUNCTION ax_particles, simage=simage, file=file, axis_on=axis_on 
 
COMMON ax_part, image, particles, nx, ny, p, np, p_now, npix 
 
@axis_com
 
; --------- read-in data ---------------- 
if NOT keyword_set(simage) then begin 
	if NOT keyword_set(file) then begin 
	   file = pickfile2(/read, filter='*.axb') 
	   a = axb_load(file=file) 
	endif else a = axb_load(file=file) 
endif 
if keyword_set(simage) then begin 
    a = simage 
    file = 'c:\'+ Label(CurBuf) + '.txt' 
endif 
 
if n_tags(a) EQ 0 then return, 0 
image=fix(a.d) 
image_min = min(image, max=Image_max) 
if image_min NE 0 OR image_max NE 1 then begin 
    text = 'ERROR (AX_PARTICLES): input data must be MASKED (0/1 only)' 
    if keyword_set(axis_on) then widget_control, uprompt, set_value = text 
    print, text 
    return, 0 
endif 
 
nx = n_elements(a.x) 
ny = n_elements(a.y) 
dx = a.x(1) - a.x(0) 
dy = a.y(1) - a.y(0) 
dr = (dx + dy)/2. 
if abs(dx-dy) GT 0.05*dx then begin 
   text = 'WARNING: non-square pixels: dx= ' + string(dx) + '  dy= ' + string(dy) 
    if keyword_set(axis_on) then widget_control, uprompt, Set_value = text 
    print, text 
endif 
; ----------- plot image to be analysed ----------------- 
ax_color 
if keyword_set(axis_on) then begin 
   wset, MainImg 
   Widget_control,/HourGlass 
endif else begin 
   window, 0, xsize=nx+40, ysize=ny+40 
   splot2d,a 
endelse 
 
; ----------- define parameters to make pixel search more efficent 
if keyword_set(axis_on) then begin 
	a_min_d = get_num(prompt = 'Analysis - discard below (um)', val = dr, group = Axis_ID) 
endif else a_min_d = get_num(prompt = 'Analysis - discard below (um)', val = dr) 
if keyword_set(axis_on) then begin 
	a_max_d = get_num(prompt = 'Analysis - discard above (# of pixels)', val = 2, group = Axis_ID) 
endif else a_max_d = get_num(prompt = 'Analysis - discard above (# of pixels)', val = 2) 
print, string(format='("Searching for particles from",f6.3," to ",f6.3," um")',a_min_d, a_max_d) 
a_min = round(3.14*a_min_d*a_min_d/(dr*dr*4)) 
a_max = round(3.14*a_max_d*a_max_d/(dr*dr*4)) 
print, string(format='("                                  ",i6," to ",i6," pixels")',a_min, a_max) 
if keyword_set(axis_on) then begin 
	d_min_d = get_num(prompt = 'Highlight - minimum size (# of pixels)', val = a_min_d, group = Axis_ID) 
endif else d_min_d = get_num(prompt = 'Highlight - minimum size (# of pixels)', val = a_min_d) 
if keyword_set(axis_on) then begin 
	d_max_d = get_num(prompt = 'Highlight - maximum size (# of pixels)', val = a_max_d, group = Axis_ID) 
endif else d_max_d = get_num(prompt = 'Highlight - maximum size (# of pixels)', val = a_max_d) 
print, string(format='("Highlight particles from",f6.3," to ",f6.3," um")',d_min_d, d_max_d) 
d_min = round(3.14*d_min_d*a_min_d/(dr*dr*4)) 
d_max = round(3.14*d_max_d*a_max_d/(dr*dr*4)) 
print, string(format='("                                 ",i6," to ",i6," pixels")',d_min, d_max) 
 
particles = intarr(nx,ny) 
; -------- particles -------------------------- 
p = 0 & np = 0 & npt = 0 
FOR i = 0, nx-1 DO BEGIN 
  FOR j = 0, ny-1 DO BEGIN 
	if image(i,j) EQ 1 then begin 
	  IF particles(i,j) EQ 0 THEN BEGIN 
; found part of a NEW particle - mark it 
		p_now = intarr(1,2) 
		npix = 0	; restart index of particles 
		add_part,i,j 
; and find all other pixels of this particle which share an edge 
		nochange = 0 
		ti = intarr(1)  &  tj = intarr(1) 
		while nochange EQ 0 AND npix LT a_max do begin 
			np_tmp = npix 
			np_low = 0 
            check = 1 
            n_loop = 1 
;            print, 'Particle ', np,' . loop ', n_loop, '. #-pixels ', npix 
			for tt = np_low,np_tmp-1 DO BEGIN 
				ti = p_now(tt,0)  & tj = p_now(tt,1) 
				test = ax_p_adj(ti(0),tj(0)) 
				if test EQ 1 then check = 0 
			endfor 
			n_loop = n_loop + 1 
			np_low = npix-1 
			if check EQ 1 then nochange = 1 
   		endwhile 
; one final check ----------- 
		if npix LT a_max then begin 
			check = 0 
			for tt = 0, npix-1 do begin 
				ti = p_now(tt,0)  & tj = p_now(tt,1) 
				test = ax_p_adj(ti(0),tj(0)) 
				if test EQ 1 then print, np, ' particle. More found !' 
			endfor 
		endif 
 
; -------- flag if particle is INSIDE display_select limits 
		disp_color = ax_color03_index	;ax_blue_color_index 
		if npix GE d_min and npix LE d_max then begin 
			for ij = 0, npix-1 do particles(p_now(ij,0),p_now(ij,1)) = 2 
			disp_color = ax_color02_index	;ax_green_color_index 
		endif 
 
; --------------- indicate this particle on the image --------- 
		if not keyword_set(axis_on) then begin 
		    plots, a.x(p_now(0,0)), a.y(p_now(0,1)) 
			for ij = 1, npix-1 do begin 
			    plots, a.x(p_now(ij,0)), a.y(p_now(ij,1)), $ 
			        /continue, color = disp_color 
			endfor 
			if npix GE d_min and npix LE d_max then begin 
				print, '# ', np,'   size ', npix, '   green' 
			endif else print, '# ', np,'   size ', npix, '   blue' 
		endif else begin 
			text = string(format='("particle ",i4,/,"#pixels = ",i5)',np,npix) 
			widget_control, uprompt, set_value = text 
		endelse 
 
; ---- OK, what do we do with this particle ? 
		if npix GE a_min AND npix LE a_max then begin	; skip if too small 
			p(np) = npix												; or too large 
			np = np + 1  ; increment particle number 
			p = [p,1]	 ; increment  particle array 
		endif 
       endif 
	  endif 
	ENDFOR 
ENDFOR 
 
; ----- remove sub-size particles ------------- 
if keyword_set(axis_on) then begin 
    low_d = get_num(Prompt='minimum size (um)',val=a_min_d,group=axis_id) 
endif else low_d = get_num(Prompt='minimum size (um)', val=a_min_d) 
low = round(3.14*low_d*low_d/(dr*dr*4)) 
p = p(where(p GE low))   ; remove all hits below user-defined size (assume it is noise) 
np = n_elements(p) 
 
; ---------- convert area in pixels to average diameter (assume circular) 
s = 2.*sqrt(p/3.14)*dr 
max_s = max(s) 
bin = 2.*max_s/np      ; default bin size for histogram 
 
; --------- set histogram bin size (um) 
if keyword_set(axis_on) then begin 
	bin = get_num(Prompt = ' size of interval for histogram (um)', val = bin, group = axis_id) 
endif else bin = get_num(Prompt = ' size of interval for histogram (um)', val = bin) 
 
; ------- save particles image in buffer 7 
if keyword_set(axis_on) then begin 
	tmp = a 
	tmp.d = particles 
	tmp.dl = 'Particles in' + tmp.dl 
	CurBuf = 7 
   	HANDLE_VALUE, Data(CurBuf), tmp, /set 
   	Label(CurBuf) = tmp.dl 
   	PlotBuf,CurBuf 
endif else begin 
endelse 
 
; ----------- report results ------------- 
t = ax_name(file) & file_name = t(1) 
print, 'Particles in ', file_name 
print, '------------' 
if keyword_set(axis_on) then begin 
    text = string(FORMAT='(i4," particles",/, "Size  Area (#-pixels)  d(um)",/,"Min",i5,10x,f6.2,/,"Max",i5,10x,f6.2)', $ 
        np, min(p), min(s),max(p), max_s) 
    widget_control, uprompt, set_value = text 
endif 
text = string(FORMAT='(i4," particles")',np) 
print, text 
print,"Size  Area (#-pixels)  d(um)" 
text = string(FORMAT='("Min",i5,14x,f6.2)', min(p), min(s)) 
print, text 
text = string(FORMAT='("Max",i5,14x,f6.2)', max(p), max_s) 
print, text 
 
; -------------- plot histogram of particle sizes ----------- 
    print, 'bin size distribution at ', bin, ' um' 
    n_bins = fix(max_s/bin) 
    if max_s/bin - n_bins GT 0 then n_bins=n_bins+1 
    x = findgen(n_bins)*bin 
    y = histogram(s,bin=bin,min=0) 
;    print, 'x:',n_elements(x),' y:',n_elements(y) 
if keyword_set(axis_on) then begin 
	s = {t:'1d', d:reform(y), x:x, xl: 'Particle size (um)',dn:reform(y),$ 
	        yl: 'Number', dl:'Histogram of ' + file_name} 
endif else begin 
    window,1,xsize=200,ysize=200 
    plot,x,y, color=ax_plot_axes_color_index	;ax_white_color_index 
endelse 
 
return, s 
 
END 
