; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	IMG_WARP
;
;LAST CHANGED: ----------------------------------- 16-Dec-11
;
;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
; 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)
;-

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)
endelse
; -------------- set up fiducials on the first image -----------
; -------- display as pixels to full width ------
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)
Npts = get_num(val = 4, prompt='# of reference pts(>2)', group=axis_id)
xr=fltarr(npts)  & yr=fltarr(npts)
FOR i = 0, npts-1 do begin
	WIDGET_CONTROL, Uprompt, SET_VALUE='Warp: select pt. '+ string(i+1)
	cursor, tx, ty, /device, /down
	xr(i)= tx(0) & yr(i) = ty(0)
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))
;	tmp = {t:'2d',x: xn, y: yn,d: dn, xl: tmp.xl,yl: tmp.yl, dl: 'T ' + tmp.dl}
; --------- store truncated data ------------
;	CurBuf = 3
;	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)

; ----  Save first image as warped image and start list of warped images ----------
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)

; -------- Process second (and other) images -----
for img_n = 1, n_files-1 do begin
	ax_wait,/off		; correct a freeze problem - perhaps from sav_nsls ??
	if NOT keyword_set(list) then begin
		; get image to warp --------
		bt = buffer(prompt = 'Select new image', group = axis_ID)
		IF bt LT 0 then RETURN
		HANDLE_VALUE, Data(bt), ref
		IF n_tags(ref) EQ 0 THEN RETURN
		IF ref.t NE '2d' then RETURN
	endif else begin
		file = filename_list(img_n)
		ref = read_bnl(file)
	endelse

	xs=fltarr(npts)  & ys=fltarr(npts)
	scale_factor = main_size/(max([n_elements(ref.x),n_elements(ref.y)]))
	x_size = scale_factor*n_elements(ref.x)
	y_size = scale_factor*n_elements(ref.y)
	dn = congrid(ref.d,x_size,y_size)
	erase
	tvscl, ax_imscl(dn)

	FOR i = 0, npts-1 do begin
		axis_log, ' ', /clear
		text = 'Image ' + strmid(string(img_n+1),2) + ' of ' + strmid(string(n_files),2)
		WIDGET_CONTROL, Uprompt, SET_VALUE='Warp: select pt. '+ string(i+1),/append
		cursor, tx, ty, /device, /down
		xs(i)= tx(0) & ys(i) = ty(0)
	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
	dn = POLY_2D(dn, P, Q) ; , 2, x_final, y_final, cubic = 0.5)
;    	dn = warp_tri(xs,ys,xr,yr,dn)  ;, $
;    	    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)
	ref.dl = 'W '+ ref.dl
	ref = {t:'2d', x:tmp.x, y:tmp.y, d:dn, xl:ref.xl, yl:ref.yl, dl:ref.dl}
	if not keyword_Set(list) then begin
		CurBuf = 0
		HANDLE_VALUE, Data(CurBuf), ref, /set
		Label(CurBuf) = ref.dl
		PlotBuf, CurBuf
	  	t = filename_list(img_n)
		wfilename_list(img_n+1) =  'w' + strmid(t,1,strlen(t))
	endif else begin
		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 warled files
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

END
