# ifndef PS_IMAGE_H
# define PS_IMAGE_H

/** \file psImage.h
 *  \brief Basic image definitions and operations.
 *  \ingroup AstroGroup
 */

/// basic image data structure.
typedef struct psImage {
    psType type; 			///< image data type and dimension
    const int ncols, nrows;		///< size of image 
    const int x0, y0;			///< data region relative to parent 
    union {
        psS8  **S8;                     ///< Pointers to char data
        psS16 **S16;                    ///< Pointers to short-integer data
        psS32 **S32;                    ///< Pointers to integer data
        psS64 **S64;                    ///< Pointers to long-integer data
        psU8  **U8;                     ///< Pointers to unsigned-char data
        psU16 **U16;                    ///< Pointers to unsigned-short-integer data
        psU32 **U32;                    ///< Pointers to unsigned-integer data
        psU64 **U64;                    ///< Pointers to unsigned-long-integer data
        psF32 **F32;                    ///< Pointers to floating-point data
        psF64 **F64;                    ///< Pointers to double-precision data
        psC32 **C32;                    ///< Pointers to complex floating-point data
        psC64 **C64;                    ///< Pointers to complex floating-point data
    } data;
    const struct psImage *parent;	///< parent, if a subimage 
    psArray *children;			///< children of this region
} psImage;

/** Functions **************************************************************/
/** \addtogroup AstroGroup Astronomy-Specific Utilities
 *  \{
 */

/*** Image structure manipulation ***/

/// Create an image of the specified size and type.
psImage *
psImageAlloc (int width,		///< image width (ncols)
	      int height,		///< image height (nrows)
	      psElemType type)		///< image data type 
;

/// Create a subimage of the specified area.
psImage *
psImageSubset(const psImage *image,	///< parent image 
	      int nx,			///< subimage width (<= image.ncols - col0)  
	      int ny,			///< subimage height (<= image.nrows - row0)  
	      int x0,			///< subimage starting col (0 <= col0 < ncols)   
	      int y0)			///< subimage starting row (0 <= row0 < nrows)   
;

/// Create a subimage of the area specified by an image section, [x1:x2,y1:y2]
psImage *
psImageSubsection(const psImage *image,	///< Parent image
		  const char *section	///< Image section, of the form [x1:x2,y1:y2]
    );


/// Create a copy of the specified image.
psImage *
psImageCopy(psImage *output, 		///< target structure for output image data
	    const psImage *input,	///< copy this image 
	    psElemType type)		///< element type of output image
;


/// Trim an image, in-place in memory
psImage *
psImageTrim(psImage *image,		///< Input image, and output
	    int x0, int x1,		///< Keep x0 to x1
	    int y0, int y1		///< Keep y0 to y1
    );
	    
/*** various image pixel extractions ***/

/** Direction in which to cut */
typedef enum {
    PS_CUT_X_POS,                     ///< Cut in positive x direction
    PS_CUT_X_NEG,                     ///< Cut in negative x direction
    PS_CUT_Y_POS,                     ///< Cut in positive y direction
    PS_CUT_Y_NEG                      ///< Cut in negative y direction
} psImageCutDirection;

/// Extract pixels from rectlinear region to a vector.
psVector *
psImageSlice(psVector *out,		///< Vector to output, or NULL
	     const psImage *input,	///< Input image
	     const psImage *mask,	///< Ignore those pixels where mask & maskVal == 1.  May be NULL
	     unsigned int maskVal,	///< Value in mask to ignore
	     unsigned int x,            ///< starting x coord of region to slice
	     unsigned int y,            ///< starting y coord of region to slice
	     unsigned int nx,           ///< length of region in x
	     unsigned int ny,           ///< length of region in y
	     psImageCutDirection direction, ///< Direction in which to cut
	     const psStats *stats       ///< defines statistics used to find output values
    );

/// Extract pixels along a line to a vector.
psVector *
psImageCut(psVector *out,		///< Vector to output, or NULL
	   const psImage *input,	///< extract cut from this image
	   const psImage *mask,		///< Ignore those pixels where mask & maskVal == 1.  May be NULL
	   unsigned int maskVal,	///< Value in mask to ignore
	   float xs, 			///< starting x coord of cut
	   float ys, 			///< starting y coord of cut
	   float xe, 			///< ending x coord of cut
	   float ye, 			///< ending y coord of cut
	   unsigned int nSamples,	///< Number of samples along the slice
	   psImageInterpolateMode mode	///< defines statistics used to find output values
	   );

/// Extract radial annulii data to a vector.
psVector *
psImageRadialCut(psVector *out,		///< Vector to output, or NULL
		 const psImage *input,	///< extract profile from this image
		 const psImage *mask,	///< Ignore those pixels where mask & maskVal == 1.  May be NULL
		 unsigned int maskVal,	///< Value in mask to ignore
		 float x, 		///< center x coord of annulii
		 float y, 		///< center y coord of annulii
		 const psVector *radii, ///< Radii of the annuli
		 const psStats *stats)	///< defines statistics used to find output values
;

/*** various image geometry manipulation ***/

/** How to do interpolation */
typedef enum {
    PS_INTERPOLATE_FLAT, 
    PS_INTERPOLATE_BILINEAR, 
    PS_INTERPOLATE_BICUBIC, 
    PS_INTERPOLATE_SINC 
} psImageInterpolateMode mode;

/// Rebin image to new scale.
psImage *
psImageRebin(psImage *out,		///< Image to output, or NULL
	     const psImage *in,		///< rebin this image
	     const psImage *restrict mask, ///< Ignore those pixels where mask & maskVal == 1.  May be NULL
	     unsigned int maskVal,	///< Value in mask to ignore
	     int scale, 		///< rebinning scale
	     const psStats *stats)	///< defines statistics used to find output values
