; Copyright (c) 1998-2000 A.P. Hitchcock  All rights reserved 
;+ 
;NAME: 
;	ALIGN 
; 
;LAST CHANGED: ----------------------------------- 31-dec-99 
; 
; PURPOSE: 
;	This procedure is the image alignment part of STACK_ANALYZE. 
; It uses 2d fourier transform correlation techniques 
; to shift one image to the same pixel scale as a reference image. 
; Numerous options to control the alignment procedure are provided 
; 
; CATEGORY: 
;	AXIS: stack analysis 
; 
; CALLING SEQUENCE: 
; align,image1,image2,xshift,yshift,image2_shifted,$ 
;    [ sobel=sobel,roberts=roberts,edgegauss=edgegauss,$ 
;    cm=cm,maxshift=maxshift,$ 
;    meanfill=meanfill,medianfill=medianfill,$ 
;    xcorimg_win=xcorimg_win,xcorimg_zoom=xcorimg_zoom,$ 
;    xplot_win=xplot_win,yplot_win=yplot_win,$ 
;    debug=debug,help=help ] 
; 
; INPUTS: 
;	IMAGE1  - reference image 
;	IMAGE2  - image to be shifted to reference image 
; 
; KEYWORDS: 
;	SOBEL	- use sobel edge sharpening 
;	ROBERTS	- use roberts edge sharpening 
;	EDGEGAUSS = # - sigma value for gaussian edge sharpening (default = 4) 
;	CM	- align on centre-of-mass of correlation FFT (default = peak) 
;	MAXSHIFT - maximum number of pixels to shift 
;	MEANFILL - use mean for image pixels not in shifted image 
;	MEDIANFILL - use median for image pixels not in shifted image 
;			(default is to leave those pixels as zero) 
;	XCORING_WIN - number of window to display correlation image 
;	XCORIMG_ZOOM - zoom factor for alignment images 
;	XPLOT_WIN - number of window to display x-shifts 
;	YPLOT_WIN - number of window to display y-shifts 
;	DEBUG - debug code implemented 
;	HELP - display help message 
; 
; OUTPUTS: 
;	XSHIFT 	- value of x-shift (in fractional pixels) 
;	YSHIFT - value of y-shift (in fractional pixels) 
;	IMAGE2_SHIFTED - shifted second image 
; 
; COMMON BLOCKS: 
;   none 
; 
; EXAMPLE: 
;	Call from stack_align.pro : 
;            align,ref_image,this_image,this_x_shift,this_y_shift,$ 
;              shifted_image, sobel = sobel, roberts=roberts, $ 
;              maxshift=maxshift, /meanfill, $ 
;              xcorimg_win = stack_align_par.diff_image_win, $ 
;              xcorimg_zoom = img_zoom 
; 
; CAUTIONS 
;	Image alignment will often fail in cases of ambigous or low contrast 
; 
; MODIFICATION HISTORY: 
; (25-mar-1998 cjj) modified to deal with 24 bit color 
;   If you provide the white, red, and green colors, then 
;   it is assumed that the color table is preloaded. 
;   "white" is an index for a color table from 0 to 255, 
;   whereas "red" and "green" are plot color indices 
;   which might be longs on 24 bit graphics devices.  CJJ 
; (11-jul-98 cjj) Significantly modified 
; Eliminated the correl option, and made /peak the default (i.e., you 
; cannot specify it anymore).  Fit the three pixels about the 
; peak to z=a+bx*cx^2 and z=a+by+cy^2 for subpixel location 
; (30-dec-99 aph) documented as part of AXIS 
;- 
 
PRO align,image1,image2,xshift,yshift,image2_shifted,$ 
    sobel=sobel,roberts=roberts,edgegauss=edgegauss,$ 
    cm=cm,maxshift=maxshift,$ 
    meanfill=meanfill,medianfill=medianfill,$ 
    xcorimg_win=xcorimg_win,xcorimg_zoom=xcorimg_zoom,$ 
    xplot_win=xplot_win,yplot_win=yplot_win,$ 
    debug=debug,help=help 
 
IF (keyword_set(help) OR (n_params() eq 0)) THEN BEGIN 
    print,'align,image1,image2,xshift,yshift,image2_shifted' 
    print,'  Modifiers to input images: /sobel or /roberts, edgegauss_sigma=' 
    print,'    (if you simply do /edgegauss, you will get edgegauss_sigma=4.)' 
    print,'  Constraints on peak finding: /cm, maxshift=' 
    print,'  Constraints on shifted image: /meanfill or /medianfill' 
    print,'  Display the correlation: xcorimg_win=, xcorimg_zoom=' 
    return 
ENDIF 
 
IF (n_elements(sobel) EQ 0) THEN sobel = 0 
IF (n_elements(roberts) EQ 0) THEN roberts = 0 
IF (n_elements(edgegauss) EQ 0) THEN BEGIN 
    edgegauss_sigma = 0. 
ENDIF ELSE IF (float(edgegauss) LE 1.) THEN BEGIN 
    edgegauss_sigma = 4. 
ENDIF ELSE BEGIN 
    edgegauss_sigma = edgegauss 
