; Copyright (c) 1998-2013 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		AX_ANGLE_STACK
;
;LAST CHANGED: ----------------------------------- 10-Nov-13 (aph)
;
; PURPOSE:
;	This function processes a set of folders containing multi-images stacks
; for generating an angle stack of a single image or image difference (map)
; used for tomography data read-in. It can aslo combine all images at all angles into a single stack
;
; CATEGORY:
;	Input / output utilty; aXis2000 or stand alone
;
; CALLING SEQUENCE:
;	Result = AX_GEN_STACK(GROUP=GROUP, CHANNEL=CHANNEL, MAP=MAP, IMAGE_NUM=IMAGE_NUM)
;
; INPUTS: none
;  [asks for *.lst is a file giving path then a set of {folder name, angle}]
;
; KEYWORDS
;	CHANNEL			number of data channel (default = 1)
;	GROUP 			ID of calling program
;	IMAGE_NUMBER	number of *.xim image to use (requested if not supplied)
;	MAP				if set, use 2 different images to generate angle stack of a 2-E map
;	REGION			region of data
;   STACK           combine all images into a single stack where E(i)= 1000*angle + image#
;        e.g.  angle = 2.5 (rounded to nearest 0.1)   with 3 images(0,1,2) have E = 2500.0, 2501.0, 2502.0
;
; OUTPUTS: binary format (*.ncb) stack written
;
; COMMON:
;	AXIS_com
;	BSIF_com
;	stack_process_com
;   VOLUME_DATA
;
; USAGE:
;
; format of list file for tomography applications is
; path
; name1, angle1
; name2, angle2
; . . . .
; where name1  is the name of the FOLDER containing the stack measured at angle 1
;
;
; MODIFICATION HISTORY:
; (09-Jun-06 gaj) first version for use in axis and stand alone (for tomography)
; (10-jun-06 aph) linked into aXis2000
; (30-jul-09 aph) changed analcom to stack_process_com (update to stack_process)
; (16-Feb-10 aph) added @ to stack_process.com inclde !!
; (01-Jun-11 aph) add image_num  parameter to select a specific image; get working !!
; (31-Jul-11 aph) generalize to get maps from 2-energies
;                 change function name from ax_gen_stack to ax_angle_stack
; (10-Nov-13 aph) generate full stack of all images
;-

Function ax_angle_stack, group=group, map=map, channel=channel, region=region,  image_num = image_num, stack = stack

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

axis_on=0
if NOT keyword_set(group) then group = 0 else begin
	t = size(axis_ID, /type)
	if t NE 0 then begin
		group = axis_ID
		axis_on=1
	endif else begin
		axis_on=0
	endelse
endelse

; --------- read list of folders and angles ------------
list_filename=pickfile2(title = "List of files", filter='*.lst', /LPATH, DEFPATH=defpath)
get_lun,lun
openr,lun,list_filename
sl_path = ''
readf, lun, sl_path		; read first line as path
t = ax_name(sl_path)
filebase=t(0)+t(1)

widget_control,hourglass=1

; ------ determine number of files to process
got_started = 0

while ~ EOF(lun) do begin
	text = ' '
    readf, lun, text
    this_filename = strtrim(strmid(text, 0, strpos(text,',')),2)
    this_angle = strmid(text, strpos(text,',')+1)
    this_name = strtrim(string(this_filename),2)
    if (got_started eq 0) then begin
; ------ OK the first line is a path .... but is it the right one ? !!!!
       filename_list = [this_name + ax_sep() + this_name]
       angle = [float(this_angle)]
       got_started = 1
    endif else begin
        filename_list = [filename_list,this_name + ax_sep() + this_name]
        angle = [angle, float(this_angle)]
    endelse
 endwhile
nfiles=n_elements(filename_list)
; for i = 0, nfiles-1 do print, angle(i),'   ',filename_list(i)
filename_ev_msec_list = strarr(nfiles)
close,lun  &  free_lun,lun

t =ax_name(list_filename)
file_short = t(1)
axis_log, 'generating an angle-stack from ' + file_short

; ---- define data channel (if not defined in keyword)
if NOT keyword_set(channel) then begin
	if axis_on then channel = get_num(prompt = 'Data channel (1 = PMT)', val = '1', group = Axis_ID) $
	   else channel = get_num(prompt = 'Data channel (1 = PMT)', val = '1')
	channel = fix(channel(0))
endif

; ---- define region (if not defined in keyword)
if NOT keyword_set(region) then begin
	if axis_on then region = get_num(prompt = 'Data region', val = '1', group = Axis_ID) $
	   else region = get_num(prompt = 'Data region', val = '1')
	region = fix(region(0))
endif

; -------------- define the signal and background images ----------
	image_n= 1		; for first file read-in
	img_bgnd = 0
	img_signal = 1
; ----- identify parameters specific to map mode
IF keyword_set(map) then begin
	if keyword_set(group) then $
	        img_bgnd = get_num(prompt='image number for BACKGROUND', val = img_bgnd, group=group) $
	   else img_bgnd = get_num(prompt='image number for BACKGROUND', val = img_bgnd)
	if keyword_set(group) then $
	       img_signal = get_num(prompt='image number for SIGNAL', val = img_signal, group=group) $
	   else img_signal = get_num(prompt='image number for SIGNAL', val = img_signal)
	image_bs=[img_bgnd,img_signal]
