Initial import

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@1 bf6ca9ba-c028-0410-8290-897cf20841d1

1 files changed, 141 insertions, 0 deletions

diff --git a/src/mrc.c b/src/mrc.c
new file mode 100644
index 0000000..4815626
--- /dev/null
+++ b/src/mrc.c
@@ -0,0 +1,141 @@
+/*
+ * 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;
+	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;
+
+}