; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_SVD
;
;LAST CHANGED: ----------------------------------- 23-Apr-11
;
;PURPOSE:
;	This procedure executes the SVD procedure to generate component maps
; from an arbitary number of images for an arbitrary number of components.
; Will run stand alone (prompting user for AXIS format files) or from AXIS.
;
;CATEGORY:
;	STAND ALONE: image, stack analysis
;
;CALLING SEQUENCE:
;	AX_SVD[,  axis=axis,  coeff = coeff, comp_images = comp_images, $
;         help=help, images = images, stack = stack, verbose = verbose]
;
;CALLED FROM AXIS:
;	Stacks->SVD

;INPUTS:
;	All inputs are through keywords. User is prompted for missing data.
;
;KEYWORDS:
;	/AXIS - if on, then called from AXIS widget
;	COEFF = fit parameter file with names of spectral files
;	COMP_IMAGES = list of component images (*.axb)
;	/HELP - set to print how-to-use statement
;	IMAGES = file with list of images to be analysed
;   STACK = name of a binary format stack (*.ncb)
;	/VERBOSE  - print additional documentation of SVD procedure
;
;OUTPUTS:
;	No explicit outputs.   A new window is created if necessary
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;	stack_process_com
;	COMMON volume_data, image_stack
;	bsif_com
;
;PROCEDURE:
;	The images are first read in - either as a binary stack (if stack keyword set)
; or from the list of images, or from pickfile dialog.
; The absorption coefficents for each component at the energies of the images
; are then read in, either through interpolation of a reference spectrum
; (read in via pickfile dialog), or from a standard absorption coefficient file
; (see example below) or by manual input.
; The SVD procedure is then applied.
; Component maps are generated. written to files, and displayed either in
; IDL windows (if run stand alone) or in AXIS buffers.
;
;EXAMPLE:
; example of file to provide list of images (same as stack_list files)
;
;	f:\stxm\9909\10-13\			; path to images
;	a1013_01.axb
;	a1013_02.axb
;	a1013_03.axb
;	a1013_04.axb
;	a1013_05.axb
;	a1013_06.axb
;	a1013_07.axb
;	a1013_08.axb
;	a1013_09.axb
;
; example of file to provide mass absorption coefficents
;	2900  155000     8200  76700
;	3400   29900   139600  52000
;   DVB
;   EGDMA
;	SF-normalized values		; ignores ANY text after ncomp lines and labels
;	species	281 285	 288.3	305	; each with n_image entries (abs. coeff. for that species)
;	DVB   (divinylbenzene)  EGDMA (acrylate)
;
;MODIFICATION HISTORY:
;	Developed by ANDY WINESETT, NCSU , 1998, 1999
;	Modified from SGU's version by GEM 7/15/99 to do 3 image SAP stacks
;	updated 7/16/99, GEM to save data files in XRD format
;	updated March 5, 1999 (SGU)
; (29-jul-99 aph) generalised to a dialog
; (10-aug-99 aph) installed in AXIS; use reform to remove need to reload maps
; (17-aug-99 aph) auto-truncate to match size of all images NB does NOT use (x,y) scales
; (14-oct-99 aph) auto-extract intensity coeff from files if coeff set
; (17-oct-99 aph) further improvements to file handling stuff; allow use of any file format
; (23-oct-99 aph) plot coefficents for feedback to user
; (27-oct-99 aph) correct error in (x,y) scales of images
; (01-jan-00 aph) AXIS standard documentation
; (27-feb-00 aph) add groupID to get_text
; (11-apr-00 aph) add residual and chi square evaluation
; (22-aug-00 aph) force interpolation
; (19-oct-00 aph) increase to 8 components max
; (28-jul-01 aph) neaten up while developing AX_CGO; ensure stand-alone operation
; (12-aug-01 aph) correct error from 28-jul-01 change
; (07-oct-01 aph) add xl to temporary files to allow use later
; (13-dec-01 aph) fix error on component map name definition
; (29-dec-01 aph) add E to 2d structure
; (26-jan-02 aph) set correct path for par file
; (13-May-03 aph) convert coeff file to list of reference spectra files
; (04-jun-03 aph) force '.par' extension (pickfile, filter= & /write); fix comp_names for no blanks
; (04-jun-03 aph) use ax_par_load and ax_par_save
; (30-dec-03 aph) add ax_sort_mono to force stack E-scale and reference spectral scales to be monotonic
; (20-jan-04 aph) use dialog_message to control residual stack writing
;				  ensure E-value (0.0) added to component maps
; (01-apr-04 aph) fix ERROR in writing component maps (duplicate E)
; (12-dec-04 aph) reformat component spectra to allow other uses
; (01-feb-05 aph) fix ax_sort_mono to use structure (ricochet change from lox change)
; (17-Sep-08 aph) replace analcom with stack_process_com
; (30-Jul-09 aph) add stack name to component map; remove 'clip at limits' option
;                 truncate short name for component maps at the first space
;				  changed spectral plotting to stack_fit approach - keeps original labels
; (28-dec-09 aph) extend to up to 16 spectra (needed for XRF fitting)
; (23-Apr-11 aph) optional selection of limits of component maps
;-

