; Copyright (c) 1998-2018 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	IMG_ALIGN
;
;LAST CHANGED: ----------------------------------- 06 Jul 2018
;
;PURPOSE:
;	This procedure aligns one image to the spatial scale of the first file
; processed which acts as the reference for all subsequent files.
; The alignment can be a shift (1-point) or stretch-shift (2-point).
;
;CATEGORY:
;	AXIS: image analysis (AXIS - captive)
;
;CALLING SEQUENCE:
;	IMG_ALIGN, IMG, FNM, one_pt=one_pt
;
;CALLED FROM AXIS:
;	->Stacks->align->{shift, stretch-shift}
;
;INPUTS:
;	IMG  - image as a aXis2000 structure
;   FNM  - filename, including #
;
;KEYWORDS:
;	ONE_PT	if set, only align at 1-point, not 2
;
;OUTPUTS:
;	Aligned, clipped and (optionally) binned file is written as
;	s____.nc (1-point)  or a_____.nc (2-point)
;
;COMMON BLOCKS:
;	@AXIS_COM	axis standard set of common blocks
;	@BSIF_COM	nsls common
;
; PROCEDURE
; manual alignment of images
; has ability to linearly distort image to account for sample stretching !!
; if one_pt is set, then only shifts selected point to reference value
; otherwise asks user for 2 fiducials and stretch-shifts each file to MATCH first file
; generates align-shift files in buffers 5,8 (X) &  6,9 (Y), 4 (stretch)
; cuts user-selected region;
; meshes to user-defined pixelation with optimum binning
; stores processed data in netCDF format for use in STACK_ANALYZE
; saves file name to generate list for stack_analyze
;
; WARNINGS:
;		only prompts user for transformation paramaters on FIRST pass
;
;       CLEAR buffers (4,5,6),8,9 to reset image alignment
;
;MODIFICATION HISTORY:
; ( 5-Mar-98 aph) first version - worked
; (xx-apr-98 aph) changed something - did not correctly cut data - broken to 6-jun-98
; ( 6-jun-98 aph) replaced direct call to wrstx_ax with sav_nsls
; (14-jun-98 aph) axis_com
; (18-sep-98 aph) cosmetics - checking possible error
; ( 8-jun-99 aph) add group to get_num calls
; (23-aug-99 aph) fixing up 'freeze' when called from axis 1.7
;				 - CURSOR not supposed to be used !!
; (19-sep-99 aph)  complete debugging freeze-up problem
;			a) 'first' used in aaa_axis changed to 'firstimg' -> correct image display
;            b) turn-off hourglass (turned on by PlotBuf) - by a kluge
; (15-apr-00 aph) AXIS standard documentation
; (20-apr-00 aph) debug to see why does not work at fine scale
; (03-Aug-09 aph) write out files if use stretch-shift with different tag - axis_c gave wrong WritePath !
; (03-Apr-18 aph) write files as *.axb unless input was *.nc
; (06-Jul-18 aph) define energy = tmp.e at start; tmp is being corrupted !
;-

PRO img_algn, img, fnm, one_pt=one_pt
@axis_com
@bsif_com

SetGraf, 'MainImg'

;HANDLE_VALUE, Data(CurBuf), tmp
tmp = img
IF n_tags(tmp) EQ 0 THEN RETURN
IF tmp.t NE '2d' THEN RETURN
help, tmp, /struct
;	tmp = {t:'2d',x:tmp.x,y:tmp.y,d:tmp.d,e:energy, xl:tmp.xl,yl:tmp.yl,dl: tmp.dl}
energy = tmp.e
print, 'Energy (eV) ', energy
print, 'Input file name ', fnm
; ------------
 	ax_wait,/off			; make sure hourglass is turned off
 	!X.S = Xscl(CurBuf,*)   ; set up device to data conversions
 	!Y.S = Yscl(CurBuf,*)
 	!Z.S = Zscl(CurBuf,*)
 	WIDGET_CONTROL,UPrompt,SET_VALUE='Click feature 1'
;
  	cursor, x1, y1, /up, /data     ; identify FIRST fiducial
