;+
; $Id: test_crosscorr.pro,v 1.2 2010/05/11 16:00:17 nathan Exp $
; NAME:
;	CROSSCORR
;
; PURPOSE:
;	To perform the cross correlation function of
;	two images
;
; CATEGORY:
;	PICO
;
; CALLING SEQUENCE:
;	result=CROSSCORR(array1, array2[, xmax, ymax])
;
; INPUTS:
;	array1, array2: Two arrays of the same size
;		containing the image data
;
; OPTIONAL INPUTS:
;	xmax, ymax: The index in x and y direction of
;		the detected shift. If two images are
;		similar to each other but shifted, the
;		best congruence would be assured if the
;		second image was shifted by (xmax,ymax).
;		See also EXAMPLES. The xmax and ymax
;               are given in tenths of pixels (one decimal)
;
; KEYWORD PARAMETERS:
;	VERBOSE: prints the time needed out
;       GAUSS:  If set, a gaussian fit is performed to
;               get the maximum of the Correlation. For
;               this fit the correlation function must be
;               rather well peaked. GAUSS2DFIT is used.
;       MASK:   If set, array1 and array2 are multiplied with
;               a pyramidal mask in order to avoid edge
;               problems by the scroll-around technique.
;               This is of particular interest if opposite
;               edges of the arrays have highly different
;               intensities (for example a step image where
;               one edge is 0 and the other 1)
;
; OUTPUTS:
;	result: A float array containing the
;		cross correlation function about the
;		centered images.
;
; OPTIONAL OUTPUTS:
;	None
;
; COMMON BLOCKS:
;	None
;
; SIDE EFFECTS:
;	Unknown
;
; EXAMPLE:
;	a=CROSSCORR(image1,image2,xmax,ymax)
;	image2=shift(image2,xmax,ymax)
;
;	The image 2 is shifted in a way that it
;	will fit best to image1
;
; RESTRICTIONS:
;	The two arrays must have the same dimensions.
;	The cross correlation will always be centered
;	on unshifted images
;
; PROCEDURE:
;	The calculation is performed by the fourier
;	space using the inverse fourier transform
;	of the product of the two individual fourier
;	transforms, similar to the Wiener-Kintchine
;	theorem:
;	result=F^(-1)( F(array1) * CONJ(F(array2)) )
;
; MODIFICATION HISTORY:
; $Log: test_crosscorr.pro,v $
; Revision 1.2  2010/05/11 16:00:17  nathan
; added debug messages
;
;	V1.0 written by Alo Epple, Dez 1994
;	V2.0 Alo Epple, Pic Du Midi, OCT-1995: using the
;		fourier algorithm rather than the CORRELATE
;		function in nested loops.
;       V2.1 Alo Epple, MPAe Lindau, FEB-1996: Detecting
;               sub-pixel shifts by SFIT to the result.
;-

FUNCTION test_crosscorr,arr1,arr2,xmax,ymax, $
                   GAUSS=gauss, VERBOSE=verbose, $
                   MASK=kmask, FLOAT=kfloat, $
                   tope=tope

time=SYSTIME(1)

IF (TOTAL(ABS(SIZE(array1)-SIZE(array2))) GT 0) THEN BEGIN
    MESSAGE,'The arrays must have the same dimensions'
ENDIF

array1=arr1
array2=arr2
;array1=arr1-MEAN(arr1)
;array2=arr2-MEAN(arr1)

sz=SIZE(array1)
xdim=sz(1)
ydim=sz(2)

IF KEYWORD_SET(kmask) THEN BEGIN
    mask_x=FINDGEN(xdim/2)
    IF ((xdim MOD 2) EQ 0) THEN BEGIN
        mask_x=[ROTATE(mask_x,2),mask_x]
    ENDIF ELSE BEGIN
        mask_x=[ROTATE(mask_x,2)+1,0,mask_x+1]
    ENDELSE

    mask_x=1-mask_x/MAX(mask_x)

    mask_y=FINDGEN(ydim/2)
    IF ((xdim MOD 2) EQ 0) THEN BEGIN
        mask_y=[ROTATE(mask_y,2),mask_y]
    ENDIF ELSE BEGIN
        mask_y=[ROTATE(mask_y,2)+1,0,mask_y+1]
    ENDELSE
    mask_y=1-mask_y/MAX(mask_y)

    mask=mask_x#mask_y
ENDIF

IF KEYWORD_SET(kmask) THEN BEGIN
    result=FLOAT(FFT(FFT(array1*mask,-1)*CONJ(FFT(array2*mask,-1)),1))
ENDIF ELSE BEGIN
    result=FLOAT(FFT(FFT(array1,-1)*CONJ(FFT(array2,-1)),1))
ENDELSE

;result=result/SQRT(TOTAL(array1)^2*TOTAL(array2)^2)
;result=result*((FLOAT(xdim)*ydim)^2)

max1 = MAX(result)
maximum=(WHERE(result EQ max1,count))(0)

xmax=maximum MOD xdim
ymax=maximum/xdim

xshift=xdim/2-xmax
yshift=ydim/2-ymax

IF keyword_set(VERBOSE) THEN help,xshift,yshift,max1,count,xdim
IF KEYWORD_SET(gauss) THEN BEGIN
    a=GAUSS2DFIT(SHIFT(result,xshift,yshift),coeff)
    xmax=coeff(4)-xshift
    ymax=coeff(5)-yshift
    IF keyword_set(VERBOSE) THEN print,'coeff=',coeff

ENDIF ELSE IF KEYWORD_SET(kfloat) THEN BEGIN
    IF (count GT 0) THEN BEGIN

    ENDIF

; Select the npix x npix pixels around the maximum to make a surface
; fit in order to get the sub-pixel maximum

    n_pix=7
    result=SHIFT(result,n_pix/2-xmax,n_pix/2-ymax)
    maxarr=result(0:n_pix-1,0:n_pix-1)
    result=SHIFT(result,-(n_pix/2-xmax),-(n_pix/2-ymax))

    a=SFIT(CONGRID(maxarr,10*n_pix,10*n_pix),2)

    fitmaximum=WHERE(a EQ MAX(a),count)

    fitsz=SIZE(a) & fitxdim=fitsz(1) & fitydim=fitsz(2)

    IF (count GT 0) THEN BEGIN
        fitmaximum=fitmaximum(0)
        fitymax=fitmaximum/fitxdim
        fitxmax=fitmaximum MOD fitxdim
    ENDIF

    fitxmax=fitxmax/10.
    fitymax=fitymax/10.

    xmax=xmax-n_pix/2+fitxmax-.5
    ymax=ymax-n_pix/2+fitymax-.5

ENDIF

IF (xmax GT xdim/2) THEN xmax=xmax-xdim
IF (ymax GT ydim/2) THEN ymax=ymax-ydim

;;;;;;;;;;;IF KEYWORD_SET(verbose) THEN MESSAGE,/INFO,/CONT,SC(SYSTIME(1)-time,DEC=2)+' sec.'

;;;;;tope = max(a)
if keyword_set(VERBOSE) THEN BEGIN
    print,'result:'
    maxmin,result
ENDIF
RETURN,result

ende:
END