; Copyright (c) 1998-2015 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	IMG_WARP
;
;LAST CHANGED: ----------------------------------- 11-Jan-15
;
;PURPOSE:
;	This procedure warps an image to another image to get a common (x,y) scale
; by matching a set of fiducial points identifed
; by a cursor position selected by the user on first image
;
;CATEGORY:
;	AXIS:  image analysis
;
;CALLING SEQUENCE:
; for procedures:
;	IMG_WARP
;
;CALLED FROM AXIS:
;	->Images->waep
;
;INPUTS: none (uses axis buffers which must be pre-loaded)
;
;KEYWORDS: none
;
;OUTPUTS:
;	The result is written to buffer 0
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;RESTRICTIONS:
;	Occasionally does not work properly for unknown reasons
;
;PROCEDURE:
;	the images are first expanded to full axis window size [11-Jan-14 aph - removed ]
; to improve position accuracy - note the routines used are limited to integer pixel
; uses Poly_warp and poly_2d routines to match the two images
; the expanded image is then congridded back to its original dimensions
;
;MODIFICATION HISTORY:
; (08-jun-99 aph) extracted from aaa_axis - not functioning
;                add group ID to get_num call
; (24-feb-00 aph) add group ID to buffer call; pixels only
; (26-apr-00 aph) develop pixel-based CUT same; large -> small version
; (20-may-00 aph) AXIS standard documentation
; (26-Aug-05 aph) add write out of the w*.sl file (screwed up!)
; (05-Mar-06 aph) correct to write first file and save list as *.sl not *.nc !!
; (16-dec-11 aph) force plotting to be done on MainImg (in a test, it spawned a new graphics window)
; (15-Jan-15 aph) modify to allow warp-alignment of 2 images or from a list (this was never completed in 2011)
;-

PRO IMG_WARP, list=list
@axis_com
on_error,2

if NOT keyword_set(list) then begin
	HANDLE_VALUE, Data(CurBuf), tmp	; get data in current buffer
	IF n_tags(tmp) EQ 0 THEN RETURN
	IF tmp.t NE '2d' then RETURN
	n_files = 2
endif else begin
; --- read in first file from list
   	stack_readlist, list, filename_list
	n_files = n_elements(filename_list)
; -------- check if file exists - this handles subdirectory changes
	file1 = filename_list(0)
	tmp = read_bnl(file1)				; list files have to be *.nc (NSLS old files)
endelse

goto, continue   ; temporary skip of the full scale plotting

; -------- display as pixels to full width ------    WHY ?????
main_size = gr_scale*360
x_final = n_elements(tmp.x)
y_final = n_elements(tmp.y)
scale_factor = main_size/(max([n_elements(tmp.x),n_elements(tmp.y)]))
x_size = scale_factor*n_elements(tmp.x)
y_size = scale_factor*n_elements(tmp.y)
dn = congrid(tmp.d,x_size,y_size)
x_step = (tmp.x(2) - tmp.x(1))/scale_factor
xn = findgen(x_size)*x_step + tmp.x(0)
y_step = (tmp.y(2) - tmp.y(1))/scale_factor
yn = findgen(y_size)*y_step + tmp.y(0)

; ------- display first image
SetGraf, 'MainImg'
loadct,0,/silent
erase
tvscl, ax_imscl(dn)

continue:
; -------------- set up fiducials on the first image -----------
w_npts = get_num(val = 4, prompt='# of reference pts(>3)', group=axis_id)
xr=fltarr(w_npts)  & yr=fltarr(w_npts)
print, '----- cursor positions of fiducial points -------------'
FOR i = 0, w_npts-1 do begin
	WIDGET_CONTROL, Uprompt, SET_VALUE='Warp: select pt. '+ string(i+1)
	cursor, tx, ty, /device, /down
	xr(i)= fix(tx(0)) & yr(i) = fix(ty(0))
	print, i, xr(i), yr(i)
endfor

; ------ select a range that will encompass the warped data ------
	cll(0) = get_num(prompt='X-min',val=cll(0), group = axis_id)
	cur(0) = get_num(prompt='X-max',val=cur(0), group = axis_id)
	cll(1) = get_num(prompt='Y-min',val=cll(1), group = axis_id)
	cur(1) = get_num(prompt='Y-max',val=cur(1), group = axis_id)
	Ilo = DIndex(cll(0),cll(1),tmp)
	Ihi = DIndex(cur(0),cur(1),tmp)
;	print, 'Ilo', Ilo, ' data:',tmp.x(0), tmp.y(0)
;	print, 'Ihi', Ihi, ' data:',tmp.x(n_elements(tmp.x)-1), tmp.y(n_elements(tmp.y)-1)
  	xn = tmp.x(Ilo(0):Ihi(0))	; cut out the desired data from tmp
	yn = tmp.y(Ilo(1):Ihi(1))
	dn = tmp.d(Ilo(0):Ihi(0),Ilo(1):Ihi(1))
	test = where(tag_names(tmp) EQ 'E',count)
	if count EQ 1 then energy = tmp.e else energy = 0.
	tmp = {t:'2d',x: xn, y: yn,d: dn, e: energy, xl: tmp.xl,yl: tmp.yl, dl: tmp.dl + ' warp-ref'}
