diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-06 12:05:22 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-06 12:05:22 +0000 |
commit | 26736a29bfa67118143faf3ac840c7b89ce934e1 (patch) | |
tree | d5c64f9ce054159c62cbb34395ea63d4bd687e5c /src/reproject.c | |
parent | 787209d148adb0d7fd2aecfba641d9232bba6f52 (diff) |
Debodge reprojection and ImageDisplay
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@123 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/reproject.c')
-rw-r--r-- | src/reproject.c | 55 |
1 files changed, 38 insertions, 17 deletions
diff --git a/src/reproject.c b/src/reproject.c index 78cced6..06611cd 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -23,11 +23,12 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect ImageReflection *refl; Reflection *reflection; int i; - double smax = 0.5e9; + double smax = 0.25e9; double tilt, omega; double nx, ny, nz, nxt, nyt, nzt; /* "normal" vector (and calculation intermediates) */ double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */ double ux, uy, uz, uxt, uyt, uzt; /* "up" vector (and calculation intermediates) */ + //int done_debug = 0; tilt = deg2rad(image.tilt); omega = deg2rad(image.omega); @@ -38,7 +39,8 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect reflection = rctx->reflections; - /* Calculate the (normalised) incident electron wavevector */ + /* Calculate the (normalised) incident electron wavevector + * n is rotated "into" the reconstruction, so only one omega step. */ nxt = 0.0; nyt = 0.0; nzt = 1.0; nx = nxt; ny = cos(tilt)*nyt + sin(tilt)*nzt; nz = -sin(tilt)*nyt + cos(tilt)*nzt; nxt = nx; nyt = ny; nzt = nz; @@ -48,14 +50,13 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect kz = nz / image.lambda; /* This is the centre of the Ewald sphere */ reflection_add(ctx->reflectionctx, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); - /* Determine where "up" is */ + /* Determine where "up" is + * See above. */ uxt = 0.0; uyt = 1.0; uzt = 0.0; ux = uxt; uy = cos(tilt)*uyt + sin(tilt)*uzt; uz = -sin(tilt)*uyt + cos(tilt)*uzt; uxt = ux; uyt = uy; uzt = uz; ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt; - reflection_add(ctx->reflectionctx, ux*50, uy*50, uz*50, 1, REFLECTION_VECTOR_MARKER_2); - - printf("angle between k and u is %f\n", rad2deg(angle_between(kx, ky, kz, ux, uy, uz))); + reflection_add(ctx->reflectionctx, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); do { @@ -106,29 +107,49 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER); /* 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 */ + cx = yi*nz-zi*ny; cy = nx*zi-nz*xi; cz = ny*xi-nx*yi; /* c = i x n */ 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] */ + rx = cy*nz-cz*ny; ry = nx*cz-nz*cx; rz = ny*cx-nx*cy; /* r = [i x n] x n */ 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 */ + // if ( (i==20) && !done_debug ) { + // reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3); + // reflection_add(ctx->reflectionctx, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4); + // reflection_add(ctx->reflectionctx, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4); + // printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc)); + // } + if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) { + psi -= M_PI_2; /* Case 1 */ + // if ( (i==20) && !done_debug ) printf("case 1\n"); + } else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) { + psi = 3*M_PI_2 - psi; /* Case 2 */ + // if ( (i==20) && !done_debug ) printf("case 2\n"); + } else if ( (psi < M_PI_2) && (disc < M_PI_2) ) { + psi = 3*M_PI_2 - psi; /* Case 3 */ + // if ( (i==20) && !done_debug ) printf("case 3\n"); + } else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) { + psi = 3*M_PI_2 + psi; /* Case 4 */ + // if ( (i==20) && !done_debug ) printf("case 4\n"); + } + + // if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi)); psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ + // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi)); psi += omega; + // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi)); + // if ( (i==20) && !done_debug ) done_debug = 1; /* 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.camera_len*tan(theta)*cos(psi); + y = image.camera_len*tan(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 = tan(theta)*cos(psi) / image.lambda; + y = tan(theta)*sin(psi) / image.lambda; x /= image.pixel_size; y /= image.pixel_size; } else { @@ -139,7 +160,7 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect /* Adjust centre */ x += image.x_centre; y += image.y_centre; - x = image.width - 1 - x; /* Debodge me */ + /* Sanity check */ if ( (x>=0) && (x<image.width) && (y>=0) && (y<image.height) ) { |