; Copyright (c) 1998-2014  A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_IMG_AVG
;
;LAST CHANGED: ----------------------------------- 29-Oct-14
;
; PURPOSE:
;	This procedure computes the mean and standard deviation
;   of all pixels in a user-selected region.
;
; CATEGORY:
;	AXIS: image analysis
;
; CALLING SEQUENCE:
;	AX_IMG_AVG [, /NOZEROS]
;
;CALLED FROM AXIS:
;	->Images->Average by region->{all pixels, no zeros}

; INPUTS:
;	All input parameters are passed as keywords or through AXIS_COM
;
; KEYWORDS:
;	NOZEROS - if set, report average of only non-zero pizels in selected region
;	WHOLE	- if set, average all pixels in the image
;
; OUTPUTS:
;	Results are printed to log window and the AXIS com window.
;
; COMMON BLOCKS:
;	AXIS_COM	standard set of common blocks
;
; PROCEDURE:
;	Uses DefROI widget for felxible region of interest definition.
;
; EXAMPLE:
;	only executable from AXIS
;
; MODIFICATION HISTORY:
; (31-dec-99 aph) isolated from AXIS_C and merged with /nozeros (corected !)
;                 AXIS standard documentation
; (03-feb-01 aph) add keyword WHOLE to report average of full image
; (14-sep-04 aph) extend format field for axis log; extend to whole without zeros
;                  extend variable size to accomodate large images; append to log
; (08-mar-05 aph) incrase precision
; (24-mar-05 aph) correct error for evlauating images bigger than pixel mapped display
; (30-nov-11 aph) improve speed & add an hourglass symbol
; (29-Oct-14 aph) fix error for region, ignore zeros
;-

PRO ax_img_avg, nozeros = nozeros , whole = whole
@axis_com
on_error,2

HANDLE_VALUE, Data(CurBuf), tmp
IF n_tags(tmp) EQ 0 THEN RETURN
IF tmp.t NE '2d' THEN RETURN

ax_wait
stand_dev = 0.

pixarea = abs((tmp.x(1)-tmp.x(0)) * (tmp.y(1)-tmp.y(0)))

IF keyword_set(whole) then begin
	if keyword_Set(nozeros) then begin	; ------------KEEP ONLY non-zero values
		nz = where(tmp.d NE 0, n_pix)
		print, 'number of non-zero pixels = ', n_pix
		print, 'number of zero pixels = ', n_elements(tmp.d)-n_pix
		print, 'total number of pixels = ', n_elements(tmp.d)
		dat_tmp = reform(tmp.d(nz), n_pix)
	endif else begin
		dat_tmp = tmp.d
		n_pix = n_elements(tmp.d)
	endelse

	t = moment(dat_tmp,/nan, sdev=stand_dev,/double)
	average = t(0)


ENDIF ELSE BEGIN
	ncol = n_elements(tmp.x) & nrow = n_elements(tmp.y)
	img_size = max([ncol, nrow])
	img_zoom =  floor(360*gr_scale/img_size)
	if img_zoom EQ 0 then img_zoom = 1 ;for images larger than 360*gr_scale
	if img_size GT 360*gr_scale then print, 'Image larger than display area; rebinned prior to display'
	SetGraf, 'MainImg'
	erase
	loadct,ColTbl,/silent
	tvscl, rebin(tmp.d, img_zoom*ncol,img_zoom*nrow)
	roi = cw_defroi(MainWid, image_size = [ncol,nrow], zoom =[img_zoom, img_zoom])
;     	print,'Zoom = ', img_zoom, '.   ROI limits: ',min(roi),max(roi)
	svec = size(roi)
	IF (svec(0) NE 0) THEN BEGIN	; proceed only if roi is real region
;		if keyword_Set(nozeros) then begin	; KEEP ONLY non-zero values
;			nz = where(tmp.d NE 0, n_pix)			; THIS DuPLICATES THE ABOVE
;			if n_pix NE 0 then dat_tmp = reform(tmp.d(nz), n_pix)
;		endif else begin
;			dat_tmp = tmp.d(roi)
;			n_pix = n_elements(roi)
;		endelse

		dat_tmp = tmp.d(roi)
		n_pix = n_elements(roi)
		if keyword_Set(nozeros) then begin	; KEEP ONLY non-zero values
			nz = where(dat_tmp NE 0, n_pix)
			print, 'number of non-zero pixels = ', n_pix
			print, 'number of zero pixels = ', n_elements(dat_tmp)-n_pix
			print, 'total number of pixels = ', n_elements(dat_tmp)
			if n_pix GT 0 then dat_tmp = dat_tmp(nz)
		endif
		if n_elements(dat_tmp) GT 1 then begin
			dat_tmp = dat_tmp(1:n_elements(dat_tmp)-1)
			t = moment(dat_tmp,/nan, sdev=stand_dev,/double)
			average = t(0)
		endif else average = dat_tmp
     	ax_color
   		PlotBuf, CurBuf
	ENDIF else begin
		axis_log, 'No region selected'
 		return
	endelse
ENDELSE

; report average value of pixels evaluated (all or non-zero)
if keyword_set(whole) then text = 'Whole image ' else text = 'Region '
if keyword_Set(nozeros) then text = text + '- only non-zero pixels' else text = text + '- all pixels'
axis_log, text
axis_log, string(n_pix,format='("Area: ",i8," pixels")')
axis_log, string(n_elements(dat_tmp)*pixarea, format='("       ", g11.2," um^2")')
axis_log, string(average, stand_dev, format='("Mean= ",g11.3," +/- ",g11.3)')

ax_wait, /off

END
