; Copyright (c) 1998-2010 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	RD_PEEM
;
; last changed ----------------------------- 21-Feb-10 (aph)
;
; PURPOSE
;   This function reads in a PEEM-2 image file in either the as-recorded 12-bit
;  or the converted 16-bit unsigned tif formats.
;  Dependending on keywords, it then corrects for gain, background effects;
;  applies median smoothing, and assigns the x,y scales and the photon energy
;  The data is returned as an AXIS internal format 2d structure
;
; CATEGORY:
;	Image processing.
;
; CALLING SEQUENCE:
;	Result = RD_PEEM(file [, bits12 = bits12, col12 = col12, Energy = energy, AXIS = AXIS, $
;         CCD_bgnd = CCD_bgnd, scale = scale, smth = smth, CCD_gain = CCD_gain, $
;         aoi_box = aoi_box, aoi_data = aoi_data, aoi_gain=aoi_gain, aoi_dark=aoi_dark, $
;		  region = region, sphinx=sphinx, NOFILTER=nof, _extra=e])
;
; INPUTS:
;	FILE 	name of file to convert (user prompted if not supplied)
;
; KEYWORDS:
;	AOI_box 	area of interest, defined by box, [x0, y0, x1, y1] acquired by PEEM-2
;					used to extract sub-image from dark and white files
;	AOI_DATA	if set, select AOI from input image
;	AOI_GAIN	if set, select AOI from CCD-gain image
;	AOI_DARK	if set, select AOI from dark image
;	AXIS		if set, indicates called from AXIS2000 widget
;	BITS12 		if set, this turns on the 12-bit read routine
;   COL12		if set, defines number of columns to drop in 12-bit mode
;	DWELL		exposure time for image in seconds
;	ENERGY 		sets energy PEEM image was recorded at
;	CCD_BGND 	CCD background image (2-d array) that is subtracted; must be same size as image
;	GROUP		name of group leader calling RD_PEEM
;	NOFILTER	if set, then pickfile dialog does not have filter set to '*.tif'
;	REGION		a 4-number vector [x_start, x_stop, y_start, y_stop]
;				defining the region of the image to be selected
;				(if scale set, in real space units. Otherwise, in pixels)
;	SCALE 		pixel size in microns
;	SMTH 		if set, applies a 3-point median smooth
;   SPHINX      if set, switched 16-bit TIF read to signed
;	CCD_GAIN	CCD background image (2-d array) that is divided; must be same size as image
;           	(gives response of CCD to uniform illumination)
;
; COMMON
;	AXIS_COM general common blocks for AXIS
;
; WARNINGS
;	1024x1024 images from CCD camera can be slow loading with older computers
;
; MODIFICATION HISTORY:
; (14-may-99 aph) first version
; (15-may-99 aph) save scale (pix_siz), region (cur, cll) values in common
; (08-jun-99 aph) add group to get_num call
; (25-jun-99 aph) report short gile name only
; (17-jul-99 aph) CCD_bgnd used
; (01-nov-99 ads) add nobits keyword & option for reading 12bit PEEM2 images
; (26-dec-99 aph) incorporated in AXIS 1.8c
; (23-oct-00 aph) 12-bit read in from arbitrary size - use Query_Tiff to obtain dimensions
;                  dump first 5 columns
; (20-apr-01 aph) getting 1024x1024 12-bit files to read
; (06-jun-01 aph) AXIS standard header; apply full background and gain image correction
;				  implement Area of Interest for gain, image and background
; (12-jul-01 aph) continuing development of improved PEEM read-in
; (31-jul-01 aph) handle bad file names properly; add axis keyword
; (04-aug-01 aph) add group keyword; remove intermediate plotting if group set
; (07-aug-01 aph) add dwell keyword and value to x-axs label if set
; (12-aug-01 aph) added keyword parameter to set 12-bit column drop; change GAIN correction to array
; (05-may-02 aph) add sphinx keyword to correct for signed format
; (12-may-02 aph) corrct name format in header.
; (16-jul-02 aph) read in 12-bit elmitec data stored in 16-bit TIF files correctly
; (23-may-04 aph) get keyword energy to be recognized and stop requesting a value
; (12-sep-07 aph) remove action of col12 function (drop first few columns); replace higher up with ax_cimage
;                 use axis_log instead of print
; (21-Feb-10 aph) fix up group ID problem (reading ALS PEEM2 format files)
;-