ENDIF else begin
	IF NOT keyword_set(stack) then begin
; ----------  define number of the single image, unless it is a stack mode read
	if NOT keyword_set(image_num) then begin
		if axis_on then image_n = get_num(prompt = 'Image number', val = '1', group = Axis_ID) $
		          else image_n = get_num(prompt = 'Image number', val = '1')
		image_n = fix(image_n(0))
	endif else image_n = fix(image_num)
	ENDIF
ENDELSE

; ------------- read a single image to set up stack dimension
; read a single file to determine the size of image_stack
sep= ax_sep()
file = sl_path + filename_list(0) + '.hdr'
print, 'reading in image # ',image_n,'  from ', file
tmp = read_sdf(file, one_image = image_n, channel=channel, region=region, /silent, /no_interpolate, /no_display)
; help, tmp, /struct
; energy = tmp.e
; axis_log, 'Energy ' + string(energy)
nx = n_elements(tmp.x)
ny = n_elements(tmp.y)
n_cols = nx  & n_rows = ny
x_start = tmp.x(0)
x_stop = tmp.x(n_elements(tmp.x)-1)
y_start = tmp.y(0)
y_stop = tmp.y(n_elements(tmp.y)-1)
xstep = (x_stop - x_start)/n_cols
ystep = (y_stop - y_start)/n_rows
image_stack = fltarr(nx, ny, nfiles)

if keyword_set(stack) then begin
; ------- read all images from each stack at each angle ------------
; determine how many images are in each stack and generate full array to contain all images
	tmp = read_sdf(file, channel=channel, region=region, /silent, /no_interpolate, /no_display, /no_save)
	t = size(image_stack)
	n_E = t(3)
	image_stack_all = fltarr(nx, ny, n_E*nfiles)
	filename_ev_msec_list_all = strarr(n_E*nfiles)

	angle_all = fltarr(n_E*nfiles)
	for i = 0, nfiles-1 do begin
		file = sl_path +  filename_list(i) + '.hdr'
		axis_log, 'reading in stack ' + filename_list(i)
		tmp = read_sdf(file, channel=fix(channel), region=region, /silent, /no_interpolate, /no_display, /no_save)
		ibase=fix(i*n_E)
		for j = 0, n_E-1 do begin		; loop through each stack
			k = fix(ibase+j)
; -------- generate stack axis parameter:  1000*angle + image#
			angle_all(k) = 100*(floor(10*angle(i))+j)
			filename_ev_msec_list_all(k) = filename_ev_msec_list(j) + '   ' + string(angle(i),format='(F7.1)')
;			print, k  , '   ', filename_ev_msec_list_all(k)
			image_stack_all(*,*,k) = image_Stack(*,*,j)
		endfor
	endfor
;  ---------- transfer all_stacks entities to normal stack entities
	image_stack = image_stack_all
	filename_ev_msec_list = filename_ev_msec_list_all

; ---------- write file for assisting later reconstruction of aligned single-angle stacks
	t = ax_name(list_filename)
	out_file = t(0) + t(1) + '_angle_E.lst'
	openw, iunit, out_file, /get_lun
	printf, iunit, list_filename
	printf, iunit, 'Filename      Angle        Energy'
	for i = 0, nfiles-1 do begin
		for j = 0, n_E-1 do begin
			sep= ax_sep()
			t=strpos(filename_list(i), sep)
			fn  = strmid(filename_list(i),t+1)
			printf, iunit, fn, angle(i), ev(j)
		endfor
	endfor
	print, 'filenames, angles, energies stored in ', out_file
	close, iunit & free_lun, iunit

endif else begin
; ----------- read  an image or map from the set of angles
	for i = 0, nfiles-1 do begin
		file = sl_path +  filename_list(i) + '.hdr'
		print, 'reading in image # ',image_n,'  from ', file
		IF keyword_set(map) then begin
			tmp=read_sdf(file, channel=fix(channel), /map , image_num=image_bs,region=region, /silent, /no_interpolate)
			help, /struct, tmp
		ENDIF ELSE BEGIN
			tmp=read_sdf(file, channel=fix(channel), one_image = image_n, region=region,  /silent, /no_interpolate)
		ENDELSE
		tmp.e = angle(i)
		t = ax_name(file)
		filename_ev_msec_list(i) = filename_list(i) + ' ' + strtrim(string(angle(i)), 2)
		axis_log, strtrim(string(i),2) + ': angle: ' + strtrim(string(tmp.e, format='(f5.1)'),2) + ' file: ' + t(1)
		image_stack(*,*,i) = tmp.d
	endfor
endelse

; -------- write stack
svec = size(image_stack)

;  ------------ set default filename
t = ax_name(list_filename)
filename = t(0) + t(1) + '.ncb'
filename = pickfile2(/write, filter='*.ncb', /Lpath, defpath = WritePath) ;, get_path=WritePath)
if filename EQ '' then return, filename		; allow user to cancel without writing

; ---- force extension
t = ax_name(filename)
filename = t(0) + t(1) + '.ncb'
if not keyword_set(stack) then ev = angle else ev=angle_all

stack_wb, filename
widget_control,hourglass=0

return, filename
END
