From e9a2b408c139a17e431dfb84384edd62a2ead7e3 Mon Sep 17 00:00:00 2001 From: taw27 Date: Thu, 30 Aug 2007 12:14:32 +0000 Subject: Merge the two reflection_add routines for different formulations git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@95 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/control.h | 6 ++-- src/itrans-lsq.c | 12 +++---- src/itrans-stat.c | 9 ++--- src/itrans-threshold.c | 27 +++----------- src/itrans-zaefferer.c | 10 ++---- src/reflections.c | 97 +++++++++++++++----------------------------------- src/reflections.h | 3 +- src/reproject.c | 17 ++++----- 8 files changed, 53 insertions(+), 128 deletions(-) (limited to 'src') diff --git a/src/control.h b/src/control.h index bae9f4b..ac3a357 100644 --- a/src/control.h +++ b/src/control.h @@ -50,8 +50,8 @@ typedef enum { typedef struct imagerecord_struct { uint16_t *image; - double tilt; - double omega; + double tilt; /* Degrees */ + double omega; /* Degrees */ FormulationMode fmode; double pixel_size; @@ -89,7 +89,7 @@ typedef struct cctx_struct { /* Basic parameters, stored here solely so they can be copied * into the ImageRecord(s) more easily */ double camera_length; - double omega; + double omega; /* Degrees */ double resolution; double lambda; double pixel_size; diff --git a/src/itrans-lsq.c b/src/itrans-lsq.c index 168a5b3..baac6b6 100644 --- a/src/itrans-lsq.c +++ b/src/itrans-lsq.c @@ -257,14 +257,12 @@ unsigned int itrans_peaksearch_lsq(ImageRecord imagerecord, ControlContext *ctx, int dx, dy; int bb_top, bb_bot, bb_lh, bb_rh; double frac; + double brightness; - printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", iter, x, y, ga, gb, gc, gd, ge, gf); - if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) { - reflection_add_from_reciprocal(ctx, (x-imagerecord.x_centre)*imagerecord.pixel_size, (y-imagerecord.y_centre)*imagerecord.pixel_size, - tilt_degrees, image[x + width*y]); - } else { - reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), tilt_degrees, image[x + width*y]); - } + printf("Fit converged after %i iterations: Centre %3i,%3i, a=%f b=%f c=%f d=%f e=%f f=%f\n", + iter, x, y, ga, gb, gc, gd, ge, gf); + brightness = image[x + width*y]; + reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), imagerecord, brightness); n_reflections++; /* Remove this peak from the image */ diff --git a/src/itrans-stat.c b/src/itrans-stat.c index f0677d4..7099c8d 100644 --- a/src/itrans-stat.c +++ b/src/itrans-stat.c @@ -495,7 +495,6 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx gsl_matrix *p; int i; double px,py; - double tilt_degrees = imagerecord.tilt; uint16_t *image = imagerecord.image; m = itrans_peaksearch_stat_createImageMatrix(image, imagerecord.width, imagerecord.height); @@ -505,12 +504,7 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx px = gsl_matrix_get(p,0,i); py = gsl_matrix_get(p,1,i); - if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) { - reflection_add_from_reciprocal(ctx, (px-imagerecord.x_centre)*imagerecord.pixel_size, (py-imagerecord.y_centre)*imagerecord.pixel_size, - tilt_degrees, 1.0); - } else { - reflection_add_from_dp(ctx, (px-imagerecord.x_centre)/imagerecord.resolution, (py-imagerecord.y_centre)/imagerecord.resolution, tilt_degrees, 1.0); - } + reflection_add_from_dp(ctx, (px-imagerecord.x_centre), (py-imagerecord.y_centre), imagerecord, 1.0); } gsl_matrix_free(m); @@ -519,3 +513,4 @@ unsigned int itrans_peaksearch_stat(ImageRecord imagerecord, ControlContext *ctx return n_reflections; } + diff --git a/src/itrans-threshold.c b/src/itrans-threshold.c index 71c39a1..519acc8 100644 --- a/src/itrans-threshold.c +++ b/src/itrans-threshold.c @@ -26,7 +26,6 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext int x, y; unsigned int n_reflections = 0; uint16_t *image = imagerecord.image; - double tilt_degrees = imagerecord.tilt; uint16_t max = 0; width = imagerecord.width; @@ -46,17 +45,8 @@ unsigned int itrans_peaksearch_threshold(ImageRecord imagerecord, ControlContext assert(y=0); assert(y>=0); - if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) { - reflection_add_from_reciprocal(ctx, - (signed)(x-imagerecord.x_centre)*imagerecord.pixel_size, - (signed)(y-imagerecord.y_centre)*imagerecord.pixel_size, - tilt_degrees, image[x + width*y]); - } else { - reflection_add_from_dp(ctx, - (signed)(x-imagerecord.x_centre)/imagerecord.resolution, - (signed)(y-imagerecord.y_centre)/imagerecord.resolution, - tilt_degrees, image[x + width*y]); - } + reflection_add_from_dp(ctx, (x-imagerecord.x_centre), (y-imagerecord.y_centre), + imagerecord, image[x + width*y]); n_reflections++; } } @@ -72,7 +62,6 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr int width, height; unsigned int n_reflections = 0; uint16_t *image; - double tilt_degrees = imagerecord.tilt; uint16_t max; int x, y; @@ -111,16 +100,8 @@ unsigned int itrans_peaksearch_adaptive_threshold(ImageRecord imagerecord, Contr assert(max_y=0); assert(max_y>=0); - if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) { - reflection_add_from_reciprocal(ctx, - (signed)(max_x-imagerecord.x_centre)*imagerecord.pixel_size, - (signed)(max_y-imagerecord.y_centre)*imagerecord.pixel_size, - tilt_degrees, image[max_x + width*max_y]); - } else { - reflection_add_from_dp(ctx, (signed)(max_x-imagerecord.x_centre)/imagerecord.resolution, - (signed)(max_y-imagerecord.y_centre)/imagerecord.resolution, - tilt_degrees, image[max_x + width*max_y]); - } + reflection_add_from_dp(ctx, (max_x-imagerecord.x_centre), (max_y-imagerecord.y_centre), + imagerecord, image[x + width*y]); n_reflections++; /* Remove it and its surroundings */ diff --git a/src/itrans-zaefferer.c b/src/itrans-zaefferer.c index 1454fd4..c6e7b39 100644 --- a/src/itrans-zaefferer.c +++ b/src/itrans-zaefferer.c @@ -82,14 +82,8 @@ unsigned int itrans_peaksearch_zaefferer(ImageRecord imagerecord, ControlContext } if ( !did_something ) { - if ( imagerecord.fmode == FORMULATION_PIXELSIZE ) { - reflection_add_from_reciprocal(ctx, (signed)(mask_x-imagerecord.x_centre)*imagerecord.pixel_size, - (signed)(mask_y-imagerecord.y_centre)*imagerecord.pixel_size, - tilt_degrees, image[mask_x + width*mask_y]); - } else { - reflection_add_from_dp(ctx, (signed)(mask_x-imagerecord.x_centre)/imagerecord.resolution, (signed)(mask_y-imagerecord.y_centre)/imagerecord.resolution, - tilt_degrees, image[mask_x + width*mask_y]); - } + reflection_add_from_dp(ctx, ((double)mask_x-imagerecord.x_centre), ((double)mask_y-imagerecord.y_centre), + imagerecord, image[mask_x + width*mask_y]); n_reflections++; } diff --git a/src/reflections.c b/src/reflections.c index d53af10..a1a8267 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -93,85 +93,46 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, } -/* x and y in metres (in real space!) */ -void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { +/* x and y in pixels, measured from centre of image */ +void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity) { - double nx, ny, nz; - - /* Real space */ - double L; + /* "Input" space */ double d; - - /* Shared */ - double theta; - double psi; - - /* Reciprocal space */ - double r; - 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); - psi = atan2(y, x); - - r = 1/ctx->lambda; - 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); - ny = x_temp*sin(omega) + y_temp*cos(omega); - nz = z_temp; - - /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way. - This is because the crystal is rotated in the experiment, not the Ewald sphere. */ - x_temp = nx; y_temp = ny; z_temp = nz; - nx = x_temp; - ny = cos(tilt)*y_temp - sin(tilt)*z_temp; - nz = -sin(tilt)*y_temp - cos(tilt)*z_temp; - - reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); - -} - -/* Alternative formulation for when the input data is already in reciprocal space units - x and y in metres^-1 (in reciprocal space) */ -void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { - - double nx, ny, nz; - - /* Input (reciprocal space) */ - double dr; - - /* Shared */ - double theta; - double psi; - double r; + /* Angular description of reflection */ + double theta, psi, k; /* Reciprocal space */ double tilt; double omega; - double x_temp, y_temp, z_temp; - tilt = M_PI*(tilt_degrees/180); /* Convert to Radians */ - omega = M_PI*(ctx->omega/180); /* Likewise */ + double x_temp, y_temp, z_temp; + double nx, ny, nz; - dr = sqrt((x*x) + (y*y)); - r = 1/ctx->lambda; - theta = atan2(dr, r); + tilt = 2*M_PI*(imagerecord.tilt/360); /* Convert to Radians */ + omega = 2*M_PI*(imagerecord.omega/360); /* Likewise */ + k = 1/imagerecord.lambda; + + /* Calculate an angular description of the reflection */ + if ( ctx->fmode == FORMULATION_CLEN ) { + x /= imagerecord.resolution; + y /= imagerecord.resolution; /* Convert pixels to metres */ + d = sqrt((x*x) + (y*y)); + theta = atan2(d, imagerecord.camera_len); + } else if ( ctx->fmode == FORMULATION_PIXELSIZE ) { + x *= imagerecord.pixel_size; + y *= imagerecord.pixel_size; /* Convert pixels to metres^-1 */ + d = sqrt((x*x) + (y*y)); + theta = atan2(d, k); + } else { + fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n"); + return; + } psi = atan2(y, x); - 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); + x_temp = k*sin(theta)*cos(psi); + y_temp = -k*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ + z_temp = k- k*cos(theta); /* Apply the rotations... First: rotate image clockwise until tilt axis is aligned horizontally. */ diff --git a/src/reflections.h b/src/reflections.h index 571f754..ae19362 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -76,8 +76,7 @@ extern Reflection *reflection_add(ReflectionContext *reflectionctx, double x, do #include "control.h" -extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); -extern void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity); +extern void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord imagerecord, double intensity); extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r); #endif /* REFLECTION_H */ diff --git a/src/reproject.c b/src/reproject.c index 1deb044..faa9433 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -26,8 +26,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect double smax = 0.5e9; double tilt, omega; - tilt = 2*M_PI*(image.tilt/360); - omega = 2*M_PI*(image.omega/360); /* Convert to Radians */ + tilt = deg2rad(image.tilt); + omega = deg2rad(image.omega); refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection)); @@ -88,8 +88,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect gy = yi - ky; gz = zi - kz; /* This is the vector from the centre of the sphere to the intersection */ theta = angle_between(-kx, -ky, -kz, gx, gy, gz); - printf("theta=%f ", theta); - theta = deg2rad(theta); cx = nz*gy - ny*gz; cy = nz*gx - nx*gz; cz = ny*gx - nx*gy; @@ -99,10 +97,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect ux = uyt*-sin(omega); uy = uyt*cos(omega); uz = uzt; - psi = angle_between(cx, cy, cz, ux, uy, uz); - psi -= 90; - printf("psi=%f\n", psi); - psi = deg2rad(psi); + psi = angle_between(cx, cy, cz, ux, uy, uz) - M_PI_2; if ( image.fmode == FORMULATION_CLEN ) { x = image.camera_len*sin(theta)*cos(psi); @@ -110,8 +105,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect x *= image.resolution; y *= image.resolution; } else if ( image.fmode == FORMULATION_PIXELSIZE ) { - x = 0; - y = 0; + x = sin(theta)*cos(psi) / image.lambda; + y = sin(theta)*sin(psi) / image.lambda; + x /= image.pixel_size; + y /= image.pixel_size; } else { fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n"); return NULL; -- cgit v1.2.3