#ifndef _H_I3_IMAGE
#define _H_I3_IMAGE
#include "i3_global.h"

/* 
The complex header has to come before the fftw3 header,
since then fftw3 uses the C99 complex type rather than defining its own version.

So then we can do maths on the numbers easily.
*/
#include "complex.h"
#include "fftw3.h"


/**
i3_image.h
Created by JAZ.  27/10/10.

Defines the i3_image struct which contains an image or chunk thereof.

For quick access the image is stored as a contiguous chunk of memory (image->data), but row pointers
are set up to the start of each row so that there is a friendly access mode too:
image->row[i][j]

Defines a constructor and destructor for the type, and a few other useful tools.
*/


/**
* A basic image type.
* Can be accessed either via the main data buffer or row-by-row.
* nx and ny are the dimensions.
*/
typedef struct i3_image {
	i3_flt * data; /**< The main pointer to the image data*/
	i3_flt **row;/**< An array of pointers to each row of the image*/
	unsigned long n; /**< The total number of pixels in the image*/
	unsigned long nx; /**< The number of pixels per row (horizontal image dimensions)*/
	unsigned long ny;/**< The number of pixels per column (vertical image dimensions)*/
} i3_image;

/**
* A type to hold the fourier transform of an image
* Like the i3_image type, it can be accessed either via the main data buffer or row-by-row.
* It stores the dimension of both the original image and the fourier transform..
*/
typedef struct i3_fourier{
	i3_cpx * data; /**< The main pointer to the fourier transformed data.  fftw_complex has two eight byte parts (real and imaginary)*/
	i3_cpx ** row; /**< An array of pointers to each row of the FT.  */
	unsigned long nx_real; /**< The x dimension size of the image that the i3_fourier instance came from.*/
	unsigned long ny_real; /**< The y dimension size of the image that the i3_fourier instance came from.*/
	unsigned long n_real; /**< The total size of the image that the i3_fourier instance came from.*/
	unsigned long nx_f; /**< The x dimension size of the fourier transform data.  Because of the packing of 
							 the negative frequency components of the FFT this is not the same as nx_real; 
							 it is nx_real/2+1.  The fact that nx_real=2j and nx_real=2j+1 yield the same nx_f 
							is why we save all info about the original image. */
	unsigned long ny_f; /**< The y dimension of the fourier transform data.  This is the same as the real y dimension. */
	unsigned long n_f; /**< The toal size of the fourier transform data */
	unsigned long x_padding_real; /**< [Not yet used] The x padding applied to the image before transforming*/
	unsigned long y_padding_real; /**< [Not yet used] The y padding applied to the image before transforming*/
} i3_fourier;


/**
* A type representing the moments of an image
*/
typedef struct i3_image_moments {
	i3_flt sum; /**< The zeroth moment (sum) of the image.*/

	i3_flt x0; /**< The x first moment of the image - the x coordinate of the image center*/
	i3_flt y0; /**< The y first moment of the image - the y coordinate of the image center*/
	       
	i3_flt qxx; /**< The xx second moment of the image */
	i3_flt qxy; /**< The xy second moment of the image*/
	i3_flt qyy; /**< The yy second moment of the image*/
	
	i3_flt e1; /**< The x-y shear of the image (computed from qij) */
	i3_flt e2; /**< The 45 degree shear of the image (computed from qij) */
	
	} i3_image_moments;


/**
 * Create an i3_image.  
 *
 * Allocates memory and sets up row pointers.  Use i3_image_zero after creation to get a zero image.
 * 
 * \param nx The x dimension of the array to be created
 * \param ny The y dimension of the array to be created
 * \return An allocated but uninitialized image.
**/
i3_image * i3_image_create(unsigned long nx, unsigned long ny);

/**
 * Create an image with the same dimensions as another image.
 * \param image The template image to use as a basis for the new one.
 * \return An uninitialized image with the same dimension as the input.
**/
i3_image * i3_image_like(i3_image * image);


/**
 * Destroy an i3_image.
 * 
 * Deallocates all memory associated with the image.
 * \param image Pointer to the image to destroy. 
**/
void i3_image_destroy(i3_image * image);

