/*
 * intensities.c
 *
 * Extract integrated intensities by relrod estimation
 *
 * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
 *
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

#include "control.h"
#include "reflections.h"
#include "image.h"
#include "reproject.h"
#include "displaywindow.h"
#include "utils.h"
#include "basis.h"

/* Add a Friedel equivalent if missing, otherwise set the
 *  intensities of both equivalents to the larger intensity */
static void intensity_friedelise(Reflection *ref, ReflectionList *list,
				 int *n_fried)
{
	Reflection *equiv;

	equiv = reflectionlist_find(list, -ref->h, -ref->k, -ref->l);

	/* Friedel equivalent exists? */
	if ( equiv == NULL ) {

		/* No, add it */
		equiv = reflection_add(list, -ref->x, -ref->y, -ref->z,
					ref->intensity, REFLECTION_GENERATED);

		equiv->h = -ref->h;
		equiv->k = -ref->k;
		equiv->l = -ref->l;

		(*n_fried)++;

	} else {

		/* Yes, use largest intensity */
		if ( equiv->intensity > ref->intensity )
			ref->intensity = equiv->intensity;
		else
			equiv->intensity = ref->intensity;

	}
}

static void intensity_measure(signed int h, signed int k, signed int l,
				ImageFeature *feature, ReflectionList *list,
				double *max,
				int *n_meas, int *n_dupl, int *n_fried)
{
	if ( (h!=0) || (k!=0) || (l!=0) ) {

		double intensity;
		Reflection *ref;

		intensity = feature->partner->intensity;

		ref = reflectionlist_find(list, h, k, l);

		/* Do we already have a value for this reflection? */
		if ( ref == NULL ) {

			/* No, add it */
			ref = reflection_add(list,
					     feature->reflection->x,
					     feature->reflection->y,
					     feature->reflection->z,
					     intensity, REFLECTION_GENERATED);

			ref->h = h;
			ref->k = k;
			ref->l = l;

			if ( intensity > *max ) *max = intensity;

			(*n_meas)++;

		} else {

			/* Yes, record only if the new value is larger */
			if ( intensity > ref->intensity ) {

				ref->x = feature->reflection->x;
				ref->y = feature->reflection->y;
				ref->z = feature->reflection->z;
				ref->intensity = intensity;

			}

			(*n_dupl)++;

		}

		intensity_friedelise(ref, list, n_fried);

	}
}

/* Extract integrated reflection intensities by estimating the spike function
 * based on the observed intensity and the calculated excitation error from
 * the lattice refinement.  Easy.	*/
void intensities_extract(ControlContext *ctx)
{
	int i, j;
	int n_meas, n_dupl, n_notf, n_fried;
	double max;
	Reflection *reflection;

	/* Free previous analysis if required */
	if ( ctx->integrated != NULL ) {
		reflectionlist_free(ctx->integrated);
	}
	ctx->integrated = reflectionlist_new();

	n_meas = 0;
	n_dupl = 0;
	n_notf = 0;
	n_fried = 0;
	max = 0;
	for ( i=0; i<ctx->images->n_images; i++ ) {

		ImageRecord *image;

		/* Reproject this pattern, if necessary */
		image = &ctx->images->images[i];
		if ( image->rflist == NULL )
			image->rflist = reproject_get_reflections(image,
							ctx->cell_lattice);

		/* For each reprojected feature */
		for ( j=0; j<image->rflist->n_features; j++ ) {

			ImageFeature *feature;
			signed int h, k, l;

			feature = &image->rflist->features[j];

			h = feature->reflection->h;
			k = feature->reflection->k;
			l = feature->reflection->l;

			/* Corresponds to a measured reflection? */
			if ( feature->partner != NULL )
				intensity_measure(h, k, l, feature,
						  ctx->integrated, &max,
						  &n_meas, &n_dupl, &n_fried);
			else
				n_notf++;
		}

	}

	/* Normalise all reflections to the most intense reflection */
	reflection = ctx->integrated->reflections;
	while ( reflection ) {
		reflection->intensity /= max;
		reflection = reflection->next;
	}

	printf("IN: %5i intensities measured\n", n_meas);
	printf("IN: %5i duplicated measurements\n", n_dupl);
	printf("IN: %5i reflections generated by Friedel's law\n", n_fried);
	printf("IN: %5i predicted reflections not found\n", n_notf);
}

static int intensities_do_save(ReflectionList *integrated, Basis *cell,
				const char *filename)
{
	FILE *fh;
	Reflection *reflection;
	UnitCell rcell;

	fh = fopen(filename, "w");

	rcell = basis_get_cell(cell);
	fprintf(fh, "#INT This file contains reflection intensities\n");
	fprintf(fh, "a %12.8f\n", rcell.a*1e9);
	fprintf(fh, "b %12.8f\n", rcell.b*1e9);
	fprintf(fh, "c %12.8f\n", rcell.c*1e9);
	fprintf(fh, "alpha %12.8f\n", rad2deg(rcell.alpha));
	fprintf(fh, "beta %12.8f\n", rad2deg(rcell.beta));
	fprintf(fh, "gamma %12.8f\n", rad2deg(rcell.gamma));

	reflection = integrated->reflections;
	while ( reflection ) {
		fprintf(fh, "%3i %3i %3i %12.8f\n",
			reflection->h, reflection->k, reflection->l,
			reflection->intensity);
		reflection = reflection->next;
	}
	fclose(fh);

	return 0;
}

static gint intensities_save_response(GtkWidget *widget, gint response,
					ControlContext *ctx)
{
	if ( response == GTK_RESPONSE_ACCEPT ) {
		char *filename;
		filename = gtk_file_chooser_get_filename(
						GTK_FILE_CHOOSER(widget));
		if ( intensities_do_save(ctx->integrated,
						ctx->cell, filename) ) {
			displaywindow_error("Failed to save cache file.",
						ctx->dw);
		}
		g_free(filename);
	}

	gtk_widget_destroy(widget);

	return 0;
}

void intensities_save(ControlContext *ctx)
{
	GtkWidget *save;

	save = gtk_file_chooser_dialog_new("Save Reflections to File",
					GTK_WINDOW(ctx->dw->window),
					GTK_FILE_CHOOSER_ACTION_SAVE,
					GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
					GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
					NULL);
	g_signal_connect(G_OBJECT(save), "response",
				G_CALLBACK(intensities_save_response), ctx);
	gtk_widget_show_all(save);
}