#include "i3_image.h"
#include "stdlib.h"
#include "string.h"
#include "float.h"
#include "i3_math.h"
#include "i3_image_fits.h"
#include "i3_load_data.h"
#ifdef _OPENMP
#include <omp.h>
#endif

#ifdef I3_USE_DOUBLE
const int i3UseDouble = 1;
#else
const int i3UseDouble = 0;
#endif

const char* i3ListOfModels=I3_LIST_OF_MODELS;

i3_image * i3_image_like(i3_image * image){
	return i3_image_create(image->nx, image->ny);
}

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) 
{
	*x = i + (i_sub+0.5)/n_s;
	*y = j + (j_sub+0.5)/n_s;
}

void i3_image_xy_from_ij(i3_image * image, int i, int j, i3_flt * x, i3_flt *y)
{
	*x = 0.5 + i;
	*y = 0.5 + j;
}

void i3_image_ij_from_xy(i3_image * image, i3_flt x, i3_flt y, int * i, int *j)
{
	*i = (int) x;
	*j = (int) y;
}

void i3_image_uxuy_from_ij(int i, int j, int nx, int ny, i3_flt * ux, i3_flt * uy)
{
	if (i >= ((i3_flt)nx)/2) *ux = ((i3_flt)i)/((i3_flt)nx) - 1.0;
    else *ux = ((i3_flt)i)/((i3_flt)nx);
    if (j >= ((i3_flt)ny)/2) *uy = ((i3_flt)j)/((i3_flt)ny) - 1.0;
    else *uy = ((i3_flt)j)/((i3_flt)ny);
}

void i3_image_ij_from_uxuy(i3_flt ux, i3_flt uy, int nx, int ny, int * i, int *j)
{
	if (ux >= 0) *i = (int)(ux * nx);
    else *i = (int)((ux + 1) * nx);
    if (uy >= 0) *j = (int)(uy * ny);
    else *j = (int)((uy + 1.) * ny);
}

// Not finished
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)
{
    if ((down_i < ((i3_flt)down_nx)/2) && (down_j < ((i3_flt)down_ny)/2)){
    	*up_i = down_i;
    	*up_j = down_j;
    } 
    if ((down_i >= ((i3_flt)down_nx)/2) && (down_j < ((i3_flt)down_ny)/2)){
    	*up_i = up_nx - (down_nx - down_i);
    	*up_j = down_j;
    }
    if ((down_i < ((i3_flt)down_nx)/2) && (down_j >= ((i3_flt)down_ny)/2)){
    	*up_i = down_i;
    	*up_j = up_ny - (down_ny - down_j);
    }	
    if ((down_i >= ((i3_flt)down_nx)/2) && (down_j >= ((i3_flt)down_ny)/2)){
    	*up_i = up_nx - (down_nx - down_i);
        *up_j = up_ny - (down_ny - down_j);
    }
}

i3_image * i3_image_create(unsigned long nx, unsigned long ny){
	/*
	Create space for an image.  Set up its main data chunk and row pointers.
	Nothing in the image is intialized.
	*/
	unsigned long n = nx*ny;
	i3_image * image = (i3_image*)malloc(sizeof(i3_image));
	/* Check for allocation error. This amount is so small that it is unlikely to be an out-of-mem error.*/
	if (!image) return NULL;


	/* Set up main data buffer and dimensions */
	image->n=n;
	image->nx=nx;
	image->ny=ny;
	image->data = (i3_flt*)malloc(sizeof(i3_flt)*n);
	
	/* Check for mem error.  Out-of-memory error is possible here.*/
	if (!image->data){
		free(image);
		return NULL;
	}

	/* Set up row pointer array*/
	image->row = (i3_flt**)malloc(ny*sizeof(i3_flt*));
	if (!image->row){
		/* Check for mem error.  Out-of-memory error is unlikely here.*/
		free(image);
		free(image->data);
		return NULL;
	}
	/* Set up row pointers so they each point to the correct row */
	for(unsigned long i=0;i<ny;i++) {image->row[i]=image->data + (i*nx);}
	return image;
}

/*  Allocate space for a fourier transform. */
i3_fourier * i3_fourier_create(unsigned long nx_real, unsigned long ny_real){
	unsigned long n_real = nx_real*ny_real;
	i3_fourier * fourier = (i3_fourier*)malloc(sizeof(i3_fourier));
	/* Check for allocation error. This amount is so small that it is unlikely to be an out-of-mem error.*/
	if (!fourier) return NULL;


	/* Set up main data buffer and dimensions */
	fourier->n_real=n_real;
	fourier->nx_real=nx_real;
	fourier->ny_real=ny_real;
	
	unsigned long nx_f = nx_real/2+1;
	unsigned long ny_f = ny_real;
	unsigned long n_f =  nx_f * ny_f;
	
	fourier->n_f=n_f;
	fourier->nx_f=nx_f;
	fourier->ny_f=ny_f;
	
	fourier->data = (i3_cpx*)malloc(sizeof(i3_cpx)*n_f);
	
	/* Check for mem error.  Out-of-memory error is possible here.*/
	if (!fourier->data){
		free(fourier);
		return NULL;
	}

	/* Set up row pointer array*/
	fourier->row = (i3_cpx**)malloc(ny_f*sizeof(i3_cpx*));
	if (!fourier->row){
		/* Check for mem error.  Out-of-memory error is unlikely here.*/
		free(fourier);
		free(fourier->data);
		return NULL;
	}
	/* Set up row pointers so they each point to the correct row */
	for(unsigned long i=0;i<ny_f;i++) {fourier->row[i]=fourier->data + (i*nx_f);}
	//for(unsigned long i=0;i<nx_f;i++) {fourier->row[i]=fourier->data + (i*nx_f);}
	for(unsigned long i=0;i<n_f;i++) fourier->data[i]=0;
	
	return fourier;
}

/*
* Free the memory associated with an i3_image. 
*/
void i3_image_destroy(i3_image * image){
	if (!image) return;
	if (image->row) free(image->row);
	if (image->data) free(image->data);
	free(image);
}

void i3_fourier_destroy(i3_fourier * fourier){
	if (!fourier) return;
	if (fourier->row) free(fourier->row);
	if (fourier->data) free(fourier->data);
	free(fourier);
}

// Copy the whole of the image into part of a new one.
int i3_image_copy_into_part(i3_image * image, long x_start, long y_start, i3_image * superimage){
	long width = image->nx;
	long height = image->ny;
	if (x_start+width>superimage->nx ||  y_start+height>superimage->ny) 
		return I3_IMAGE_FAIL;
	for(int row=0;row<height;row++) memcpy(superimage->row[y_start+row]+x_start, image->row[row], width*sizeof(i3_flt));
	return I3_OK;
}

int i3_image_copy_part_into(i3_image * image, long x_start, long y_start, i3_image * subimage){
	unsigned long width = subimage->nx;
	unsigned long height = subimage->ny;
	
	/* Check that the sub-image does not  exceed the bounds of the image*/
	if (x_start<0 || y_start<0) return I3_IMAGE_FAIL;
	if (x_start+width>image->nx ||  y_start+height>image->ny) return I3_IMAGE_FAIL;
	/* Copy the data row by row */
	for(int row=0;row<height;row++) memcpy(subimage->row[row],image->row[y_start+row]+x_start,width*sizeof(i3_flt));
	return I3_OK;
}




i3_image * i3_image_copy_part(i3_image * image, long x_start, long y_start, unsigned long width, unsigned long height){
	i3_image * subimage = i3_image_create(width,height);
	int fail = i3_image_copy_part_into(image,x_start, y_start, subimage);
	if (fail){
		i3_image_destroy(subimage);
		i3_print(i3_verb_noisy, "Failed to copy part of an image - full image:%ldx%ld  sub image: %ldx%ld,  start: %ld, %ld", image->nx, image->ny, width, height, x_start, y_start);
		return NULL;
	}
	return subimage;
}

i3_image * i3_image_copy(i3_image * image){
	return i3_image_copy_part(image,0,0,image->nx,image->ny);
}

int i3_image_copy_into(i3_image * image, i3_image * clone){
	return i3_image_copy_part_into(image,0,0,clone);
}

void i3_image_copy_from_pointer(i3_image * image, i3_flt * data){
	size_t nbytes = image->n*sizeof(i3_flt);
	memcpy(image->data, data, nbytes);
}


int i3_fourier_copy_into( i3_fourier * ft_src, i3_fourier * ft_dst ){
		
	/* Check that the sub-image does not  exceed the bounds of the image*/
	if ( ft_src->nx_f !=  ft_dst->nx_f ||  ft_src->ny_f !=  ft_dst->ny_f ) I3_FATAL("source and destination fourier have to be of same size",0);

	for(int i=0;i<ft_dst->n_f;i++) ft_dst->data[i] = ft_src->data[i];
	//memcpy(ft_dst->data[0],ft_src->data[0],ft_src->n_f*sizeof(i3_cpx));
	return I3_OK;
	

}

/*
Compute the maximum pixel value in the image.
*/
i3_flt i3_image_max(i3_image *image){
	i3_flt max_value = FLT_MIN;
	for (int p=0;p<image->n;p++) if (image->data[p]>max_value) max_value=image->data[p];
	return max_value;
}