pro ax_svd,  axis=axis,  coeff = coeff, comp_images = comp_images, $
    help=help, images = images, stack = stack, verbose = verbose

@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

IF keyword_set(help) THEN BEGIN
    print,'AX_SVD'
    print,'Executes SVD 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, '	COEFF  = fit parameter file (list of filenames of reference spectra'
	print, '	HELP   = print this how-to-use statement'
	print, '	COMP_IMAGES = list of component images'
	print, '    VERBOSE = print additional documentation of SVD procedure'
    return
ENDIF

; determine if AXIS is running (therefore may have called ax_svd )
; 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, ' Singular value decomposition 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'
		if axis_on then	widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text $
		   else print, text
		file = pickfile2(/read, filter='*.sl', title = 'File with list of images', /LPATH, DEFPATH=defpath)
		if strlen(file) EQ 0 then return
	endelse
	stack_readlist, file, image_list
	nimg = n_elements(image_list)
	svd_e = fltarr(nimg)
ENDIF

; ------------------- READ IN IMAGE DATA --------
if keyword_set(stack) then begin  ;read in a binary format stack---------
	t = size(stack) 	; use type identifier to see if stack is a string
	if t(1) NE 7 then begin
		stack = pickfile2(/read, title = ' Select stack file', filter='*.ncb', /LPATH, DEFPATH=defpath)
		if strlen(stack(0)) EQ 0 then return
	endif else begin
		check = findfile(stack)
		if strlen(check(0)) EQ 0 then return
	endelse
	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)
	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)
	stackname_short=t(1)
endif else begin
; --------- do not have an image list or stack, so try for image-by-image read in
	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 LE 0 then return
	svd_e = fltarr(nimg)
	image_list = strarr(nimg)
	file = pickfile2(/read, title = ' OD image 1', filter='*.axb', /LPATH, DEFPATH=defpath)
	if strlen(file) EQ 0 then return
	image_list(0) = file
endelse

; -------- read in individual files if stack was not read in--------
if NOT keyword_Set(stack) then begin
	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 AX_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))
			file = pickfile2(/read, title = text, filter=last_image_ext, /LPATH, DEFPATH=defpath)
			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 AX_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'
				if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text $
				    else print, 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))
			if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text $
				    else print, 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
endif

; -------------- define coefficents at energies of images -------------
print, ' Input data - intensity coefficients'
t = ' '
IF NOT keyword_set(coeff) THEN BEGIN
    t = dialog_message('Read fit parameter file ?', /question)
	if t(0) EQ 'Yes' then begin
	    par_file = pickfile2(title = ' Select fit parameter file', $
	           filter='*.par', /LPATH, DEFPATH=defpath )
	    if strlen(par_file(0)) NE 0 then coeff = par_file
	endif
ENDIF
IF keyword_set(coeff) THEN BEGIN
; ------- read spectral file info from a file
	check = findfile(coeff)
	if strlen(check(0)) NE 0 then pars = ax_par_load(coeff) else goto, get_user_info
	if n_tags(pars) EQ 0 then begin
		goto, get_user_info
	endif else begin
		ncomp = pars.n
		comp_names = pars.names
; ---------- check if any have spaces - truncate at first space
		for i=0, n_elements(comp_names)-1 do begin
			test= strpos(comp_names(i),' ')
			if test GT 0 then comp_names(i) = strmid(comp_names(i), 0, test)
			comp_names(i) = strtrim(comp_names(i),2)
		endfor
		comp_files = pars.files
	endelse

; ------ if fails then ask user to locate reference files one-by-one
ENDIF ELSE BEGIN
 get_user_info:
	if axis_on then ncomp = get_num(prompt = 'Number of components',val=ncomp, group=axis_ID) $
	   else ncomp = get_num(prompt = 'Number of components',val=ncomp)
	if ncomp LE 0 then return
	comp_names = strarr(ncomp)
	comp_files = strarr(ncomp)