/**
 * Copy a rectangular chunk of an i3_image to a new i3_image.
 * \param image The original full-size image
 * \param x_start The x-coordinate of the top-left corner of the full-size image to start from
 * \param y_start The y-coordinate of the top-left corner of the full-size image to start from
 * \param width The x-size in pixels of the new image.
 * \param height The y-size in pixels of the new image.
 * \return A new copy of part of the image, or NULL if you tried to copy something outside the image bounds.
**/
i3_image * i3_image_copy_part(i3_image * image, long x_start, long y_start, unsigned long width, unsigned long height);

/**
 * Copy a rectangular chunk of an i3_image into an existing i3_image.
 * \param image The original full-size image
 * \param x_start The x-coordinate of the top-left corner of the full-size image to start from
 * \param y_start The y-coordinate of the top-left corner of the full-size image to start from
 * \param subimage The image into which to put the chunk.  The size of the chunk is determined by the image dimensions.
 * \return A status code, zero for success.
**/
int i3_image_copy_part_into(i3_image * image, long x_start, long y_start, i3_image * subimage);


/** 
 * Copy a complete image.
 * \param image The image to copy.
 * \return A new copy of the image
**/
i3_image * i3_image_copy(i3_image * image);


int i3_image_copy_into_part(i3_image * subimage, long x_start, long y_start, i3_image * image);

/** 
 * Copy a complete image into an existing one.
 * \param image The image to copy.
 * \param clone The image into which to put the copy.  Must have same size as image.
 * \return A status code, zero for success.
**/
int i3_image_copy_into(i3_image * image, i3_image * clone);

int i3_fourier_copy_into( i3_fourier * ft_src, i3_fourier * ft_dst );

void i3_image_copy_from_pointer(i3_image * image, i3_flt * data);

/**
 * Get the maximum value of an image.
 * \param image The image to get the value from.
 * \return The maximum pixel in the image.  NaNs will be ignored.
**/
i3_flt i3_image_max(i3_image *image);
i3_flt i3_image_min(i3_image *image);
i3_flt i3_image_max_masked(i3_image *image, i3_image * mask);
i3_flt i3_image_min_masked(i3_image *image, i3_image * mask);



/**
 * Get the L2 norm of an image treated as a 1-D vector.
 * \param image The image to get the value from.
 * \return The norm = RMS = i3_sqrt(sum pixel**2) of the image
**/
i3_flt i3_image_norm(i3_image * image);

/**
 * Get the square of L2 norm of an image treated as a 1-D vector.
 * \param image The image to get the value from.
 * \return The norm = RMS = sum pixel**2 of the image
**/
i3_flt i3_image_norm_sq(i3_image * image);


/**
 * Save an image to a text file, with ny columns and nx rows per column.
 * \param image The image to save
 * \param filename The name to the file to save the image to.
 * \return A status code, zero for success.
**/
int i3_image_save_text(i3_image *image, char * filename);

/**
 * Save an image to a flat binary file; just a straight dump of the data.
 * \param image The image to save
 * \param filename The name to the file to save the image to.
 * \return A status code, zero for success.
**/
int i3_image_save_binary(i3_image *image, char * filename);

/**
 * Load an image from a text file with the format produced by i3_image_save_text
 * \param filename The name to the file to save the image to.
 * \return The loaded image of NULL if not found or wrong format.
**/
i3_image * i3_image_load_text(char * filename);

/**
 * Load an image from a fits file.
 * \param filename The name to the file to save the image to.
 * \return The loaded image of NULL if not found or wrong format.
**/
i3_image * i3_image_load_fits(char * filename);

/** 
 Convolution and other image processing functions
*/

/**
 * Convolve two images together to make a new one.  If one of the images is used repeatedly you should use i3_image_fourier instead to precompute that.
 * \param image1 The first image.
 * \param image2 The second image.
 * \return The convolution of the two images.
**/
i3_image * i3_image_convolve(i3_image * image1, i3_image * image2);

/** 
 *  Print first element of each row. This function is of no use and used for debugging purposes only.
**/
void i3_image_print_first_element_of_each_row(i3_image * image);

/**
 * Compute the fourier transform of an image.
 * DEPRECATED - use i3_image_fourier instead.
 * \param image The image to transform.
 * \return Pointer to the fourier transform array.
**/
i3_cpx * i3_image_fft(i3_image *image);

/**
 * Convolve an image with another in-place, using the other's precomputed fourier transform
 * DEPRECATED - use i3_image_convolve_fourier instead.
 * \param image The image to convolve
 * \param PSF_FT The precomputed fourier transform of the convolving function.
**/
void i3_image_convolve_fft(i3_image * image, i3_cpx * PSF_FT);


