; Copyright (c) 1998-2016 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	READ_NSRRC_PEEM.PRO
;
;LAST CHANGED: ----------------------------------- 23-Feb-16
;
;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)
; (10-Feb-16 aph) read (ofset, nx,ny,n_img) info from the header
;                  use uint (unsigned 32-bit) not int (signed 32-bit)
; (23-Feb-16 aph) change read-in filename filter to find *.Pem ans well as *.pem files
;-

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   [ from email of 1 Sep 2013 )
;  Position
;  char TD[6];              // PEEM3D
;: Offset is unsigned int ->  8    F09B0000=>39920 (high bytes at back) // total bytes from the beginning of file to the Data
; image number->12                E2000000=>226 images
; nx->16                                   9202=>658
; ny->18                                   0502=>517
;
; From the actual (xxx) code that writes the *.pem files
;typedef struct
;{
;  char TD[6];              // PEEM3D
;  uns32 Offset;            // in bytes, the Distance from the File Header to 3d Data
;  uns32 PageNumber;        // Scan Number
;  uns16 nX;                // GetWindowSizeX
;  uns16 nY;                // GetWindowSizeY
;  uns16 Ori;               // The Orientation of the Graph. There are probably two kind of orientation of the Plot
;                                 //  One is the definition of the LABWINDOW plot intensityand the other
;                                 // Orientation is difined by the ordinaray BMP files,
;                                 / Here we define the 0 is the BMP     (TOP->down)
;                                 //                   1 is the Labwindow (Bottom->up)
;  double Ver;
;  int TopLayer;
;  int BottomLayer;
;  int Reserve[6];
;
; [APH] Do you record any other types of data files with the PEEM ?
;  e.g. a TEY or other spectral signal to act as an Io ?
;
; [GC, 1-sep-2013] Yes, I have record the Io and sample current in this file.
; It is in the structure I send you last time, but for now, it is difficult
; to get the data because the position is changed from file to file
; since the number is different from each file.
; I attached the file structure for PEEM at end of this letter
; for your reference. We actually saved Energy, sample current,
; Io mesh current for each images, and for the file structure
; itself, the energy can be nonlinear,
; however, I did not implement non-linear scan
; in the control program so I guess it is still linear for now
; ----------------------------------------------

if n_elements(file) eq 0 then begin  ;popup file dialog box
   filter =['*.Pem','*.pem; *.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

; ----------- read number of images and nx, ny from first header (NB each image has an associated header)
openr,unit,file,/get_lun
a=strarr(8)					; looks like the space reserved for the label is actually 8 characters not 6 ? [=2 32-bit words]
;					without the 2 extra character reads, the subsequent variables are gibberish
t=' '
for i = 0, 7 do begin
   readu, unit, t
   a(i) = t
endfor
print, a(0:5)
offset = ulong(1)
readu, unit, offset
print, 'offset = ', offset
n_img=ulong(1)
readu, unit, n_img
print, '# of images = ', n_img
nx = uint(1) & ny = uint(1)
readu, unit, nx
print, 'nx = ', nx
readu, unit, ny
print, 'ny = ', ny
close, unit, /all		; start again
free_lun, unit

; --------- skip header info ----------
openr,unit,file,/get_lun
if n_elements(offset) gt 0 then t =  offset else t = 39920  ; not clear if the 39920 value is same for all files
tmp =bytarr(t)
readu, unit, tmp

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

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

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
	avg_img = avg_img + float(img)/float(n_img)
endfor

close, unit, /all
free_lun, unit

; -------- get photon energy parameters
; try to guess from the file name
;  assumes just prior to the '.pem' extension is first-last  where first is first E, and last is last E

dot = strpos(/reverse_search, file, '.')
dash = strpos(/reverse_search,file,'-')
if dash gt 0 and dash lt dot then begin
	t_nd= dot-dash
	t = strmid(file, dash+1,t_nd)     &  print, t
	PEEM_E_last  = float(t)
	t = strmid(file, dash-t_nd+1, t_nd) & PEEM_E_first = float(t)
endif

if n_elements(PEEM_E_first) EQ 0 then PEEM_E_first = 0.
	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 = 1000.
	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 pixel size
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
	ax_wait		; put on egg timer
  	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)
	ax_wait, /off		; turn off egg timer
 endif

;----- write stack as aXis format binary (*.ncb)
t=ax_name(file)
outfile=fileshort+'_' + t(2) + '.ncb'		; t(0) + ax_sep() +
outfile=PICKFILE2(title='Name of stack file',/write, file=outfile, FILTER='*.ncb', /LPATH, DEFPATH=defpath)

; 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
