/*
 * mapping.c
 *
 * 3D Mapping
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

#include "reflections.h"
#include "control.h"
#include "itrans.h"
#include "image.h"

static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {

	/* "Input" space */
	double d, x, y;
	ImageRecord *imagerecord;
	
	/* Angular description of reflection */
	double theta, psi, k;
	
	/* Reciprocal space */
	double tilt;
	double omega;
	
	double x_temp, y_temp, z_temp;
	double nx, ny, nz;
	
	imagerecord = refl->parent;
	x = refl->x;	y = refl->y;
	tilt = 2*M_PI*(imagerecord->tilt/360);	/* Convert to Radians */
	omega = 2*M_PI*(imagerecord->omega/360);	/* Likewise */
	k = 1/imagerecord->lambda;
	
	/* Calculate an angular description of the reflection */
	if ( imagerecord->fmode == FORMULATION_CLEN ) {
		x /= imagerecord->resolution;
		y /= imagerecord->resolution;	/* Convert pixels to metres */
		d = sqrt((x*x) + (y*y));
		theta = atan2(d, imagerecord->camera_len);
	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
		x *= imagerecord->pixel_size;
		y *= imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
		d = sqrt((x*x) + (y*y));
		theta = atan2(d, k);
	} else {
		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
		return -1;
	}
	psi = atan2(y, x);
	
	x_temp = k*sin(theta)*cos(psi);
	y_temp = k*sin(theta)*sin(psi);
	z_temp = k - k*cos(theta);
	
	/* Apply the rotations...
		First: rotate image clockwise until tilt axis is aligned horizontally. */
	nx = x_temp*cos(omega) + y_temp*sin(omega);
	ny = -x_temp*sin(omega) + y_temp*cos(omega);
	nz = z_temp;
	
	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
	x_temp = nx; y_temp = ny; z_temp = nz;	
	nx = x_temp;
	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
	
	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
	x_temp = nx; y_temp = ny; z_temp = nz;	
	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
	nz = z_temp;
	
	*ddx = nx;
	*ddy = ny;
	*ddz = nz;
	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
	
	return 0;

}

ReflectionList *mapping_create(ControlContext *ctx) {

	int i;
	
	/* Pass all images through itrans */
	printf("MP: Analysing images..."); fflush(stdout);
	for ( i=0; i<ctx->images->n_images; i++ ) {
		ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
	}
	printf("done.\n");
	
	/* Create reflection list for measured reflections */
	ctx->reflectionlist = reflectionlist_new();
	printf("MP: Mapping to 3D..."); fflush(stdout);
	for ( i=0; i<ctx->images->n_images; i++ ) {
	
		int j;
		
		/* Iterate over the features in this image */
		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
		
			double nx, ny, nz, twotheta;
			
			if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
				reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
			}

		}
		
	}
	printf("done.\n");
	
	return ctx->reflectionlist;

}