/** 
 * Shift the image so its zero is at new point.  Useful for point spread functions
 * where you want to shift the zero point to the middle of the image.
 * \param image The image to shift
 * \param center_x The x-coordinate of the pixel which is moved to (0,0)
 * \param center_y The y-coordinate of the pixel which is moved to (0,0)
**/
void i3_image_shift_center(i3_image * image, int center_x, int center_y);

/**
 * Stack a number of images together into one average image.
 * \param images An array of image to stack.
 * \param number_images The number of images in the stack.
 * \return The average of all the images.
**/
i3_image * i3_image_stack(i3_image ** images, int number_images);

/**
 * Compute the moments of an image.  See i3_image_moments for details of what is calculated.
 * \param image The image whose moments to get.
 * \param moments The output moments
**/
void i3_image_compute_moments(i3_image *image, i3_image_moments * moments);

/**
 * Compute the moments of an image after applying a gaussian weighting.  See i3_image_moments for details of what is calculated.
 * This method is iterative because the center of the weighting is set to the computed weighted image center at each iteration.
 * \param image The image whose moments to get.
 * \param weight_radius The one-sigma radius of the gaussian used to weight the image before computing the transforms  Should be a little larger than the galaxy image.
 * \param iterations The number of iterations used.  Five seems enough.
 * \param moments The output moments

**/
void i3_image_compute_weighted_moments(i3_image *image, i3_flt weight_radius, int iterations, i3_image_moments * moments);


/**
 * INTERNAL FUNCTION
 * Computes the quadrupole parts of an image.  This is for internal use by i3_image_compute_moments or i3_image_compute_weighted_moments
 * \param image The image whose moments to get.
 * \param moments The input/output moments - expect x0,y0,sum to be already filled.  Fill in qij, e12 as well.
**/
void i3_image_compute_quadrupole(i3_image *image, i3_image_moments * moments);

/**
 * Fill in a weighted version of an image, with a given center and radius.
 * \param image The image that is to be weighted (will not be changed)
 * \param weighted The output weighted image (should be pre-allocated)
 * \param x0 The x-center of the weight functions
 * \param y0 The y-center of the weight functions
 * \param sigma The one-sigma width of the weight funtion radius.
**/
void i3_weight_image(i3_image * image, i3_image * weighted, i3_flt x0, i3_flt y0, i3_flt sigma);


/**
 * Fill an image with a constant
 * \param image The image to fill.
 * \param fill_value The value to fill the image with.
**/
void i3_image_fill(i3_image * image, i3_flt fill_value);

/**
 * Fill an ft of an image with a constant value, both real and imaginary parts
 * \param image The image to fill.
 * \param fill_value The value to fill the image with.
**/
void i3_fourier_fill(i3_fourier * ft, i3_cpx fill_value);


/**
 * Add white noise to the image - gaussian distributed noise with mean zero and fixed, uncorrelated variance.
 * \param image The image to add noise to. Will be modified.
 * \param noise_sigma The standard deviation of the noise to add.
**/
void i3_image_add_white_noise(i3_image * image, i3_flt noise_sigma);

/**
 * Add two images together into first image space.
 * \param image The image to add the second image to. Will be modified.
 * \param image2 The second image.
**/
void i3_image_add_image_into(i3_image * image, i3_image * image2);

/**
 * Add two images together into third image space.
 * \param image1 The image to add the second image to. Will be modified.
 * \param image2 The second image.
 * \param image_result result image that will contain the sum of image1 and image2
**/
void i3_image_add_images_into(i3_image * image1, i3_image * image2, i3_image * image_result);


/**
 * Add two images together into first image space, weighted by some coefficients.
 * \param image The image to add the second image to. Will be modified.
 * \param A1 The factor by which to multiply the first image
 * \param image2 The second image.
 * \param A2 The factor by which to multiply the second image
**/
void i3_image_weighted_add_image_into(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2);

/**
 * Add two images together into third image space, each weighted by some coefficients.
 * \param image1 The first image scaled with A1.
 * \param image2 The second image scaled with A2.
 * \param image3 The third image that will contain the weighted sum of image 1 and image 2.
**/
void i3_image_weighted_add_images_into(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2, i3_image * image_result);