i3_flt i3_image_min(i3_image *image){
	i3_flt min_value = FLT_MAX;
	for (unsigned int p=0;p<image->n;p++) if (image->data[p]<min_value) min_value=image->data[p];
	return min_value;
}

i3_flt i3_image_max_masked(i3_image *image, i3_image * mask){
	if (mask==NULL) return i3_image_max(image);
	// If mask is zero, ignore value
	i3_flt max_value = FLT_MIN;
	for (int p=0;p<image->n;p++) if ((mask->data[p]!=0) && (image->data[p]>max_value)) max_value=image->data[p];
	return max_value;
}

i3_flt i3_image_min_masked(i3_image *image, i3_image * mask){
	if (mask==NULL) return i3_image_min(image);
	i3_flt min_value = FLT_MAX;
	for (unsigned int p=0;p<image->n;p++) if ((mask->data[p]!=0) && (image->data[p]<min_value)) min_value=image->data[p];
	return min_value;
}

int i3_image_value_location(i3_image * image, i3_flt value, int * nx, int * ny)
{
	*nx = 0;
	*ny = 0;
	
	for (int j=0; j<image->ny; j++){
		for (int i=0; i<image->nx; i++){
			if (image->row[j][i]==value){
				*nx = i;
				*ny = j;
				return 1;
			}
		}
	}
	return 0;
}

int i3_image_maxloc(i3_image * image, int * nx, int * ny)
{
	i3_flt max_value = i3_image_max(image);
	return i3_image_value_location(image, max_value, nx, ny);
}

int i3_image_minloc(i3_image * image, int * nx, int * ny)
{
	i3_flt min_value = i3_image_min(image);
	return i3_image_value_location(image, min_value, nx, ny);
}




i3_flt i3_image_max_index(i3_image *image, int * x, int * y){
	i3_flt max_value = FLT_MIN;
	int max_index = -1;
	for (unsigned int p=0;p<image->n;p++) if (image->data[p]>max_value) {
		max_value=image->data[p];
		max_index = p;
	}
	*x = max_index/image->nx;
	*y = max_index%image->ny;
	return max_value;
}


i3_flt i3_image_sum_masked(i3_image * image, i3_image * mask){
	if (mask==NULL) return i3_image_sum(image);
	i3_flt sum = 0.0;
	for (unsigned int p=0;p<image->n;p++) if (mask->data[p]!=0) sum += image->data[p];
	return sum;
}


i3_flt i3_image_sum(i3_image * image){
	i3_flt sum = 0.0;
	for (unsigned int p=0;p<image->n;p++) sum += image->data[p];
	return sum;
}

int i3_image_nonzero_count(i3_image * image)
{	
	int count = 0;
	for (unsigned int p=0;p<image->n;p++) if (image->data[p]!=0) count++;
		return count;
}

i3_flt i3_image_mean(i3_image * image){
	return i3_image_sum(image)/image->n;
}

i3_flt i3_image_mean_masked(i3_image * image, i3_image * mask){
	return i3_image_sum_masked(image, mask)/i3_image_nonzero_count(mask);
}


i3_flt i3_image_stdev(i3_image * image){
	i3_flt mu = i3_image_mean(image);
	i3_flt std=0.0;
	for (unsigned int p=0;p<image->n;p++) std += i3_pow(image->data[p]-mu,2);
	std/=image->n;
	return std;
}

i3_flt i3_image_norm(i3_image * image){
	i3_flt norm = 0.0;
	for (unsigned int p=0;p<image->n;p++) norm += i3_pow(image->data[p],2);
	return i3_sqrt(norm);
}

i3_flt i3_image_norm_sq(i3_image * image){
	i3_flt norm = 0.0;
	for (unsigned int p=0;p<image->n;p++) norm += i3_pow(image->data[p],2);
	return norm;
}

void i3_image_extract_into(i3_image * image, int i0, int j0, int di, int dj, i3_image * into){
	for (int j=0;j<into->ny;j++){
		for (int i=0;i<into->nx;i++){
			if (i0+i*di<0 || i0+i*di>=image->nx) I3_FATAL("Image extraction bad (x)!",(int) i0+i*di);
			if (j0+j*dj<0 || j0+j*dj>=image->ny) I3_FATAL("Image extraction bad (y)!",(int) j0+j*dj);
			into->row[j][i]=image->row[j0+j*dj][i0+i*di];
		}
	}
}

int i3_image_save_text(i3_image *image, char * filename){
	FILE * output_file = fopen(filename,"w");
	if (!output_file) return I3_PATH_ERROR;
	for(int j=0;j<image->ny;j++){
		for(int i=0;i<image->nx;i++) fprintf(output_file,"%e\t",image->row[j][i]);
		fprintf(output_file,"\n");
	}
	fclose(output_file);
	return I3_OK;
}

int i3_count_file_lines(char * filename){
	FILE * f = fopen(filename,"r");
	if (!f) return 0;
	int count=0;
	while(!feof(f)){
		int c = fgetc(f);
		if (c=='\n') count++;
	}
	fclose(f);
	return count;
}

int i3_line_length(char * filename){
	FILE * f = fopen(filename,"r");
	if (!f) return 0;
	int count=0;
	int c=0;
	while(c!='\n'){
		c = fgetc(f);
		count++;
	}
	fclose(f);
	return count;	
}

int i3_count_line_words(char * filename){
	int line_chars = i3_line_length(filename)+100; //extra just in case
	FILE * f = fopen(filename,"r");
	if (!f) return 0;
	char * line = malloc(line_chars);
	fgets(line,line_chars,f);
  	char * pch = strtok (line," \t");
	int count=0;
	while (pch != NULL){
		count++;
    	pch=strtok(NULL, " \t");
  }
	//count--;
	fclose(f);
	return count;

}

i3_image * i3_image_load_text(char * filename){
	
	int nx = i3_count_line_words(filename);
	if (!nx) return NULL;
	int ny = i3_count_file_lines(filename);
	if (!ny) return NULL;

	FILE * input_file = fopen(filename,"r");
	if (!input_file) return NULL;

	i3_image * image = i3_image_create(nx,ny);
	for (int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
#ifdef I3_USE_DOUBLE
			fscanf(input_file,"%lf",image->row[j]+i);
#else
			fscanf(input_file,"%f",image->row[j]+i);
#endif
		}
	}
	fclose(input_file);
	return image;
}

i3_image * i3_image_load_fits(char * filename){
	return i3_read_fits_image(filename);
}

int i3_image_save_binary(i3_image *image, char * filename){
	FILE * output_file = fopen(filename,"wb");
	if (!output_file) return I3_PATH_ERROR;
	fwrite(image->data,sizeof(i3_flt),image->n,output_file);
	fclose(output_file);
	return I3_OK;
}

/*
Convolve two images together using Fourier transforms.
If you want to convolve one image repeatedly with many others (e.g. with a point-spread function)
Then for speed you should use i3_convolve_fourier

*/
i3_image * i3_image_convolve(i3_image * image1, i3_image * image2){
	/*
	Check that images have the same dimension and are square.
	*/
	
	int n = image1->nx;
	int ncomplex = n/2+1;

	i3_cpx * complex_space1 = malloc(n*ncomplex*sizeof(i3_cpx));
	i3_cpx * complex_space2 = malloc(n*ncomplex*sizeof(i3_cpx));

	i3_image * convolved_image = i3_image_create(image1->nx, image1->ny);

#ifdef I3_FFTW_USE_DOUBLE

	double * real_space = malloc(n*n*sizeof(double));
	fftw_plan plan1 =  fftw_plan_dft_r2c_2d(n,n,real_space, complex_space1,FFTW_MEASURE);
	fftw_plan plan2 =  fftw_plan_dft_r2c_2d(n,n,real_space, complex_space2,FFTW_MEASURE);
	fftw_plan plan3 =  fftw_plan_dft_c2r_2d(n,n,complex_space1, real_space,FFTW_MEASURE);

	for(int i=0;i<n*n;i++) real_space[i]=image1->data[i];
	fftw_execute(plan1);
	for(int i=0;i<n*n;i++) real_space[i]=image2->data[i];
	fftw_execute(plan2);
	for(int i=0;i<n*ncomplex;i++) complex_space1[i]*=complex_space2[i];
	fftw_execute(plan3);
	for(int i=0;i<n*n;i++) convolved_image->data[i] = real_space[i];
	
	i3_flt factor = 1.0 / (1.0* image1->nx * image1->ny);
	i3_image_scale(convolved_image,factor);
	
	fftw_destroy_plan(plan1);
	fftw_destroy_plan(plan2);
	free(real_space);

#else

	fftwf_plan plan1 =  fftwf_plan_dft_r2c_2d(n,n,image1->data, complex_space1,FFTW_MEASURE);
	fftwf_plan plan2 =  fftwf_plan_dft_r2c_2d(n,n,image2->data, complex_space2,FFTW_MEASURE);
	fftwf_plan plan3 =  fftwf_plan_dft_c2r_2d(n,n,complex_space1,convolved_image->data,FFTW_MEASURE);

	fftwf_execute(plan1);
	fftwf_execute(plan2);
	for(int i=0;i<n*ncomplex;i++) complex_space1[i]*=complex_space2[i];
	fftwf_execute(plan3);
	
	i3_flt factor = 1.0 / (1.0* image1->nx * image1->ny);
	i3_image_scale(convolved_image,factor);
	
	fftwf_destroy_plan(plan1);
	fftwf_destroy_plan(plan2);

#endif

	free(complex_space1);
	free(complex_space2);
	return convolved_image;
	
	
}

