FUNCTION AutoCorr, imageP, XLAGMAX = xLagMax, $
	YLAGMAX = yLagMax, $
	NORMALIZED = normalized,  $
	SQUARESUM = squareSum, FAST = fast

;+
; PURPOSE: 
;	Computes the autocorrelation of an image or a vector.
; CALLING SEQUENCE:
;	result = AutoCorr( image )
; INPUT:
;	image: a 2D array (or 1D array), real
; RESULT:
;	the autocorrelated image
; KEYWORD:
;	/NORMALIZED: if set, the autocorrelation is normalized.
;	XLAGMAX: the max. lag used in the "square sum" method, 
;		 the limit of the array returned in the FFT method.
;		 default:  nx-1 (nx: nb elements in x)
;	YLAGMAX: the max. lag used in the "square sum" method, 
;		 the limit of the array returned in the FFT method.
;		 default: ny-1 (ny: nb elements in y)
;	SQUARESUM: output of the sum of the elements squared.
;	FAST: if set the Wiener-Khinchin theorem is used. 
;		Otherwise it is calculated from its definition.
;		
; PROCEDURE:
;	With the Wiener-Khintchin Theorem: the
;	autocorrelation of an image is the backward fourrier
;	transform of the power spectrum.
;	With the definition: see the manual.
;	With FFT: In the two-dimensional case, the 0 lag
;	is shifted to the middle of the picture (i.e. the
;	image is symmetrical. In the FFT 1dimensional case,
;	however, only one half of the autocorrelation
;	is shown (and therefore not shifted).
; RESTRICTION:
;	The Fourier method is much faster, but much less
;	accurate than the method based on the algorithm.
;       CAUTION: the FFT method assumes periodic boundaries,
;	i.e. the function is repeated infinitely.
; MODIFICATION HISTORY:
;	Created in August 1991 by A.Csillaghy
;		Inst. of Astronomy, ETH Zurich
;	FFT Case: 1d no shift in July 1995, ACs.
;-

  Print, 'Working ... '

  image = imageP - Avg( imageP )

  nx = N_Elements( image( *, 0 ) )
  ny = N_Elements( image( 0, * ) )

  squareSum = Total ( image^2 )

  IF Keyword_Set( XLAGMAX ) EQ 0 THEN xLagMax = nx
  IF Keyword_Set( YLAGMAX ) EQ 0 THEN yLagMax = ny

  xLagMax = xLagMax < (nx -1) > 0
  yLagMax = yLagMax < (ny -1) > 0

  IF Keyword_Set( FAST ) THEN BEGIN
    nEls = nx*ny
    ac =  Float( FFT( Abs( FFT( image, -1 )  )^2 , + 1 ))/nEls
    IF Keyword_Set( NORMALIZED ) THEN ac = ac/ac(0,0) 
    ac = ac(0:xLagMax, 0:yLagMax)
    IF ny GT 1 THEN  RETURN, Shift( ac, xLagMax/2, yLagMax/2 ) $
    ELSE RETURN, ac(0:(xLagMax/2))
  ENDIF

  ac = FltArr( xLagMax+1, yLagMax+1 )
  FOR j = 0, yLagMax DO $
    FOR i = 0, xLagMax  DO  BEGIN
        array1 = image( 0:nx-i-1, 0:ny-j-1)
        array2 = image( i: nx-1, j:ny-1)
        ac( i, j ) = Total( array1*array2 )
    ENDFOR

  IF Keyword_Set( NORMALIZED ) THEN BEGIN
    RETURN,  ac / squareSum 
  ENDIF ELSE RETURN, ac

END