/*
 * 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 "reflections.h"
#include "control.h"
#include "imagedisplay.h"
#include "itrans.h"

ReflectionContext *mapping_create(ControlContext *ctx) {

	ReflectionContext *rctx;
	int i;
	
	/* Create reflection context */
	rctx = reflection_init();
	
	if ( !ctx->prealign ) {
	
		int max_x, max_y;
		uint16_t max_val;
		int twidth, theight;
		int x, y;
		uint16_t *image_total;
		ImageDisplay *sum_id;
		
		/* Find centre of image stack
		 * Determine maximum size of image to accommodate, and allocate memory */
		twidth = 0;  theight = 0;
		for ( i=0; i<ctx->n_images; i++ ) {
			if ( ctx->images[i].width > twidth ) twidth = ctx->images[i].width;
			if ( ctx->images[i].height > theight ) theight = ctx->images[i].height;
		}
		image_total = malloc(twidth * theight * sizeof(uint16_t));
		memset(image_total, 0, twidth * theight * sizeof(uint16_t));
		
		/* Add the image stack together */
		for ( i=0; i<ctx->n_images; i++ ) {
			int xoffs, yoffs;
			xoffs = (twidth - ctx->images[i].width)/2;
			yoffs = (theight - ctx->images[i].height)/2;
			for ( y=0; y<theight; y++ ) {
				for ( x=0; x<twidth; x++ ) {
					image_total[x + twidth*y] += ctx->images[i].image[x + twidth*y]/ctx->n_images;
				}
			}
		}
		
		/* Locate the highest point */
		max_val = 0;	max_x = 0;	max_y = 0;
		for ( y=0; y<theight; y++ ) {
			for ( x=0; x<twidth; x++ ) {
				if ( image_total[x + twidth*y] > max_val ) {
					max_val = image_total[x + twidth*y];
					max_x = x;  max_y = y;
				}
			}
		}
		
		/* Record this measurement on all images */
		for ( i=0; i<ctx->n_images; i++ ) {
			ctx->images[i].x_centre = max_x;
			ctx->images[i].y_centre = max_y;
		}
		
		/* Display */
		sum_id = imagedisplay_open(image_total, twidth, theight, "Sum of All Images");
		imagedisplay_mark_point(sum_id, max_x, max_y);
		imagedisplay_add_tilt_axis(sum_id, max_x, max_y, ctx->omega);
		
	}
	
	/* Pass all images through itrans
	 *   (let itrans add the reflections to rctx for now) */
	ctx->reflectionctx = rctx;
	for ( i=0; i<ctx->n_images; i++ ) {
		itrans_process_image(ctx->images[i], ctx);
	}
	
	return rctx;

}