/*
 * intensities.c
 *
 * Extract integrated intensities by relrod estimation
 *
 * (c) 2007 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"

/* 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 *new;
					
					/* Perform relrod calculation of doom here.
					 * TODO: Figure out if this is even possible. */
					intensity = feature->partner->intensity;
					
					new = reflection_add(ctx->integrated,
							feature->reflection->x, feature->reflection->y, feature->reflection->z,
							intensity, REFLECTION_GENERATED);
					
					if ( new != NULL ) {
						new->h = h;
						new->k = k;
						new->l = l;
						//printf("IN: Adding %3i %3i %3i, intensity=%f\n", h, k, l, intensity);
						if ( intensity > max ) max = intensity;
						n_meas++;
					} else {
						printf("IN: Duplicate measurement for %3i %3i %3i\n", h, k, l);
						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);

}