;
; NB IDL says should NOT use cursor in Widgets - BAD BOY !!
;
	curval = CONVERT_COORD(X1, Y1, /DEVICE, /TO_DATA)
  	CurInd1 = Dindex(x1(0), y1(0), tmp)
  	if keyword_set(one_pt) THEN begin
  		text=string(12398./sd.wavelength, CurInd1(0), (xlock(0) - x1(0)),  CurInd1(1), (ylock(0) - y1(0)), $
  		format="('E = ',f8.3,' eV',/,'   pixel   del(um)',/, 'x  ',i4,5x,f6.3,/,'y  ',i4,5x,f6.3)")
  		WIDGET_CONTROL,Uprompt, Set_value = text
  		print, text
  	endif
    if NOT keyword_set(one_pt) THEN BEGIN                                                                   ; stretch/shift
  		WIDGET_CONTROL,UPrompt,SET_VALUE='Click feature 2'                                                  ; stretch/shift
  		cursor, x2, y2, /data, wait=4     ; identify SECOND fiducial                                         ; stretch/shift
  		CurInd2 = Dindex(x2(0), y2(0), tmp)                                                                 ; stretch/shift
  		text=string(12398./sd.wavelength, CurInd1(0), CurInd2(0),  CurInd1(1), CurInd2(1), $
  		format="('E = ',f8.3,' eV',/,'   pt-1    pt-2',/,'x  ',i4,5x,i4,/,'y  ',i4,5x,i4)")
  		WIDGET_CONTROL,Uprompt, Set_value = text
  		print, text
    endif
  	lock = 0		; LOCK = 0 means define conversion parameters first pass
  					; LOCK = 1 means conversion parameters LOCKED
; save positions of first point in buffers 5, 6
  	algn = {t:'1d', x:fltarr(2), d:fltarr(2), xl:'E', dn:fltarr(2), dl:'X-shift-1'}
  	HANDLE_VALUE, Data(5),tmp1
  	first_pass=0
  	IF n_tags(tmp1) EQ 0 THEN BEGIN
  		lock = 1                                         ; FIRST PASS  KLUGE !!
  	    algn.x(0) = 12398./sd.wavelength - .001          ; add 2nd pt 1 meV away
  	    algn.x(1) = 12398./sd.wavelength
  		algn.d(0) = x1(0)
 	    algn.d(1) = x1(0)
  	    algn.dn = algn.d
  	    xlock(0) = x1(0)                                  ; save as fiducial x1
  	    pix_siz = fix(1000.*abs(tmp.x(1) - tmp.x(0) ))    ; save for default pixel size
  	    HANDLE_VALUE, Data(5), algn, /set   ; store first X-1 value
  	    f_store = {n:0,f:""}
		HANDLE_VALUE, Data(10), f_store, /set
  	    first_pass=1
  	    PlotBuf,5
  	ENDIF ELSE BEGIN
  		x=[tmp1.x,1]
  		d=[tmp1.d,1]
		npt = n_elements(tmp1.x)
  		x(npt) = 12398./sd.wavelength
  		d(npt) = x1(0)				 ; save x-data value
  		xs=UNIQ(x,sort(x))           ; force the data into regular order
  		xt=x(xs)  &  dt = d(xs)
  		algn = {t:'1d', x:xt, d:dt, xl:'E', dn:d, dl:'X-shift-1'}
  		HANDLE_VALUE, Data(5), algn, /set    ; store updated X-1 values
  	ENDELSE

  	HANDLE_VALUE, Data(6),tmp1       ; get any existing alignment data
  	IF n_tags(tmp1) EQ 0 THEN begin
  	    algn.d(0) = y1(0)               ; FIRST point add a marker  (Ref=0)
  	    algn.d(1) = y1(0)
  	    algn.dn = algn.d
  	    ylock(0)= y1(0)              ; save as fiducial y1
  	    algn.dl = 'Y-shift-1'
  	    HANDLE_VALUE, Data(6), algn, /set   ; store first Y-1 value
  	    PlotBuf,6
	ENDIF ELSE BEGIN
 		x=[tmp1.x,1]
  		d=[tmp1.d,1]
		npt = n_elements(tmp1.x)
  		x(npt) = 12398./sd.wavelength
  		d(npt) = y1(0)   ; save y-data value
  		xs=UNIQ(x,sort(x))
  		xt=x(xs)  &  dt = d(xs)
  		algn = {t:'1d', x:xt, d:dt, xl:'E', dn:d, dl:'Y-shift=1'}
  		HANDLE_VALUE, Data(6), algn, /set   ; store updated Y-1 values
  	ENDELSE

