; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		AX_TOMO_READ
;
;LAST CHANGED: -----------------------------------  30-Oct-11 (aph)
;
; PURPOSE:
;	This function processes a set of folders of files containing multi-images stacks for tomography read-in
;
; CATEGORY:
;	Input / output utilty; stand alone
;
; CALLING SEQUENCE:
;	Result = AX_TOMO_READ(FILE=FILE, ONE_IMAGE=ONE_IMAGE, MAP=MAP, STACK=STACK, AXB=AXB, GROUP=GROUP)
;
; INPUTS: none
;
; KEYWORDS
;	CHANNEL		data channel to read
; 	FILE		list of folders to process
;	ONE_IMAGE  	read in one image at each angle
;	MAP        	read in map (diff of OD of 2 images) at each angle
;	STACK		read in stack of stacks (not yet implemented)
;   TIF			read in set of tif files (from SVD) (ABORTED 14-Jul-07 due to no color)
;	AXB			read in set of axb files (from SVD)
;
; OUTPUTS:
;	folder with selected files
;
; COMMON:
;	AXIS
;	BSIF
;	STACK_PROCESS_COM
;   VOLUME_DATA
;
; MODIFICATION HISTORY:
; (15-dec-05 aph) first version fur use  in axis and stand alone (for tomography)
; (17-may-06 gaj) added /no_interpolate keyword to call of read_sdf() to aviod rounding error of number of pixels
; (10-jul-06 gaj) fixed some bugs with Stack~Tomography~map
; (14-jul-07 aph) read set of '*.axb' files (component maps) and convert to *.ncb angle sequence
; (30-Jul-09 aph) replace analcom with stack_process_com
; (30-Oct-11 aph) try to solve why it does not run in VM version (no change to this code)
;-

Function ax_tomo_read, file=file, one_image=one_image, map=map, stack=stack, axb=axb, $ ; tif=tif, $
              channel=channel, group=group

@axis_com
@bsif_com
@stack_process_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 tomography folders/files", filter='*.lst', /LPATH, DEFPATH=defpath)
if strlen(list_filename) EQ 0 then return, 'no file'
get_lun,lun
openr,lun,list_filename
sl_path = ''
readf, lun, sl_path		; check if first line is only a path
t = ax_name(sl_path)
filebase=t(0)+t(1)

widget_control,hourglass=1

; ------ determine number of files to process
got_started = 0
;if strlen(t(1)) EQ 0 or  strlen(t(2)) EQ 0 then begin
while ~ EOF(lun) do begin
	text = ' '
    readf, lun, text
    this_filename = strmid(text, 0, strpos(text,','))
    this_angle = strmid(text, strpos(text,',')+1)
    if (got_started eq 0) then begin
; ------ OK the first line is a path .... but is it the right one ? !!!!
	   t=ax_name(this_filename)
	   this_filename = t(1)        ; only save main part of name in case extensions included
       filename_list = [string(this_filename)]
       angle = [float(this_angle)]
;       print, this_filename, this_angle
       got_started = 1
    endif else begin
	   t=ax_name(this_filename)
	   this_filename = t(1)        ; only save main part of name in case extensions included
       filename_list = [filename_list,string(this_filename)]
       angle = [angle, float(this_angle)]
;       print, this_filename, this_angle
    endelse
 endwhile
nfiles=n_elements(filename_list)
filename_ev_msec_list = strarr(nfiles)
close,lun  &  free_lun,lun

; select item to read if not defined by a keyword
if keyword_Set(one_image) then item = 'image'
if keyword_Set(map) then item = 'map'
if keyword_Set(tif) then item = 'tif'
if keyword_Set(axb) then item = 'axb'
if keyword_Set(stack) then item = 'stack'
if NOT keyword_Set(one_image) AND $
	NOT keyword_Set(map) AND $
	NOT keyword_Set(axb) AND $
	NOT keyword_Set(stack) then begin
		t=dialog_message('single image? (Yes; No=map)', /question, title='what to read')
		if t EQ 'Yes' then item = 'image' else item = 'map'
endif