/**
 * Add two images together into a new space, weighted by some coefficients.
 * \param image The image to add the second image to. 
 * \param A1 The factor by which to multiply the first image
 * \param image2 The second image.
 * \param A2 The factor by which to multiply the second image
 * \return The weighted image sum
**/
i3_image * i3_image_weighted_add_image(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2);

/**
 * Create a new version of the given image with the center shifted from the corner.
 * This is useful following Fourier Transforms.
 * \param image The input image to be centered (unmodified).
 * \return The re-centered image.
**/
i3_image * i3_image_centered(i3_image * image);

/**
 * Fills the second image with first image but with the center shifted from the corner.
 * This is useful following Fourier Transforms.
 * \param image The input image to be centered (unmodified).
 * \param result The output image, with the new center.
**/
void i3_image_centered_into(i3_image * image, i3_image * result);


/**
 * Create a new version of the image, rotated by 90 degrees.
 * I think only square images are currently used, but this is easy to remedy if required.
**/
i3_image * i3_image_rotate90(i3_image * image);

/**
 * Copy the first image into the second, but rotated by 90 degrees.
 * I think only square images are currently used, but this is easy to remedy if required.
**/
void i3_image_rotate90_into(i3_image * image, i3_image * result);

/**
 * Copy the first image into the second, but flipped left to right.
 * I think only square images are currently used, but this is easy to remedy if required.
**/
void i3_image_fliplr_into(i3_image * image, i3_image * result);

/**
Transposes the square image into an existing image.
*/
void i3_image_transpose_into(i3_image * image, i3_image * result);


/**
 * Scale the complete image by a constant factor.
**/
void i3_image_scale(i3_image * image, i3_flt x);

/**
 * Multiply one image by a second (mask) image
**/
 int i3_image_multiply_mask(i3_image * image, i3_image * mask);

/**
 * Add to the complete image a constant factor.
**/
void i3_image_addconst(i3_image * image, i3_flt x);

/**
 * Create a downsampled version of a high res image into a low res version.
**/
void i3_image_downsample_into(i3_image * high_resolution, i3_image * low_resolution);

/**
 * Create a downsampled version of a high res fourier image into a low res version.
**/
void i3_image_fourier_downsample_into(i3_fourier * high_resolution, i3_fourier * low_resolution);

/**
 * Create a downsampled version of a high res image in a new file with the given dimensions.
**/
i3_image * i3_image_downsample(i3_image * high_resolution, int nx, int ny);

/**
 * Create a padded version of the image.
**/
	i3_image * i3_image_padded_central(i3_image * image, int npad);
i3_image * i3_image_padded_corner(i3_image * image, int npad);

/**
 * Create a sheared copy of an existing image using a bicubic interpolation scheme.  This scheme should be improved soon.
**/
i3_image  * i3_image_shear(i3_image * image, i3_flt x0, i3_flt y0, i3_flt e1, i3_flt e2);

/**
 * Create a sheared version of an existing image into a waiting image space.  This uses the same bicubic scheme as i3_image_shear.
**/
void i3_image_shear_into(i3_image * image, i3_image * sheared, i3_flt x0, i3_flt y0, i3_flt e1, i3_flt e2);

/**
 * Zero an image; set all its pixels to zero.
**/
void i3_image_zero(i3_image * image);

/**
 * Zero a fourier image; set all its pixels to zero.
**/
void i3_fourier_zero(i3_fourier * fourier_image);

/**
 * Cut off an image at a given circular radius from a central point.  Set it to zero beyond that radius.
**/
void i3_image_truncate_circular(i3_image * image,i3_flt x0,i3_flt y0,i3_flt radius);


int i3_image_nonzero_count(i3_image * image);

/**
* Compute the sum of all the pixels in an image.
*/
i3_flt i3_image_sum(i3_image * image);
i3_flt i3_image_sum_masked(i3_image * image, i3_image * mask);

/**
 * Compute the mean of all the pixels in an image.
**/
i3_flt i3_image_mean(i3_image * image);
i3_flt i3_image_mean_masked(i3_image * image, i3_image * mask);

/**
 * Compute the standard deviation of all the pixels in an image.
**/
i3_flt i3_image_stdev(i3_image * image);

/**
 *  Extract a (sparse) regular grid of pixels from an image into another existing image.
**/
void i3_image_extract_into(i3_image * image, int x0, int y0, int dx, int dy, i3_image * into);