void i3_image_shift_center(i3_image * image, int center_x, int center_y){
	/* 
	Shift an image so that its center moves to (0,0)
	*/
	i3_image * buffer = i3_image_like(image);
	int i2,j2;
	for(int j=0;j<image->ny;j++){
		j2 = (j+center_y)%image->ny;
		for(int i=0;i<image->nx;i++){
			i2 = (i+center_x)%image->nx;
			buffer->row[j][i]=image->row[j2][i2];
			
		}
	}

	for(int p=0;p<image->n;p++)image->data[p]=buffer->data[p];
	i3_image_destroy(buffer);

}

void i3_image_print_first_element_of_each_row(i3_image * image){
	/* 
	Print first element of each row. This function is of no use and used for debugging purposes only.
	*/
	for(int j=0;j<image->ny;j++){
	  fprintf(stdout, "Row[%d] %.16f.\n",j,image->row[j][0]);
	}
}

i3_cpx * i3_image_fft(i3_image *image){


	int nx = image->nx;
	int ny = image->ny;
	int ncomplex = ny/2+1;
	i3_cpx * complex_space = malloc(nx*ncomplex*sizeof(i3_cpx));

#ifdef I3_FFTW_USE_DOUBLE

	double * real_space = malloc(nx*ny*sizeof(double));
	fftw_plan fwd_transform =  fftw_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_PATIENT);
	for(unsigned long i=0;i<image->n;i++) real_space[i]=(double)(image->data[i]);
	fftw_execute(fwd_transform);
	fftw_destroy_plan(fwd_transform);	
	free(real_space);

#else

	fftwf_plan fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_WISDOM_ONLY);
	if (!fwd_transform){
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_EXHAUSTIVE);
		free(real_space);
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_WISDOM_ONLY);
		if (!fwd_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	fftwf_execute(fwd_transform);
	fftwf_destroy_plan(fwd_transform);	
#endif
	return complex_space;
}

i3_fourier * i3_image_fourier(i3_image *image){


	int nx = image->nx;
	int ny = image->ny;

#ifdef I3_FFTW_USE_DOUBLE
	double * real_space = malloc(nx*ny*sizeof(double));
	i3_fourier * fourier = i3_fourier_create(nx,ny);
	fftw_plan fwd_transform =  fftw_plan_dft_r2c_2d(nx,ny,real_space, fourier->data, FFTW_PATIENT);
	for(unsigned long i=0;i<image->n;i++) real_space[i]=(double)(image->data[i]);
	fftw_execute(fwd_transform);
	fftw_destroy_plan(fwd_transform);
	free(real_space);
#else
	//printf("we are NOT using I3_FFTW_USE_DOUBLE \n");
	//float * real_space = malloc(nx*ny*sizeof(float)); for(unsigned long i=0;i<image->n;i++) real_space[i]=(float)(image->data[i]);
	i3_fourier * fourier = i3_fourier_create(nx,ny);
	fftwf_plan fwd_transform;
	 
	fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, fourier->data, FFTW_WISDOM_ONLY);
	
	if (!fwd_transform){
		printf("i3_image_fourier: FT is learning the fwd transfrom of size %d x %d.\n",nx,ny);
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, fourier->data,FFTW_EXHAUSTIVE);
		fftwf_destroy_plan(fwd_transform);
		free(real_space);
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, fourier->data,FFTW_WISDOM_ONLY);
		if (!fwd_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	
	//fftwf_plan fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, fourier->data, FFTW_PATIENT);
	//for(unsigned long i=0;i<image->n;i++) real_space[i]=(float)(image->data[i]);
	fftwf_execute(fwd_transform);

	fftwf_destroy_plan(fwd_transform);

	//free(real_space);
#endif
	return fourier;
}



