; Copyright (c) 1998-2008 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_PC
;
;LAST CHANGED: ----------------------------------- 17-sep-08
;
;PURPOSE:
;	This procedure executes the PCOMP principle component analysis
; on an arbitary number of images for arbitrary number of components
;
;CATEGORY:
;	STAND ALONE: image analysis
;
;CALLING SEQUENCE:
;	AX_PC [,stack = ncb_file, standardize = standardize, axis=1]
;
;CALLED FROM AXIS:
;	->Stacks->Principle Components
;
; INPUTS:
;	all inputs are either from files identified by keywords or by a pickfile dialog
;
; KEYWORD PARAMETERS:
;	AXIS = if on, then called from AXIS widget
;	HELP = print how-to-use statement
;	IMAGES = file with list of images to be analysed
;   STACK = name of a binary format stack (*.ncb)
;	STANDARDIZE = standardize input (mean of zero and variance of one).

; OUTPUTS:
;	all outputs are either to files via dialog or to AXIS buffers

; COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;	stack_process_com
;	COMMON volume_data, image_stack
;	bsif_com
;
;PROCEDURE:
;  This procedure uses the IDL code (PCOMP.SAV - nb no source provided)
;  to compute the eigenvalues (strengths) and eigenvectors (principle components)
;  from a set of spectromicroscopy images
;
;MODIFICATION HISTORY:
; (28-oct-99 aph) first written (adapted ax_svd to use stack read_in)
; (31-dec-99 aph) AXIS standard documentation
; (27-feb-00 aph) add groupID to get_text
; (14-may-00 aph) make cancel abort first time
; (17-Sep-08 aph) replace analcom with stack_process_com
;-

pro ax_pc, stack=stack, images = images, standardize = standardize, $
           axis = axis, help=help

@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

IF keyword_set(help) THEN BEGIN
    print,'AX_SVD'
    print,'Executes PCOMP procedure to generate component thickness maps'
    print, 'Uses AXIS format image files (*.axb) as input/output'
    print, ' KEYWORDS: '
    print, '	AXIS   = if on, then called from AXIS widget'
    print, '    STACK = name of a binary format stack (*.ncb)'
	print, '	IMAGES = file with list of images to be analysed'
	print, '	HELP   = print this how-to-use statement'
	print, '    VERBOSE = print additional documentation of SVD procedure'
    return
ENDIF

; determine if AXIS is running (therefore may have called ax_pc )
; either when AXIS keyword is supplied or if any widget active
if  keyword_set(axis) then axis_on = 1 else axis_on = widget_info(/active)
print, ' '
print, ' Principal component image analysis '
print, ' -------------------------------------------'
print, ' Input data - images'
; assumes all image files have similar SIZE; and ignores (x,y) scales (no interpolation)

IF keyword_set(images) THEN BEGIN
; ------- read image information from file
	check = findfile(images)
	if check(0) EQ images then begin
	 	file = images
	endif else begin
		text = 'Cannot find '+ images + '. Please select image list file'
		print, text
		widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
		file = pickfile2(/read, filter='*.sl')
		if strlen(file) EQ 0 then return
	endelse
	stack_readlist, file, image_list
	nimg = n_elements(image_list)
	svd_e = fltarr(nimg)
ENDIF

if NOT keyword_set(images) AND NOT keyword_set(stack) then begin
	if axis_on then begin
	     nimg = get_num(Prompt = ' Number of images', val=nimg, group=axis_ID)
	endif else  nimg = get_num(Prompt = ' Number of images', val=nimg)
	if nimg EQ 0 then return
	svd_e = fltarr(nimg)
	image_list = strarr(nimg)
	text =  ' Select OD image 1'
	print, text
	widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
	file = pickfile2(/read, filter='*.axb')
	if strlen(file) EQ 0 then return
	image_list(0) = file
ENDIF

; ------------------- READ IN IMAGE DATA --------
if keyword_set(stack) then begin  ;read in a binary format stack---------
		stack_rb,stack
		test = size(image_stack)
		nimg = test(3) & nx = test(1) & ny = test(2)
		text = string(format='("Stack: ",i4," images. ",/,"       ",i4," x ",i4)', nimg, nx, ny)

		if keyword_set(axis) then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
		x_step = (x_stop - x_start)/n_cols
		y_step = (y_stop - y_start)/n_rows
		xd = x_start + findgen(n_cols)*x_step
		yd = y_start + findgen(n_rows)*x_step
		svd_e = ev
		od = fltarr(nimg,nx,ny)
		for i = 0, nimg-1 do od(i,*,*) = image_stack(*,*,i)
		t = ax_name(stack)
		data_path = t(0)
