PRO SplineImage, image, axis, n, direction, first, step, nbEls

;+
; NAME:
;	Image Spline Interpolation
; PURPOSE:
;	Assuming the image is composed of irregular measurements
;	in the X or Y axis, this procedure recomputes the
;	image values at regular steps, using the spline interpolation
;	technique.
; CALLING SEQUENCE:
;	SplineImage, image, axis, n [ , direction 
;		[ , first, step, nbEls ] ]
; INPUT PARAMETERS:
;	image: a 2d real (or int) array
;	axis: the axis to regularize (X or Y)
;	n: the number of elements in axis
;	direction: 'X' or 'Y'. If not given, 'Y' is default.
;	first: the first element of the new (regular) axis. 
;		Default: axis(0)
;	step: the step between two axis values.
;		Default: (axis(n-1)-axis(0))/n
;	nbels: the number of elements of the new axis.
;		Default: n
; OUTPUT PARAMETERS:
;	image: the redimensioned and computed image
;	axis: the redimensioned and computed axis
; SIDE EFFECT:
;	Can be very cpu-time consuming. You can run this 
;	procedure with less priority when you start wave with
;	the "nohup" command or stop the process and run it in the
; 	background with th "bg" command.
; SEE ALSO:
;	Spline
; MODIFICATION HISTORY
;	Created: A.Csillaghy, ETHZ, Jan 93.
;-


  IF N_Params() LT 3 THEN BEGIN
    Message, 'Usage: SplineImage,image,axis,n ' + $
	'[,direction [,first,step,nbEls]]', /INFO, /CONT
    RETURN
  END

  IF N_Elements( direction ) EQ 0 THEN direction = 'Y' 
  IF N_Elements( first ) EQ 0 THEN first = Min(axis)
  IF N_Elements( nbEls ) EQ 0 THEN nbEls = n
  IF N_Elements( step ) EQ 0 THEN $
	step = (Double( axis(n-1) ) - axis(0) ) / (nbEls-1)

  newAxis = DIndGen( nbEls ) * Abs(step) + first
  IF step LT 0 THEN newAxis = Reverse( newAxis )
  siz = Size( image )

  IF siz(0) EQ 1 THEN nbSplines = 1 $
  ELSE BEGIN
    IF direction NE 'X' THEN image = Transpose( image )
    nbSplines = N_Elements( image( 0, * ) )
  ENDELSE

  newImage = FltArr( nbEls, nbSplines )

  IF nbSplines GT 10 THEN $
    Print, 'Computing spline interpolations.'

  CASE 1 OF
    axis(0) GT axis(n-1) AND newAxis(0) GT newAxis(nbEls-1): $    
      BEGIN
        axisRev = Reverse( axis )
        newAxisRev = Reverse( newAxis )
        DisplayRestTime, /INIT
        FOR i = 0,nbSplines-1 DO BEGIN
          newImage(*,i) = Reverse( Spline( axisRev, $
	Reverse(image(*,i)), newAxisRev  ) )
          IF i EQ 10 THEN DisplayRestTime, Double(i)/nbSplines
        ENDFOR
      END
    axis(0) GT axis(n-1): $
      BEGIN
        axisRev = Reverse( axis )
        DisplayRestTime, /INIT
        FOR i = 0,nbSplines-1 DO BEGIN
          newImage(*,i) =  Spline( axisRev , $
	Reverse(image(*,i)),  newAxis ) 
          IF i EQ 10 THEN DisplayRestTime,Double(i)/nbSplines
        ENDFOR
      END
    ELSE: BEGIN
        DisplayRestTime, /INIT
        FOR i = 0,nbSplines-1 DO BEGIN
          newImage(*,i) =  Spline( axis, image(*,i),  newAxis ) 
          IF i EQ 10 THEN DisplayRestTime, Double(i)/nbSplines
        ENDFOR
      END
  ENDCASE

  axis = newAxis
  IF direction NE 'X' THEN image = Transpose( newImage ) $
  ELSE image = newImage

END ; Image spline interpolation.