/*
 * renderer.c
 *
 * Render Fourier Transform output array into display array
 *
 * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2D - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <fftw3.h>

#include "renderer.h"

double renderer_width(int width, int height, double gamma, int nx, int ny) {
	return nx*width + ny*height*fabs(cos(gamma));
}

double renderer_height(int width, int height, double gamma, int nx, int ny) {
	return height*sin(gamma)*ny;
}

/* Provide the coordinate in the image (ie xn,yn) given x,y in the Fourier array */
double renderer_map_x(double x, double y, int width, int height, double gamma, int nx, int ny) {

	double xn, offs;

	if ( gamma > M_PI_2 ) {
		offs = (double)ny*height*cos(M_PI-gamma);
	} else {
		offs = 0.0;
	}
	xn = x + offs + y*cos(gamma);
	
	return xn;

}

/* Same, but return y */
double renderer_map_y(double x, double y, int width, int height, double gamma, int nx, int ny) {
	return y*sin(gamma);
}

static double renderer_map_reverse_x(double xn, double yn, int width, int height, double gamma, int nx, int ny) {

	double offs, xd, yd;

	if ( gamma > M_PI_2 ) {
		offs = (double)ny*height*cos(M_PI-gamma);
	} else {
		offs = 0.0;
	}
	
	yd = (double)yn / sin(gamma);
	xd = (double)xn - offs - yd*cos(gamma);
	
	while ( xd < 0 ) xd += width;
	xd = fmod(xd, width);
	
	return xd;
	
}

static double renderer_map_reverse_y(double xn, double yn, int width, int height, double gamma, int nx, int ny) {

	double yd;

	yd = (double)yn / sin(gamma);
	while ( yd < 0 ) yd += height;
	yd = fmod(yd, height);
	
	return yd;
	
}

static double renderer_interpolate_linear_re(double xd, double yd, fftw_complex *out, int width, int height) {

	double frac, re, re1, re2;
	int x, y;

	/* Get the left-hand point value */
	x = floor(xd);  y = floor(yd);
	if ( x >= width ) x -= width;
	if ( y >= height ) y -= height;
	re1 = out[(height-1-y)+height*(width-1-x)][0];
	
	/* Get the right-hand point value */
	x++;
	if ( x >= width ) x -= width;
	re2 = out[(height-1-y)+height*(width-1-x)][0];
	
	frac = fmod(xd, 1);
	re = (1-frac)*re1 + frac*re2;
	
	return re;

}

static double renderer_interpolate_linear_im(double xd, double yd, fftw_complex *out, int width, int height) {

	double frac, im, im1, im2;
	int x, y;

	/* Get the left-hand point value */
	x = floor(xd);  y = floor(yd);
	if ( x >= width ) x -= width;
	if ( y >= height ) y -= height;
	im1 = out[(height-1-y)+height*(width-1-x)][1];
	
	/* Get the right-hand point value */
	x++;
	if ( x >= width ) x -= width;
	im2 = out[(height-1-y)+height*(width-1-x)][1];
	
	frac = fmod(xd, 1);
	im = (1-frac)*im1 + frac*im2;
	
	return im;

}

static double renderer_interpolate_bilinear_re(double xd, double yd, fftw_complex *out, int width, int height) {

	double frac, re, re1, re2;
	
	/* Get the lower interpolated value */
	re1 = renderer_interpolate_linear_re(xd, yd, out, width, height);
	
	/* Get the upper interpolated value */
	yd++;
	if ( yd >= height ) yd -= height;
	re2 = renderer_interpolate_linear_re(xd, yd, out, width, height);
	
	frac = fmod(yd, 1);
	re = (1-frac)*re1 + frac*re2;
	
	return re;

}

static double renderer_interpolate_bilinear_im(double xd, double yd, fftw_complex *out, int width, int height) {

	double frac, im, im1, im2;
	
	/* Get the lower interpolated value */
	im1 = renderer_interpolate_linear_im(xd, yd, out, width, height);
	
	/* Get the upper interpolated value */
	yd++;
	if ( yd >= height ) yd -= height;
	im2 = renderer_interpolate_linear_im(xd, yd, out, width, height);
	
	frac = fmod(yd, 1);
	im = (1-frac)*im1 + frac*im2;
	
	return im;

}

ComplexArray renderer_draw(fftw_complex *out, int width, int height, double gamma, int nx, int ny) {

	ComplexArray cxar;
	int width_n, height_n;
	int xn, yn;
	
	width_n = (int)renderer_width(width, height, gamma, nx, ny);
	height_n = (int)renderer_height(width, height, gamma, nx, ny);
	
	cxar.re = malloc(height_n*width_n*sizeof(double));
	cxar.im = malloc(height_n*width_n*sizeof(double));
	
	for ( yn=0; yn<height_n; yn++ ) {
		for ( xn=0; xn<width_n; xn++ ) {
	
			double xd, yd;
			double re, im;
			
			/* Map this pixel onto the FFTW output array */
			xd = renderer_map_reverse_x(xn, yn, width, height, gamma, nx, ny);
			yd = renderer_map_reverse_y(xn, yn, width, height, gamma, nx, ny);
			
			/* Extract (interpolatively) the value */
			re = renderer_interpolate_bilinear_re(xd, yd, out, width, height);
			im = renderer_interpolate_bilinear_im(xd, yd, out, width, height);
			
			/* Store the values */
			cxar.re[xn+width_n*yn] = re;
			cxar.im[xn+width_n*yn] = im;
			
		}
	}
	
	return cxar;

}