function Rd_PEEM, file, bits12=bits12, col12 = col12, energy=energy, axis = axis, $
         CCD_bgnd=CCD_bgnd, sphinx=sphinx, scale=scale, smth=smth, CCD_gain=CCD_gain, dwell = dwell, $
         aoi_box = aoi_box, aoi_data = aoi_data, aoi_gain=aoi_gain, aoi_dark=aoi_dark, $
         region = region, group = group, NOFILTER=nof, _extra=e
on_error,2
@axis_com

; determine if AXIS is running (therefore may have called ax_read_peem)
; either when AXIS parameter in ax_peem_com is set or if any widget active

if keyword_set(axis) then begin
	axis_on = 1
	if not keyword_set(group) then group = axis_ID  ; cover cases where called from axis but group not set
endif else axis_on = 0
 print, 'AXIS_ON set to ', axis_on

if n_elements(file) eq 0 then begin  ;popup file dialog box
   fltr='*.tif'
   if keyword_set(nof) then fltr='*'
   file=PICKFILE2(/Read, FILTER=fltr, /LPATH, DEFPATH=defpath)
endif
s = 0
IF strlen(file) EQ 0 THEN BEGIN
	axis_log, 'Non-existent file '+ file
	return, s
ENDIF

; define area of interest if aoi_box is set
if keyword_set(aoi_box) then begin
	aoi = aoi_box	; rename the pxiel limits for convenience
	aoi_width = aoi(1) - aoi(0) + 1
	aoi_height = aoi(3) - aoi(2) + 1
	print, 'Area of Interest (x0, y0, width, height)',aoi(0), aoi(2), aoi_width, aoi_height
endif

; --------------------------------- read in image as 12-bit or 16-bit TIFF
WIDGET_CONTROL, /Hourglass
name_bits = ax_name(file)
fileshort = name_bits(1)
IF keyword_set(bits12) then begin
;
; routine to read illegal 12-bit tif files from ALS PEEM2 system
; structure is 2 x 12 bit pixels split over 3x 8bit bytes
; files recorded on a PC which is 'little endian', byte order needs
; to be swapped if using a 'big endian' system - e.g. Mac or Sun.
;
; seem to be missing the first byte of each line - therefore skip first two pixels.
; (Added ADS 12/11/99)
;
    axis_log, 'RD_PEEM: reading in 12-bit image: ' + fileshort
	test = query_tiff(file,tif_info)
    img_x = tif_info.dimensions(0)
    img_y = tif_info.dimensions(1)
    print, 'size: ', img_x, ' x ', img_y
; -------
; there are 2 cases for mapping 2 *12-bits into 3 * 8-bits
	test = float(img_x/2.)-floor(img_x/2)
; CASE 1:  img_x  is EVEN		e.g. 2 12-bits takes 3 8-bits
    IF test EQ 0 then rawline = 3*floor(img_x/2) else $
; CASE 2:  img_x  is ODD		e.g. 3 12-bits takes 5 8-bits
    rawline = 3*floor(img_x/2) + 2
;	print, ' each line will have ', rawline, ' bytes'
	rawimg=bytarr(rawline,img_y)
    openr, filelun, file, /GET_LUN , /SWAP_IF_BIG_ENDIAN
    readu, filelun, rawimg
    close, filelun, /all
    a = uintarr(img_x,img_y)
    b = uintarr(img_x,img_y)
    aone = uintarr(img_x,img_y)
; read image by extracting 3 12-bit pixels from every pair of 16-bit pixels  (* = each line)
    FOR i=2,rawline-3,3 DO BEGIN
; 		if i LT 8 then print, '3 bytes: ', uint(rawimg(i,0)), uint(rawimg(i+1,0)), uint(rawimg(i+2,0))
        ki=(i*2)/3
        aone(ki,*)=ishft(uint(rawimg(i,*)),4) + ishft(uint(rawimg(i+1,*) AND 240),-4)
        ki=ki+1
        aone(ki,*)=ishft(uint(rawimg(i+1,*) AND 15),8) + uint(rawimg(i+2,*))
;        if i LT 8 then print, aone(ki-1,0), aone(ki,1)
    ENDFOR
;  --   amplify signal to occupy top 12bits, packing bottom 4 bits with zero
    aone = ishft(aone,4)		; THIS CONVERTS 12-bit to 16-bit data !!

; ------- have data - mask on area of interest if needed
    IF keyword_set(aoi_data) then begin
		b = aone(aoi(0):aoi(1), aoi(2):aoi(3))
		aone = b
	ENDIF

