; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	AX_IMAGE_APPEND
;
;LAST CHANGED:-------------------------- 10-Dec-11
;
;PURPOSE
;  This procedure appends the currently displayed image
;  to a user-selected second image with averaging in the common region.
;  The final image is interpolated to the mesh of the more finely sampled
;   image, unless this will result in an array with any dimension > 3000
;
;CATEGORY:
;	AXIS: image analysis
;
;CALLING SEQUENCE:
; 	ax_image_append
;
;INPUT: none (uses AXIS_COM common blocks for input/output)
;
;KEYWORDS: none
;
;OUTPUT:
;	(Image in buffer 0) = (finer pixel image) appended to (coarser pixel image)
;
;COMMON BLOCKS:
;	AXIS_COM	standard set of common blocks
;
;MODIFICATION HISTORY:
; (29-jun-00 aph) first version
; (03-feb-11 aph) ensure energy is defined on output; default to energy of first image
; (11-Dec-11 aph) make more reliable by changing algorithm & introducing meshing to same pixel size
;-

PRO AX_IMAGE_APPEND
@axis_com
on_error,2

HANDLE_VALUE, Data(CurBuf), tmp
IF n_tags(tmp) EQ 0 THEN RETURN
IF tmp.t EQ '2d' then begin
    bt = buffer(group = AXIS_ID)
    IF bt EQ -1 THEN RETURN				  ; abort procedure if no other image selected
  	HANDLE_VALUE, Data(bt), div
  	IF n_tags(div) EQ 0 THEN RETURN
	IF (div.t EQ '2d') THEN BEGIN     ; for images only
			ax_wait
			xmin1 = min(tmp.x, max=xmax1)
			xmin2 = min(div.x, max=xmax2)
			xmin = min([xmin1,xmin2])
			xmax = max([xmax1,xmax2])
			ymin1 = min(tmp.y, max=ymax1)
			ymin2 = min(div.y,max=ymax2)
			ymin = min([ymin1,ymin2])
			ymax = max([ymax1,ymax2])
; ------- check step size and adjust each image to finest step size
			xstep = tmp.x(2)-tmp.x(1)		; ASSUMES STEP IS THE SAME !
			ystep = tmp.y(2)-tmp.y(1)
			xstep2 = div.x(2)-div.x(1)
			ystep2 = div.y(2)-div.y(1)
			step = float(min([xstep, xstep2, ystep, ystep2]))

; ------ set up template image that is larger that either image; force maxmum of 3000
			nx = ceil(2 + (xmax - xmin)/step)
			if nx GT 3000 then begin
				axis_log, 'appended image limited to 3000 x-pixels'
				step = (xmax - xmin)/3000.
				nx = ceil(2 + (xmax - xmin)/step)
			endif
			ny = ceil(2 + (ymax - ymin)/step)
			if ny GT 3000 then begin
				axis_log, 'appended image limited to 3000 y-pixels'
				step = (ymax - ymin)/3000.
				nx = ceil(2 + (xmax - xmin)/step)
				ny = ceil(2 + (ymax - ymin)/step)
			endif
			x = findgen(nx)*step + xmin-step
			y = findgen(ny)*step + ymin-step
			print, 'Append image size is ', nx, ny
			d_avg = (mean(tmp.d) + mean(div.d))/2.
			d = make_array(nx, ny, /float, value=d_avg)
			dl = string(format='(A," AP ",A)',tmp.dl,div.dl)
;-------- set energy to energy of first image
			if n_elements(tmp.e) GT 0 then energy = tmp.e else energy = 0.
	 	   	s = {t:'2d', x:x, y:y, d:d, e: energy, xl:tmp.xl, yl:tmp.yl, dl:dl}

			print, 'AX_IMAGE_APPEND: mesh both images to finest step (um): ', step
			tmp_mesh = intarr(2)
			tmp_mesh(0)=fix((tmp.x(n_elements(tmp.x)-1)-tmp.x(0))/step)
			tmp_mesh(1)=fix((tmp.y(n_elements(tmp.y)-1)-tmp.y(0))/step)
			div_mesh = intarr(2)
			div_mesh(0)=fix((div.x(n_elements(div.x)-1)-div.x(0))/step)
			div_mesh(1)=fix((div.y(n_elements(div.y)-1)-div.y(0))/step)