; ---- define data channel (if not defined in keyword)

if not keyword_Set(axb) then begin
	if not keyword_set(channel) then begin
		channel = 1
		if keyword_set(group) then $
		       channel=get_num(prompt='channel', val = channel, group=group) $
		  else channel=get_num(prompt='channel', val = channel)
	endif
endif

; ---- defne path separator (operating system dependent)
sep= ax_sep()

; ------------- read a single image to set up stack dimension
; read a single file to determine the size of image_stack
;file = sl_path + filename_list(0) + sep + filename_list(0) + '.hdr'
;tmp=read_sdf(file, channel=1, one_image=one_image, /silent, /no_interpolate )
;energy = tmp.e
;nx = n_elements(tmp.x)
;ny = n_elements(tmp.y)
;image_stack = fltarr(nx, ny, nfiles)

print, 'Reading tomography data set as defined in ', list_filename

CASE item OF

'image': BEGIN
; ------------- read a single image for each angle

if keyword_set(group) then $
	one_image=fix(get_num(prompt='image number', val = one_image, group=group)) $
	else one_image=fix(get_num(prompt='image number', val = one_image))

file = sl_path + filename_list(0) + sep + filename_list(0) + '.hdr'
;print, 'first read in file: ', file
tmp = read_sdf(file, channel=fix(channel), one_image=one_image, /silent, /no_interpolate) ; read a single image to set up stack dimension
energy = tmp.e
nx = n_elements(tmp.x)
ny = n_elements(tmp.y)
image_stack = fltarr(nx, ny, nfiles)
axis_log, 'Energy ' + string(energy)
for i = 0, nfiles-1 do begin
	file = sl_path + filename_list(i) + sep + filename_list(i) + '.hdr'
	print, file
	tmp=read_sdf(file, channel=fix(channel), one_image=one_image, /silent, /no_interpolate)
	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
END


'map': BEGIN

; define parameters
img_bgnd = 0
img_signal = 1
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)
   file = sl_path + filename_list(0) + sep + filename_list(0) + '.hdr'
tmp=read_sdf(file, channel=fix(channel), one_image=img_bgnd, /silent, /no_interpolate )
energy = tmp.e
nx = n_elements(tmp.x)
ny = n_elements(tmp.y)
image_stack = fltarr(nx, ny, nfiles)
axis_log,'Background Energy ' + string(tmp.e)

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)
tmp=read_sdf(file, channel=fix(channel), one_image=img_signal, /silent, /no_interpolate )
axis_log, 'Signal Energy ' + string(tmp.e)

image_num=[img_bgnd,img_signal]
for i = 0, nfiles-1 do begin
	file = sl_path + filename_list(i) + sep + filename_list(i) + '.hdr'
	tmp=read_sdf(file, channel=1, /map , image_num=image_num, /silent, /no_interpolate)
	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

END

'stack': BEGIN

n_ang = nfiles
dim = indgen(3,nfiles)



for i = 0, nfiles-1 do begin
	file = sl_path + filename_list(i) + sep + filename_list(i) + '.ncb'
	print, file
	stack_rb, file
	im_size = size(image_stack)
	if im_size(0) EQ 2 then n_img = 1 else n_img = im_size(3)
	dim(0,i) = im_size(1)
	dim(1,i) = im_size(2)
	dim(2,i) = n_img
endfor

nx = max(dim(0,*))
ny = max(dim(1,*))
n_energies = max(dim(2,*))

print, nx, ny, n_energies, n_ang


image_stack_4 = fltarr(nx, ny, n_energies, n_ang)

ev_2 = fltarr(nfiles)
filename_ev_msec_list_2 = strarr(n_energies,nfiles)
shft = fltarr(2,nfiles)




sel_nrj = 0				; tilt angle offset
sel_nrj = get_num(val=sel_nrj,Prompt=('Select energy '+ ' 0,1,...') , group = Axis_ID)