;  --   invert image (so's right way up.)
     a = reverse(aone,2)
; ---- remove last 5 columns - header information - seems to be variable in length
; ---- NB if can find length, could just skip prior to read-in !!
;	 if NOT keyword_set(col12) then col_drop = 5 else col_drop = col12
;	 if col_drop LT 0 then col_drop = 0
;	 print, '12-bit: columns to drop = ', col_drop
;     t = size(a) &	b=a(col_drop:(t(1)-1),0:(t(2)-1)) & a = b
;    t = size(a) & print,' converted to ', t(1), ' x ', t(2)

ENDIF ELSE BEGIN
;
; ------ read in 16-bit tifs ----------------
; -------  mask on area of interest if needed
	axis_log, ' RD_PEEM: reading 16-bit TIF from ' + fileshort
	IF keyword_set(sphinx) then sign_mode=0 else sign_mode=1
	IF keyword_set(aoi_data) then begin
		a = read_tiff(file, unsigned=sign_mode, sub_rect =[aoi(0), aoi(2), aoi_width, aoi_height])
	ENDIF ELSE BEGIN
		a = read_tiff(file, unsigned=sign_mode)
	ENDELSE
; --- subtract 2^15 when read in signed
	if min(a) LT -16000 then a = a + 2^15
;	reverse y-order
	a = reverse(a,2)
	t = size(a)
	print, 'size ', t(1), ' x ', t(2)
ENDELSE


; Then - process background and GAIN (white, or CCD pattern noise) corrections

; ------------- BACKGROUND  / dark correction ------------
; CCD_bgnd is the dark image - this should be on the same pixel scale, OR
; be such that the aoi paramters can be applied to extract the dark
; NB much more efficient in stack processing to cut out aoi_dark and
; simply supply the correct dark-image
IF keyword_set(CCD_bgnd) THEN begin
; ----- extract area of interest from  CCD_bgnd if AOI_dark set ------
	IF keyword_set(aoi_dark) then CCD_bgnd = CCD_Bgnd[aoi(0):aoi(1),aoi(2):aoi(3)]
; check that CCD_bgnd and data are the same size
	tg = size(CCD_bgnd)  &  ta = size(a)
	if tg(1) NE ta(1) OR tg(2) NE ta(2) then begin
		axis_log, ' CCD backgrond (DARK) does not match image size'
		test = DIALOG_MESSAGE('apply Area of Interest ?',/question)
		if test EQ 'Yes' then begin
			aoifile = PICKFILE2(/Read, title = 'Area of Interest file', FILTER='*.aoi')
			if strlen(aoifile) NE 0 then begin
				openr, lun, aoifile, /get_lun
				txt = '' &	aoi = intarr(4)
				readf, lun, txt
				close, lun & free_lun, lun
				start = strpos(txt,'Box1')
				txt = strmid (txt, start+5)
				reads, txt, aoi
				t = aoi(2) & aoi(2) = aoi(1) & aoi(1) = t
;				print, 'Area of interest: X: ', aoi(0), ' to ',aoi(1), '. Y: ', $
;						aoi(2), ' to ', aoi(3)
			endif else CCD_Bgnd = 0.
		endif else CCD_Bgnd = 0.
	endif
	a = a - CCD_Bgnd	; if CCD_Bgnd is an image, must be same dimension as a !!
ENDIF

; --------- load GAIN correction image ------ **** Adapted to Area of Interest ******
IF keyword_set(CCD_gain) THEN BEGIN
	gain = CCD_gain
; ----- correct gain for area of interest if needed
	IF keyword_set(aoi_gain) then gain = gain[aoi(0):aoi(1),aoi(2):aoi(3)]
; check that gain and data are the same size
	tg = size(gain)  &  ta = size(a)
	if tg(1) NE ta(1) OR tg(2) NE ta(2) then begin
		axis_log, ' gain array does not match image size'
		test = DIALOG_MESSAGE('apply Area of Interest ?',/question)
		if test EQ 'Yes' then begin
			aoifile = PICKFILE2(/Read, title = 'Area of Interest file', FILTER='*.aoi')
			if strlen(aoifile) NE 0 then begin
				openr, lun, aoifile, /get_lun
				txt = '' &	aoi = intarr(4)
				readf, lun, txt
				close, lun & free_lun, lun
				start = strpos(txt,'Box1')
				txt = strmid (txt, start+5)
				reads, txt, aoi
				t = aoi(2) & aoi(2) = aoi(1) & aoi(1) = t
;					print, 'Area of interest: X: ', aoi(0), ' to ',aoi(1), '. Y: ', $
;							aoi(2), ' to ', aoi(3)
			endif else gain = 1.
		endif else gain = 1.
	endif
ENDIF ELSE gain = 1.

; ------------- execute GAIN correction
a = float(a)/float(gain)

; ------------- replace NaN and inf values with adjacent pixel values ---
ntests=1
test = a(where(finite(a) EQ 1)) & check = n_elements(a) - n_elements(test)
WHILE check NE 0 AND ntests LT 5 DO BEGIN
;	print, 'Ratio generated ', check, ' math errors. Attempt to correct ', ntests
	cor = where(finite(a) EQ 0) +1      ;take next point
  ;	cor(where(cor EQ n_elements(a))) = n_elements(a) - 2  ; unless at edge !
	a(where(finite(a) EQ 0)) = a(cor)
	test = a(where(finite(a) EQ 1)) & check = n_elements(a) - n_elements(test)
	ntests = ntests + 1
ENDWHILE
test = a(where(finite(a) EQ 1)) & check = n_elements(a) - n_elements(test)
IF check GT 0 THEN BEGIN
	print, 'Residual math errors - replace with average'
	a(where(finite(a) EQ 0)) = mean(a(where(finite(a) EQ 1)))
ENDIF

; -- keyword set SMTH option
IF keyword_set(smth) THEN a = median(a,3)   ; median filter to remove hot spots

; -- keyword set ENERGY - if not specified in call then prompt for it
IF not keyword_set(energy) THEN BEGIN
	IF keyword_set(group) OR axis_on THEN $
		Energy = get_num(Prompt = 'Energy', Val = Energy, group = group) $
			ELSE Energy = get_num(Prompt = 'Energy', Val = Energy)
ENDIF

; -- keyword set SCALE option - set (x,y) scale if requested ----------
IF keyword_set(scale) THEN pix_siz = scale
IF not keyword_set(scale) THEN BEGIN
	IF keyword_set(group) OR axis_on THEN $
		pix_siz = get_num(Prompt = 'um / pixel', Val = pix_siz, group = axis_ID) $
	 	 ELSE pix_siz = get_num(Prompt = 'um / pixel', Val = pix_siz)
ENDIF
if keyword_Set(dwell) then begin
	xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f5.2," s")', $
      energy, dwell)