;			print, 'tmp_mesh ', tmp_mesh
;			print, 'div_mesh ', div_mesh
			tmp = ax_interp(tmp, mesh=tmp_mesh)
			div = ax_interp(div, mesh=div_mesh)

; ------ insert tmp.d and div.d into the template (d.d) placing the finest one last
			if xstep GT xstep2 then begin
				Ilo = DIndex(xmin1,ymin1,s)
				d(Ilo(0):Ilo(0)+n_elements(tmp.x)-1,Ilo(1):Ilo(1)+n_elements(tmp.y)-1) = tmp.d
				Ilo = DIndex(xmin2,ymin2,s)
				d(Ilo(0):Ilo(0)+n_elements(div.x)-1,Ilo(1):Ilo(1)+n_elements(div.y)-1) = div.d
			endif else begin
				Ilo = DIndex(xmin2,ymin2,s)
				d(Ilo(0):Ilo(0)+n_elements(div.x)-1,Ilo(1):Ilo(1)+n_elements(div.y)-1) = div.d
				Ilo = DIndex(xmin1,ymin1,s)
				d(Ilo(0):Ilo(0)+n_elements(tmp.x)-1,Ilo(1):Ilo(1)+n_elements(tmp.y)-1) = tmp.d
			endelse

; ---------- remove first and last rows & columns (made template 2 larger to handle rounding errors)
			tx = x(1:nx-2)
			ty = y(1:ny-2)
			td = d(1:nx-2,1:ny-2)
			s = {t:'2d', x:tx, y:ty, d:td, e: energy, xl:tmp.xl, yl:tmp.yl, dl:dl}

; ---------- store appended image in buffer 0
     	   	CurBuf = 0
   	   		HANDLE_VALUE, Data(CurBuf), s, /set
   	   		Label(CurBuf) = s.dl
   	   		PlotBuf,CurBuf
   	   		ax_wait, /off
   	 ENDIF ELSE WIDGET_CONTROL, Uprompt, SET_VALUE='Append buffers: only for 2 images'
ENDIF ELSE WIDGET_CONTROL, Uprompt, SET_VALUE='Append buffers: only for images'
END

; *********************************************************
; --------- OLD CODE (pre-Dec-11) - save in case there is a fatal flaw with new way
; *********************************************************

;			goto, skip_mesh

;			if (xstep-xstep2) NE 0 then begin
;			    axis_log, 'Different X-step sizes ' + string(xstep) + string(xstep2)
; ------------ mesh larger pixel to smaller pixel size
;				if xstep LT xstep2 then begin	; tmp.x finer mesh than div.x - remesh div.x
;					mesh = intarr(2)
;					mesh(0)=round(float(n_elements(div.x))*xstep/xstep2)
;					if mesh(0) GT 2500 then mesh(0) = 2500	; cap size to 2500 in case if large mismatch
;					mesh(1)=n_elements(div.y)
;					div = ax_interp(div, mesh=mesh)
;				endif else begin				; div.x finer mesh than tmp.x - remesh tmp.x
;					mesh = intarr(2)
;					mesh(0)=round(float(n_elements(tmp.x))*xstep2/xstep)
;;					if mesh(0) GT 2500 then mesh(0) = 2500	; cap size to 2500 in case if large mismatch
;					mesh(1)=n_elements(tmp.y)
;					tmp = ax_interp(tmp, mesh=mesh)
;				endelse
;			endif
;			ystep = tmp.y(2)-tmp.y(1)
;			ystep2 = div.y(2)-div.y(1)
;			if (ystep-ystep2) NE 0 then begin	 ; tmp.y finer mesh than div.y - remesh div.y
;				axis_log, 'Different Y-step sizes ' + string(ystep) + string(ystep2)
; ------------ mesh larger pixel to smaller pixel size
;				if ystep LT ystep2 then begin
;					mesh = intarr(2)
;					mesh(0)=n_elements(div.x)
;					mesh(1)=round(float(n_elements(div.y))*ystep/ystep2)
;					if mesh(1) GT 2500 then mesh(1) = 2500	; cap size to 2500 in case if large mismatch
;					div = ax_interp(div, mesh=mesh)
;				endif else begin		; div.y finer mesh than tmp.y - remesh tmp.y
;					mesh = intarr(2)
;					mesh(0)= n_elements(tmp.x)
;					mesh(1)=round(float(n_elements(tmp.y))*ystep2/ystep)
;					if mesh(1) GT 2500 then mesh(1) = 2500	; cap size to 2500 in case if large mismatch
;					tmp = ax_interp(tmp, mesh=mesh)
;				endelse
;			endif
;		skip_mesh:


; ------------- old way (pre Dec-11)
; identify indices in d which correspond to tmp and div
;			Ilo = DIndex(xmin1,ymin1,s)
;			Ihi = DIndex(xmax1,ymax1,s)
;			print, 'Append tmp.d into d matrix index values of:  X:', Ilo(0), Ihi(0),' Y:',Ilo(1), Ihi(1)
;			help, d
;			ta_s = size(tmp.d)
;			if Ihi(0)-Ilo(0) GE ta_s(1) then ta_x =(Ihi(0)-Ilo(0)) - ta_s(1)+1  else ta_x=0
;			if Ihi(1)-Ilo(1) GE ta_s(2) then ta_y =(Ihi(1)-Ilo(1)) - ta_s(2)+1  else ta_y=0
;			print, 'tmp matrix index range: X: 0 ', Ihi(0)-Ilo(0)-ta_x,' Y: 0 ', Ihi(1)-Ilo(1)-ta_y
;			ta = tmp.d(0:(Ihi(0)-Ilo(0)-ta_x),0:(Ihi(1)-Ilo(1)-ta_y))
;			help, ta
;			d(Ilo(0):Ihi(0)-ta_x,Ilo(1):Ihi(1)-ta_y) = ta

;			Ilo = DIndex(xmin2,ymin2,s)
;			Ihi = DIndex(xmax2,ymax2,s)
;			ta_s = size(div.d)
;			if Ihi(0)-Ilo(0) GE ta_s(1) then ta_x =(Ihi(0)-Ilo(0)) - ta_s(1)+1  else ta_x=0
;			if Ihi(1)-Ilo(1) GE ta_s(2) then ta_y =(Ihi(1)-Ilo(1)) - ta_s(2)+1  else ta_y=0
;			print, 'Append div.d into:index values of  X:', Ilo(0), Ihi(0),' Y:',Ilo(1), Ihi(1)
;			print, 'div matrix index range: X: 0 ', Ihi(0)-Ilo(0)-ta_x,' Y: 0 ', Ihi(1)-Ilo(1)-ta_y
;			ta = tmp.d(0:(Ihi(0)-Ilo(0)-ta_x),0:(Ihi(1)-Ilo(1)-ta_y))
;			d(Ilo(0):Ihi(0)-ta_x, Ilo(1):Ihi(1)-ta_y-1) = ta

; use average of the common section if needed
;	  		t = Ax_Mesh(tmp,div,2)  	; NB tmp and div are returned with common area only
;			if t EQ 1 then begin
;				dn = (tmp.d + div.d)/2.
;				xmin = min(tmp.x, max=xmax)
;				ymin = min(tmp.y, max=ymax)
;				Ilo = DIndex(xmin,ymin,s)
;				d(Ilo(0):Ilo(0)+n_elements(tmp.x)-1,Ilo(1):Ilo(1)+n_elements(tmp.y)-1) = dn
;			endif
; ------------- old way (pre Dec-11)
;				Ihi = DIndex(xmax,ymax,s)
;				if Ihi(0)-Ilo(0) GE ta_s(1) then ta_x =(Ihi(0)-Ilo(0)) - ta_s(1)+1  else ta_x=0
;				if Ihi(1)-Ilo(1) GE ta_s(2) then ta_y =(Ihi(1)-Ilo(1)) - ta_s(2)+1  else ta_y=0
;				print, 'Append common into: X:', Ilo(0), Ihi(0),' Y:',Ilo(1), Ihi(1)
;				ta = tmp.d(0:(Ihi(0)-Ilo(0)-ta_x),0:(Ihi(1)-Ilo(1)-ta_y))
;				d(Ilo(0):Ihi(0)-ta_x,Ilo(1):Ihi(1)-ta_y) = ta
