From 74d58e5f81c7b010aa4db3627f119f9df47c10fa Mon Sep 17 00:00:00 2001 From: taw27 Date: Thu, 23 Aug 2007 10:17:19 +0000 Subject: Implement fast reprojection git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@72 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/control.c | 32 +++++++--------- src/control.h | 12 ++++-- src/itrans-zaefferer.c | 6 --- src/main.c | 2 +- src/mapping.c | 2 + src/mrc.c | 2 +- src/qdrp.c | 22 ----------- src/readpng.c | 2 +- src/reflections.c | 6 +-- src/reproject.c | 99 +++++++++++++++++++++++++++++++++++++++++++++++++- src/reproject.h | 2 +- 11 files changed, 128 insertions(+), 59 deletions(-) diff --git a/src/control.c b/src/control.c index b626d7c..e0318c0 100644 --- a/src/control.c +++ b/src/control.c @@ -13,31 +13,25 @@ #include "control.h" -int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size) { +int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt) { ctx->images[ctx->n_images].tilt = tilt; - ctx->images[ctx->n_images].omega = omega; - ctx->images[ctx->n_images].pixel_size = pixel_size; - ctx->images[ctx->n_images].camera_len = 0; + ctx->images[ctx->n_images].omega = ctx->omega; ctx->images[ctx->n_images].image = image; ctx->images[ctx->n_images].width = width; ctx->images[ctx->n_images].height = height; + ctx->images[ctx->n_images].lambda = ctx->lambda; + ctx->images[ctx->n_images].fmode = ctx->fmode; - ctx->n_images++; - - return ctx->n_images - 1; - -} - -int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len) { - - ctx->images[ctx->n_images].tilt = tilt; - ctx->images[ctx->n_images].omega = omega; - ctx->images[ctx->n_images].pixel_size = 0; - ctx->images[ctx->n_images].camera_len = camera_len; - ctx->images[ctx->n_images].image = image; - ctx->images[ctx->n_images].width = width; - ctx->images[ctx->n_images].height = height; + if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + ctx->images[ctx->n_images].pixel_size = ctx->pixel_size; + ctx->images[ctx->n_images].camera_len = 0; + ctx->images[ctx->n_images].resolution = ctx->resolution; + } else if ( ctx->fmode == FORMULATION_CLEN ) { + ctx->images[ctx->n_images].pixel_size = 0; + ctx->images[ctx->n_images].camera_len = ctx->camera_length; + ctx->images[ctx->n_images].resolution = 0; + } ctx->n_images++; diff --git a/src/control.h b/src/control.h index 1952dab..83ba0fb 100644 --- a/src/control.h +++ b/src/control.h @@ -52,10 +52,17 @@ typedef struct imagerecord_struct { int16_t *image; double tilt; double omega; + + FormulationMode fmode; double pixel_size; double camera_len; + double lambda; + double resolution; + int width; int height; + int x_centre; + int y_centre; } ImageRecord; @@ -87,12 +94,10 @@ typedef struct cctx_struct { unsigned int started; unsigned int camera_length_set; unsigned int omega_set; - unsigned int max_d_set; unsigned int resolution_set; unsigned int lambda_set; /* Miscellaneous modes and operations */ - unsigned int max_d; unsigned int first_image; /* Size of input images - assumed the same throughout. */ @@ -114,8 +119,7 @@ typedef struct cctx_struct { } ControlContext; -extern int control_add_image_ps(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double pixel_size); -extern int control_add_image_cl(ControlContext *ctx, int16_t *image, int width, int height, double tilt, double omega, double camera_len); +extern int control_add_image(ControlContext *ctx, int16_t *image, int width, int height, double tilt); #endif /* CONTROL_H */ diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c index 8ac845b..aaa8110 100644 --- a/src/itrans-zaefferer.c +++ b/src/itrans-zaefferer.c @@ -34,12 +34,6 @@ unsigned int itrans_peaksearch_zaefferer(int16_t *image, ControlContext *ctx, do double dx1, dx2, dy1, dy2; double dxs, dys; double grad; - - if ( ctx->max_d ) { - if ( (x-ctx->x_centre)*(x-ctx->x_centre) + (y-ctx->y_centre)*(y-ctx->y_centre) > ctx->max_d*ctx->max_d ) { - continue; - } - } /* Get gradients */ dx1 = image[x+ctx->width*y] - image[(x+1)+ctx->width*y]; diff --git a/src/main.c b/src/main.c index 4b42809..4e5723a 100644 --- a/src/main.c +++ b/src/main.c @@ -208,7 +208,7 @@ int main(int argc, char *argv[]) { #endif ctx->filename = strdup(argv[1]); - ctx->max_d = 0; + ctx->n_images = 0; if ( stat(argv[1], &stat_buffer) == -1 ) { fprintf(stderr, "File '%s' not found\n", argv[1]); diff --git a/src/mapping.c b/src/mapping.c index 0a0dcdf..a8d50b3 100644 --- a/src/mapping.c +++ b/src/mapping.c @@ -75,6 +75,8 @@ ReflectionContext *mapping_create(ControlContext *ctx) { * (let itrans add the reflections to rctx for now) */ ctx->reflectionctx = rctx; for ( i=0; in_images; i++ ) { + ctx->images[i].x_centre = ctx->x_centre; + ctx->images[i].y_centre = ctx->y_centre; itrans_process_image(ctx->images[i].image, ctx, ctx->images[i].tilt); } diff --git a/src/mrc.c b/src/mrc.c index 26a31de..abc096c 100644 --- a/src/mrc.c +++ b/src/mrc.c @@ -85,7 +85,7 @@ int mrc_read(ControlContext *ctx) { fseek(fh, mrc.next + sizeof(MRCHeader) + mrc.nx*mrc.ny*2*i, SEEK_SET); fread(image, mrc.nx*mrc.ny*2, 1, fh); - control_add_image_ps(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt, ctx->omega, ctx->pixel_size); + control_add_image(ctx, image, mrc.nx, mrc.ny, ext[i].a_tilt); } diff --git a/src/qdrp.c b/src/qdrp.c index 6cc9296..96adfda 100644 --- a/src/qdrp.c +++ b/src/qdrp.c @@ -166,26 +166,6 @@ static int qdrp_parseline(ControlContext *ctx, const char *line) { } - if ( strncasecmp(line+pos, "max_d", 5) == 0 ) { - - char *max_d_s; - - skip_chars; - skip_whitespace; - mark = pos; - skip_chars; - max_d_s = strndup(line+mark, pos-mark); - ctx->max_d = strtod(max_d_s, NULL); - free(max_d_s); - ctx->max_d_set = 1; - if ( ctx->max_d == 0 ) { - printf("QD: max_d=0 is Not Allowed!\n"); - return 1; - } - printf("QD: max_d = %i pixels\n", ctx->max_d); - - } - if ( strncasecmp(line+pos, "omega", 5) == 0 ) { char *omega_s; @@ -216,10 +196,8 @@ int qdrp_read(ControlContext *ctx) { ctx->camera_length_set = 0; ctx->omega_set = 0; - ctx->max_d_set = 0; ctx->resolution_set = 0; ctx->lambda_set = 0; - ctx->max_d = 0; ctx->started = 0; ctx->first_image = 1; diff --git a/src/readpng.c b/src/readpng.c index 68a92c3..20fc182 100644 --- a/src/readpng.c +++ b/src/readpng.c @@ -141,7 +141,7 @@ int readpng_read(const char *filename, double tilt, ControlContext *ctx) { png_destroy_read_struct(&png_ptr, &info_ptr, &end_info); fclose(fh); - control_add_image_cl(ctx, image, width, height, tilt, ctx->omega, ctx->camera_length); + control_add_image(ctx, image, width, height, tilt); return 0; diff --git a/src/reflections.c b/src/reflections.c index 03d901b..e16bcff 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -111,10 +111,10 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt double tilt; double omega; double x_temp, y_temp, z_temp; - + tilt = 2*M_PI*(tilt_degrees/360); /* Convert to Radians */ omega = 2*M_PI*(ctx->omega/360); /* Likewise */ - + d = sqrt((x*x) + (y*y)); L = ctx->camera_length; theta = atan2(d, L); @@ -124,7 +124,7 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt x_temp = r*sin(theta)*cos(psi); y_temp = -r*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ z_temp = r- r*cos(theta); - + /* Apply the rotations... First: rotate image clockwise until tilt axis is aligned horizontally. */ nx = x_temp*cos(omega) + y_temp*-sin(omega); diff --git a/src/reproject.c b/src/reproject.c index f62e8f3..c9023af 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -10,21 +10,118 @@ */ #include +#include #include "control.h" #include "reflections.h" #define MAX_IMAGE_REFLECTIONS 1024 -ImageReflection *reproject_get_reflections(ImageRecord *image, size_t *n, ReflectionContext *ref) { +ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx) { ImageReflection *refl; + Reflection *reflection; int i; + double smax = 1e9; + double tilt, omega; + + tilt = 2*M_PI*(image.tilt/360); + omega = 2*M_PI*(image.omega/360); /* Convert to Radians */ refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection)); i = 0; + reflection = rctx->reflections; + do { + + double xl, yl, zl; + double nx, ny, nz; + double xt, yt, zt; + double a, b, c; + double s1, s2, s; + + /* Get the coordinates of the reciprocal lattice point */ + xl = reflection->x; + yl = reflection->y; + zl = reflection->z; + + /* Now calculate the (normalised) incident electron wavevector */ + xt = 0; + yt = - sin(tilt); + zt = - cos(tilt); + nx = xt*cos(omega) + yt*-sin(omega); + ny = xt*sin(omega) + yt*cos(omega); + nz = zt; + + /* Next, solve the reprojection equation to calculate the excitation error */ + a = nx*nx + ny*ny + nz*nz; + b = 2*(xl*nx + yl*ny + zl*nz - nz/image.lambda); + c = xl*xl + yl*yl + zl*zl - 2*zl/image.lambda; + s1 = (-b + sqrt(b*b-4*a*c))/(2*a); + s2 = (-b - sqrt(b*b-4*a*c))/(2*a); + if ( s1 < s2 ) s = s1; else s = s2; + + /* Skip this reflection if s is large */ + if ( s <= smax ) { + + double xddd, yddd, zddd; + double xdd, ydd, zdd; + double xd, yd, zd; + double theta, psi; + double x, y; + + /* Determine the intersection point */ + xddd = xl + s*nx; yddd = yl + s*ny; zddd = zl + s*nz; + + /* Invert the image->3D mapping to get the image coordinates */ + xdd = xddd; + ydd = (yddd/cos(tilt) - zddd*tan(tilt)/cos(tilt))/(1+tan(tilt)*tan(tilt)); + zdd = (-zddd-ydd*sin(tilt))/cos(tilt); + + yd = (ydd-xdd*tan(omega))/(sin(omega)*tan(omega)+cos(omega)); + xd = (xdd+yd*sin(omega))/cos(omega); + zd = zdd; + + if ( image.fmode == FORMULATION_CLEN ) { + psi = atan2(-yd, xd); + theta = acos(1-zd*image.lambda); + x = image.camera_len*sin(theta)*cos(psi); + y = image.camera_len*sin(theta)*sin(psi); + x *= image.resolution; + } else if ( image.fmode == FORMULATION_PIXELSIZE ) { + x = xd / image.pixel_size; + y = yd / image.pixel_size; + } else { + fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n"); + return NULL; + } + + /* Adjust centre */ + x += image.x_centre; + y += image.y_centre; + + /* Sanity check */ + if ( (x>=0) && (x=0) && (y MAX_IMAGE_REFLECTIONS ) break; + + printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y); + + } else { + fprintf(stderr, "Reflection failed sanity test\n"); + } + + } + + reflection = reflection->next; + + } while ( reflection ); *n = i; return refl; diff --git a/src/reproject.h b/src/reproject.h index c65549c..ac7b240 100644 --- a/src/reproject.h +++ b/src/reproject.h @@ -19,7 +19,7 @@ #include "control.h" #include "reflections.h" -extern ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *ref); +extern ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, ReflectionContext *rctx); #endif /* REPROJECT_H */ -- cgit v1.2.3