i3_image * i3_fourier_image(i3_fourier *fourier){


	int nx = fourier->nx_real;
	int ny = fourier->ny_real;

#ifdef I3_FFTW_USE_DOUBLE
	i3_image * image = i3_image_create(nx,ny);
	double * real_space = malloc(nx*ny*sizeof(double));
	fftw_plan rev_transform =  fftw_plan_dft_c2r_2d(nx,ny,fourier->data,real_space,FFTW_WISDOM_ONLY);
	if (!rev_transform){
		printf("i3_image_fourier: FT is learning inverse fwd transfrom of size %d x %d.\n",nx,ny);
		fftw_complex * complex_space = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fourier->n_f);
		rev_transform =  fftw_plan_dft_c2r_2d(nx,ny,complex_space,real_space,FFTW_EXHAUSTIVE);
		fftw_destroy_plan(rev_transform);
		rev_transform =  fftw_plan_dft_c2r_2d(nx,ny,fourier->data,real_space,FFTW_WISDOM_ONLY);
		if (!rev_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	
	fftw_execute(rev_transform);
	
	for(unsigned long i=0;i<image->n;i++) image->data[i] = real_space[i];
	fftw_destroy_plan(rev_transform);
	free(real_space);
#else
	i3_image * image = i3_image_create(nx,ny);
	fftwf_plan rev_transform =  fftwf_plan_dft_c2r_2d(nx,ny,fourier->data, image->data, FFTW_WISDOM_ONLY);
	
	if (!rev_transform){
		printf("i3_image_fourier: FT is learning the inverse transfrom of size %d x %d.\n",nx,ny);
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		rev_transform =  fftwf_plan_dft_c2r_2d(nx,ny,fourier->data,real_space,FFTW_EXHAUSTIVE);
		free(real_space);
		fftwf_destroy_plan(rev_transform);
		fftwf_plan rev_transform =  fftwf_plan_dft_c2r_2d(nx,ny,fourier->data, image->data, FFTW_WISDOM_ONLY);
		if (!rev_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	
	//fftwf_plan fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, fourier->data, FFTW_PATIENT);
	//for(unsigned long i=0;i<image->n;i++) real_space[i]=(float)(image->data[i]);
	fftwf_execute(rev_transform);
	
	fftwf_destroy_plan(rev_transform);
	
	//free(real_space);
#endif
	return image;
}





void i3_image_convolve_fft(i3_image * image, i3_cpx * PSF_FT){
	
	int nx = image->nx;
	int ny = image->ny;
	int ncomplex = ny/2+1;
	i3_cpx * complex_space = malloc(nx*ncomplex*sizeof(i3_cpx));

#ifdef I3_FFTW_USE_DOUBLE

	// set up
	double * real_space = malloc(nx*ny*sizeof(double));
	fftw_plan fwd_transform =  fftw_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_MEASURE);
	fftw_plan inv_transform =  fftw_plan_dft_c2r_2d(nx,ny,complex_space, real_space,FFTW_MEASURE);
	for(int i=0;i<nx*ny;i++) real_space[i]=image->data[i];

	// fwd
	fftw_execute(fwd_transform);

	// convolution
	for(int i=0;i<nx*ncomplex;i++) complex_space[i]*=PSF_FT[i];

	// inv fft	
	fftw_execute(inv_transform);
	i3_flt scaling = (1./image->nx)/image->ny;

	for(int i=0;i<nx*ny;i++) image->data[i] = (i3_flt)(real_space[i]*scaling);

	fftw_destroy_plan(fwd_transform);
	fftw_destroy_plan(inv_transform);
	free(real_space);


#else
	

	// set up
	fftwf_plan fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_MEASURE);
	fftwf_plan inv_transform =  fftwf_plan_dft_c2r_2d(nx,ny,complex_space, image->data,FFTW_MEASURE);
	
	// fwd
// 	fftwf_execute(fwd_transform);

	// convolution
	for(int i=0;i<nx*ncomplex;i++) complex_space[i]*=PSF_FT[i];

	// inv fft	
	fftwf_execute(inv_transform);
	i3_flt scaling = (1./image->nx)/image->ny;

	for(int i=0;i<nx*ny;i++) image->data[i] *= scaling;

	fftwf_destroy_plan(fwd_transform);
	fftwf_destroy_plan(inv_transform);

#endif

	free(complex_space);
	
}

void i3_image_convolve_fourier(i3_image * image, i3_fourier * kernel){
	if (image->nx!=kernel->nx_real){
		printf("\nImage: %ld x %ld\nKernel:%ld x %ld\n",image->nx,image->ny,kernel->nx_real,kernel->ny_real);
		I3_FATAL("Tried to convolve image with incompatibly sized kernel (x dimension)",1);
	}
	if (image->ny!=kernel->ny_real) {
		printf("\nImage: %ld x %ld\nKernel:%ld x %ld\n",image->nx,image->ny,kernel->nx_real,kernel->ny_real);
		I3_FATAL("Tried to convolve image with incompatibly sized kernel (y dimension)",1);
	}

	int nx = image->nx;
	int ny = image->ny;
	int nxf = kernel->nx_f;
	int nyf = kernel->ny_f;
	i3_cpx * complex_space = malloc(nxf*nyf*sizeof(i3_cpx));

#ifdef I3_FFTW_USE_DOUBLE

	// setup
	double * real_space = malloc(nx*ny*sizeof(double));
	fftw_plan fwd_transform;
	fftw_plan inv_transform;
	
	fwd_transform =  fftw_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_ESTIMATE);
	inv_transform =  fftw_plan_dft_c2r_2d(nx,ny,complex_space, real_space,FFTW_ESTIMATE);
	
	
	for(int i=0;i<nx*ny;i++) real_space[i]=image->data[i];
	
	// fwd transform
	fftw_execute(fwd_transform);
	
	//convolve
	for(int i=0;i<nxf*nyf;i++) complex_space[i]*=kernel->data[i];

	// inv fft and scale
	fftw_execute(inv_transform);
	i3_flt scaling = (1.f/image->nx)/image->ny;

	for(int i=0;i<nx*ny;i++) image->data[i] = (i3_flt)(real_space[i]*scaling);
	
	
	fftw_destroy_plan(fwd_transform);
	fftw_destroy_plan(inv_transform);
	
	
	
	free(real_space);
#else
	
	// setup
	fftwf_plan fwd_transform;
	fftwf_plan inv_transform;  
	
	
	fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_WISDOM_ONLY);
	inv_transform =  fftwf_plan_dft_c2r_2d(nx,ny,complex_space, image->data,FFTW_WISDOM_ONLY);
	

	if (!fwd_transform){
		printf("i3_image_fourier: FT is learning the fwd transfrom of size %d x %d.\n",nx,ny);
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_PATIENT);
		free(real_space);
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_WISDOM_ONLY);
		if (!fwd_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	if(!inv_transform){
		printf("i3_image_fourier: FT is learning the inverse transfrom of size %d x %d.\n",nx,ny);
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		inv_transform =  fftwf_plan_dft_c2r_2d(nx,ny,complex_space,real_space,FFTW_PATIENT);
		free(real_space);
		inv_transform =  fftwf_plan_dft_c2r_2d(nx,ny,complex_space,image->data,FFTW_WISDOM_ONLY);
		if (!inv_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}

		
	// fwd transform
	fftwf_execute(fwd_transform);
	
	//convolve
	for(int i=0;i<nxf*nyf;i++) complex_space[i]*=kernel->data[i];

	// inv fft and scale
	fftwf_execute(inv_transform);
	i3_flt scaling = (1.f/image->nx)/image->ny;
	i3_image_scale(image,scaling);

	
	fftwf_destroy_plan(fwd_transform);
	fftwf_destroy_plan(inv_transform);
	
#endif

	free(complex_space);

/*
	fftwf_plan fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_PATIENT|FFTW_WISDOM_ONLY);
	if (!fwd_transform){
		i3_flt * real_space = malloc(nx*ny*sizeof(i3_flt));
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,real_space, complex_space,FFTW_PATIENT);
		free(real_space);
		fwd_transform =  fftwf_plan_dft_r2c_2d(nx,ny,image->data, complex_space,FFTW_PATIENT|FFTW_WISDOM_ONLY);
		if (!fwd_transform) I3_FATAL("Joe and Tomek do not understand FFTW - complain to them.",1);
	}
	fftwf_execute(fwd_transform);
	fftwf_destroy_plan(fwd_transform);	*/

}





i3_image * i3_fourier_real_part(i3_fourier * fourier){
	i3_image * image = i3_image_create(fourier->nx_f, fourier->ny_f);
	for (int p=0;p<fourier->n_f;p++) image->data[p]=creal(fourier->data[p]);	
	return image;
}

i3_image * i3_fourier_imaginary_part(i3_fourier * fourier){
	i3_image * image = i3_image_create(fourier->nx_f, fourier->ny_f);
	for (int p=0;p<fourier->n_f;p++) image->data[p]=cimag(fourier->data[p]);
	return image;
}

i3_image * i3_fourier_absolute(i3_fourier * fourier){
	i3_image * image = i3_image_create(fourier->nx_f, fourier->ny_f);

	for (int p=0;p<fourier->n_f;p++) image->data[p]=i3_cabs(fourier->data[p]);

	return image;
}

i3_image * i3_image_stack(i3_image ** images, int number_images){
	/* Create output image from dimensions o first image. */
	i3_image * stacked_image = i3_image_create(images[0]->nx,images[0]->ny);
	
	/* Loop through images adding in data. */
	for(int i=0;i<number_images;i++){
		i3_image * image = images[i];
		for(int p=0;p<stacked_image->n;p++) stacked_image->data[p]+=image->data[p];
	}
	/* Divide to get the average.*/
	i3_image_scale(stacked_image,1.0/number_images);

	return stacked_image;
	
}

void i3_image_compute_center(i3_image * image, i3_image_moments * moments){
	/* Compute the zeroth moments of an image - the x0 and y0 means.*/
	i3_flt v;
	moments->x0=0.0;
	moments->y0=0.0;
	moments->sum=0.0;
	i3_flt x,y;
	for(int j=0;j<image->ny;j++){
		for(int i=0;i<image->nx;i++){
			v = image->row[j][i];
			i3_image_xy_from_ij(image, i, j, &x, &y);
			moments->sum += v;
			moments->x0  += x*v;
			moments->y0  += y*v;
		}
	}
	moments->x0/=moments->sum;
	moments->y0/=moments->sum;
}

void i3_weight_image(i3_image * image, i3_image * weighted, i3_flt x0, i3_flt y0, i3_flt sigma){
	/* Weight an image by a gaussian*/
	i3_flt x,y;
	i3_flt tau = -0.5/(sigma*sigma);
	for(int j=0;j<image->ny;j++){
		for(int i=0;i<image->nx;i++){
			i3_image_xy_from_ij(image, i, j, &x, &y);
			i3_flt dx2 = i3_pow(x-x0,2);
			i3_flt dy2 = i3_pow(y-y0,2);
			weighted->row[j][i] = image->row[j][i] * i3_exp(tau*(dx2+dy2));
		}
	}
	
}

void i3_image_compute_weighted_moments(i3_image *image, i3_flt weight_radius, int iterations, i3_image_moments * moments){
	i3_image * weighted_image = i3_image_copy(image);
	moments->x0 = image->nx/2.0;
	moments->y0 = image->ny/2.0;
	for (int i=0;i<iterations;i++){
		i3_weight_image(image,weighted_image,moments->x0,moments->y0,weight_radius);
		i3_image_compute_center(weighted_image,moments);
	}
	//i3_image_compute_quadrupole(image,moments);
	// SLB: I assume it is supposed to be the below? That's what I'd expect from a function with this title anyway
	i3_image_compute_quadrupole(weighted_image,moments);
	i3_image_destroy(weighted_image);
}

/*
Compute the quadrupole moment and ellipticity of an image.
NB - this function assumes that the first moment is already calculated.
*/
void i3_image_compute_quadrupole(i3_image *image, i3_image_moments * moments){
	i3_flt v;
	i3_flt qxx,qxy,qyy;
	
	qxx=0.0;
	qxy=0.0;
	qyy=0.0;
	i3_flt dx,dy,x,y;
	for(int j=0;j<image->ny;j++){
		for(int i=0;i<image->nx;i++){
			i3_image_xy_from_ij(image, i, j, &x, &y);
			v = image->row[j][i];
			dx=x-moments->x0;
			dy=y-moments->y0;
			qxx += v*dx*dx;
			qxy += v*dx*dy;
			qyy += v*dy*dy;
		}
	}
	qxx/=moments->sum;
	qxy/=moments->sum;
	qyy/=moments->sum;
	i3_flt ellipticity_denominator = qxx+qyy+2*i3_sqrt(qxx*qyy-qxy*qxy);

	moments->qxx=qxx;
	moments->qxy=qxy;
	moments->qyy=qyy;
	moments->e1 = (qxx-qyy)/ellipticity_denominator;
	moments->e2 = (2*qxy)/ellipticity_denominator;
	
	
}

void i3_image_compute_moments(i3_image *image, i3_image_moments * moments){
	i3_image_compute_center(image, moments);
	i3_image_compute_quadrupole(image, moments);
}

void i3_image_fill(i3_image * image, i3_flt fill_value){
	for(int p=0;p<image->n;p++)	image->data[p]=fill_value;
}

void i3_fourier_fill(i3_fourier * ft, i3_cpx fill_value){

	for(int i=0;i<ft->n_f;i++){
		ft->data[i] = fill_value;	
	}
}

void i3_image_zero(i3_image * image){
	memset(image->data,0,sizeof(i3_flt)*image->n);
}

void i3_fourier_zero(i3_fourier * fourier_image){
	memset(fourier_image->data,0,sizeof(i3_cpx)*fourier_image->n_f);
//	i3_image_fill(image,0.0);
}



void i3_image_add_white_noise(i3_image * image, i3_flt noise_sigma){
	for (int p=0;p<image->n;p++) image->data[p]+=i3_random_normal() * noise_sigma;
}

void i3_image_add_image_into(i3_image * image1, i3_image * image2){
	for (int p=0;p<image1->n;p++) image1->data[p]+=image2->data[p];
}

void i3_image_add_images_into(i3_image * image1, i3_image * image2, i3_image * image_result){
	for (int p=0;p<image1->n;p++) image_result->data[p] = image1->data[p]+image2->data[p];
}

void i3_image_weighted_add_image_into(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2){
	for (int p=0;p<image1->n;p++){
		image1->data[p]= A1*image1->data[p] + A2*image2->data[p];
	}
}

void i3_image_weighted_add_images_into(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2, i3_image * image3){
	for (int p=0;p<image1->n;p++){
		image3->data[p]= A1*image1->data[p] + A2*image2->data[p];
	}
}

i3_image * i3_image_weighted_add_image(i3_image * image1, i3_flt A1, i3_image * image2, i3_flt A2){
	i3_image * result = i3_image_create(image1->nx,image1->ny);
	for (int p=0;p<image1->n;p++){
		result->data[p]= A1*image1->data[p] + A2*image2->data[p];
	}
	return result;
}

i3_image * i3_image_centered(i3_image * image){
	i3_image * result = i3_image_create(image->nx,image->ny);
	i3_image_centered_into(image,result);
	return result;
}

void i3_image_centered_into(i3_image * image, i3_image * result){
	int nx = image->nx;
	int ny = image->ny;
	int i_new, j_new;
	int i_center, j_center;

	if(nx%2==1){ 
	 	i_center = nx/2+1;
		j_center = ny/2+1;
	}
	else{ 
	 	i_center = nx/2;
		j_center = ny/2;
	}

	for(int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			i_new = i_center + i;
			j_new = j_center + j;
			if (i_new>=nx) i_new-=nx;
			if (j_new>=ny) j_new-=ny;
			result->row[j_new][i_new]=image->row[j][i];
		}
	}
}