ENDELSE 
IF (n_elements(cm) EQ 0) THEN cm = 0 
IF (n_elements(meanfill) EQ 0) THEN meanfill = 0 
IF (n_elements(medianfill) EQ 0) THEN medianfill = 0 
IF (n_elements(xcorimg_zoom) EQ 0) THEN xcorimg_zoom = 1 
IF (xcorimg_zoom LT 1) THEN xcorimg_zoom = 1 
 
svec=size(image1) 
IF (svec(0) NE 2) THEN BEGIN 
    print,'image1 has '+strtrim(string(svec(0)),2)+' rather than 2 dimensions' 
    return 
ENDIF 
nx1=svec(1) 
ny1=svec(2) 
 
IF (n_elements(maxshift) EQ 0) THEN maxshift = fix(min(0.8*[nx1/2,ny1/2])) 
 
svec=size(image2) 
IF (svec(0) NE 2) THEN BEGIN 
    print,'image2 has '+strtrim(string(svec(0)),2)+' rather than 2 dimensions' 
    return 
ENDIF 
nx2=svec(1) 
ny2=svec(2) 
 
IF ((nx1 NE nx2) OR (ny1 NE ny2)) THEN BEGIN 
    print,'Images must match in size' 
    return 
ENDIF 
nx=nx1 
ny=ny1 
xcenter = 0 
ycenter = 0 
 
;; We make our own copy of the input image and call it fft in anticipation 
;; of what we will do to it later. 
IF (sobel NE 0) THEN BEGIN 
    fft1 = sobel(image1) 
    fft2 = sobel(image2) 
ENDIF ELSE IF (roberts NE 0) THEN BEGIN 
    fft1 = roberts(image1) 
    fft2 = roberts(image2) 
ENDIF ELSE BEGIN 
    fft1 = image1 
    fft2 = image2 
ENDELSE 
 
IF (edgegauss_sigma NE 0.) THEN BEGIN 
    edgegauss,fft1,edgegauss_sigma,dc1,/zero_edge 
    edgegauss,fft2,edgegauss_sigma,dc2,/zero_edge 
ENDIF 
 
fft1 = shift(temporary(fft1),nx/2,ny/2) 
fft1 = fft(fft1,-1,/overwrite) 
fft1 = shift(temporary(fft1),nx/2,ny/2) 
 
fft2 = shift(temporary(fft2),nx/2,ny/2) 
fft2 = fft(fft2,-1,/overwrite) 
fft2 = shift(temporary(fft2),nx/2,ny/2) 
 
fft1 = conj(temporary(fft1)) 
fft2 = fft1*temporary(fft2) 
fft1 = 0 
 
fft2 = shift(temporary(fft2),nx/2,ny/2) 
fft2 = fft(fft2,1,/overwrite) 
fft2 = shift(temporary(fft2),nx/2,ny/2) 
fft2 = abs(temporary(fft2)) 
 
;; Find the maximum 
xleft = (nx/2-maxshift)>0 
xright = (xleft+2*maxshift)<(nx-1) 
ybot = (ny/2-maxshift)>0 
ytop = (ybot+2*maxshift)<(ny-1) 
dummy = max(fft2(xleft:xright,ybot:ytop),index) 
ycenter = ybot+(index / (xright-xleft+1)) 
xcenter = xleft+(index mod (xright-xleft+1)) 
x_integer_center = xcenter 
y_integer_center = ycenter 
 
IF keyword_set(debug) THEN BEGIN 
    print,'Maximum value is at ['+$ 
      strtrim(string(xcenter),2)+','+$ 
      strtrim(string(ycenter),2)+']' 
ENDIF 
 
IF (cm NE 0) THEN BEGIN 
    max_fft2 = max(fft2,min=min_fft2) 
    cm_scale=1./(max_fft2-min_fft2) 
    cm_threshold=0.5 
    cm_image=fft2(xleft:xright,ybot:ytop) 
    cm_indices=where(((cm_image-min_fft2)*cm_scale) GE cm_threshold,n_cm) 
    IF (n_cm GE 4) THEN BEGIN 
        xpositions=xleft+(cm_indices mod (xright-xleft+1)) 
        ypositions=ybot+(cm_indices/(xright-xleft+1)) 
        inverse_mass_total=1./total(cm_image(cm_indices)) 
        xcenter=total(xpositions*cm_image(cm_indices))*inverse_mass_total 
        ycenter=total(ypositions*cm_image(cm_indices))*inverse_mass_total 
        IF keyword_set(debug) THEN BEGIN 
            print,'Center of mass is at ['+$ 
              strtrim(string(xcenter,format='(f10.2)'),2)+','+$ 
              strtrim(string(ycenter,format='(f10.2)'),2)+']' 
        ENDIF 
    ENDIF 