; --------- store truncated data ------------
	CurBuf = 4
	tmp.dl = tmp.dl
 	HANDLE_VALUE, Data(CurBuf), tmp, /set    ; store converted image in Buffer 4
  	PlotBuf,CurBuf

; --- define final (x,y) range as that of first file (truncated if done)
x_final = n_elements(tmp.x)
y_final = n_elements(tmp.y)


;start processing the rest of the list of  images to be warped ----------
if keyword_set(list) then begin
	t = ax_name(file1)
	file = t(0) + 'w' + strmid(t(1),1,strlen(t(1))) + '.' + t(2)
	file = sav_nsls(tmp, file = file)
	t = ax_name(file)
	wfilename_list = strarr(n_files + 1)
	wfilename_list(0) =  t(0)
	wfilename_list(1) =  t(1) + '.' + t(2)
endif

; -------- Process second (and other) images -----
xs=fltarr(w_npts)  & ys=fltarr(w_npts)
for img_n = 1, n_files-1 do begin
	if NOT keyword_set(list) then begin
		bt = buffer(prompt = 'Select new image', group = axis_ID)	; get next  image to warp --------
		CurBuf = bt(0)
		IF bt LT 0 then RETURN
		HANDLE_VALUE, Data(bt), tmp
		IF n_tags(tmp) EQ 0 THEN RETURN
		IF tmp.t NE '2d' then RETURN
		PlotBuf, CurBuf
	endif else begin
		ax_wait,/off					; correct a freeze problem - perhaps from sav_nsls ??
		file = filename_list(img_n)
		tmp = read_bnl(file)			; list files have to be *.nc (NSLS old files)
	endelse

	dn = tmp.d		; routine manipulates the dn array

;	scale_factor = main_size/(max([n_elements(tmp.x),n_elements(tmp.y)]))
;	x_size = scale_factor*n_elements(tmp.x)
;	y_size = scale_factor*n_elements(tmp.y)
;	dn = congrid(tmp.d,x_size,y_size)
;	erase
;	tvscl, ax_imscl(dn)

	text = 'Image ' + strmid(string(img_n+1),2) + ' of ' + strmid(string(n_files),2)
	axis_log,text
	print, '----- positions in image ', img_n, ' ---------'
	FOR i = 0, w_npts-1 do begin
		WIDGET_CONTROL, Uprompt, SET_VALUE='Warp: select pt. '+ string(i+1),/append
		cursor, tx, ty, /device, /down
		xs(i)= fix(tx(0)) & ys(i) = fix(ty(0))
		print, i, xs(i), ys(i)
	endfor
;	print, 'Warp: X-points -  (data, ref, diff)'
;	print, fix(xs), fix(xr), fix(xs-xr)
;	print, 'Warp: Y-points - (data, ref, diff)'
;	print, fix(ys), fix(yr), fix(ys-yr)
;	polywarp, xs, ys, xr, yr, 1, p, q, /double
;	dn = POLY_2D(dn, p, q) ; , 2, x_final, y_final, cubic = 0.5)
    	dn = warp_tri(xs,ys,xr,yr,dn, /tps)  ;, $
;    	    output_size=[n_elements(tmp.x),n_elements(tmp.y)], /extrapolate)
;		tvscl, ax_imscl(dn)
    	help, dn
;    	print, 'Match to the mesh of first file ', x_final, y_final
	dn = congrid(dn,x_final, y_final)
;	tmp.dl =tmp.dl + ' warped'
	test = where(tag_names(tmp) EQ 'E',count)
	if count EQ 1 then energy = tmp.e else energy = 0.		; make sure that the energy for the file is defined
	ref = {t:'2d', x:tmp.x, y:tmp.y, d:dn, e: energy, xl:tmp.xl, yl:tmp.yl, dl:tmp.dl}
	if not keyword_Set(list) then begin
		CurBuf = 5								; assume images to warp-align are in Buffers 1,2
		HANDLE_VALUE, Data(CurBuf), ref, /set	; output images in Buffers 4,5
		Label(CurBuf) = ref.dl
		PlotBuf, CurBuf
	endif else begin
		t = filename_list(img_n)
		wfilename_list(img_n+1) =  'w' + strmid(t,1,strlen(t))
		t = ax_name(file)
		fileshort = 'w' + strmid(t(1),1,strlen(t(1))) + '.' + t(2)
		file = t(0) + fileshort
		file = sav_nsls(ref, file = file)
		wfilename_list(img_n+1) = fileshort
	endelse
ENDFOR

; --------- write out stack_list for warped files
if keyword_set(list) then begin
	t = ax_name(list)
	wlist_short = 'w_' + t(1)
	wlist = t(0) + wlist_short + '.sl'
	axis_log, 'list of warped files = ' + wlist
	stack_writelist, wlist, wfilename_list
endif

END
