diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-02-05 21:12:57 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-02-05 21:12:57 +0000 |
commit | 05b5d261682b9136fb46476a64eab6980b0dba64 (patch) | |
tree | d7faa450b69cf2104ffff7fc89a95914284d23af /src/mrc.c |
Initial import
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@1 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/mrc.c')
-rw-r--r-- | src/mrc.c | 141 |
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; + +} |