i3_image * i3_image_rotate90(i3_image * image){
	i3_image * result = i3_image_create(image->nx,image->ny);
	i3_image_rotate90_into(image, result);
	return result;
}

void i3_image_rotate90_into(i3_image * image, i3_image * result){
	I3_WARNING("NOT CERTAIN THAT i3_image_rotate90_into DEFINITELY WORKS RIGHT WAY - PLEASE CHECK");
	int nx = image->nx;
	int ny = image->ny;
	if (nx!=ny) {fprintf(stderr,"ERROR - tried to rotate non-square image - not coded yet.\n");}
	for(int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			result->row[j][i]=image->row[ny-i-1][j];
		}
	}
}

void i3_image_fliplr_into(i3_image * image, i3_image * result){
	I3_WARNING("NOT CERTAIN THAT i3_image_fliplr_into DEFINITELY WORKS RIGHT WAY - PLEASE CHECK");
	int nx = image->nx;
	int ny = image->ny;
	if (nx!=ny) {fprintf(stderr,"ERROR - tried to rotate non-square image - not coded yet.\n");}
	for(int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			result->row[j][i]=image->row[j][nx-i-1];
		}
	}
}

void i3_image_transpose_into(i3_image * image, i3_image * result){
	int nx = image->nx;
	int ny = image->ny;
	if (nx!=ny) {fprintf(stderr,"ERROR - tried to transpose non-square image - not coded yet.\n");}
	for(int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			result->row[j][i]=image->row[i][j];
		}
	}
	
}

void i3_image_scale(i3_image * image, i3_flt x){
	for (int p=0;p<image->n;p++) image->data[p]*=x;
}

int i3_image_multiply_mask(i3_image * image, i3_image * mask)
{
	if (image->nx!=mask->nx || image->ny!=mask->ny) return 1;
	for (int p=0;p<image->n;p++) image->data[p] *= mask->data[p];
	return 0;
}


void i3_image_addconst(i3_image * image, i3_flt x){
	for (int p=0;p<image->n;p++) image->data[p]+=x;
}

i3_image * i3_image_downsample(i3_image * high_resolution, int nx_low, int ny_low){
	i3_image * low_resolution = i3_image_create(nx_low,ny_low);
	i3_image_downsample_into(high_resolution,low_resolution);
	return low_resolution;
}

void i3_image_downsample_into(i3_image * high_resolution, i3_image * low_resolution){
	int nx_high = high_resolution->nx;
	int ny_high = high_resolution->ny;
	int nx_low = low_resolution->nx;
	int ny_low = low_resolution->ny;
	
	i3_flt fx = (nx_high*1.0)/nx_low;
	i3_flt fy = (ny_high*1.0)/ny_low;
	
	//I think that this always works and there are no floating point surprises.

	i3_image_zero(low_resolution);

	if (round(fx)==fx && round(fy)==fy){
		//Easy case; the number of pixels in the high-res image is an integer multiple of the number in the low-res.
		int rx = (int) fx;
		int ry = (int) fy;

		for (int j=0;j<ny_high;j++){
			for (int i=0;i<nx_high;i++){
				low_resolution->row[j/ry][i/rx] += high_resolution->row[j][i];
			}
		}
		i3_image_scale(low_resolution,1.0/(rx*ry));
	}
	else{
		I3_FATAL("Have not yet coded non-integer down-samplings.  See function i3_image_downsample_into", 1);
	}	
}

// not finished.
// not sure this is correct so far! 
void i3_image_fourier_downsample_into(i3_fourier * high_resolution, i3_fourier * low_resolution){
 
	int nx_f_low = low_resolution->nx_f; 
	int ny_f_low = low_resolution->ny_f;
	int nx_high = high_resolution->nx_real; 
	int ny_high = high_resolution->ny_real; 
	int nx_low = low_resolution->nx_real; 
	int ny_low = low_resolution->ny_real;

	i3_fourier_zero(low_resolution);

	for (int j_low=0; j_low<ny_f_low; j_low++){
		for (int i_low=0; i_low<nx_f_low; i_low++){
            int i_hi;
            int j_hi;
            i3_image_fourier_upsampled_ij_from__downsampled_ij(i_low, j_low, nx_high, ny_high, nx_low, ny_low, &i_hi, &j_hi);
            low_resolution->row[j_low][i_low] = high_resolution->row[j_hi][i_hi];

		}
	} 

}

i3_image * i3_image_padded_central(i3_image * image, int npad){
	i3_image * padded = i3_image_create(image->nx+2*npad, image->ny+2*npad);
	i3_image_fill(padded,0.0);
	i3_image_copy_into_part(image, npad, npad, padded);
	return padded;
}

i3_image * i3_image_padded_corner(i3_image * image, int npad){
	i3_image * padded = i3_image_create(image->nx+npad, image->ny+npad);
	i3_image_zero(padded);
	i3_image_copy_into_part(image, 0, 0, padded);
	return padded;
	
}

i3_image * i3_image_shear(i3_image * image, i3_flt x0, i3_flt y0, i3_flt e1, i3_flt e2){
	i3_image * sheared = i3_image_like(image);
	i3_image_shear_into(image,sheared,x0,y0,e1,e2);
	return sheared;
}


i3_flt cubicInterpolate (i3_flt p[4], i3_flt x) {
	return p[1] + (-0.5*p[0] + 0.5*p[2])*x + (p[0] - 2.5*p[1] + 2.0*p[2] - 0.5*p[3])*x*x + (-0.5*p[0] + 1.5*p[1] - 1.5*p[2] + 0.5*p[3])*x*x*x;
}

i3_flt bicubicInterpolate (i3_flt p[4][4], i3_flt x, i3_flt y) {
	i3_flt arr[4];
	arr[0] = cubicInterpolate(p[0], y);
	arr[1] = cubicInterpolate(p[1], y);
	arr[2] = cubicInterpolate(p[2], y);
	arr[3] = cubicInterpolate(p[3], y);
	return cubicInterpolate(arr, x);
}

i3_flt i3_image_bicubic_interpolate(i3_image * image, i3_flt x, i3_flt y){
	i3_flt p[4][4];
	int xl = floor(x)-1;
	int yl = floor(y)-1;
	i3_flt xi = x-xl+1;
	i3_flt yi = y-yl+1;
	
	for (int i=0;i<4;i++) 
		for(int j=0;j<4;j++)
		p[i][j]=image->row[xl-1+i][yl-1+j];
	return bicubicInterpolate(p,xi,yi);
}

i3_flt i3_image_bilinear_interpolate(i3_image * image, i3_flt x, i3_flt y){
	int xlow=floor(x);
	int xhigh=xlow+1;
	int ylow=floor(y);
	int yhigh=ylow+1;
	
	i3_flt rxlow = (x-xlow);
	i3_flt rylow = (y-ylow);
	i3_flt rxhigh = (xhigh-x);
	i3_flt ryhigh = (yhigh-y);
	printf("%f  %f     -  %d,  %d - %f  %f  %f  %f\n",x,y,xlow,ylow,rxlow,rxhigh, rylow,ryhigh);
	i3_flt value_xlow = rylow*image->row[xlow][ylow] + ryhigh*image->row[xlow][yhigh];
	i3_flt value_xhigh = rylow*image->row[xhigh][ylow] + ryhigh*image->row[xhigh][yhigh];
	return value_xlow * rxlow + value_xhigh * rxhigh;
}

void i3_image_truncate_circular(i3_image * image,i3_flt x0,i3_flt y0,i3_flt radius){
	int nx = image->nx;
	int ny = image->ny;
	i3_flt t2 = radius*radius;
	i3_flt x,y;
	for (int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			i3_image_xy_from_ij(image,i,j,&x,&y);
			i3_flt r2 = i3_pow(x-x0,2) + i3_pow(y-y0,2);
			if (r2>t2) image->row[j][i] = 0.0;
		}
	}
	
}

