; Copyright (c) 1998-2012  A.P. Hitchcock  All rights reserved
;+
;NAME:
;  STACK_ROTATE
;
;LAST CHANGED: ----------------------------------- 26-Nov-12
;
; 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
;-

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)

; --- generate standard 2d Axis structure to allow mimic of ax_rotate
pix_siz = (x_stop - x_start)/n_cols  ; ASSUMES square pixels !!
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)

; ----------- determine pixel size (square)
if NOT keyword_set(pixelsize) then begin
	if axis_on then pix_siz = get_num(Prompt = 'pixel size (um)', val = pix_siz, group = axis_ID) $
		else pix_siz = get_num(Prompt = 'pixel size (um)', val = pix_siz)
endif else pix_siz = pixelsize

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

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

