PRO rigid_al,datcub,datal,edges,CUMUL=cumul
;
; RIGID ALIGNMENT, where ALIGN_C2 is used to calculate shifts
; and single-frame reference or cumulative method can be elected.
;
; INPUTS: datcub = input data cube
;		   correlation window vertices (in the program)
;
; OUTPUTS: datal = aligned data cube
;          edges = array [X_dep,X_fin,Y_dep,Y_fin] - final FOV limits
;                   (edges lost due to image shifts may be removed)
;
; KEWORDS: cumul - if set, cumulative method is used;
;                  otherwise (default), the middle frame of the series
;                  is taken as the reference
;
; CALLED PROGS: acre, align_c2, fivepoin, maxloc, mean, morphing_image,
;               correl_optimize and its subroutines
;
;Le 26 Mars 2005
;auteurs Thierry Roudier OMP/CNRS, Michal Sobotka, AsU Ondrejov
;modification 18 Jan 2010, Michal
;
on_error,1

si=size(datcub)

;----------------------------------------------------------------------
; correlation window vertices to be set here:
x1=2   ;10
x2=si(1)-3   ;265
y1=2
y2=si(2)-3

;----------------------------------------------------------------------
; rigid alignment parameters
wx=x2-x1+1-130	; 66: this allows the shifts up to +/- 32 pixels
wy=y2-y1+1-130	;  (see comments in ALIGN_C2)
;----------------------------------------------------------------------

print,'now calculating rigid-alignment shifts'

dx=fltarr(si(3))      ; arrays to store shifts
dy=dx                 ;  (subpixel accuracy)

if keyword_set(cumul) then begin

; cumulative method -------------------------------------------------
cdx=0.                ; accumulators
cdy=0.

for j=1,si(3)-1 do begin

  a=datcub(x1:x2,y1:y2,j-1)     ; reference subimage
  b=datcub(x1:x2,y1:y2,j)       ; correlated subimage

  offset=align_c2(a,b,flg,wx,wy)	; subpixel accuracy
  xoff=offset(0)
  yoff=offset(1)

;  when align_c2 fails (flg=1) we call correl_optimize...
   if (flg eq 1) then begin
       print,j,'    correl_optimize called...'
       correl_optimize,a,b, XOFF_INIT=0 ,YOFF_INIT=0, $
       xoff,yoff, MAGN=6,/NUMPIX
   endif

  cdx=cdx+xoff                  ; cumulative shift in x
  cdy=cdy+yoff                  ; cumulative shift in y
  dx(j)=cdx
  dy(j)=cdy

endfor

endif else begin

; single-frame reference method --------------------------------------
nr=si(3)/2                   ; middle of the series
a=datcub(x1:x2,y1:y2,nr)     ; reference subimage

for j=0,si(3)-1 do begin

;  print,j
  b=datcub(x1:x2,y1:y2,j)       ; correlated subimage

  offset=align_c2(a,b,flg,wx,wy)	; subpixel accuracy
  xoff=offset(0)
  yoff=offset(1)

;  when align_c2 fails (flg=1) we call correl_optimize...
   if (flg eq 1) then begin
       print,j,'    correl_optimize called...'
       correl_optimize,a,b, XOFF_INIT=0 ,YOFF_INIT=0, $
       xoff,yoff, MAGN=6,/NUMPIX
   endif

  dx(j)=xoff
  dy(j)=yoff

endfor

endelse   ; -----------------------------------------------------------

print,'now shifting images'

datal=datcub*0.

for j=0,si(3)-1 do begin
;  datal(0,0,j)=shift(datcub(*,*,j),round(dx(j)),round(dy(j)))
  morphing_image,datcub(*,*,j),im,dx(j),dy(j)	;subpixel shift
  datal(0,0,j)=im
endfor


; Allows to extract the common field of view after alignment.
; inputs are cumulative dx and dy

nx=si(1)                         ; size of the full FOV
ny=si(2)

midx=min(dx)<0                   ; maximum/minimum shifts
madx=max(dx)>0
midy=min(dy)<0
mady=max(dy)>0

edges=intarr(4)

edges(0)=round(madx)            ; edges of the common FOV
edges(1)=round(nx-1+midx)
edges(2)=round(mady)
edges(3)=round(ny-1+midy)
print,'edges x1,x2,y1,y2:',edges

END