void i3_image_shear_into(i3_image * image, i3_image * sheared, i3_flt x0, i3_flt y0, i3_flt e1, i3_flt e2){
	I3_WARNING("SHEARING CODE COULD BE BETTER - currently uses bicubic interpolation - UNTESTED");
	int nx = image->nx;
	int ny = image->ny;
	int npad = (nx>ny)?nx:ny;
	i3_flt x,y;
	i3_image * padded = i3_image_padded_central(image,npad);
	for (int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			i3_image_xy_from_ij(image,i,j,&x,&y);
			
			i3_flt dx = x-x0;
			i3_flt dy = y-y0;
			i3_flt xi = (1-e1)*dx - e2*dy + x0;
			i3_flt yi = e2*dx + (1+e1)*dy + y0;
// #warning Check the signs of these.			
			// i3_flt xi = ((1+e1)*dx+e2*dy)/det + x0;
			// i3_flt yi = (e2*dx+(1-e1)*dy)/det + y0;
			sheared->row[j][i] = i3_image_bicubic_interpolate(padded,xi+npad,yi+npad); 
		}
	}
	
	i3_image_destroy(padded);
}


void i3_fftw_save_wisdom(){
	FILE * wisdom_file = fopen(I3_FFTW_WISDOM_FILENAME,"w");
	if (!wisdom_file) {I3_WARNING("Could not save fftw wisdom - could not open I3_FFTW_WISDOM_FILENAME"); return;};
#ifdef I3_FFTW_USE_DOUBLE
	fftw_export_wisdom_to_file(wisdom_file);
#else
	fftwf_export_wisdom_to_file(wisdom_file);
#endif	
}

void i3_fftw_load_wisdom(){
	FILE * wisdom_file = fopen(I3_FFTW_WISDOM_FILENAME,"r");
	if (!wisdom_file) {I3_WARNING("Could not load fftw wisdom - could not open I3_FFTW_WISDOM_FILENAME"); return;};
#ifdef I3_FFTW_USE_DOUBLE
	fftw_import_wisdom_from_file(wisdom_file);
#else
	fftwf_import_wisdom_from_file(wisdom_file);
#endif

}

void i3_image_convolve_downsample_into(i3_image * high_res, i3_fourier * kernel, i3_image * low_res){
	int upsampling = high_res->nx/low_res->nx;
	i3_image * image = i3_image_padded_corner(high_res, high_res->nx);
	i3_image_convolve_fourier(image, kernel);
	i3_image_extract_into(image, image->nx/2+1, image->nx/2+1, upsampling, upsampling, low_res);
	i3_image_destroy(image);
}

void i3_image_wiener_filter(i3_image * image, i3_image * model, i3_flt noise){
	for (int p=0;p<image->n;p++){
		image->data[p] *= model->data[p]/(model->data[p] + noise);
	}

}

void i3_image_dsample_into( i3_image * img_src, i3_image * img_dst ){
// TK 20110320 downsampling here means just selecting each nth element
// for now works only with odd number of subpixels and square images	
	int n_sub = img_src->nx / img_dst->nx ;

	// if(n_sub%2 == 0){
	// 	I3_FATAL("abort: please use odd number of subpixels for now",1);
	// }
	if(n_sub == 1){
		i3_image_copy_into( img_src, img_dst );	return;
	}
//	int N = img_dst->n;
	int nx = img_dst->nx;
	int ny = img_dst->ny;
	int dn_sub = n_sub/2;
	for (int j=0;j<ny;j++){
		for (int i=0;i<nx;i++){
			img_dst->row[j][i] = img_src->row[j*n_sub + dn_sub][i*n_sub + dn_sub];  
			
		}
	}
}

void i3_image_dsample_cut_into( i3_image * img_src, i3_image * img_dst, int i_start, int j_start, int n_sub ){
// TK 20110320 downsampling here means just selecting each nth element
// for now works only with odd number of subpixels and square images	
	//int n_sub = img_src->nx / img_dst->nx ;

	// if(n_sub%2 == 0){
	// 	// ifI3_FATAL("abort: please use odd number of subpixels or zero padding for now",1);
	// }
	if(n_sub == 1){
		i3_image_copy_into( img_src, img_dst );	return;
	}

	for (int j_dst=0; j_dst<img_dst->ny; j_dst++){
		for (int i_dst=0; i_dst<img_dst->nx; i_dst++){
			int i_src = i_dst*n_sub + i_start;
			int j_src = j_dst*n_sub + j_start;
			if (i_src>=img_src->nx || j_src>=img_src->nx) I3_FATAL("Bugger.",i_dst+j_dst);
			img_dst->row[j_dst][i_dst] = img_src->row[j_src][i_src];
		}
	}
}

i3_image * i3_image_dsample( i3_image * img_src, int n_sub){
// TK 20110320 downsampling here means just selecting each nth element
// for now works only with odd number of subpixels and square images
	int nx_low = img_src->nx/n_sub;
	int ny_low = img_src->ny/n_sub;		
	i3_image * img_dst = i3_image_create(nx_low,ny_low);
	i3_image_dsample_into( img_src,img_dst );
	return img_dst;
	 
}

void i3_fourier_dsample_into( i3_fourier * ft_src,  i3_fourier * ft_dst){

// TK 20110320 downsampling here means just selecting each nth element
// for now works only with odd number of subpixels and square images

// this still doesn't work properly, because the ft is not being aliased
	
	int n_sub = ft_src->ny_f / ft_dst->ny_f ;

	printf("%d\n",n_sub);

	if(n_sub%2 == 0){
		I3_FATAL("abort: please use odd number of subpixels for now",1);
	}
	if(n_sub == 1){
		i3_fourier_copy_into( ft_src, ft_dst );	return;
	}

	//int dn_sub = n_sub/2;
	for (int j=0;j<ft_dst->ny_f;j++){
		for (int i=0;i<ft_dst->nx_f;i++){
			ft_dst->row[j][i] = ft_src->row[j*n_sub][i*n_sub];  
			
		}
	}

}


void i3_convolve_real_same_into( i3_image * image, i3_image * kernel, i3_image * result){

// works with kernel which is of even number of pixels

int nx = image->nx;
int ny = image->ny; 
int ki = kernel->nx;
int kj = kernel->ny;

int half_i = ki/2;
int half_j = kj/2;

int xi = 0;
int yj = 0;

for(int x = 0; x < nx; x++){
	for(int y = 0; y < ny; y++){
		for(int i = 0; i<ki; i++){
			for(int j = 0; j<kj; j++){

				xi = x - half_i + i;
				yj = y - half_j + j;
  	
				if( xi >= 0 && xi < nx && yj >= 0 &&  yj < ny )
				{
					//printf("x = %d y = %d i = %d j = %d xi = %d yj = %d \n", x,y,i,j,xi,yj );
					result->row[x][y] += image->row[xi][yj] * kernel->row[i][j];
				}
				else
				{
					//printf("rejected x = %d y = %d i = %d j = %d xi = %d yj = %d \n", x,y,i,j,xi,yj );
				}  		

			}
		}			
		
	}
}



}

void i3_multiply_fourier_fourier_into(i3_fourier * f1, i3_fourier * f2, i3_fourier * result){
	if((!f1) || (!f2) || (!result)) return ;

	if (f1->nx_real!=f2->nx_real) {
		fprintf(stderr, "%ld, %ld, %ld, %ld\n", f1->nx_real, f2->nx_real, f1->ny_real, f2->ny_real);
		I3_FATAL("Incompatible x dimension sizes in i3_multiply_fourier_fourier_into",1);
	}
	if (f1->ny_real!=f2->ny_real) I3_FATAL("Incompatible y dimension sizes in i3_multiply_fourier_fourier_into",1);
	if (f1->nx_real!=result->nx_real) I3_FATAL("Incompatible x dimension sizes (result) in i3_multiply_fourier_fourier_into",1);
	if (f1->ny_real!=result->ny_real) I3_FATAL("Incompatible y dimension sizes (result) in i3_multiply_fourier_fourier_into",1);
	
	for(int i=0;i<f1->n_f;i++) result->data[i] = f1->data[i]*f2->data[i]; 
}


i3_image * i3_image_cut_out_stamp(i3_image * image, int stamp_size){
	//Figure out the amount of image to cut
	int cutx = image->nx - stamp_size;
	int cuty = image->ny - stamp_size;

	//If the image is smaller than the stamp then fail.
	if (cutx<0 || cuty<0) return NULL;
	
	//If the image is already at the stamp size then just return a copy of it.
	if (cutx==0 && cuty==0) return i3_image_copy(image);
	i3_image * output = i3_image_copy_part(image,cutx/2,cuty/2,stamp_size,stamp_size);
	if (output->nx!=stamp_size || output->ny!=stamp_size) I3_FATAL("Failed to cut out correct size of stamp (Joe's fault).",5);
	return output;
}

void i3_e12_to_etheta(i3_flt e1, i3_flt e2, i3_flt * e, i3_flt * theta)
{
  // BARNEY QUERY: SHOULD THERE BE FACTORS OF 2/0.5 HERE, AREN'T WE IN STANDARD POLAR COORDS RATHER THAN SPIN-2 PSEUDO VECTOR SPACE?
	*e = i3_sqrt(e1*e1+e2*e2);
	*theta = i3_atan2(e2,e1);
}