ENDIF ELSE BEGIN 
    xpts=xcenter+[-1,0,1] 
    ypts=fft2((xcenter-1):(xcenter+1),ycenter) 
    tri_fit,xpts,ypts,xfit,xpeak 
    xpts=ycenter+[-1,0,1] 
    ypts=fft2(xcenter,(ycenter-1):(ycenter+1)) 
    tri_fit,xpts,ypts,yfit,ypeak 
    xcenter=xpeak 
    ycenter=ypeak 
    IF keyword_set(debug) THEN BEGIN 
        print,'Linear+parabolic center is at ['+$ 
          strtrim(string(xcenter,format='(f10.2)'),2)+','+$ 
          strtrim(string(ycenter,format='(f10.2)'),2)+']' 
    ENDIF 
ENDELSE 
 
IF (n_elements(xplot_win) NE 0) THEN BEGIN 
    wset,xplot_win 
    wshow,xplot_win 
    xpts=x_integer_center+findgen(13)-6 
    zpts=fft2((x_integer_center-6):(x_integer_center+6),y_integer_center) 
    plot,xpts,zpts,xtitle='X center',psym=2,ystyle=16 
    oplot,xcenter+[0.,0.],!y.crange,linestyle=1 
ENDIF 
 
IF (n_elements(yplot_win) NE 0) THEN BEGIN 
    wset,yplot_win 
    wshow,yplot_win 
    ypts=y_integer_center+findgen(13)-6 
    zpts=fft2(x_integer_center,(y_integer_center-6):(y_integer_center+6)) 
    plot,ypts,zpts,xtitle='Y center',psym=2,ystyle=16 
    oplot,ycenter+[0.,0.],!y.crange,linestyle=1 
ENDIF 
 
IF (n_elements(xcorimg_win) NE 0) THEN BEGIN 
    wset,xcorimg_win 
    wshow,xcorimg_win 
 
    red_index = !d.table_size-1 
    green_index = !d.table_size-2 
    white_index = !d.table_size-3 
    r = byte(0.5+(255.*findgen(!d.table_size)/ $ 
                  float(!d.table_size-3))<255.) 
    g = r 
    b = r 
    r(green_index:red_index) = byte([0,255]) 
    g(green_index:red_index) = byte([255,0]) 
    b(green_index:red_index) = byte([0,0]) 
    tvlct,r,g,b 
 
    IF (long(!d.n_colors) EQ 255L) THEN BEGIN 
        red_plot_index = red_index 
        green_plot_index = green_index 
        white_plot_index = white_index 
    ENDIF ELSE IF (long(!d.n_colors) EQ 16777216L) THEN BEGIN 
        ; For 24 bit color, use (r+256L*(g+256L*b)) 
        red_plot_index = 255L 
        green_plot_index = 255L*256L 
        white_plot_index = long(white_index)+$ 
            256L*(long(white_index)+256L*long(white_index)) 
    ENDIF ELSE BEGIN 
        print,'Unknown display depth' 
        red_plot_index = red_index 
        green_plot_index = green_index 
        white_plot_index = white_index 
    ENDELSE 
 
    IF (xcorimg_zoom GT 1) THEN BEGIN 
        tv,rebin(bytscl(alog(fft2),top=white_index),$ 
            (xcorimg_zoom*nx),(xcorimg_zoom*ny),/sample) 
    ENDIF ELSE BEGIN 
        tv,bytscl(alog(fft2),top=white_index) 
    ENDELSE 
    plots,(nx/2)*xcorimg_zoom+[-8,-4],$ 
              (ny/2)*xcorimg_zoom+[0,0],$ 
              /device,color=green_plot_index 
    plots,(nx/2)*xcorimg_zoom+[4,8],$ 
              (ny/2)*xcorimg_zoom+[0,0],$ 
              /device,color=green_plot_index 
    plots,(nx/2)*xcorimg_zoom+[0,0],$ 
              (ny/2)*xcorimg_zoom+[-8,-4],$ 
              /device,color=green_plot_index 
    plots,(nx/2)*xcorimg_zoom+[0,0],$ 
              (ny/2)*xcorimg_zoom+[4,8],$ 
              /device,color=green_plot_index 
    plots,[fix(0.5+xcenter*xcorimg_zoom)],$ 
            [fix(0.5+ycenter*xcorimg_zoom)],/device,$ 
            color=red_plot_index,psym=6 
ENDIF 
 
fft2 = 0 
xshift = xcenter - double(nx/2) 
yshift = ycenter - double(ny/2) 
 
IF ((abs(xshift) GT 0.2) OR (abs(yshift) GT 0.2)) THEN BEGIN 
    p = [xshift,0.,1.,0.] 
    q = [yshift,1.,0.,0.] 
    IF (meanfill NE 0) THEN BEGIN 
        missing = total(image2)/float(nx*ny) 
        image2_shifted = poly_2d(image2,p,q,1,missing=missing,cubic=-0.5) 
    ENDIF ELSE IF (medianfill NE 0) THEN BEGIN 
        missing = median(image2) 
        image2_shifted = poly_2d(image2,p,q,1,missing=missing,cubic=-0.5) 
    ENDIF ELSE BEGIN 
        image2_shifted = poly_2d(image2,p,q,1,cubic=-0.5) 
    ENDELSE 
ENDIF ELSE BEGIN 
    image2_shifted = image2 
ENDELSE 
 
return 
END 
 