; ------- read coefficent information from spectra files
	for i = 0, ncomp-1 do begin
	 	text = 'Spectral file for component ' + strcompress(string(fix(i+1)))
		if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
		comp_files(i) = pickfile2(title = ' Spectrum of component '+ strtrim(string(i),2), $
	           filter='*.txt', /LPATH, DEFPATH=defpath )
	    if comp_files(i) EQ '' then return
		tmp = spc_load(file = comp_files(i))
;; ------ define name - keep short as used in map name
		text = 'short name of component ' + strcompress(string(fix(i+1)))
		if axis_on then comp_names(i) = get_text(prompt = text, val = tmp.dl, group = axis_ID) $
		    else comp_names(i) = get_text(prompt = text, val = tmp.dl)
		 test= strpos(comp_names(i),' ')
		if test GT 0 then comp_names(i) = strmid(comp_names(i), 0, test)
		comp_names(i) = strtrim(comp_names(i),2)
	endfor

; ------ save coefficient array for later use (optional)
	par_file = pickfile2(filter='*.par', /write, path = DefPath, $
	           title = 'Name of fit parameter file')
	if strlen(par_file(0)) NE 0 then ax_par_save,par_file,ncomp,comp_names,comp_files
ENDELSE

; check that the stack data is monotonic - if not, force monotonic
; ----- (aph; added 30-dec-03 - interpol failures may be cause of 'baulky' fits) -------------
;   ------ converted to structure format 1-feb-05 aph) --------------------
d={x:ev,d:ev}
d = ax_sort_mono(d, /axis)
ev = d.x

; read in files and interpolate to the energy scale of the image set
	A=fltarr(ncomp, n_elements(svd_e))

; interpolate and store in AXIS data buffer
model = fltarr(ncomp, n_elements(ev))
for i = 0, ncomp-1 do begin
	tmp = spc_load(file=comp_files(i))
	tmp = ax_sort_mono(tmp,/axis)			; FORCE MONOTONIC !
	model(i,*)=interpol(tmp.d,tmp.x,ev)
	tmp = {t:'1d', d: model(i,*), x: ev, xl: 'X', dn:model(i,*), dl: 'interp '+ tmp.dl}
	a(i,*) = interpol(tmp.d, tmp.x, svd_e)
	CurBuf = i+1
	if curBuf GT 8 then curBuf = 1		; for 9-buffer aXis need to write-over previous maps
	if  n_tags(tmp) NE 0 AND keyword_set(AXIS) THEN BEGIN
		HANDLE_VALUE, Data(Curbuf), tmp, /SET
		PlotBuf, Curbuf
	endif
endfor

; ------- document parameters
if keyword_set(verbose) then begin
	print,'OD values of components'
	for i = 0, ncomp-1 do begin
		print, 'absorption coefficent for ', comp_names(i)
		print, a(i,*)
	endfor
endif

; ------------ execute SVD procedure --------------------
widget_control,/hourglass
SVDC,A,W,U,V
if min(w) LT 1E-30 then begin
	print, 'Warning!: small W values - output may not be valid'
	widget_control, Bad_ID=Bad_ID, Uprompt, $
	     SET_VALUE= 'Warning!: small W values - output may not be valid'
endif
comp_map=fltarr(ncomp,nx,ny)

for x=0,nx-1 do begin
	text = string(format='("Line ",i4," of ",i4)',x,nx)
	if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt,SET_VALUE= text
	for y=0,ny-1 do begin
		B=[od(0,x,y)]
		for i = 1, nimg-1 do B= [B,od(i,x,y)]
		temp=SVSOL(U,W,V,B)
		for i = 0,ncomp-1 do comp_map(i,x,y) = temp(i)
	endfor
endfor

; ------- apply clip if desired ----
t = dialog_message('Clip at limits ?', /question, /default_no)
; t='No'		; remove but leave in the code as may want later
if t(0) EQ 'Yes' then begin
	if n_elements(limits) NE 2 then limits = fltarr(2)
	if axis_on THEN limits(0) = get_num(prompt='lower limit', val = limits(0), group = axis_ID) else $
		limits(0) = get_num(prompt='lower limit',  val = limits(0))
	if axis_on THEN limits(1) = get_num(prompt='upper limit',  val = limits(1), group = axis_ID) else $
		limits(1) = get_num(prompt='upper limit',  val = limits(1))
	low = where(comp_map LT limits(0),count)
	if count GT 0 then comp_map(low) = limits(0)
	hi = where(comp_map GT limits(1),count)
	if count GT 0 then comp_map(hi) = limits(1)
endif

