diff options
Diffstat (limited to 'src/reproject.c')
| -rw-r--r-- | src/reproject.c | 46 |
1 files changed, 25 insertions, 21 deletions
diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..3bfc20a 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -45,7 +45,7 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { } -ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) { +ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) { ImageFeatureList *flist; Reflection *reflection; @@ -54,7 +54,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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; + int done_debug = 0; + int i = 0; flist = image_feature_list_new(); @@ -70,7 +71,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * kx = nx / image->lambda; ky = ny / image->lambda; kz = nz / image->lambda; /* This is the centre of the Ewald sphere */ - //reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); + reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1); /* Determine where "up" is * See above. */ @@ -78,7 +79,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); + reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); reflection = reflectionlist->reflections; do { @@ -125,8 +126,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * double psi, disc; - //reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED); - //reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER); + reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED); + reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER); /* Calculate azimuth of point in image (clockwise from "up", will be changed later) */ cx = yi*nz-zi*ny; cy = nx*zi-nz*xi; cz = ny*xi-nx*yi; /* c = i x n */ @@ -135,33 +136,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * /* Finally, resolve the +/- Pi ambiguity from the previous step */ 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 ( (i==20) && !done_debug ) { - // reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3); - // reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4); - // reflection_add(ctx->reflectionlist, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4); - // printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc)); - // } + if ( (i==20) && !done_debug ) { + reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3); + reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4); + reflection_add(ctx->reflectionlist, 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"); + 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"); + 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"); + 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("case 4\n"); } - // if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi)); + 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)); + 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; + 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 ) { @@ -197,7 +201,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * } /* else this is the central beam so don't worry about it */ } - + i++; reflection = reflection->next; } while ( reflection ); |