void i3_etheta_to_e12(i3_flt e, i3_flt theta, i3_flt * e1, i3_flt * e2)
{
 // BARNEY QUERY: SHOULD THERE BE FACTORS OF 2/0.5 HERE, AREN'T WE IN STANDARD POLAR COORDS RATHER THAN SPIN-2 PSEUDO VECTOR SPACE?
	*e1 = e * i3_cos(theta);
	*e2 = e * i3_sin(theta);
}

void i3_unit_shear_matrix(i3_flt e1, i3_flt e2, i3_flt * m1, i3_flt * m2, i3_flt * m3)
{
	i3_flt e_abs = (e1*e1+e2*e2);

	*m1 = (1 + e_abs - 2*e1)/(1-e_abs);
	*m2 = -2*e2/(1-e_abs);
	*m3 = (1 + e_abs + 2*e1)/(1-e_abs);
}
// Compare to the real space shear matrix the fourier version
// swaps the axis lengths "a" and "b".  This is equivalent to
// flipping the signs of e1 and e2.
// (Rowe, private communication)
void i3_unit_shear_matrix_fourier(i3_flt e1, i3_flt e2, i3_flt * m1, i3_flt * m2, i3_flt * m3)
{
	i3_unit_shear_matrix(-e1, -e2, m1, m2, m3);
}

void i3_unit_shear_matrix_old(i3_flt e1, i3_flt e2, i3_flt * m1, i3_flt * m2, i3_flt * m3)
{
	i3_flt e = i3_sqrt(e1*e1+e2*e2);
	i3_flt theta = 0.5*i3_atan2(e2,e1);
	i3_flt f = (1-e)/(1+e);
	
	// i3_flt r1=i3_sqrt(1.0*f);	/* The ellipse axes.  Sorry for the notation confusion. */
	// i3_flt r2=i3_sqrt(1.0/f);
	// i3_flt b = i3_pow(r1,-2);  /* Inverses of the ellipse axes, squared */
	// i3_flt a = i3_pow(r2,-2);
	
	i3_flt a = f;
	i3_flt b = 1./f;
	
        // SORRY FLIPPED a and b TEMPORARILY!

	i3_flt c = i3_cos(theta);
	i3_flt s = i3_sin(theta);
	
	/* The elements of the weight matrix ((m1,m2/2),(m2/2,m3))  */
	/* Put the minus signs in here so we do not need to later. */
	i3_flt c2 = c*c;
	i3_flt s2 = s*s;
	*m1 = a*c2+b*s2;	/* Note that a is actually 1/the usual a^2 */
	*m2 = (a-b)*s*c;
	*m3 = a*s2+b*c2;
}

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)
{
  /** 
   *  Used Mathematica (http://www.calc101.com/webMathematica/derivatives.jsp#topdoit)
   *  for computing the derivatives. Some code bits useful to have for copying
   *  (1-Sqrt[x^2+y^2])/(1+Sqrt[x^2+y^2])
   *  cos[0.5*arctan[x/y]]^2
   **/

  i3_flt e = e1*e1+e2*e2;
  i3_flt denom = (1-e)*(1-e);
  
  *dm1_e1 = -2*(1+e1*e1-e2*e2-2*e1)/denom;
  *dm2_e1 = -4*e1*e2/denom;
  *dm3_e1 =  2*(1+e1*e1-e2*e2+2*e1)/denom;
  
  *dm1_e2 =  4*(1-e1)*e2/denom;
  *dm2_e2 = -2*(1-e1*e1+e2*e2)/denom;
  *dm3_e2 =  4*(1+e1)*e2/denom;

}

void i3_unit_shear_matrix_jac_old(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)
{
  /** 
   *  Used Mathematica (http://www.calc101.com/webMathematica/derivatives.jsp#topdoit)
   *  for computing the derivatives. Some code bits useful to have for copying
   *  (1-Sqrt[x^2+y^2])/(1+Sqrt[x^2+y^2])
   *  cos[0.5*arctan[x/y]]^2
   **/

  if( (e1==0) && (e2 == 0) ){
        *dm1_e1 = -2.0;
	*dm2_e1 =  0.0;
	*dm3_e1 =  2.0;

	*dm1_e2 =  0.0;
	*dm2_e2 = -2.0;
	*dm3_e2 =  0.0;
  }
  else{
	i3_flt e = i3_sqrt(e1*e1+e2*e2);
	i3_flt theta = 0.5*i3_atan2(e2,e1);
	i3_flt f = (1-e)/(1+e);
	
	i3_flt a = f;
	i3_flt b = 1./f;
	
	i3_flt c = i3_cos(theta);
	i3_flt s = i3_sin(theta);

	i3_flt c2 = c*c;
	i3_flt s2 = s*s;
	
	// Note: all following derivatives are taken wrt e1 || e2

	i3_flt da = -2.0/(e*(e+1)*(e+1)); // * e1 || * e2
	i3_flt db =  2.0/(e*(e-1)*(e-1)); // * e1 || * e1

	i3_flt dc2 =  0.5*i3_sin(2*theta)/(e*e); // * e2 || * (-e1)
	i3_flt ds2 = -dc2; // * e2 || * (-e1)

	i3_flt dcs = 0.5*(s2-c2)/(e*e); // * e2 || * (-e1)

	*dm1_e1 = e1 * (da*c2 + db*s2) + e2 * (a*dc2 + b*ds2);
	*dm2_e1 = e1 * (da-db) * c * s + e2 * (a-b) * dcs;
	*dm3_e1 = e1 * (da*s2 + db*c2) + e2 * (a*ds2 + b*dc2);

	*dm1_e2 = e2 * (da*c2 + db*s2) - e1 * (a*dc2 + b*ds2);
	*dm2_e2 = e2 * (da-db) * c * s - e1 * (a-b) * dcs;
	*dm3_e2 = e2 * (da*s2 + db*c2) - e1 * (a*ds2 + b*dc2);
  }
}


void i3_e12_to_e12beermat(i3_flt e1in, i3_flt e2in, i3_flt emax, i3_flt * e1out, i3_flt * e2out)
{
  // Convert e1, e2 into polar coordinates (e, theta)
  //JAZ could switch this to the function above now
  i3_flt ein = i3_sqrt(e1in * e1in + e2in * e2in);
  i3_flt tin = i3_atan2(e2in, e1in); // NOTE NO FACTOR OF 0.5 HERE, THESE ARE STANDARD POLAR COORDS!

  // Define switches for locating what quadrant on the beer mat on which to place the e1 e2
  i3_flt de = (ein / emax) - (i3_flt) ((int) (ein / emax)) ; // Homemade % for floats
  de*=emax;
  int eswitch = ((int) (ein / emax)) % 2;
  int tswitch = ((int) (0.5 * ein / emax)) % 2;

  // Use these to build the beermat-ed location in e, theta polar coords
  i3_flt eout = eswitch * emax + pow(-1., eswitch) * de;
  i3_flt tout = tswitch * M_PI + tin;

  // Output in e1, e2
  *e1out = eout * i3_cos(tout);
  *e2out = eout * i3_sin(tout); // NOTE NO FACTOR OF 2 HERE, THESE ARE STANDARD POLAR COORDS!
}


void i3_r_to_rbeermat(i3_flt rin, i3_flt rparam, i3_flt * rout){
  if(i3_fabs(rin) < rparam){
    *rout = .5 * (rparam + rin * rin / rparam); 
  }else{
    *rout = i3_fabs(rin);
  }
}

i3_flt i3_image_get_fwhm(i3_image * image, i3_flt h, i3_flt x0){

	int nx = image->nx;
	
	// (x1,f1) and (x2,f2) are the indices of the pixels and pixel values, between which the fwhm (x3,f3), lies
	i3_flt f1 = 0;
	int x1 = 0;
	i3_flt f2 = 0;
	int x2 = 0;

	int s = nx/2+1;

	i3_flt * profile = malloc(sizeof(i3_flt)*s);
	i3_flt * diff = malloc(sizeof(i3_flt)*s);
	for(int i=0;i<s;i++) profile[i] = image->row[nx/2][i];	

	int max_ind = nx/2;
	i3_flt max_val = profile[max_ind];
	i3_flt f3 = max_val*h;

	for(int i=0;i<nx/2;i++)	diff[i] = i3_fabs(profile[i] - f3); 

	//printf("max_val fwxm f3 %2.4f %2.4f %2.4f \n",max_val,fwxm,f3);

	x1=i3_array_min_index( diff, nx/2 );
	f1=profile[x1];
	if( f1 < f3 ) 	x2 = x1+1;
	else 		x2 = x1-1;
 	f2 = profile[x2];



	i3_flt a = (f1-f2)/((i3_flt)x1 - (i3_flt)x2);
	i3_flt b = f1 - a*(i3_flt)x1; 
	i3_flt x3 = (f3 - b)/a;

	i3_flt fwm = 2.*i3_fabs(max_ind-x3);
	
	free(profile);
	free(diff);

	//printf("max_val f1 f2 f3 %2.6f %2.6f %2.6f %2.6f \n",max_val,f1,f2,f3);
	//printf("max_ind x1 x2 x3 %d %d %d %2.6f\n",max_ind,x1,x2,x3);
	//printf("xy0 fwm %2.4f %2.4f \n",x0,fwm);


	return fwm;

}