; ---- output component thickness maps
print, ' Output - component maps'
IF keyword_set(comp_images)  THEN BEGIN
; ------- output component images according to information from file
	check = findfile(comp_images)
	if strlen(check) EQ 0 then 	begin
		file = pickfile2(/read, filter='*.sl', /LPATH, DEFPATH=defpath)
	endif else file = comp_images
	stack_readlist, file, comp_image_list
ENDIF ELSE BEGIN
	t = STRARR(3)
	IF keyword_set(images) THEN t = ax_name(images)
	IF keyword_set(stack) THEN t = ax_name(stack)
	rootname = t(1)
	text = ' Root name for maps '
	if axis_on then	rootname = get_text(Prompt = text, val = rootname, group = axis_ID) $
	   else rootname = get_text(Prompt = text, val = rootname)
; -------- if no path included, use data_path
	rootname = t(0) + rootname
	t = ax_name(rootname)
	if t(0) EQ '' then rootname = data_path + t(1)

ENDELSE
overwrite_all = 0    ; reset the overwrite flag
;print, 'DEBUG'
;print, 'Names for component maps'
;for i = 0, ncomp-1 do print, comp_names(i)
for i = 0, ncomp-1 do begin
	widget_control,/hourglass
	tmp = comp_map(i,*,*) & tmp = reform(tmp)
	s = {t:'2d', x:xd, y:yd, e: 0.0, d:tmp, xl:'X', yl:'Y', dl: comp_names(i)+ ' SVD '+ stackname_short}
	s = ax_order(s)
	file = rootname + '_' + strtrim(strcompress(comp_names(i)),2) + '.axb'
;	print, 'Writing component map ',string(i+1),' to file ', file
	test = findfile(file)
	if test(0) EQ file AND overwrite_all NE 1 then begin
	    if axis_on then t = file_overwrite(group = axis_ID, val = 1) $
	        else t = file_overwrite(val = 1)
	    if t EQ 0 then file = pickfile2(/write, LFILE = file, /LPATH, DEFPATH=defpath)
	    if t EQ 2 then overwrite_all = 1
	endif
	file = axb_save(s, file=file)

; ------- store and plot the component maps in successive buffers
	if axis_on then begin
		if ncomp EQ 4 then offset = 5 else offset = 4
		if ncomp GT 4 then offset = 1
		CurBuf = i + offset
		if curBuf GT 8 then curBuf = 1		; for 9-buffer aXis need to write-over previous maps
	   	HANDLE_value, Data(CurBuf),s,/set
		Label(CurBuf) = s.dl
		Plotbuf, CurBuf
	endif
endfor

; ------- generate residuals stack and map of chi square
; generate image_stack from OD data (NB does not have to come from a stack!)
widget_control,/hourglass
image_stack=fltarr(nx, ny, nimg)
fit_stack = image_stack
for i = 0, nimg-1 do begin
	image_stack(*,*,i)= od(i,*,*)
	for j = 0, ncomp-1 do fit_stack(*,*,i) = fit_stack(*,*,i) + a(j,i)*comp_map(j,*,*)
endfor
image_stack = image_stack - fit_stack	; compute residual
; --------- query user to see if want to save residuals stack
if dialog_message('save residuals stack?',/question,/default_no, title = 'save residuals stack?') EQ 'Yes' then begin
	text = 'save residuals stack'
	if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text else print, text
	file = pickfile2(/WRITE, FILTER='*.ncb', title = 'File for residual stack', /LPATH, DEFPATH=defpath)
	if strlen(file) GT 0 then stack_wb, file
endif
; -------- compute standard deviation
chisq = fltarr(nx,ny)
for i = 0, nimg-1 do $
   chisq =  chisq + sqrt(image_stack(*,*,i) * image_stack(*,*,i))
chisq = chisq/(nimg-1)
s = {t:'2d', x:xd, y:yd, d:chisq, E: 0.0, xl:'X', yl:'Y', dl: 'residuals' + ' SVD ' + stackname_short}
s = ax_order(s)
file = rootname + '_chi.axb'
if test(0) EQ file AND overwrite_all NE 1 then begin
    if axis_on then t = file_overwrite(group = axis_ID, val = 1) $
        else t = file_overwrite(val = 1)
    if t EQ 0 then file = pickfile2(/write, LFILE = file, /LPATH, DEFPATH=defpath)
    if t EQ 2 then overwrite_all = 1
endif
file = axb_save(s, file=file)
if axis_on then begin
	CurBuf = 9
   	HANDLE_value, Data(CurBuf),s,/set
	Label(CurBuf) = s.dl
	Plotbuf, CurBuf
endif

return

end