;  SAVE positions of second point unless one_pt flag is set
if NOT keyword_set(one_pt) THEN BEGIN                                                                   ; stretch/shift
  	HANDLE_VALUE, Data(8),tmp1                                                                          ; stretch/shift
  	IF n_tags(tmp1) EQ 0 THEN BEGIN                                                                     ; stretch/shift
  		algn.d(0) = x2(0)                                                                                   ; stretch/shift
 	    algn.d(1) = x2(0)                                                                                   ; stretch/shift
  	    algn.dn = algn.d                                                                                ; stretch/shift
  	    xlock(1)= x2(0)                     ; save as fiducial x2                                       ; stretch/shift
  	    algn.dl = 'X-shift-2'                                                                           ; stretch/shift
  	    HANDLE_VALUE, Data(8), algn, /set   ; store first x-2 value                                     ; stretch/shift
  	    Plotbuf,8                                                                                       ; stretch/shift
  	ENDIF ELSE BEGIN                                                                                    ; stretch/shift
  		x=[tmp1.x,1]                                                                                    ; stretch/shift
  		d=[tmp1.d,1]                                                                                    ; stretch/shift
		npt = n_elements(tmp1.x)                                                                        ; stretch/shift
  		x(npt) = 12398./sd.wavelength                                                                   ; stretch/shift
  		d(npt) = x2(0)              ; save x2 data                                                      ; stretch/shift
  		xs=UNIQ(x,sort(x))                                                                              ; stretch/shift
  		xt=x(xs)  &  dt = d(xs)                                                                         ; stretch/shift
  		algn = {t:'1d', x:xt, d:dt, xl:'E', dn:d, dl:'X-shift-2'}                                       ; stretch/shift
  		HANDLE_VALUE, Data(8), algn, /set    ; store updated x-2 values                                 ; stretch/shift
  	ENDELSE                                                                                             ; stretch/shift

  	HANDLE_VALUE, Data(9),tmp1       ; get any existing alignment data                                  ; stretch/shift
  	IF n_tags(tmp1) EQ 0 THEN begin                                                                     ; stretch/shift
 	    algn.d(0) = y2(0)                                                                                   ; stretch/shift
	    algn.d(1) = y2(0)                                                                                   ; stretch/shift
  	    algn.dn = algn.d                                                                                ; stretch/shift
  	    ylock(1) =  y2(0)                    ; save as fiducial y2                                       ; stretch/shift
  	    algn.dl = 'Y-shift-2'                                                                           ; stretch/shift
 	    HANDLE_VALUE, Data(9), algn, /set   ; store first y-2 value                                     ; stretch/shift
 	    PlotBuf,9                                                                                       ; stretch/shift
	ENDIF ELSE BEGIN                                                                                    ; stretch/shift
 		x=[tmp1.x,1]                                                                                    ; stretch/shift
  		d=[tmp1.d,1]                                                                                    ; stretch/shift
		npt = n_elements(tmp1.x)                                                                        ; stretch/shift
  		x(npt) = 12398./sd.wavelength                                                                   ; stretch/shift
  		d(npt) = y2(0)                                                                             ; stretch/shift
  		xs=UNIQ(x,sort(x))                                                                              ; stretch/shift
  		xt=x(xs)  &  dt = d(xs)                                                                         ; stretch/shift
  		algn = {t:'1d', x:xt, d:dt, xl:'E', dn:d, dl:'Y-shift-2'}                                       ; stretch/shift
  		HANDLE_VALUE, Data(9), algn, /set   ; store updated y-2 values                                  ; stretch/shift
  	ENDELSE                                                                                             ; stretch/shift
                                                                                                        ; stretch/shift
; Convert x, y scales of image to that of first image (stored as xlock(2), ylock(2))
; FIRST VERSION OF THIS was wrong !!! - MUST STRETCH X, Y INDEPENDENTLY !!                              ; stretch/shift
; BUT KEEP THIS "AVERAGE STRETCH" CALCULATION AS A FAST WAY TO TRACK EFFECT
    dnorm = sqrt( (xlock(1)- xlock(0))^2 + (ylock(1) - ylock(0))^2 )                                    ; stretch/shift
    dist  = sqrt((x2(0)- x1(0))^2 + (y2(0) - y1(0))^2 )                                                   ; stretch/shift
	slope = dnorm/dist                                                                                  ; stretch/shift
