; Copyright (c) 1998-2009  A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_CGO
;
;LAST CHANGED: ----------------------------------- 308-Jul-09
;
;PURPOSE:
;	This procedure generates component maps by using
; pixel-by-pixel conjugate gradient optimization curve fitting. It handles
; an arbitary number of images fit to up to 8 components.
; Except for individual image by image read in mode, AX_CGO
; assumes all image files have similar SIZE.
; AX_CGO ignores (x,y) scales (no interpolation)
;
;CATEGORY:
;	STAND ALONE: image, stack analysis
;
;CALLING SEQUENCE:
;	AX_CGO [,  axis=axis, comp_images = comp_images, coeff = coeff, $
;         help=help, images = images, stack = stack, verbose = verbose]
;
;CALLED FROM AXIS:
;	->Stacks->maps ->CGO curve fit

;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 = file with names of files to store 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:
;	Component maps, residuals etc are output to disk.
;
;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 file of absorption coefficients
; in a  standard format (see example below, written by AX_CGO); or by manual input.
; The CGO_OPTIMIZE curve fit procedure is then applied on a pixel-by-pixel basis.
; Component maps are generated. written to files, and displayed either in
; IDL windows (if run stand alone) or in AXIS buffers.
;
; EXAMPLE:
; Absorption coefficient data file format (version as of Jul-01)
; # components
;  (line per energy; E, sig-1  sig-2  sig-3   sig-4)
;  (labels for the components)
; -----------------------
;      4.00000
;      281.824  0.000687768  0.000191430  0.000572100  0.000390513
;      282.199  0.000678105  0.000192197  0.000618820  0.000361705
;      282.574  0.000753450  0.000189890  0.000628327  0.000574156
;      282.949  0.000803872  0.000187579  0.000724339  0.000398248
; ....
;  mtx
;  san
;  pipa
;  fg
;
;MODIFICATION HISTORY:
; (27-jul-01 aph) first developed from ax_svd
; (06-oct-01 aph) add xl to temp plots to allow re-use
; (14-may-03 aph) modify paramter files to store names only
; (28-may-03 aph) force '.par' extension (pickfile, filter= & /write); comp_names
; (04-jun-03 aph) use ax_par_save & ax_par_load
; (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
; (21-jan-05 aph) increase default convergence to 1e-10
; (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 label
;                 truncate short name for component maps at the first space
;				  changed spectral plotting to stack_fit approach - keeps original labels
;-

pro ax_cgo,  axis=axis, comp_images = comp_images, coeff = coeff, $
    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_CGO'
    print,'Executes conjugate gradient optimization 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, '	COEFF  = fit parameter file with names of spectral files'
	print, '	COMP_IMAGES = list of component images'
	print, '	HELP   = print this how-to-use statement'
	print, '	IMAGES = file with list of images to be analysed'
    print, '	STACK = name of a binary format stack (*.ncb)'
	print, '	VERBOSE = print additional documentation of SVD procedure'
    return
ENDIF

; determine if AXIS is running (therefore may have called ax_cgo)
; 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, ' Image sequence analysis by conjugate gradient optimization (CGO) fits'
print, ' --------------------------------------------------------------------------------'
print, ' Input data - images'

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')
		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')
		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')
	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_CGO'
		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)
			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_CGO'
			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_CGO'
				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) -------------
d={x:ev,d:ev}
d = ax_sort_mono(d, /axis)
ev = d.x

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

; -------------- define starting coefficents of fit (all 0 to begin with)
xp = fltarr(ncomp)
xp_start = xp

; -------- let user choose all positive coefficients if desired
t = dialog_message('Force only positive values ?', /question, /default_no)
if t(0) EQ 'Yes' then begin
  	pos_only = 1
  	text = 'forcing only positive coefficients'
  	if axis_on then widget_control, Bad_ID=Bad_ID, Uprompt, SET_VALUE=text
  	print,text
endif else pos_only = 0

; ---- optional limits
IF NOT keyword_set(limits) AND pos_only NE 1 then begin
	t = dialog_message('Set limits ?', /question, /default_no)
	if t(0) EQ 'Yes' then begin
		limits = fltarr(2)
		if axis_on THEN limits(0) = get_num(prompt='lower limit', group = axis_ID) else $
			limits(0) = get_num(prompt='lower limit')
		if axis_on THEN limits(1) = get_num(prompt='upper limit', group = axis_ID) else $
			limits(1) = get_num(prompt='upper limit')
		if limits(0) GT limits(1) then begin		; force lower/upper
			tt = limits(0) & limits(0) = limits(1) & limits(1)= tt
		endif
		tt = dialog_message('Swap if hit limits ?', /question)
		if tt(0) EQ 'Yes' then swap = 1 else swap = 0
	endif
ENDIF ELSE begin
	swap = 0 & limits = 0
ENDELSE

; ----- set tolerance
   tol = 1.e-010
   if axis_on then begin
	     tol = get_num(Prompt = ' Tolerance', val=tol, group=axis_ID)
	endif else  tol = get_num(Prompt = ' Tolerance', val=tol)

; ------ set maximum number of iterations
	numit = 100

; ------------ execute pixel-by-pixel CGO procedure --------------------
widget_control,/hourglass
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
		xp = xp_start		; begin fit from standard guess of coefficients
		xp = CG_OPTIMIZE(xp,'cgex_func', 'cgex_dfunc',  LINMIN='cgex_linmin', $
		        tolerance = tol, pos_only = pos_only, $
		        limits = limits, swap = swap, $
	            num_iter=numit, OBJECTIVE=obj, $
	     		A = A, b = od(*,x,y), _extra=e)
		for i = 0,ncomp-1 do comp_map(i,x,y) = xp(i)
	endfor
endfor

; ------- apply clip if desired ----
t = dialog_message('Clip at limits ?', /question, /default_no)
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 	file = pickfile2(/read, filter='*.sl') 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
for i = 0, ncomp-1 do begin
	widget_control,/hourglass
	tmp = comp_map(i,*,*) & tmp = reform(tmp)
	s = {t:'2d', x:xd, y:yd, d:tmp, E: 0.0, xl:'X', yl:'Y', dl: comp_names(i) + ' CGO '+ stackname_short}
	s = ax_order(s)
	filename = rootname + '_' + comp_names(i) + '.axb'
	t = axb_save(s, file=filename)

; ------- 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
	   	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
	filename = pickfile2(/WRITE, FILTER='*.ncb', title = 'File for residual stack', /LPATH)
	if strlen(file) GT 0 then stack_wb, filename
endif

; ------------  compute residuals
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' + ' CGO '+ stackname_short}
s = ax_order(s)
filename = rootname + '_chi.axb'
t = axb_save(s, file=filename)
if axis_on then begin
	CurBuf = 9
   	HANDLE_value, Data(CurBuf),s,/set
	Label(CurBuf) = s.dl
	Plotbuf, CurBuf
endif

return

end
