/*
 * intensities.c
 *
 * Extract intensities from patterns
 *
 * (c) 2006-2010 Thomas White <taw@physics.org>
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


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

#include "image.h"
#include "intensities.h"
#include "cell.h"
#include "sfac.h"
#include "diffraction.h"


#define MAX_HITS (1024)


struct reflhit {
	signed int h;
	signed int k;
	signed int l;
	double min_distance;
	int x;
	int y;
};


static double sum_nearby_points(float *data, int width, int x, int y)
{
	int dx, dy;
	double intensity = 0;

	for ( dx=-3; dx<=3; dx++ ) {
	for ( dy=-3; dy<=3; dy++ ) {
		intensity += data[(x+dx) + width*(y+dy)];
	}
	}

	return intensity;
}


void output_intensities(struct image *image, UnitCell *cell)
{
	int x, y;
	double ax, ay, az;
	double bx, by, bz;
	double cx, cy, cz;
	struct reflhit hits[MAX_HITS];
	int n_hits = 0;
	int i;

	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);

	for ( x=0; x<image->width; x++ ) {
	for ( y=0; y<image->height; y++ ) {

		double hd, kd, ld;  /* Indices with decimal places */
		double dh, dk, dl;  /* Distances in h,k,l directions */
		signed int h, k, l;
		struct rvec q;
		double dist;
		int found = 0;
		int j;

		q = get_q(image, x, y, 1, NULL, 1.0/image->lambda);

		hd = q.u * ax + q.v * ay + q.w * az;
		kd = q.u * bx + q.v * by + q.w * bz;
		ld = q.u * cx + q.v * cy + q.w * cz;

		h = (signed int)rint(hd);
		k = (signed int)rint(kd);
		l = (signed int)rint(ld);

		dh = hd - h;
		dk = kd - k;
		dl = ld - l;
		dist = sqrt(pow(dh, 2.0) + pow(dk, 2.0) + pow(dl, 2.0));
		if ( dist > 0.1 ) continue;

		for ( j=0; j<n_hits; j++ ) {
			if ( (hits[j].h == h) && (hits[j].k == k)
			                      && (hits[j].l == l) ) {

				if ( dist < hits[j].min_distance ) {
					hits[j].min_distance = dist;
					hits[j].x = x;
					hits[j].y = y;
				}

				found = 1;

			}
		}

		if ( !found ) {
			hits[n_hits].min_distance = dist;
			hits[n_hits].x = x;
			hits[n_hits].y = y;
			hits[n_hits].h = h;
			hits[n_hits].k = k;
			hits[n_hits].l = l;
			n_hits++;
			assert(n_hits < MAX_HITS);
		}

	}
	}

	STATUS("Found %i reflections\n", n_hits);

	/* Explicit printf() used here (not normally allowed) because
	 * we really want to output to stdout */
	printf("New pattern: %7.5f %7.5f %7.5f %7.5f\n",
	       image->orientation.w, image->orientation.x,
	       image->orientation.y, image->orientation.z);
	for ( i=0; i<n_hits; i++ ) {

		double intensity;

		/* Bounds check */
		if ( hits[i].x + 3 >= image->width ) continue;
		if ( hits[i].x - 3 < 0 ) continue;
		if ( hits[i].y + 3 >= image->height ) continue;
		if ( hits[i].y - 3 < 0 ) continue;

		intensity = sum_nearby_points(image->data, image->width,
		                              hits[i].x, hits[i].y);

		printf("%3i %3i %3i %6f (at %i,%i)\n",
		       hits[i].h, hits[i].k, hits[i].l, intensity,
		       hits[i].x, hits[i].y);

	}

	/* Blank line at end */
	printf("\n");
}