endif else begin
; -------- read in individual files --------
	file = image_list(0)
	; --------------------- make the default output path be the same as file locations
	t = ax_name(file)
	data_path = t(0)
	s = axis_read_image(file=file)
	if n_tags(s) EQ 0 then begin
		print, ' Image file ', file, ' not found. Terminating SVD'
		return
	endif
	svd_e(0) = 12398./sd.wavelength
	nx = fix(n_elements(s.x)) & ny = fix(n_elements(s.y))
	od = fltarr(nimg,nx,ny)
	od(0,*,*) = s.d
	xd = s.x   &   yd = s.y      ; save for constructing component maps
	for i = 1, nimg-1 do begin
		if keyword_set(images) then file = image_list(i) else begin
			text = ' Select OD image ' + strcompress(string(i+1))
			print, text
			widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
			file = pickfile2(/read, filter=last_image_ext)
			if strlen(file) EQ 0 then return
			image_list(i) = file
		endelse
		s = axis_read_image(file=file)
		if n_tags(s) EQ 0 then begin
			print, ' Image file ', file, ' not found. Terminating SVD'
			return
		endif
		svd_e(i) = 12398./sd.wavelength
		if n_elements(s.x) EQ nx AND n_elements(s.y) EQ ny then begin
			od(i,*,*) = s.d
		endif else begin
			nxt = fix(n_elements(s.x))
			nyt = fix(n_elements(s.y))
			if abs(nxt-nx) GT 5 or abs(nyt-ny) GT 5 then begin
				text = 'Image sizes differ by more than 5 pixels. Terminating AX_SVD'
				print, text
				widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
				return
			endif
			if nxt LT nx then nx = nxt
			if nyt LT ny then ny = nyt
			text = 'WARNING: image size reduced to ' + strcompress(string(nx)) + ' x ' + strcompress(string(ny))
			print, text
			widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
			tmp = fltarr(nimg,nx,ny)
			for j = 0, i-1 do tmp(j,0:nx-1,0:ny-1) = od(j,0:nx-1, 0:ny-1)
			tmp(i,0:nx-1,0:ny-1) = s.d(0:nx-1,0:ny-1)
			od = tmp
			xd = xd(0:nx-1)
			yd = yd(0:ny-1)
		endelse
	endfor
endelse

;----------- execute PCOMP procedure
widget_control,/hourglass
ncomp = nx*ny
pc_eigen = make_array(ncomp,/float,value=1.)
comp_spec = make_array(nimg,nimg,/float,value=1.)
t = reform(od,nimg,ncomp)
if keyword_set(standardize) then $
   pc_out = pcomp(t,eigenvalue=pc_eigen, /standardize, coefficients=comp_spec) $
   else pc_out = pcomp(t,eigenvalue=pc_eigen, coefficients=comp_spec)


; ---- output components
print, ' Output - principle components'
IF keyword_set(images) THEN t = ax_name(images)
IF keyword_set(stack) THEN t = ax_name(stack)
rootname = t(0)+t(1)
text = ' Root name for principle components'
rootname = get_text(Prompt = text, val = rootname, group = axis_ID)
t = ax_name(rootname)
if t(0) EQ '' then rootname = data_path + t(1)

; ---------- display eigenvalue map
x = indgen(nimg) & x = x + 1
s = {t:'1d', x: x , d:pc_eigen, dn: pc_eigen,  xl:'Number', dl:'eigenvalues'}
file = rootname + '_eigen.txt'
file = spc_save(s, file=file)
if axis_on then begin
	CurBuf = 1
   	HANDLE_value, Data(CurBuf),s,/set
	Label(CurBuf) = s.dl
	Plotbuf, CurBuf
endif

; ----------- fill AXIS buffers with 8 most significant components
for i = 1, 8 do begin
	widget_control,/hourglass
	tmp = comp_spec(*,i-1) & tmp = reform(tmp)
	s = {t:'1d', x:svd_E, d:tmp, dn:tmp, xl:'E (eV)', dl: 'pc ' + strcompress(string(i))}
	if axis_on then begin
		CurBuf = i+1
	   	HANDLE_value, Data(CurBuf),s,/set
		Label(CurBuf) = s.dl
		Plotbuf, CurBuf
	endif
; ----- save component signal
	file = rootname + '_' + strcompress(string(i)) + '.txt'
	file = spc_save(s, file=file)
endfor

return
end