/** 
 * FFTW saves information about how to quickly perform transforms of a given size as "wisdom"
 * which we save to this file by default.
**/
#define I3_FFTW_WISDOM_FILENAME ".i3_fftw_wisdom.dat"

/**
 * Save any accumulated wisdom about how to quickly perform FFTs to the standard file.
 * You can use the function atexit(i3_fftw_save_wisdom); just after loading your wisdom to ensure that
 * any new wisdom is saved to file when the program ends.
**/
void i3_fftw_save_wisdom();

/**
 Load saved FFTW wisdom from the standard file.  This speeds up planning the first Fourier transform you do by using saved information about what method is fastest.
*/
void i3_fftw_load_wisdom();

/**
 * Allocate an instance of an i3_fourier, based on real pixel dimensions 
*/
i3_fourier * i3_fourier_create(unsigned long nx_real, unsigned long ny_real);

/**
 * Compute the fourier transform of an image, as an i3_image instance.
**/
i3_fourier * i3_image_fourier(i3_image *image);

/**
 * Destroy and deallocate an i3_fourier.
**/
void i3_fourier_destroy(i3_fourier * fourier);

/**
 * In-place convolve an image with a fourier kernel (for example, a saved PSF fourier transform.)
**/
void i3_image_convolve_fourier(i3_image * image, i3_fourier * kernel);

/**
 * Get an image representing the real part of an i3_fourier instance.
**/
i3_image * i3_fourier_real_part(i3_fourier * fourier);

/**
 * Get an image representing the imaginary part of an i3_fourier instance.
**/
i3_image * i3_fourier_imaginary_part(i3_fourier * fourier);

/**
 * Get an image representing the absolute value i3_fourier instance.
**/
i3_image * i3_fourier_absolute(i3_fourier * fourier);

/**
 * Perform a combination convolution and downsampling of a high-resolution image into a low-res one.
 * The i3_fourier kernel should me the product of a 
**/
void i3_image_convolve_downsample_into(i3_image * high_res, i3_fourier * kernel, i3_image * low_res);

/**
 * Downsample a fourier a image. Downsampling in fourier space is equivalent to cutting out high
 * ux and uy values. In the fft convention, These values are towards the centre.
 * e.g. to downsample from an 4x4 to 2x2 grid:
 *
 *  0,3 1,3 2,3 3,3
 *  0,2 1,2 2,2 3,2   ---\  0,3 3,3       here i,j denotes F(i,j)
 *  0,1 1,1 2,1 3,1   ---/  0,0 3,0
 *  0,0 1,0 2,0 3,0
 *
**/
void i3_image_fourier_downsample_into(i3_fourier * high_resolution, i3_fourier * low_resolution);

/**
 * Apply a wiener filter to an image based on some model image and the assumed noise standard deviation.
**/
void i3_image_wiener_filter(i3_image * image, i3_image * model, i3_flt noise);


/**
 * Downsample a source image into a destination image.
**/
void i3_image_dsample_into( i3_image * img_src, i3_image * img_dst );
void i3_image_dsample_cut_into( i3_image * img_src, i3_image * img_dst, int x_start, int y_start, int n_sub );

i3_image * i3_image_dsample( i3_image * img_src, int n_sub);
void i3_fourier_dsample_into( i3_fourier * ft_src,  i3_fourier * ft_dst);


/**
broken
**/
void i3_convolve_real_same_into( i3_image * image, i3_image * kernel, i3_image * result);


/**
 * Compute the coordinates and value of the brightest pixel
**/
i3_flt i3_image_max_index(i3_image *image, int * x, int * y);

void i3_multiply_fourier_fourier_into(i3_fourier * f1, i3_fourier * f2,i3_fourier * result);


i3_image * i3_image_cut_stamp(i3_image * image, int stamp_size);

i3_image * i3_fourier_image(i3_fourier *fourier);


/**
 * Convert image coordinates into pixel coordinates.  
 * Depending on ordering and flipping conventions there can be multiple ways to 
 * go from an image x/y position to a pixel i/j or vice versa.  These conventions
 * are collected here in this function, which at the moment has a fixed simple form.
 * \param image The image whose coordinates are desired
 * \param i The i coordinate
 * \param j The j coordinate
 * \param x (output) The x coordinate
 * \param y (output) The y coordinate
**/
void i3_image_xy_from_ij(i3_image * image, int i, int j, i3_flt * x, i3_flt *y);

