; Copyright (c) 1998-2013 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	READ_NSRRC_PEEM.PRO
;
;LAST CHANGED: ----------------------------------- 27-Aug-13
;
;PURPOSE:
;	This reads an NSRCC PEEM stack from an NSRRC data file (extension: *,pem)
; and returns it as an aXis2000 format stack binary file (*.cb)
; with options to bin
;
;CATEGORY:
;	STAND ALONE: utilities
;
;CALLING SEQUENCE:
;	Result = READ_NSRRC_PEEM(file, BIN=BIN, GROUP=GROUP, NOFILTER=nof, VERBOSE=VERBOSE, _extra=e )
;
;INPUTS:
;	FILE - NSRRC PEEM image or stack (multi-image) data file
;
;KEYWORD PARAMETERS:
;	FILE	 file name
;	GROUP	 name of group calling this program
;	NOFILTER - use *.* as filter
; 	VERBOSE - if set, print additional documentation
;	_EXTRA - pass on other call line parameters
;
; COMMON BLOCKS:
;   AXIS_COM 			set of common blocks for AXIS
;	BSIF_COMMON 		for output stack parameters
;	NSRRC_SPEM_COM   	parameters specific to NSRRC SPEM read-in
;	stack_process_com   for output stack parameters
;	VOLUME_DATA			output stack
;
;MODIFICATION HISTORY:
; (26-Aug-13 aph) first version (start from read_nsrrc_spem.pro (NSRRC SPEM)
;-

function READ_NSRRC_PEEM, file=file, group=group, NOFILTER=nof, verbose = verbose,  _extra=e

on_error,2
@axis_com
@nsrrc_peem_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

; ---- determine if aXis2000 is running (if so, assume this was called from it)
axis_on  = 0
if widget_info(axis_ID, /active) NE 0 then axis_on = 1		; this works even if axis_ID is undefined

; information from Gung Chiang <gcyin@nsrrc.org.tw> on NSRRC PEEM format
; -------------
; single LARGE file:   *.pem
; [PEEM_HEADER]   (header, size of 39920 bytes at 32 bit system)
; [BODY]  	  	   (Energy stacks)*pixel_num*sizeof(unsigned short) (2 8-bit bytes for each element)
; ------------------
;  by trial and error I found, skip a byte array 39920 long then read each image as a 16-bit integer array of
;  size nx * ny
;

; in the header file , there is the following ???
; ----------------------------------------------
; no information provided about the header

if n_elements(file) eq 0 then begin  ;popup file dialog box
   filter ='*.pem'
   if keyword_set(nof) then filter='*'
   file=PICKFILE2(/Read, FILTER=filter, /LPATH, DEFPATH=defpath)
endif

s = 0
IF strlen(file) LE 0 THEN return, s

; ------- SHOULD READ USING IDL's binary file reader

openr,unit,file,/get_lun

; file has an ascii text header, "PEEM3D" - is that part of the 39920 byte header ?
; see what can be read using guesses
;check = 'Yes'
;a = 0. ; '                      '
;while check eq 'Yes' do begin
;	readu, unit, a
;	print, a
;	check = dialog_message(/question, 'Continue ? (default = Yes)')
;endwhile
; POINT_LUN, unit, 0	; return file to start

; --------- skip header info ----------
t =  39920
tmp =bytarr(t)		; test file (2005-Jun_23-Cu(110)+6MLMn+6MlCo-628-673_film area PhaseP9p9_01.Pem)
						; is 150.2 Mb or 150.16 Mb without header  ==> 220 or 221 images  (actually 226 !!)
readu, unit, tmp

; ------ get the image dimensions and number of images ---------
if n_elements(PEEM_nx) EQ 0 then PEEM_nx = 658
	if axis_on EQ 0 then $
		PEEM_nx = get_num(prompt='number of columns (nx)', val =PEEM_nx)  $
	else PEEM_nx = get_num(prompt='number of columns (nx)', val =PEEM_nx, group = axis_ID)

if n_elements(PEEM_ny) EQ 0 then PEEM_ny = 517
	if axis_on EQ 0 then $
		 PEEM_ny = get_num(prompt='number of rows (ny)', val =PEEM_ny) $
	else PEEM_ny = get_num(prompt='number of rows (ny)', val =PEEM_ny, group = axis_ID)

if n_elements(PEEM_nimg) EQ 0 then PEEM_nimg = 226
	if axis_on EQ 0 then $
		 PEEM_nimg = get_num(prompt='number of images', val =PEEM_nimg) $
	else PEEM_nimg = get_num(prompt='number of images', val =PEEM_nimg, group = axis_ID)

if n_elements(PEEM_dark) EQ 0 then PEEM_dark = 200
	if axis_on EQ 0 then $
		 PEEM_dark = get_num(prompt='dark level (ususally 200)', val =PEEM_dark) $
	else PEEM_dark = get_num(prompt='dark level (ususally 200)', val =PEEM_dark, group = axis_ID)


; generate local variables
nx = PEEM_nx & ny = PEEM_ny & n_img = PEEM_nimg

; test data nx = 658, ny = 517,  nimg = 226	; 340,186 pixels  = 0.68 Mb/image

img = intarr(nx, ny)

; ---------- get the set of images --------
t = ax_name(file) & fileshort = t(1)
print, 'reading ', n_img, ' NSRRC PEEM image(s) from file ' + fileshort

;readu, unit, img
;close, unit, /all
;free_lun, unit

; ------------  read full stack
avg_img = fltarr(nx, ny)
image_stack = intarr(nx,ny,n_img)

for i = 0, n_img-1 do begin
	axis_log, strtrim(string(i+1),2) + ' of ' + strtrim(string(fix(n_img)),2)
	readu, unit, img
	image_stack(*,*, i) = img-PEEM_dark
	avg_img = avg_img + float(img-PEEM_dark)/float(n_img)
endfor

close, unit, /all
free_lun, unit

; ------ test stack has 226 images from 628 to 673 eV; assumed to be equal spaced
; get photon energy parameters

if n_elements(PEEM_E_first) EQ 0 then PEEM_E_first = 628.
	if axis_on EQ 0 then $
		 PEEM_E_first = get_num(prompt='first energy', val = PEEM_E_first) $
	else PEEM_E_first = get_num(prompt='first energy', val = PEEM_E_first, group = axis_ID)

if n_elements(PEEM_E_last) EQ 0 then PEEM_E_last = 673.
	if axis_on EQ 0 then $
		 PEEM_E_last = get_num(prompt='last energy', val = PEEM_E_last) $
	else PEEM_E_last = get_num(prompt='last energy', val = PEEM_E_last, group = axis_ID)

; -- prepare energy (ev) array
ev=findgen(n_img)
E_step = (PEEM_E_last - PEEM_E_first)/float(n_img-1)
ev = PEEM_E_first + eV*E_step

; - get dwell time / image
if n_elements(PEEM_dwell) EQ 0 then PEEM_dwell = 1000.		; default is 1000 msec dwell
	if axis_on EQ 0 then $
		 PEEM_dwell = get_num(prompt='exposure time/image (milli-seconds)', val = PEEM_dwell) $
	else PEEM_dwell = get_num(prompt='exposure time/image (milli-seconds)', val = PEEM_dwell, group = axis_ID)

; -- prepare arrays for filename_ev_ms list
filename_list = strarr(n_img)
filename_ev_list = strarr(n_img)
for i =0, n_img-1 do begin
	filename_list(i) = string(FORMAT="(A,'_',i3)",fileshort,0)
	filename_ev_list(i)=string(FORMAT="(A,'_',i3,' ',F8.3,'  ',F8.2)", fileshort,i, ev(i), PEEM_dwell(0))
endfor
filename_ev_msec_list = filename_ev_list

continue:

; - get dwell time / image
if n_elements(PEEM_pixel_size) EQ 0 then PEEM_pixel_size = 0.02		; default is 20 nm
	if axis_on EQ 0 then $
		 PEEM_pixel_size = get_num(prompt='pixel size (microns)', val = PEEM_pixel_size) $
	else PEEM_pixel_size = get_num(prompt='pixel size (microns)', val = PEEM_pixel_size, group = axis_ID)

; ----- define parameters needed for stack write out
x=findgen(nx) & y=findgen(ny)
x = x*PEEM_pixel_size
y = y*PEEM_pixel_size

n_cols = nx  & n_rows = ny
n_data = n_img
svec = size(image_stack)
x_start = x(0)
x_stop = x(nx-1)
y_start = y(0)
y_stop = y(ny-1)
x_step = (x_stop - x_start)/n_cols
y_step = (y_stop - y_start)/n_rows


; --------- optional binning of the stack ---------------------
if n_elements(PEEM_bin) EQ 0 then PEEM_bin = 2
if axis_on EQ 0 then $
	 PEEM_bin = get_num(prompt='bin (1 = no binning)', val =PEEM_bin) $
else PEEM_bin = get_num(prompt='bin (1 = no binning)', val =PEEM_bin, group = axis_ID)
PEEM_bin = floor(PEEM_bin)
bin=PEEM_bin
if bin GT 1 then begin
  	if (float(n_cols)/float(bin))-fix(n_cols/bin) GT 0 $
      		OR (float(n_rows)/float(bin))-fix(n_rows/bin) then begin
  	    nct =fix(n_cols/bin)*bin  &  nrt = fix(n_rows/bin)*bin
  	    xtra_c = n_cols - nct + 1
  	    xtra_r = n_rows - nrt + 1
  	    print, 'truncating image to ', fix(nct), fix(nrt)
  	    x_stop = x_stop*float(nct)/float(n_cols)  ; correct axis length for truncation
  	    y_stop = y_stop*float(nrt)/float(n_rows)
        tmp = image_stack                      ;force size to integer mutiple of bin
    	image_stack = intarr(nct,nrt,n_data)
		image_stack = tmp(0:n_cols-xtra_c,0:n_rows-xtra_r,0:n_data-1)
;		avg_img = avg_img(0:n_cols-xtra_c,0:n_rows-xtra_r)
  	 endif
     n_cols = fix(n_cols/bin)
     n_rows = fix(n_rows/bin)
  	 x_step = (x_stop - x_start)/n_cols		; re-generate X, Y arrays for binned stack
	 y_step = (y_stop - y_start)/n_rows
;	 x = findgen(n_cols) & x = x_start + x*x_step
;	 y = findgen(n_rows) & x = y_start + y*y_step
     print, 'bin image ', fix(bin), ' times,  to ',n_cols,'  x', n_rows
     image_stack = rebin(image_stack,n_cols,n_rows,n_data)
; ------- also bin the average image image
;     avg_img = rebin(avg_img, n_cols, n_rows)
 endif

;----- write stack as aXis format binary (*.ncb)
t=ax_name(file)
outfile=t(0) + ax_sep() + fileshort+'_' + t(2) + '.ncb'
; print, 'sending data to write_sb. name ', outfile
ax_wait		; put on egg timer
stack_wb, outfile

; ------ generate and return the average image
xl = 'X (um)  dwell (ms) = ' + string(PEEM_dwell(0), format='(F8.1)')
yl = 'NSRRC PEEM AVERAGE image   Y (um)'

s = {t:'2d', x:x, y:y, d: avg_img, $
         e: mean(ev), xl: xl, yl:yl, dl: fileshort}

return, s

end
