diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-04 11:06:39 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-04 11:06:39 +0000 |
commit | 49d079959cac6df46d3d66fd8391a7dc67b3d721 (patch) | |
tree | 59d6e6393b2bc1d92cf777070a07043aa7102434 | |
parent | 5778a311d57a9f8b2347f616e8dba81c97f0a018 (diff) |
Resolve +/-Pi ambiguity
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@107 bf6ca9ba-c028-0410-8290-897cf20841d1
-rw-r--r-- | src/displaywindow.c | 35 | ||||
-rw-r--r-- | src/ipr.c | 2 | ||||
-rw-r--r-- | src/reflections.h | 5 | ||||
-rw-r--r-- | src/reproject.c | 121 |
4 files changed, 112 insertions, 51 deletions
diff --git a/src/displaywindow.c b/src/displaywindow.c index d087bf4..45dbc7f 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -448,6 +448,8 @@ static void displaywindow_gl_create_list(ControlContext *ctx) { reflection = ctx->reflectionctx->reflections; quadric = gluNewQuadric(); do { + + double vmul;; if ( reflection->type == REFLECTION_CENTRAL ) { @@ -459,8 +461,39 @@ static void displaywindow_gl_create_list(ControlContext *ctx) { gluSphere(quadric, 0.2, 32, 32); glPopMatrix(); - break; /* Assume there's only one. */ + } + + if ( reflection->type == REFLECTION_VECTOR_MARKER_1 ) { + + vmul = 100; + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); + glBegin(GL_LINES); + glVertex3f(0.0, 0.0, 0.0); + glVertex3f(reflection->x*vmul, reflection->y*vmul, reflection->z*vmul); + glEnd(); + + } + + if ( reflection->type == REFLECTION_VECTOR_MARKER_2 ) { + + vmul = 1/1e9; + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, yellow); + glBegin(GL_LINES); + glVertex3f(0.0, 0.0, 0.0); + glVertex3f(reflection->x*vmul, reflection->y*vmul, reflection->z*vmul); + glEnd(); + + } + + if ( reflection->type == REFLECTION_VECTOR_MARKER_3 ) { + vmul = 1/1e9; + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue); + glBegin(GL_LINES); + glVertex3f(0.0, 0.0, 0.0); + glVertex3f(reflection->x*vmul, reflection->y*vmul, reflection->z*vmul); + glEnd(); + } reflection = reflection->next; @@ -263,7 +263,7 @@ int ipr_refine(ControlContext *ctx) { /* Select an image */ i = ipr_random_image(ctx); - i = 45; /* Temporary! */ + i = 58; /* Temporary! */ cur = imagedisplay_open(ctx->images[i], "Current Image", IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS); refl = reproject_get_reflections(ctx->images[i], &n, lat, ctx); diff --git a/src/reflections.h b/src/reflections.h index ae19362..5c0a1d7 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -20,7 +20,10 @@ typedef enum { REFLECTION_NORMAL, /* Measured - expressed as x, y, z position */ REFLECTION_CENTRAL, /* The central beam */ REFLECTION_GENERATED, /* Generated and intensity-measured - expressed as h, k, l index */ - REFLECTION_MARKER + REFLECTION_MARKER, + REFLECTION_VECTOR_MARKER_1, + REFLECTION_VECTOR_MARKER_2, + REFLECTION_VECTOR_MARKER_3 } ReflectionType; typedef struct reflection_struct { diff --git a/src/reproject.c b/src/reproject.c index faa9433..5594451 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -49,8 +49,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect /* Now calculate the (normalised) incident electron wavevector */ xt = 0; - yt = - sin(tilt); - zt = - cos(tilt); + yt = sin(tilt); + zt = cos(tilt); nx = xt*cos(omega) + yt*-sin(omega); ny = xt*sin(omega) + yt*cos(omega); nz = zt; @@ -71,16 +71,15 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect double kx, ky, kz; double cx, cy, cz; double ux, uy, uz; - double uxt, uyt, uzt; - double theta, psi; + double uyt, uzt; + double theta; double x, y; + double rx, ry, rz; /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; - reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_CENTRAL); - reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER); - /* Calculate Bragg angle and azimuth of point in image */ + /* Calculate Bragg angle */ kx = nx / image.lambda; ky = ny / image.lambda; kz = nz / image.lambda; /* This is the centre of the Ewald sphere */ @@ -88,54 +87,80 @@ 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); - cx = nz*gy - ny*gz; - cy = nz*gx - nx*gz; - cz = ny*gx - nx*gy; - uxt = 0; - uyt = cos(tilt); - uzt = -sin(tilt); - ux = uyt*-sin(omega); - uy = uyt*cos(omega); - uz = uzt; - 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); - y = image.camera_len*sin(theta)*sin(psi); - x *= image.resolution; - y *= image.resolution; - } else if ( image.fmode == FORMULATION_PIXELSIZE ) { - 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; - } + if ( theta > 0.0 ) { - /* Adjust centre */ - x += image.x_centre; - y += image.y_centre; + double psi, disc; - /* Sanity check */ - if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) { - - /* Record the reflection */ - refl[i].x = x; - refl[i].y = y; - i++; + /* Determine where "up" is */ + uyt = cos(tilt); uzt = -sin(tilt); /* tilt it (uxt not needed) */ + ux = uyt*-sin(omega); uy = uyt*cos(omega); uz = uzt; /* rotate it */ + + /* Calculate azimuth of point in image (clockwise from "up", will be changed later) */ + cx = nz*yi - ny*zi; cy = nx*zi - nz*xi; cz = ny*xi - nx*yi; /* c = (-n) x i */ + psi = angle_between(cx, cy, cz, ux, uy, uz); + + /* Finally, resolve the +/- Pi ambiguity from the previous step */ + rx = nz*cy - ny*cz; ry = nx*cz - nz*cx; rz = ny*cx - nx*cy; /* r = (-n) x [(-n) x i] */ + disc = angle_between(rx, ry, rz, ux, uy, uz); + if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) psi -= M_PI_2; /* Case 1 */ + else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) psi = 3*M_PI_2 - psi; /* Case 2 */ + else if ( (psi < M_PI_2) && (disc < M_PI_2) ) psi = 3*M_PI_2 - psi; /* Case 3 */ + else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) psi = 3*M_PI_2 + psi; /* Case 4 */ + psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ + if ( (reflection->h == 2) && (reflection->k == 2) && (reflection->l == 0) ) { + printf("c = %f, %f, %f\n", cx, cy, cz); + reflection_add(ctx->reflectionctx, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_2); + printf("r = %f, %f, %f\n", rx, ry, rz); + reflection_add(ctx->reflectionctx, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_3); + reflection->type = REFLECTION_MARKER; + printf("n = %f, %f, %f\n", kx, ky, kz); + reflection_add(ctx->reflectionctx, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); + } else { + reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_GENERATED); + reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER); + } - if ( i > MAX_IMAGE_REFLECTIONS ) { - fprintf(stderr, "Too many reflections\n"); - break; + /* Calculate image coordinates from polar representation */ + if ( image.fmode == FORMULATION_CLEN ) { + x = image.camera_len*sin(theta)*cos(psi); + y = image.camera_len*sin(theta)*sin(psi); + x *= image.resolution; + y *= image.resolution; + } else if ( image.fmode == FORMULATION_PIXELSIZE ) { + 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; } - //printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y); + /* Adjust centre */ + x += image.x_centre; + y += image.y_centre; + + /* Sanity check */ + if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) { + + /* Record the reflection */ + refl[i].x = x; + refl[i].y = y; + i++; + + if ( i > MAX_IMAGE_REFLECTIONS ) { + fprintf(stderr, "Too many reflections\n"); + break; + } + + //printf("Reflection %i at %i,%i\n", i, refl[i-1].x, refl[i-1].y); + + } else { + //fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y); + } - } else { - //fprintf(stderr, "Reflection failed sanity test (x=%f, y=%f)\n", x, y); - } + } /* else this is the central beam so don't worry about it */ } |