function polymap, x, c

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Computes a polynomial function where both the input value and the 
;;	coefficients can be full images. The polynomial is computed 
;;	as y = c[0] + c[1]*x + c[2]*x^2 + ... where the index in c is the
;;	index in the last dimension
;;
;; CALLING SEQUENCE:
;;	y = polyMap(x,c)
;;
;; REQUIRED INPUTS:
;;	x - Input value to the polynomial function. Can be any number of dimensions
;;	c - coefficients for the polynomial. This can either be an image cube where
;;		each plane is the same size as x or it can be a vector specifying the
;;		coefficients. The image cube option only works when applying the data
;;		to a two dimensional x image.
;;
;; OUTPUTS:
;;	y - the result of the polynomial function at all points in x
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> y = polyMap(x,c)
;;
;; PROCEDURES USED (i.e. called directly!):
;;
;; MODIFICATION HISTORY:
;;   2005-06-07  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

on_error,2

;; Get the size of the inputs
szx = size(x)
szc = size(c)

;; Determine what type of coefficients were used: image cube or vector
IF szc[0] EQ 1 THEN BEGIN
	;; A vector of coefficients was used
	return, poly(x,c)
ENDIF ELSE IF szc[0] EQ 3 THEN BEGIN
	;; Check inputs
	IF max(abs(szx[1:2]-szc[1:2])) NE 0 THEN $
		message, 'Coefficient cube must be the same size as the input image'
	IF szx[0] NE 2 THEN message, 'Input must be a two dimensional image'

	;; Perform the polynomial operation
	out = c[*,*,szc[3]-1]
	FOR i=szc[3]-2, 0, -1 DO out = x*out + c[*,*,i]
	return, out
ENDIF ELSE message, 'Coefficients can only be a vector or an image cube'

END