diff options
Diffstat (limited to 'src/reflections.c')
-rw-r--r-- | src/reflections.c | 102 |
1 files changed, 86 insertions, 16 deletions
diff --git a/src/reflections.c b/src/reflections.c index 1a80f71..63eaa54 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -13,11 +13,14 @@ #include <config.h> #endif +#define _GNU_SOURCE #include <stdlib.h> #include <math.h> #include <stdio.h> +#include <string.h> #include "reflections.h" +#include "utils.h" static void reflection_addfirst(ReflectionContext *reflectionctx) { /* Create first items on lists - saves faffing later. Corresponds to a central marker. @@ -124,11 +127,11 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, } -/* 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) { +int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta) { /* "Input" space */ - double d; + double d, x, y; + ImageRecord *imagerecord; /* Angular description of reflection */ double theta, psi, k; @@ -140,24 +143,26 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord double x_temp, y_temp, z_temp; double nx, ny, nz; - tilt = 2*M_PI*(imagerecord.tilt/360); /* Convert to Radians */ - omega = 2*M_PI*(imagerecord.omega/360); /* Likewise */ - k = 1/imagerecord.lambda; + imagerecord = refl->parent; + x = refl->x; y = refl->y; + 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 */ + if ( imagerecord->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 */ + theta = atan2(d, imagerecord->camera_len); + } else if (imagerecord->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; + return -1; } psi = atan2(y, x); @@ -177,14 +182,36 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord nx = x_temp; ny = cos(tilt)*y_temp + sin(tilt)*z_temp; nz = -sin(tilt)*y_temp + cos(tilt)*z_temp; - + /* Finally, reverse the omega rotation to restore the location of the image in 3D space */ x_temp = nx; y_temp = ny; z_temp = nz; nx = x_temp*cos(-omega) + y_temp*sin(-omega); ny = -x_temp*sin(-omega) + y_temp*cos(-omega); nz = z_temp; - reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + *ddx = nx; + *ddy = ny; + *ddz = nz; + *twotheta = theta; /* Radians. I've used the "wrong" nomenclature above */ + + return 0; + +} + +/* 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, twotheta; + ImageReflection refl; + + refl.parent = imagerecord; + refl.x = x; + refl.y = y; + refl.intensity = intensity; + + if ( !reflection_map_to_space(&refl, &nx, &ny, &nz, &twotheta) ) { + reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + } } @@ -194,3 +221,46 @@ void reflection_add_from_reflection(ReflectionContext *rctx, Reflection *r) { rctx->last_reflection = r; } +double reflection_largest_g(ReflectionContext *reflectionctx) { + + double max = 0.0; + Reflection *reflection; + + reflection = reflectionctx->reflections; + while ( reflection ) { + if ( reflection->type == REFLECTION_NORMAL ) { + double mod; + mod = modulus(reflection->x, reflection->y, reflection->z); + if ( mod > max ) max = mod; + } + reflection = reflection->next; + }; + + return max; + +} + +Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z) { + + double max = +INFINITY; + Reflection *reflection; + Reflection *best = NULL; + + reflection = reflectionctx->reflections; + while ( reflection ) { + if ( reflection->type == REFLECTION_NORMAL ) { + double mod; + mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z); + if ( mod < max ) { + max = mod; + best = reflection; + } + } + reflection = reflection->next; + }; + + return best; + + +} + |