;	print, 'standard distance, this image, ratio', dnorm,dist, slope                                    ; stretch/shift
                                                                                                        ; stretch/shift
  	HANDLE_VALUE, Data(7),tmp1       ; get any existing stretch data                                    ; stretch/shift
  	IF n_tags(tmp1) EQ 0 THEN begin                                                                     ; stretch/shift
 	    algn.d(0) = slope - 0.001                                                                         ; stretch/shift
	    algn.d(1) = slope                                                                               ; stretch/shift
  	    algn.dn = algn.d                                                                                ; stretch/shift
  	    algn.dl = 'Stretch'                                                                             ; stretch/shift
	    HANDLE_VALUE, Data(7), algn, /set   ; store first slope value                                   ; stretch/shift
	    PlotBuf,7                                                                                       ; stretch/shift
	ENDIF ELSE BEGIN                                                                                    ; stretch/shift
 		x=[tmp1.x,1]                                                                                    ; stretch/shift
  		d=[tmp1.d,1]                                                                                    ; stretch/shift
		npt = n_elements(tmp1.x)                                                                        ; stretch/shift
  		x(npt) = 12398./sd.wavelength                                                                   ; stretch/shift
  		d(npt) = slope                                                                                  ; stretch/shift
  		xs=UNIQ(x,sort(x))                                                                              ; stretch/shift
  		xt=x(xs)  &  dt = d(xs)                                                                         ; stretch/shift
  		algn = {t:'1d', x:xt, d:dt, xl:'E', dn:dt, dl:tmp1.dl}                                          ; stretch/shift
  		HANDLE_VALUE, Data(7), algn, /set   ; store more slope values                                   ; stretch/shift
  	ENDELSE                                                                                             ; stretch/shift
; Report to user what has happened	; stretch/shift
	text = string(slope, xlock(0) - x1(0), ylock(0) - y1(0), $
	     format = "('XY scale:',/,'stretched by ',f7.3,/,'X-shift: ',f7.3,/,'Y-shift ',f7.3)")
   	WIDGET_CONTROL, Uprompt, SET_VALUE=text
 	print, text
;   	    + strmid(strcompress(string(slope)),0,6) + ' X-shift    ' $                                     ; stretch/shift
;   	     +  strmid(strcompress(string(xlock(0) - x1(0))),0,6) $                                         ; stretch/shift
;   	     + ' Y-shift    ' + strmid(strcompress(string(ylock(0) - y1(0))),0,6)                           ; stretch/shift
ENDIF		; end of storage of (X2,Y2, stretch)                                                        ; stretch/shift


; NOW generate the corrected (X,Y) scales based on
; difference of THIS files (X1,Y1,X2,Y2) values and those from the first file (in xlock, ylock)
	xslope = 1.0	& yslope = 1.0						; for one-pt conversions, no stretch
    if NOT keyword_set(one_pt) THEN BEGIN                                                               ; stretch/shift
		xslope = (X1(0) - X2(0)) /(Xlock(0) - Xlock(1))         ; generate slopes for stretching        ; stretch/shift
		yslope = (Y1(0) - Y2(0)) /(Ylock(0) - Ylock(1))                                                 ; stretch/shift
  	endif		                                                                                        ; stretch/shift
  	xint = Xlock(0) - X1(0)
  	yint = Ylock(0) - Y1(0)
; convert (x,y) axes for BOTH shift (s___) and stretch-shift (a____)
	tmp.x = xslope*tmp.x + xint
  	tmp.y = yslope*tmp.y + yint

; DEBUG added to understand where positional error arise (aph 19-may-99)

; --------- store calibrated data ------------
	CurBuf = 2
	tmp.dl = 'C ' + tmp.dl
  	HANDLE_VALUE, Data(CurBuf), tmp, /set    ; store converted image in Buffer 4
  	PlotBuf,CurBuf

; Cut-out constant range
; ------------ truncate to common range *** MUST come before meshing !! ****
	if lock EQ 1 then begin    ; NB ONLY GET THIS INFORMATION FIRST-TIME - identified when
   	    CurBuf = 0            ;  user CLEAR buffers 4,5,6,8,9 to reset image alignment
   		tmp.dl = 'RS ' + tmp.dl
   		HANDLE_VALUE, Data(CurBuf), tmp, /set
   		Label(CurBuf) = tmp.dl
   		PlotBuf,CurBuf			; plot to let user see
		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)
	endif
;Set the RANGE of (x,y)-space to be extracted
	!x.range(0) = cll(0) & !x.range(1) = cur(0)
	!y.range(0) = cll(1) & !y.range(1) = cur(1)
	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 = create_struct('t','2d','x',xn,'y',yn,'d',dn, 'e', tmp.e, $
			'xl',tmp.xl,'yl',tmp.yl,'dl','c' + tmp.dl)

; --------- store truncated data ------------
	CurBuf = 3
	tmp.dl = tmp.dl + ' T '
  	HANDLE_VALUE, Data(CurBuf), tmp, /set    ; store converted image in Buffer 4
  	PlotBuf,CurBuf

