/*
 * 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"

/* 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;
	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;
	max = 0;
	for ( i=0; i<ctx->images->n_images; i++ ) {

		ImageRecord *image;

		image = &ctx->images->images[i];
		if ( image->rflist == NULL )
			image->rflist = reproject_get_reflections(image,
							ctx->cell_lattice);

		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;

			if ( feature->partner != NULL ) {

				if ( (h!=0) || (k!=0) || (l!=0) ) {

					double intensity;
					Reflection *ref;

					intensity = feature->partner->intensity;

					ref = reflectionlist_find(
						      ctx->integrated, h, k, l);

					if ( ref == NULL ) {

						Reflection *new;

						new = reflection_add(
							ctx->integrated,
							feature->reflection->x,
							feature->reflection->y,
							feature->reflection->z,
							intensity,
							REFLECTION_GENERATED);

						new->h = h;
						new->k = k;
						new->l = l;

						if ( intensity > max )
							max = intensity;

						n_meas++;

					} else {

						if ( intensity > ref->intensity ) {

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

						}

						n_dupl++;

					}

				}

			} else {
				//printf("IN: %3i %3i %3i not found\n", h, k, l);
				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 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, "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);
}