1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
/*
* mapping.c
*
* 3D Mapping
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "reflections.h"
#include "control.h"
#include "itrans.h"
#include "image.h"
static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
/* "Input" space */
double d, x, y;
ImageRecord *imagerecord;
/* Angular description of reflection */
double theta, psi, k;
/* Reciprocal space */
double tilt;
double omega;
double x_temp, y_temp, z_temp;
double nx, ny, nz;
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 ( 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 (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 -1;
}
psi = atan2(y, x);
x_temp = k*sin(theta)*cos(psi);
y_temp = k*sin(theta)*sin(psi);
z_temp = k - k*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;
/* 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;
*ddx = nx;
*ddy = ny;
*ddz = nz;
*twotheta = theta; /* Radians. I've used the "wrong" nomenclature above */
return 0;
}
ReflectionList *mapping_create(ControlContext *ctx) {
int i;
/* Pass all images through itrans */
printf("MP: Analysing images..."); fflush(stdout);
for ( i=0; i<ctx->images->n_images; i++ ) {
ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
}
printf("done.\n");
/* Create reflection list for measured reflections */
ctx->reflectionlist = reflectionlist_new();
printf("MP: Mapping to 3D..."); fflush(stdout);
for ( i=0; i<ctx->images->n_images; i++ ) {
int j;
/* Iterate over the features in this image */
for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
double nx, ny, nz, twotheta;
if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
}
}
}
printf("done.\n");
return ctx->reflectionlist;
}
|
