; Copyright (c) 1998-2015  A.P. Hitchcock  All rights reserved
;+
;NAME:
;  STACK_ROTATE
;
;LAST CHANGED: ----------------------------------- 22-Mar-15
;
; PURPOSE:
;	This procedure rotates each image of a stack
;
; CATEGORY:
;	stack processing; stand alone operation
;
; CALLING SEQUENCE:
;	STACK_ROTATE, FILE=FILE
;
; INPUTS: none required
;
; KEYWORDS:
; 	FILE	name of the output file
;
; OUTPUTS: 	differentiated stack
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; stack_process_com
; BSIF_com
; volume_data, image_stack
;
; SIDE EFFECTS: none
;
; RESTRICTIONS: none
;
; MODIFICATION HISTORY:
; (12-Jan-09 aph) first version
; (22-Mar-15 aph) -  (x,y) axes / scales of rotated stack written correctly
;                    alert user to stacks with non-square pixels; use avg. value
;-

PRO stack_rotate, file=file, angle = angle, pixelsize = pixelsize, maxDist = maxDist, minval = minval
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com

on_error,2

; 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)

fltr = '*.ncb'

IF not keyword_set(file) then begin
	file = pickfile2(/READ, TITLE = 'STACK 1', FILTER=fltr, /LPATH, DEFPATH=defpath)
	if strlen(file) LE 0 THEN return  ; bail-out if no filename
ENDIF
t=ax_name(file)
file_short=t(1)

stack_rb, file
stack = image_stack
t1 = size(image_stack)
n_cols = t1(1)
n_rows = t1(2)
n_imgs = t1(3)
print, 'stack parameters: n_cols, n_rows, n_imgs ', n_cols, n_rows, n_imgs


ax_wait		; turn on hourglass- this could be slow

; generate average of stack to define the rotation parameters from average  ; NOT YET WRITTEN
;for i=0, n_imgs-1 do begin
;	tmp = moment(image_stack(*,*,i))
;endfor

img_0 = image_stack(*,*,0)

pix_siz_x = (x_stop - x_start)/n_cols
print, pix_siz_x
pix_siz_y = (y_stop - y_start)/n_rows
pix_avg = (pix_siz_x + pix_siz_y)/2.
axis_log, 'x-pixel_size (nm) = ' + string(1000.*pix_siz_x, format='(F6.2)')
axis_log, 'y-pixel_size (nm) = ' + string(1000.*pix_siz_y, format='(F6.2)')
axis_log, 'RATIO x/y = ' + string(pix_siz_x/pix_siz_y, format='(F6.2)')
axis_log, '  average pixel_size = ' + string(1000*pix_avg, format='(F6.2)')

; --- generate standard 2d Axis structure to allow mimic of ax_rotate
; pix_siz = (x_stop - x_start)/n_cols  ; ASSUMES square pixels !!

; ----------- determine pixel size (square)
if NOT keyword_set(pixelsize) then begin
	t_size = 1000*pix_avg
	if axis_on then t_siz = get_num(Prompt = 'pixel size (nm)', val = t_size, group = axis_ID) $
		else t_siz = get_num(Prompt = 'pixel size (nm)', val = t_size)
	pix_siz = t_size/1000.
endif else pix_siz = pixelsize

x = findgen(n_cols)*pix_siz
y = findgen(n_rows)*pix_siz
tmp = {t:'2d', x:x, y:y, d:img_0, e: ev(0), dl: file_short}
print, 'energy = ', tmp.e


; ------------ determine rotation angle
if NOT keyword_set(angle) then begin
	if axis_on then angle = get_num(Prompt = 'Angle (deg)', Val = ax_angle, group = AXIS_ID) $
		else  angle = get_num(Prompt = 'Angle (deg)', Val = ax_angle)
endif
ax_angle = angle	; save for later use

; ----- force rotation about centre of the stack
xc = (x_stop - x_start) / 2
yc = (y_stop - y_start) / 2
centre = [xc,yc]
ax_centre = [xc,yc]		; store for future use
ind = Dindex(xc,yc,tmp)


; ----------  extend image to map max possible field in which to rotate
if NOT keyword_set(maxDist) then begin
	img_ll = [tmp.x(0),tmp.y(0)]
	img_lr = [tmp.x(n_elements(tmp.x)-1),tmp.y(0)]
	img_ul = [tmp.x(0), tmp.y(n_elements(tmp.y)-1)]
	img_ur = [tmp.x(n_elements(tmp.x)-1), tmp.y(n_elements(tmp.y)-1)]
	del_ll = distance(centre,img_ll)
	del_lr = distance(centre,img_lr)
	del_ul = distance(centre,img_ul)
	del_ur = distance(centre,img_ur)
	maxDist = max([del_ll, del_lr, del_ul, del_ur])
	if axis_on then  maxDist = get_num(Prompt = 'maxDist from rotation centre', val = maxDist, group = axis_ID) $
	 	else maxDist = get_num(Prompt = 'maxDist from rotation centre', val = maxDist)
endif


npix = 2*maxDist/pix_siz+1


; ----------- determine image value for out-of-bound regions
if not keyword_set(minval) then val_extend = min(image_stack) else val_extend = minval

print, 'Center (x,y) = ', centre
print, 'maxDist = ', maxDist
print, 'Pixel size (nm) = ', 1000.*pix_siz
print, 'fill boundaries with ', val_extend

; rotate each image in stack
image_Stack2 = make_array(npix, npix, n_imgs, /float, value=val_extend)
img_x = -1.*maxDist + pix_siz*findgen(npix)
img_y = -1.*maxDist + pix_siz*findgen(npix)

for i = 0, n_imgs-1 do begin
	img_d = make_array(npix, npix, /float, value=val_extend)
	n_xoff = abs(fix((maxDist - abs(xc - tmp.x(0)))/pix_siz))
	n_yoff = abs(fix((maxDist - abs(yc - tmp.y(0)))/pix_siz))
	img_d(n_xoff: n_xoff+n_elements(tmp.x)-1, n_yoff:n_yoff+n_elements(tmp.y)-1) = image_stack(*,*,i)

	R_img = ROT(img_d, angle,1.0, ind(0)+n_xoff, ind(1)+n_yoff, /INTERP, $
	     missing = val_extend, cubic = -0.5)    ; angle in DEGREE CLOCKWISE

	image_Stack2(*,*,i) = R_img
endfor

image_stack = image_stack2

ax_wait,/off		; turn off hourglass (sometimes sticks on)

; --------  change the values defining the X and Y_x( axis lengths, since have expanded the space
x_start = img_x(0)     &   y_start=img_y(0)
x_stop = img_x(npix-1)   &   y_stop = img_y(npix-1)
n_cols = npix  & nrows = npix

; ----- write out result
file = pickfile2(/READ, TITLE = 'name of rotated stack', FILTER=fltr, /LPATH, DEFPATH=defpath)
if strlen(file) NE 0 then begin
	t = ax_name(file)
	filename = string(t(0)) + string(t(1)) + '.ncb'
	stack_wb, filename
endif

return
end