endif else xl = string(FORMAT='("x (um)     E = ",f8.3," eV")', energy)

t = size(a)
x = pix_siz*indgen(t(1))  & y = pix_siz*indgen(t(2))
;
; -- keyword set REGION option - extracts subset of image ---------------
IF keyword_set(region) THEN BEGIN
	t = size(region)
	IF t(1) NE 4 THEN BEGIN  ; ask user for xmin, xmax, ymin, ymax
		xlow = min(x, max=xhi) & ylow = min(y, max=yhi)
		IF keyword_set(group) OR axis_on THEN BEGIN
			xlow = get_num(Prompt = 'X_min',Val=xlow, group = group)
			xhi  = get_num(Prompt = 'X_max',Val=xhi, group =  group)
			ylow = get_num(Prompt = 'Y_min',Val=ylow, group = group)
			yhi  = get_num(Prompt = 'Y_max',Val=yhi, group =  group)
		ENDIF ELSE BEGIN
			xlow = get_num(Prompt = 'X_min',Val=xlow)
			xhi  = get_num(Prompt = 'X_max',Val=xhi)
			ylow = get_num(Prompt = 'Y_min',Val=ylow)
			yhi  = get_num(Prompt = 'Y_max',Val=yhi)
		ENDELSE
	ENDIF ELSE BEGIN
		xlow = region(0)  &  xhi = region(1)
		ylow = region(2)  &  yhi = region(3)
	ENDELSE
; -------------- save region for future use -------------
	cll(0) = xlow & cll(1) = ylow
	cur(0) = xhi  & cur(1) = yhi
	xind=where(x GE xlow AND x LE xhi)
	yind=where(y GE ylow AND y LE yhi)
	d = a(xind(0):xind(n_elements(xind)-1), yind(0):yind(n_elements(yind)-1))
	s = {t:'2d', x: x(xind), y: y(yind), d: d, e: energy, xl:xl, yl: 'PEEM   ', dl: fileshort}
ENDIF ELSE BEGIN
	s = {t:'2d', x: x, y: y, d: a, e: energy, xl:xl, yl: 'PEEM   ', dl: fileshort}
ENDELSE

return, s
END