for i = 0, nfiles-1 do begin
	file = sl_path + filename_list(i) + sep + filename_list(i) + '.ncb'
	print, file
	stack_rb, file
	sz = size(image_stack)
	image_stack_4(nx-sz(1):*,ny-sz(2):*,n_energies-sz(3):*,i) = image_stack(*,*,*)
	ev_2(i) = angle(i)
	filename_ev_msec_list_2(*,i) = filename_ev_msec_list(*)

	pr = where(total(image_stack(*,*,0),2) GT 150)
	mi = min(pr)
	ma = max(pr)
	wd = ma-mi
	hwd = wd/2
	ct = mi + hwd
	ofst = ct - 80
	print, ofst
	shft(1,i) = ofst

endfor



openw, iunit, 'out.aln', /get_lun

printf, iunit,'! Pixel shifts after alignment'
printf, iunit,'PLOTIT(' + string(nfiles)

printf, iunit,'0,' + string(shft(1,*)) + ',' + string(shft(0,*))
close, iunit
free_lun, iunit

image_stack = fltarr(nx, ny, nfiles)
image_stack(*,*,*) = image_stack_4(*,*,sel_nrj,*)

filename_ev_msec_list = strarr(nfiles)
filename_ev_msec_list(*) = filename_ev_msec_list_2(sel_nrj,*) + ' ' + string(ev_2(*))

print, size(image_stack)
n_energies = nfiles
n_img = nfiles
ev = ev_2



pr = where(total(image_stack(*,*,0),2) GT 150)
mi = min(pr)
ma = max(pr)
wd = ma-mi
hwd = wd/2
ct = mi + hwd
ofst = 80 - ct
print, ofst


END


'tif': BEGIN
; ------------------ read in set of tif files (e.g. from SVD of stacks) ----------------
; since the alignment and subsequent clipping ends up with different image sizes
; this routine pads each tif with 0's to a size 10% larger than the size of the first tif
; then truncates to the size of the largest (nx,ny) in the set of tifs
;
; PROBLEMS:  image_stack is only set up for gray scale
;            tif is rescaled to 0-255, does not retain oriignal intensity scale
; aborted development, but jsut copy 'axb' read in and replace axb_load with img_load (/tif)
END

'axb': BEGIN
; ------------------ read in set of '*.axb' files (e.g. component maps from SVD of stacks) ----------------
; since the alignment and subsequent clipping ends up with different image sizes
; this routine pads each image with 0's to a size 10% larger than the size of the first image
; then truncates to the size of the largest (nx,ny) in the set of tifs
;
   file = sl_path +  filename_list(0) + '.axb'
   print, 'Reading ', file
   tmp = axb_load(file=file)
   nx = n_elements(tmp.x)
   nxg = nx
   nx = nx +fix(0.1*nx)
   ny = n_elements(tmp.y)
   nyg = ny
   ny = ny +fix(0.1*ny)
   image_stack=fltarr(nx,ny,nfiles)

for i = 0, nfiles-1 do begin
	file = sl_path + filename_list(i) + '.axb'
	print, 'Reading ', file
	tmp = axb_load(file=file, /silent)
	nxt = n_elements(tmp.x)
	nyt = n_elements(tmp.y)
	if nxt GT nxg then nxg = nxt
	if nyt GT nyg then nyg = nyt
;	extend to (nx, ny)
	tarr = fltarr(nx,ny)
	tarr(0:nxt-1, 0:nyt-1) = tmp.d

	t = ax_name(file)
	filename_ev_msec_list(i) = filename_list(i) + ' ' + strtrim(string(angle(i)), 2)
	axis_log, strtrim(string(i),2) + ':  ' + strtrim(string(angle(i), format='(f5.1)'),2) + '  -  ' + t(1)
	image_stack(*,*,i) = tarr
endfor


END


ELSE:

ENDCASE


; -------- write stack
if  item eq 'stack' then begin

	svec = size(image_stack)
	n_cols = nx
	n_rows = ny

endif else begin
	svec = size(image_stack)
	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
endelse


;  ------------ 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)
; ---- force extension
t = ax_name(filename)
filename = t(0) + t(1) + '.ncb'
ev = angle
stack_wb, filename
widget_control,hourglass=0

return, filename
END