/**
 * Convert image subpixel indices into coordinates.  
 * Depending on ordering and flipping conventions there can be multiple ways to 
 * go from an image x/y position to a pixel i/j or vice versa.  These conventions
 * are collected here in this function, which at the moment has a fixed simple form.
 * \param image The image whose coordinates are desired
 * \param i The i coordinate
 * \param j The j coordinate
 * \param i_sub The subpixel i coordinate
 * \param j_sub The subpixel j coordinate
 * \param x (output) The x coordinate
 * \param y (output) The y coordinate
**/
void i3_image_xy_from_ij_subpix(i3_image * image, int i, int j, int i_sub, int j_sub, int n_s, i3_flt * x, i3_flt *y);

/**
 * Inverse of i3_image_xy_from_ij
**/
void i3_image_ij_from_xy(i3_image * image, i3_flt x, i3_flt y, int * i, int *j);

/**
 * Get ux and uy values from i and j. Note the convention is that ux and uy go up to (
 * (but not equal to) 0.5 and then become negative. 
 * e.g. in 1d, for i = [0,1,2,3,4], kx = [0,0.2,0.4,-0.4,-0.2]
 **/
void i3_image_uxuy_from_ij(int i, int j, int nx, int ny, i3_flt * ux, i3_flt * uy);

/**
 * Reverse of i3_image_uxuy_from_ij
 **/
void i3_image_ij_from_uxuy(i3_flt ux, i3_flt uy, int nx, int ny, int * i, int *j);

/**
 * Get the corresponding i and j positions in an upsampled fourier image from i and j
 * in the downsampled fourier image. Need this for downsampling a fourier image.
**/
void i3_image_fourier_upsampled_ij_from__downsampled_ij(int down_i, int down_j, int up_nx, int up_ny, int down_nx, int down_ny, int * up_i, int * up_j);

/**
 * Generate the unit shear matrix components (M11,M12,M13)
 * Computes the matrix elements that can be applied to an (x,y) position vector
 * to generate the position in the sheared frame.  The matrix is symmetric,
 * so the M21=M12 component is not used.
 * \param e1 The shear in the x-y directions
 * \param e2 The shear in the diagonal directions
 * \param m1 (output) The M11 matrix element
 * \param m2 (output) The M12 matrix element
 * \param m3 (output) The M22 matrix element
**/
void i3_unit_shear_matrix(i3_flt e1, i3_flt e2, i3_flt * m1, i3_flt * m2, i3_flt * m3);

/**
 * Generate the unit shear matrix components (M11,M12,M13) IN FOURIER SPACE
 * Computes the matrix elements that can be applied to a (ux,uy) fourier position vector
 * to generate the position in the sheared frame.  The matrix is symmetric,
 * so the M21=M12 component is not used.
 * \param e1 The shear in the x-y directions
 * \param e2 The shear in the diagonal directions
 * \param m1 (output) The M11 matrix element
 * \param m2 (output) The M12 matrix element
 * \param m3 (output) The M22 matrix element
**/
void i3_unit_shear_matrix_fourier(i3_flt e1, i3_flt e2, i3_flt * m1, i3_flt * m2, i3_flt * m3);


void i3_unit_shear_matrix_jac(i3_flt e1, i3_flt e2, i3_flt * dm1_e1, i3_flt * dm2_e1, i3_flt * dm3_e1, i3_flt * dm1_e2, i3_flt * dm2_e2, i3_flt * dm3_e2);

/**
 * Compute the mean location of the image flux
 * The x mean is defined by x_mu = Sum_p flux_p x_p
 * and the y mean similarly.  The sum is over the image pixels.
 * The moments fill in only the x0 and y0 components of the structure passed in.
 * \param image The image to use
 * \param moments A moments structure to fill in
**/
void i3_image_compute_center(i3_image * image, i3_image_moments * moments);

/**
 * Cut out a postage stamp from a larger image (will fail if image is smaller).
 * If the image is the same size as the postage stamp returns a copy.
 * Always tries to centre the output postage stamp on the input image centre.
 **/
i3_image * i3_image_cut_out_stamp(i3_image * image, int stamp_size);


void i3_image_circular_mask(i3_image * image);


