/*
 * mrc.c
 *
 * Read the MRC tomography format
 *
 * (c) 2007 Thomas White <taw27@cam.ac.uk>
 *  dtr - Diffraction Tomography Reconstruction
 *
 */

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

#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>

#include "mrc.h"
#include "control.h"
#include "imagedisplay.h"
#include "itrans.h"
#include "reflections.h"

int mrc_read(ControlContext *ctx) {

	FILE *fh;
	MRCHeader mrc;
	MRCExtHeader ext[1024];
	unsigned int i;
	unsigned int extsize;
	double sx, sy, sz;
	unsigned int nx, ny, nz;
	double pixel_size;
	int16_t *image_total;
	unsigned int x, y;
	int max_x, max_y;
	int16_t max_val;
	ImageDisplay *sum_id;
	
	fh = fopen(ctx->filename, "rb");

	/* Read primary header */
	fread(&mrc, sizeof(MRCHeader), 1, fh);
	printf("%i images in series\n", mrc.nz);
	if ( mrc.mode != 1 ) {
		fprintf(stderr, "MR: Unknown MRC image mode\n");
		fclose(fh);
		return -1;
	}
	if ( (mrc.nxstart != 0) || (mrc.nystart != 0) ) {
		fclose(fh);
		fprintf(stderr, "MR: Image origin must be at zero: found at %i,%i\n", mrc.nxstart, mrc.nystart);
		return -1;
	}
	
	/* Read all extended headers, one by one */
	extsize = 4*mrc.numintegers + 4*mrc.numfloats;
	if ( extsize > sizeof(MRCExtHeader) ) {
		fclose(fh);
		fprintf(stderr, "MR: MRC extended header is too big - tweak mrc.h\n");
		return -1;
	}
	for ( i=0; i<mrc.nz; i++ ) {
		fread(&ext[i], extsize, 1, fh);
	}
	
	pixel_size = ext[0].pixel_size;
	printf("pixel_size=%f\n", pixel_size);
	nx = mrc.nx;
	ny = mrc.ny;
	nz = (mrc.nx > mrc.ny)?mrc.nx:mrc.ny;
	sx = nx * pixel_size;
	sy = ny * pixel_size;
	sz = nz * pixel_size;
	ctx->reflectionctx = reflection_init();
	ctx->fmode = FORMULATION_PIXELSIZE;
	ctx->first_image = 1;
	ctx->width = mrc.nx;
	ctx->height = mrc.ny;
	
	printf("Judging centre...\n");
	image_total = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
	for ( y=0; y<mrc.ny; y++ ) {
		for ( x=0; x<mrc.nx; x++ ) {
			image_total[x + mrc.nx*y] = 0;
		}
	}
	for ( i=0; i<mrc.nz; i++ ) {
		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
		fread(image, mrc.nx*mrc.ny*2, 1, fh);
		for ( y=0; y<mrc.ny; y++ ) {
			for ( x=0; x<mrc.nx; x++ ) {
				image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz;
			}
		}
	}
	sum_id = imagedisplay_open(image_total, mrc.nx, mrc.ny, "Sum of All Images");
	/* Locate the highest point */
	max_val = 0;	max_x = 0;	max_y = 0;
	for ( y=0; y<ctx->height; y++ ) {
		for ( x=0; x<ctx->width; x++ ) {
			if ( image_total[x + ctx->width*y] > max_val ) {
				max_val = image_total[x + ctx->width*y];
				max_x = x;  max_y = y;
			}
		}
	}
	imagedisplay_mark_point(sum_id, max_x, max_y);
	ctx->x_centre = max_x;
	ctx->y_centre = max_y;
	
	for ( i=0; i<mrc.nz; i++ ) {
	
		int16_t *image = malloc(mrc.ny * mrc.nx * sizeof(int16_t));
		
		printf("Image #%3i: tilt=%f omega=%f L=%f\t", i, ext[i].a_tilt, ext[i].tilt_axis, ext[i].magnification);
		ctx->camera_length = ext[i].magnification;
		ctx->lambda = 2.51e-12;	/* 200kV. Fudged until Max puts the HT voltage in the MRC headers */
		ctx->omega = ext[i].tilt_axis;
		ctx->pixel_size = ext[i].pixel_size;
		
		fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET);
		fread(image, mrc.nx*mrc.ny*2, 1, fh);

		for ( y=0; y<mrc.ny; y++ ) {
			for ( x=0; x<mrc.nx; x++ ) {
				image_total[x + mrc.nx*y] += image[x + mrc.nx*y]/mrc.nz;
			}
		}

		itrans_process_image(image, ctx, ext[i].a_tilt);
		free(image);
		
	}

	fclose(fh);
	
	return 0;

}