; interpolate to specified pixel size
	if lock EQ 1 then begin
		print, 'cut out ',n_Elements(xn), ' x ',n_Elements(yn),' pixels'
		pix_siz = get_num(Prompt='Pixels (nm)',val=fix(pix_siz), group = axis_id)
	endif
	nx = fix( abs( 1000.*(cur(0) - cll(0))) /pix_siz )
	ny = fix( abs( 1000.*(cur(1) - cll(1))) /pix_siz )
	print, 'Congrid data to ', nx, ny
	xn = congrid(xn,nx,/interp)
	yn = congrid(yn,ny,/interp)
	dn = congrid(dn,nx,ny,/interp)
	tmp.dl = tmp.dl + ' cut '
	tmp = create_struct('t','2d','x',xn,'y',yn,'d',dn, 'e', energy, $
			'xl',tmp.xl,'yl',tmp.yl,'dl',tmp.dl)
	CurBuf = 0
   	HANDLE_VALUE, Data(CurBuf), tmp, /set
   	Label(CurBuf) = tmp.dl
   	PlotBuf,CurBuf
	x = tmp.x  &  y = tmp.y  & d = tmp.d
	n_cols = n_elements(x)
	n_rows = n_elements(y)

; ---------- bin to improve statistics if desired
	if lock EQ 1 THEN BEGIN
		bintst=get_num(prompt='Bin *',val=2, group = axis_id)
		factor(1) = bintst(0)
	endif
	bin = fix(factor(1))      ; use common block parameter factor to store bin amount
	IF bin NE 1 then begin
      if (float(n_cols)/float(bin))-fix(n_cols/bin) GT 0  then begin
        tmp1 = tmp.d                ;force size to integer mutiple of bin
        nct = fix(n_cols/bin)*bin  &  nrt = fix(n_rows/bin)*bin
        xtra_c = n_cols - nct + 1
        xtra_r = n_rows - nrt + 1
        if lock EQ 1 THEN print, 'truncate image to ', fix(nct), fix(nrt)
        x = intarr(nct)
        y = intarr(nrt)
        d = intarr(nct,nrt)
        x = tmp.x(0:n_cols-xtra_c)
        y = tmp.y(0:n_rows-xtra_r)
		d = tmp1(0:n_cols-xtra_c,0:n_rows-xtra_r)
      endif

   	  n_cols=fix(n_cols/bin)
      n_rows=fix(n_rows/bin)
      if lock EQ 1 THEN print, 'Data binned to ',n_cols,'  x', n_rows
	  x = congrid(x,n_cols)
	  y = congrid(y,n_rows)
      d = rebin(d,n_cols,n_rows)
      tmp.dl = tmp.dl + ' BN '
    endif

; -------- generate structure with cut data and correct X,Y scales
 	tmp = {t:'2d',x:x,y:y,d:d,e:energy, xl:tmp.xl,yl:tmp.yl,dl:tmp.dl}

; ------ OUTPUT aligned file as  *.axb  [  use *.nc format ONLY if *.nc' input ]

 	fnm = 'a' + fnm                                      ; use 1st character to define
 	if keyword_set(one_pt) THEN begin				; type of alignment  (s = 1 point (shift), t = 2 points)
 		strput, fnm,'s'
 	endif else strput, fnm,'t'
;	if first_pass EQ 1 then begin
;		file_nc=get_text(val=fnm, Prompt='Filename',group = axis_ID)
;	file=file_nc(0) + '.nc'
;	endif
	file=fnm + '.axb'
	if !VERSION.OS_FAMILY EQ 'unix' then sep = '/' else sep = '\'
	print, 'img_align: Path = ', WritePath, ' file name = ', file
		file = WritePath + sep + file            ; use WritePath  (in axis_com.mon) to define path
	t = axb_save(tmp, file=file, /silent)
;	write_file = sav_nsls(tmp,file=file)		to write *.nc

	HANDLE_VALUE, Data(10), f              ; save file name
   	f.n=f.n+1	; increment number of files
   	tmpf=[f.f,strcompress(string(file))]
   	if first_pass EQ 1 then  tmpf = [file]
   	f_new={n: f.n, f:tmpf}
   	HANDLE_VALUE, Data(10), f_new, /set
   	tmp.dl = tmp.dl + ' align'
   	CurBuf=4
  	HANDLE_VALUE, Data(CurBuf), tmp, /set    ; store converted image in Buffer 4
  	PlotBuf,CurBuf
  	close,/all
END