;

/** Resample image to new scale */
psImage *
psImageResample(psImage *out,		///< Image to output, or NULL
		const psImage *in,	///< resample this image
		int scale,		///< Factor by which to resample
		psImageInterpolateMode mode ///< Interpolation method
    );

/// Rotate image by given angle.
psImage *
psImageRotate(psImage *out,		///< Image to output, or NULL
	      const psImage *in,	///< rotate this image
	      float angle		///< rotate by this amount anti-clockwise (degrees)
	      psC64 exposed,		///< Set exposed pixels to this value
	      psImageInterpolateMode mode ///< Interpolation method
    );

/// Shift image by an arbitrary number of pixels in either direction.
/// Drop edge pixels, and set newly exposed pixels to 'exposed'.
psImage *
psImageShift(psImage *out,		///< Image to output, or NULL
	     const psImage *in,		///< shift this image
	     float dx,			///< shift by this amount in x
	     float dy,			///< shift by this amount in y
	     psC64 exposed,		///< set exposed pixels to this value
	     psImageInterpolateMode mode ///< Interpolation method
);

/// Roll image by an integer number of pixels in either direction.
/// Edge pixels wrap to the other side (no values are lost).
psImage *
psImageRoll(psImage *out,		///< Image to output, or NULL
	    const psImage *in,	///< roll this image
	    int dx, 			///< roll this amount in x
	    int dy)			///< roll this amount in y
;

/// Determine statistics for image (or subimage).
psStats *
psImageGetStats(psStats *stats,		///< defines statistics to be calculated & target
		const psImage *in, 	///< image (or subimage) to calculate stats
		const psImage *mask,	///< Ignore those pixels where mask & maskVal == 1.  May be NULL
		unsigned int maskVal	///< Value in mask to ignore
    );

/// Construct a histogram from an image (or subimage).
psHistogram *
psImageHistogram(psHistogram *hist,	///< input histogram description & target
		 const psImage *in	///< determine histogram of this image
		 const psImage *mask,	///< Ignore those pixels where mask & maskVal == 1.  May be NULL
		 unsigned int maskVal	///< Value in mask to ignore
    );

/// Fit a 2-D polynomial surface to an image.
psPolynomial2D *
psImageFitPolynomial(psPolynomial2D *coeffs, ///< coefficient structure carries in desired terms & target
		     const psImage *input) ///< image to fit
;

/// Evaluate a 2-D polynomial surface to image pixels.
psImage *
psImageEvalPolynomial(psImage *input,	///< image to produce
		      const psPolynomial2D *coeffs) ///< coefficient structure carries in desired terms
;

/** Perform interpolation of image pixel values */
psF64
psImagePixelInterpolate(const psImage *input, ///< Image on which to interpolate
			float x, float y, ///< Fractional pixel coordinates at which to interpolate
			psC64 unexposedValue, ///< Value for pixels that don't have an input value
			psImageInterpolateMode mode ///< Interpolation method
    );


/*** image input/output routines ***/
/// Read an image or subimage from a named file.
psImage *
psImageReadSection(psImage *output, 	///< place data in this structure for output 
		   int x0, 		///< starting x coord of region
		   int y0, 		///< starting y coord of region
		   int nx, 		///< x size of region (-1 for full range)
		   int ny, 		///< y size of region (-1 for full range)
		   int z, 		///< plane of interest
		   const char *extname,	///< MEF extension name ("PHU" for primary header)
		   int extnum,		///< MEF extension number (-1 for "PHU", 0 : Nextend - 1)
		   const char *filename) ///< file to read data from
;

/// Write an image section to named file (which may exist).
psImage *
psImageWriteSection(const psImage *input, ///< image to write out
		    int x, 		///< starting x coord of region		   
		    int y, 		///< starting y coord of region		   
		    int z, 		///< plane of interest			   
		    const char *extname, ///< MEF extension name ("PHU" for primary header)
		    int extnum, 	///< MEF extension number (-1 for "PHU", 0 : Nextend - 1)
		    const char *filename) ///< file to write data to		   
;

/*** basic pixel manipulations ***/
/// Clip image values outside of range to given values.  Return number of clipped pixels.
int
psImageClip(psImage *input, 		///< Image to be clipped, and output
	    double min,			///< clip pixels with values < min
	    double vmin, 		///< set min-clipped pixels to vmin
	    double max,			///< clip pixels with values > max
	    double vmax)		///< set max-clipped pixels to vmax
;

/** Clipping for complex data */
int
psImageClipComplexRegion(psImage *input, ///< Image to be clipped, and output
			 complex double min, ///< clip pixels with values < min
			 complex double vmin, ///< set min-clipped pixels to vmin
			 complex double max, ///< clip pixels with values > max
			 complex double vmax ///< set max-clipped pixels to vmax
    );

/// Clip NaN image pixels to given value.  Return number of clipped pixels.
int
psImageClipNaN(psImage *input,		///< clip this image & target
	       float value)		///< set nan pixels to this value
;

/// Overlay subregion of image with another image.  Return number of pixels replaced.
int 
psImageOverlaySection(psImage *image,	///< input image & target
		      const psImage *overlay, ///< image to overlay 
		      int x0,		///< x offset of overlay subimage 
		      int y0,		///< y offset of overlay subimage 
		      const char *op)	///< overlay operation 
;

/* \} */ // End of AstroGroup Functions

# endif
/* image overlay operations 
   PS_OVERLAY_EQUALS   = '=', 
   PS_OVERLAY_ADD      = '+', 
   PS_OVERLAY_SUBTRACT = '-', 
   PS_OVERLAY_MULTIPLY = '*', 
   PS_OVERLAY_DIVIDE   = '/', 
*/