/**
 * Convert (e1, e2) ellipticity into (|e|, \theta) coordinates.
 * Uses e1 = |e|cos(2\theta), e2 = |e|sin(2\theta)
 **/
void i3_e12_to_etheta(i3_flt e1, i3_flt e2, i3_flt * e, i3_flt * theta);

/**
 * Convert (|e|, \phi) ellipticity into (e1, e2) coordinates.
 * Uses e1 = |e|cos(2\theta), e2 = |e|sin(2\theta)
 **/
void i3_etheta_to_e12(i3_flt e, i3_flt theta, i3_flt * e1, i3_flt * e2);

/**
 * Converts (e1, e2) in BEERMAT coordinates, for a beermat of extent emax,
 * to an (e1, e2) on the unit disc.
 * These map (one way only) an unbounded 2D plane to the closed outer 
 * surface of a beermat topology (representing the bound unit disc of e1, e2 
 * space).
 **/
void i3_e12_to_e12beermat(i3_flt e1in, i3_flt e2in, i3_flt emax, i3_flt * e1out, i3_flt * e2out);

/**
 * Converts an input radius rin to a "beermat" output value rout defined as 
 *
 *  rout = { |rin|                        for |rin| >= rparam
 *         { 0.5*rparam + c * |rin|^2     for |rin| < rparam
 * 
 * where rparam / 2 is our minimum allowed value of the beermat radius, and rparam is the radius 
 * below which the parabola kicks in.  This function has (by construction)
 *
 *  d(rout)/d(rin) = 0, rout = .5 * rparam, at rin = 0;
 *
 * and 
 *
 *  d(rout)/d(rin) = 1, rout = rparam, at rin = rparam.
 *
 * The by-hand derived solution for c is in fact then simply:
 *
 * c =  1 / (2 * rparam)
 *
 *\param rin - input radius, can be any real number
 *\param rparam - beermat r parameter, sets threshold below which the parabolic function is used
 *\returns rout - the new beermat output radius
 **/
void i3_r_to_rbeermat(i3_flt rin, i3_flt rparam, i3_flt * rout);

/**
* Calculates image FWHM by creating a 1D profile crossection of the image.
* \param image - image
* \param h - use -0.5 for FWHM, for 1/4 height use 0.25, etc.
* \param x0 - x location of the center of the object.   
 **/
i3_flt i3_image_get_fwhm(i3_image * image, i3_flt h, i3_flt x0);

/**
* Reweights weighting map with segmentation mask by
* setting all pixels to zero in the weight_stamp that are
* assigned to a different object in the segmentation mask
 **/
void i3_modify_weight_map_by_segmentation_mask(gal_id identifier, i3_image * weight_stamp, i3_image * mask_stamp);

void i3_modify_weight_map_by_segmentation_mask_nearest_pixels(gal_id identifier, i3_image * weight_stamp, i3_image * mask_stamp, int stamp_size);

void i3_ellipticity_sum(i3_flt e1, i3_flt e2, i3_flt g1, i3_flt g2, i3_flt * out1, i3_flt * out2);

i3_flt i3_image_subtract_background(i3_image * galaxy);
i3_flt i3_image_subtract_weighted_background(i3_image * galaxy, i3_image * weight);

void i3_e_linear_to_quadratic(i3_flt e1_lin, i3_flt e2_lin, i3_flt * e1_quad, i3_flt * e2_quad);
void i3_e_quadratic_to_linear(i3_flt e1_quad, i3_flt e2_quad, i3_flt *e1_lin, i3_flt *e2_lin);


/** 
* Find the location where an image pixel has exactly the given value
**/
int i3_image_value_location(i3_image * image, i3_flt value, int * nx, int * ny);

/** 
* Find the location (x,y) coordinates of the maximum value of the image
**/
int i3_image_maxloc(i3_image * image, int * nx, int * ny);

/** 
* Find the location (x,y) coordinates of the minimum value of the image
**/
int i3_image_minloc(i3_image * image, int * nx, int * ny);


/** 
* Find the standard FWHM of an image - through the peak, half max.
**/
i3_flt i3_image_standard_fwhm(i3_image * image);

/**
 * Copy a collection of tiles next to each other into an image
 **/
int i3_tile_images(i3_image * target, i3_image ** tiles, int n);

i3_image * i3_create_tiled_image(int n, i3_image ** images);


#endif