i3_flt i3_image_standard_fwhm(i3_image * image)
{
	int x0, y0;
	i3_image_maxloc(image, &x0, &y0);
	return i3_image_get_fwhm(image, 0.5, x0);
}

void i3_modify_weight_map_by_segmentation_mask(gal_id  identifier, i3_image * weight_stamp, i3_image * mask_stamp){

  if (mask_stamp->n != weight_stamp->n)
    I3_FATAL("abort in [i3_image.c / i3_modify_weight_map_by_segmentation_mask]: something is wrong with size of segmentation and weight mask",1);

  // check  first whether mask stamp contains pixels with values equal to identifier
  // if this is the case, well, great! If not assume galaxy is centered and take the
  // the center value of the mask stamp to mask galaxy 
  //printf("identifier %i \n", identifier);
  //printf("float identifier %f \n", (i3_flt) identifier);
  int count = 0;
  for (int i=0; i<mask_stamp->n; i++){
    if (mask_stamp->data[i] == (i3_flt) identifier){
      count += 1;
    }
  }
  //printf("count %i\n", count);
  
  gal_id  value;
  if (count > 0)
    value = identifier;
  else
    value = mask_stamp->data[mask_stamp->n/2];
  
  //printf("value %i\n",value);		  

  for (int i=0; i<mask_stamp->n; i++){
    //if (mask_stamp->data[i] == 0 || mask_stamp->data[i] == (i3_flt) identifier)
    if (mask_stamp->data[i] == 0 || mask_stamp->data[i] == (i3_flt) value)
      continue;
    else
      weight_stamp->data[i] = 0.;
  }
}

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

  if (mask_stamp->n != weight_stamp->n)
    I3_FATAL("abort in [i3_image.c / i3_modify_weight_map_by_segmentation_mask]: something is wrong with size of segmentation and weight mask",1);

  // check  first whether mask stamp contains pixels with values equal to identifier
  // if this is the case, well, great! If not assume galaxy is centered and take the
  // the center value of the mask stamp to mask galaxy
  // Also count number of pixels with non-zero value in segmentation map - will use this to initialize array 
  // in nearest neighbour check 
  int count_target = 0;
  int count_objs = 0;
  //printf("identifier %f\n",identifier);
  for (int i=0; i<mask_stamp->n; i++){
    if (mask_stamp->data[i] == identifier){
      count_target += 1;
    }
    if (mask_stamp->data[i] != 0){
      count_objs += 1;
    }
  }
  printf("number of obj pixels %i\n", count_objs);
  
  gal_id  value;
  if (count_target > 0)
    value = identifier;
  else
    value = mask_stamp->data[mask_stamp->n/2];		  
  
  // Make an array with index of each object pixel. Note this is a 1d index, i, where
  // x = i % stamp_size, y = i / stamp_size.
  int obj_inds[count_objs];
  int count = 0; //int count_zero = 0;
  for (int i=0; i<mask_stamp->n; i++){
  	if (mask_stamp->data[i] == 0){
  	  continue;	
  	}
  	else {
      obj_inds[count] = i;
      count++;
   	}
  }

  // Now for each pixel, check distance to all object pixels, and find
  // the closest
  for (int i=0; i<mask_stamp->n; i++){

    i3_flt min_dist_sq = stamp_size * stamp_size;
    int min_dist_ind = 0;

    for (int j=0; j<count_objs; j++){

      int this_x = i%stamp_size;
      int this_y = i/stamp_size;
      int obj_x = obj_inds[j]%stamp_size;
      int obj_y = obj_inds[j]/stamp_size;
      i3_flt dist_sq = pow((this_x - obj_x),2) + pow((this_y - obj_y),2);

      if (dist_sq < min_dist_sq){
      	min_dist_sq = dist_sq;
      	min_dist_ind = obj_inds[j];
      }
      else
      	continue;
    }

    // If the the closest object pixel is from the target, leave weight map = 1,
    // otherwise, set weight map = 0.
    
    if (mask_stamp->data[min_dist_ind] == (i3_flt) value){
    	continue;
    }
    else
      weight_stamp->data[i] = 0.;
  }
}


i3_flt i3_image_subtract_background(i3_image * galaxy){

	i3_flt background_estimate = 0;
	int image_size = galaxy->nx;
	int count = 0;
	int iiy=0;
	for (int iix=0; iix<image_size;iix++){
		background_estimate += galaxy->row[iix][iiy];
		count +=1;
	}
	iiy=image_size-1;
	for (int iix=0; iix<image_size;iix++){
		background_estimate += galaxy->row[iix][iiy];
		count +=1;
	}
		int iix=0;
	for (int iiy=0; iiy<image_size;iiy++){
		background_estimate += galaxy->row[iix][iiy];
		count +=1;
	}
	iix=image_size-1;
	for (int iiy=0; iiy<image_size;iiy++){
		background_estimate += galaxy->row[iix][iiy];
		count +=1;
	}
	background_estimate=background_estimate/ count;
	i3_image_addconst(galaxy,-background_estimate);		
	return background_estimate;

}

i3_flt i3_image_subtract_weighted_background(i3_image * galaxy, i3_image * weight){
	// Average over all pixels, accounting for weight.
	// Do not subtract, just return.

	if (galaxy->n != weight->n) I3_FATAL("abort in [i3_image.c / i3_image_subtract_background_using_mask]: something is wrong with size of galaxy and segmentation mask",1);

	i3_flt sum = 0.0;
	i3_flt w_sum = 0.0;

	int j=0;
	for (int i=0; i<galaxy->nx;i++){
		sum += galaxy->row[j][i]*weight->row[j][i];
		w_sum += weight->row[j][i];
	}

	j=galaxy->ny-1;
	for (int i=0; i<galaxy->nx;i++){
		sum += galaxy->row[j][i]*weight->row[j][i];
		w_sum += weight->row[j][i];
	}

	int i=0;
	for (int j=0; j<galaxy->ny;j++){
		sum += galaxy->row[j][i]*weight->row[j][i];
		w_sum += weight->row[j][i];
	}

	i=galaxy->nx-1;
	for (int j=0; j<galaxy->ny;j++){
		sum += galaxy->row[j][i]*weight->row[j][i];
		w_sum += weight->row[j][i];
	}

	i3_flt bg = sum/w_sum;

	if (!isnan(bg)) i3_image_addconst(galaxy,-bg);		

	return bg;
}


void i3_ellipticity_sum(i3_flt e1, i3_flt e2, i3_flt g1, i3_flt g2, i3_flt * out1, i3_flt * out2){
	// (e+g)/(1+g*e) =  y/x
	// y = e+g
	// x = 1+g*e
	i3_flt x_re = 1 + g1*e1+g2*e2;
	i3_flt x_im = g1*e2-g2*e1;
	
	i3_flt y_re=e1+g1;
	i3_flt y_im=e2+g2;
	
	// y/x = (y_re + i*y_im) / (x_re + i*x_im) 
	//     = (y_re+i*y_im)(x_re-i*x_im) / [(x_re + i*x_im)(x_re-i*x_im) ]
	//     = [(y_re*x_re+y_im*x_im) + i*(y_im*x_re-x_im*y_re)] / (x_re**2 + y_re**2)
	
	i3_flt bottom = x_re*x_re + x_im*x_im;
	i3_flt top_re = x_re*y_re + x_im*y_im;
	i3_flt top_im = x_re*y_im  -x_im*y_re;
	
	*out1 = top_re/bottom;
	*out2 = top_im/bottom;
}




void i3_e_linear_to_quadratic(i3_flt e1_lin, i3_flt e2_lin, i3_flt * e1_quad, i3_flt * e2_quad)
{
	i3_flt denom = 1+e1_lin*e1_lin+e2_lin*e2_lin;
	*e1_quad = 	2*e1_lin / denom;
	*e2_quad = 	2*e2_lin / denom;
}


void i3_e_quadratic_to_linear(i3_flt e1_quad, i3_flt e2_quad, i3_flt *e1_lin, i3_flt *e2_lin)
{
	i3_flt denom = 1+ i3_sqrt(1-(e1_quad*e1_quad+e2_quad*e2_quad));
	*e1_lin = e1_quad/denom;
	*e2_lin = e2_quad/denom;
}

i3_image * i3_create_tiled_image(int n, i3_image ** images){
	int nx=0;
	int ny=0;
	for (int i=0; i<n; i++){
		nx+=images[i]->nx;
		if (images[i]->ny > ny) ny = images[i]->ny;
	}
	i3_image * stack = i3_image_create(nx, ny);
	i3_image_zero(stack);
	i3_tile_images(stack, images, n);
	return stack;
}


int i3_tile_images(i3_image * target, i3_image ** tiles, int n)
{
    int x_start = 0;
    for (int i = 0; i< n; i++){
        int status = i3_image_copy_into_part(tiles[i], x_start, 0, target);
        if (status) return status;
        x_start += tiles[i]->nx;
    }
    return 0;
}


void i3_image_circular_mask(i3_image * image)
{
	float d = image->nx/2.0;
	//Central pixel
	for(int j=0;j<image->ny;j++){
		float dy = j-d;
		for(int i=0;i<image->nx;i++){
			float dx = i-d;
			if (dx*dx+dy*dy>d*d) image->row[j][i] = 0.0;
		}
	}
}